matlab语言与应用02：基础知识与编程入门

《高等应用数学问题的MATLAB求解（第四版）》由薛定宇编写，东北大学提供相关课程资源，书中代码可在MATLAB R2016b版本运行。核心内容为现代科学运算中MATLAB语言的应用与高等应用数学问题的求解。

02 matlab语言程序设计基础

02.01 数据结构

变量，字母开头，区分大小写。

   保留常量

   eps, i, j, pi, NaN, Inf, i = sqrt(-1)

   lastwarn, lasterr

数值型，双精度 64位，8字节

   single(), 32位单精度。

   符号型, sym (A),常用于公式推导和求解析解。

   变量声明 syms var_list var_props

   显示符号变量的任何精度:vpa(A), vpa(A, n)

   默认精度–32位10进制有效数字

   显示符号变量的属性 assumptions()

   设置符号变量的属性 assume(), assumeAlso()

syms x real; assume(x >= -1); assumeAlso(x < 5);

例2-1 符号型数值与双精度数值的区别

   1/3 的双精度存储

   双精度存储 – 0.3333…33, 15位有限数字

   符号型存储 – 1/3:sym(1/3),始终存储1/3

a = 1/3; b = sym(1/3);a, b

例2-2 圆周率 𝜋  π

   求出𝜋  π 的300位有效数字。

vpa(pi, 300)

例2-3 试用符号型数据结构表示数值12345678901234567890

A = sym(12345678901234567890)

%% 显示结果:sym()内的数值被截断了
A =
12345678901234567168

应该这样:

B = sym('12345678901234567890')

%%结果正确
B =
12345678901234567890

02.02 矩阵与向量的输入方法

直接赋值语句: variable = expression

   保留变量:ans,存放最近依次无赋值变量语言的运算结果

赋值语句的末尾加一个分号可以阻止显示运算结果。

例 2-5 矩阵输入方法
𝐴=147258369  A= [    1   2   3     4   5   6     7   8   9    ]  

A = [1, 2, 3; 4, 5, 6; 7, 8, 9]

%% 动态定维
A = [1, 2, 3; 4 5, 6; 7, 8 9];
A = [[A; [1 2 3]], [1; 2; 3; 4]]
inv(A) % 求矩阵A的逆矩阵

例 2-6 复数矩阵输入
𝐵=1+𝑗94+𝑗67+𝑗32+𝑗85+𝑗58+𝑗23+𝑗76+𝑗40+𝑗1  B= [    1 + j 9   2 + j 8   3 + j 7     4 + j 6   5 + j 5   6 + j 4     7 + j 3   8 + j 2   0 + j 1    ]  

B = [1+9i 2+8i, 3+7i; 4+6i 5+5i, 6+4i; 7 + 3i , 8 + 2i, 0 + 1i]

% 注意空格，当1 +9i ,+前有空格时，会提示[串联的矩阵的维度不一致。]错误
B = [1 +9i 2+8i, 3+7i; 4+6i 5+5i, 6+4i; 7 + 3i , 8 + 2i, 0 + 1i]
% [1 + 9i]+前后空格正确
B = [1 + 9i 2+8i, 3+7i; 4+6i, 5+5i, 6+4i; 7 + 3i, 8 + 2i, 0 + 1i]

函数调用

函数调用  语句

   [returned.arguments] = function_name(input_arguments)

函数调用: [U S V] = svd(X)

   函数可以通过不同方式被调用:

   内核函数: *.m 函数

   匿名函数、inline 函数(不建议使用)

   重载函数、私有函数等

冒号表达式
𝑣=𝑠1:𝑠2:𝑠3  v= s 1 : s 2 : s 3  
开始与𝑠1,步长𝑠2,终止于𝑠3   开  始  与  s 1 , 步  长  s 2 , 终  止  于  s 3  

   默认步长为1

例2-7 冒号 表达式  生成

   用不同的步距生成 𝑡∈[0,𝜋]  t∈[0,π]  间的向量

v1 = 0: 0.2: pi % 不包含pi
v1a = linspace(0, pi, 50) % 包含pi
v2 = 0: -0.1: pi % 空的 1×0 double 行矢量
v3 = 0: pi % 0     1     2     3
v4 = pi: -1: 0 %3.1416    2.1416    1.1416    0.1416

子矩阵提取

B = A(v1, v2)

v1 表示子矩阵要保留的行号构成的向量

   v2表示子矩阵的列好构成的向量

   :,表示要提取所有行或列,取决于其位置

   end 的使用

% matlab 子矩阵获取
B = A(1:2:end, :), C = A([1 1 1 1], :)

02.03 矩阵的代数运算

矩阵表示: 矩阵A，n行m列，被称作n x m 矩阵
转置

% matlab 共轭转置(Hermite转置)
C = A'

% matlab 一般转置
C = A .'

𝐵=𝐴𝑇,𝑏𝑗𝑖=𝑎𝑖𝑗,𝑖=1,⋯,𝑛,𝑗=1,⋯,𝑚  B=AT, b  j i  = a  i j  ,i=1,⋯,n,j=1,⋯,m

加减法

   数学表示:𝑐𝑖𝑗=𝑎𝑖𝑗±𝑏𝑖𝑗,𝑖=1,⋯,𝑛,𝑗=1,⋯,𝑚   c  i j  = a  i j  ± b  i j  ,i=1,⋯,n,j=1,⋯,m

% 加法matlab表示
C = A + B    C = A - B

任一变量可为标量；如果矩阵A、B维度不同，matlab报错。

乘法

   数学表示: 𝐶=𝐴𝐵  C=AB
𝑐𝑖𝑗=∑𝑘=1𝑚𝑎𝑖𝑘𝑏𝑘𝑗,𝑖=1,⋯,𝑛,𝑗=1,⋯,𝑚   c  i j  = ∑  k = 1   m   a  i k   b  k j  ,i=1,⋯,n,j=1,⋯,m

% 乘法 matlab 表示
C = A * B

除法

   矩阵左除:

   求解线性方程组:AX = B

% matlab 求解左除
X = A \B

最小二乘解；若A为非奇异方阵，则 𝑋=𝐴−1𝐵  X=A − 1 B

   使得误差最小: min‖𝐴𝑋−𝐵‖2  min‖AX−B ‖ 2  

矩阵右除:

   求解线性方程组: XA = B

%矩阵右除
X = B / A

等效的运算: B / A = (A' \ B')'

矩阵翻转

   左右翻转: 𝑏𝑖𝑗=𝑎𝑖,𝑛+1−𝑗   b  i j  = a  i , n + 1 − j   , B = fliplr(A)

   上下翻转:𝑐𝑖𝑗=𝑎𝑚1−𝑖,𝑗   c  i j  = a   m 1  − i , j   , B = flipud(A)

   旋转90°:𝑑𝑖𝑗=𝑎𝑗,𝑛+1−𝑖   d  i j  = a  j , n + 1 − i   , D = rot90(A)

   如果旋转180°:D = rot90(A, k)

矩阵乘方

   求矩阵A的x次幂:

   数学描述: 𝐹=𝐴𝑥  F=Ax

% matlab 矩阵的乘方
F = A^x

A必须为方阵

   x为整数

   x为非整数
−1⎯⎯⎯⎯⎯√𝑘    − 1  k  开方的多解:旋转

   得出一个解，乘以k-1次标量 𝛾=𝑒2𝜋𝑗/𝑘  γ=e 2 π j  /  k  

例2-9 矩阵的三次方根
𝐴=147258360  A= [    1   2   3     4   5   6     7   8   0    ]  

% 矩阵的3次方根
A = [1, 2, 3; 4, 5, 6; 7, 8, 0];
C = A^(1/3), e = norm(A-C^3)

另两个根

j1 = exp(sqrt(-1)*2*pi/3);
A1 = C * j1, A2 = C * j1^2,
norm(A-A1^3), norm(A-A2^3)

点运算

   矩阵对应元素的直接运算

   例如:

   C = A .* B, 𝑐𝑖𝑗=𝑎𝑖𝑗𝑏𝑖𝑗   c  i j  = a  i j   b  i j  

   B = A .^ A, 𝑏𝑖𝑗=𝑎𝑎𝑖𝑗𝑖𝑗   b  i j  = a  i j    a  i j    
𝐴=147258360  A= [    1   2   3     4   5   6     7   8   0    ]  

   B = A .^ A = 114477225588336600=1256823543431251677721627466561   [    11   22   33     44   55   66     77   88   00    ] = [    1   4   27     256   3125   46656     823543   16777216   1    ]  

A = [1, 2, 3; 4, 5, 6; 7, 8, 0]; B = A .^ A

绘图:𝑦=𝑓(𝑥)=𝑥2,𝑦𝑖=𝑥2𝑖,𝑦=𝑥.2  y=f(x)=x2, y i = x i 2 ,y=x.2

x = [-2, -1, 0, 1, 2], y = x .^2; plot(x, y)

02.04 矩阵的其它运算

逻辑运算 – 相应元素间的运算

   与运算: A & B

   或运算: A | B

   非运算: B = ~A

   异或运算: xor(A, B)

   比较关系符:>, >=, <, <=, ==, ~=, find(), all(), any()

例2-10 比较运算实例

%比较运算实例
A = [1, 2, 3; 4, 5, 6; 7, 8, 0];
i = find(A >= 5)'
[i,j] = find(A >= 5)
a1 = all(A >= 5), a2 = any(A >= 5)

解析结果的化简与变换

   函数simplify()用于数学公式的化简
𝑠1   s 1   = simplify(s)

   其它常用化简函数

   numden(), collect(), expand(), factor()

例2-11 多项式处理

   化简多项式
𝑃(𝑠)=(𝑠+3)2(𝑠2+3𝑠+2)(𝑠3+12𝑠2+48𝑠+64)  P(s)=(s+3)2(s2+3s+2)(s3+12s2+48s+64)

   用不同函数求解

% 多项式化简
syms s;
P = (s+3)^2 * (s^2 + 3*s + 2) * (s^3 + 12*s^2 + 48*s + 64)
P1 = simplify(P)

% 多项式因式分解
P3 = factor(P), P4 = prod(P3)

变量替换
𝑓1=𝑠𝑢𝑏𝑠(𝑓,𝑥1,𝑥∗1)   f 1 =subs(f, x 1 , x 1 ∗ )
𝑓1=𝑠𝑢𝑏𝑠(𝑓,{𝑥1,𝑥2,⋯,𝑥𝑛},{𝑥∗1,𝑥∗2,⋯,𝑥∗𝑛})   f 1 =subs(f,{ x 1 , x 2 ,⋯, x n },{ x 1 ∗ , x 2 ∗ ,⋯, x n ∗ })
转换成LaTex表示 – 需要LaTex环境解释
𝑓1=𝑙𝑎𝑡𝑒𝑥(𝑓)   f 1 =latex(f)

例2-12 双线性变换

   试用 s = (z-1)/(z+1)对 P(s) 进行双线性变换

   多项式表达式:𝑃(𝑠)=(𝑠+3)2(𝑠2+3𝑠+2)(𝑠3+12𝑠2+48𝑠+64)  P(s)=(s+3)2(s2+3s+2)(s3+12s2+48s+64)

% matlab 变量替换代码
syms s z;
P = (s+3)^2*(s^2+3*s+2)*...
(s^3+12*s^2+48*s+64);
P1 = simplify(subs(P, s, (z-1)/(z+1)))
latex(P1)
% '\frac{8\, z\, {\left(2\, z + 1\right)}^2\, \left(3\, z + 1\right)\, {\left(5\, z + 3\right)}^3}{{\left(z + 1\right)}^7}'

基本数论运算

例 2-13 不同的取整函数

   运用各种函数，对下面的数据进行取整运算

   -0.2765, 0.5772, 1.4597, 2.1091, 1.191, -1.6187

% 各种取整
A = [-0.2765, 0.5772, 1.4597, 2.1091, 1.191, -1.6187];
v1 = floor(A), v2 = ceil(A),
v3 = round(A), v4 = fix(A)

例2-14 Hilbert矩阵

   假设3x3的Hilbert矩阵可以右hilb()定义,试对其进行有理数变换。
ℎ𝑖,𝑗=1𝑖+𝑗−1   h  i , j  =  1  i + j − 1    

% Hilbert矩阵有理数变换
A = hilb(3); [n, d] = rat(A)

例2-15 最大公约数和最小公倍数

   给定两个整数1856120和1483720。

   求其最大公约数和最小公倍数。

   求其所得出的最小公倍数的质因数分解。

m = sym(1856120); n = sym(1483720);
gcd(m, n), lcm(m, n),
factor(lcm(n, m))

例2-16 找出1000以内全部质数

% 找出1000以内所有质数
A = 1:1000; B = A(isprime(A))

例2-17 全排列计算

   全排列函数: perms()

   给出1到5的全排列.

% matlab 全排列
P = perms(1:5), size(P)
P = perms('abcde'), size(P)

02.05 matlab 流程结构

循环语句
for循环 : for i = v, loop   programs body, end

while循环:

   while (condition)

   loop structure body,

   end

例2-18 求解 ∑𝑖=1100𝑖   ∑  i = 1   100  i

% for loop
s1 = 0; for i = 1: 100, s1 = s1 + i; end, s1
% while loop
s2 = 0; i = 1;
while(i <= 100), s2 = s2 + i; i = i + 1; end, s2
% 简单sum函数
sum(1:100)

例2-19 while循环

   用循环求解最小的m,使下式成立 ∑𝑖=1𝑚𝑖>10000   ∑  i = 1   m  i>10000

s = 0; m = 0;
while (s <= 10000), m = m + 1; 
s = s + m; end, s, m

例2-20 向量化编程

   求和 𝑆=∑𝑖=1100000(12𝑖+13𝑖)  S= ∑  i = 1   100000  (  1 2i  +  1 3i  )

% 循环求解
tic, s = 0;
for i = 1: 100000, s = s + 1/2^i + 1/3^i;
end; toc
% 时间已过 0.038848 秒。
% %%%%%%%%%%%
% 向量化编程
tic, i = 1:100000;
s = sum(1 ./2 .^i + 1./ 3 .^i); toc
% 时间已过 0.019701 秒。

条件语句

if (condition 1)
    statement group 1
elseif (condition 2)
    statement group 2
else
    statement group n + 1
end

例2-21 for循环+if

   用for循环求解最大的m，试下式成立:∑𝑖=1𝑚𝑖>10000   ∑  i = 1   m  i>10000

s = 0;
for i = 1: 10000, s = s + i;
    if s > 10000, break;
    end,
end, i, s

开关语句

switch switch expression
case expression 1, statements 1
case {expression2, expression 3, ..., expression m}, statements 2
    ...
otherwise, statements n
end

try catch

try,
    statement group1,
catch,
    statement group2,
end

02.06 函数编写

例2-22 M-脚本文件实现

   问题：∑𝑖=1𝑚𝑖>1000   ∑  i = 1   m  i>1000

   m-脚本以m-文件形式存取

% 以下内容保持到文件s10000.m
s = 0; m = 0;
while (s <= 10000),
    m = m + 1,
    s = s + m;
end,
s, m
% 直接执行s10000,即调用s10000.m脚本运行

% 函数基本结构
function [return vars] = funcname(input vars)
  comments led by %
  input and output variables check
  main body of the function

例2-23 m-函数实现

% 函数实例
function [ m, s ] = findsum( k )
%findsum summary
% findsum Detailed description
  s = 0; m = 0;
  while (s <= k)
    m = m + 1;
    s = s + m;
  end;
end
% call function
[m1, s1] = findsum(145323)

例2-24 Hilbert长方形矩阵

   编写一个函数生成 n x m Hilbert 矩阵
ℎ𝑖𝑗=1𝑖+𝑗−1   h  i j  =  1  i + j − 1    

function [ output_args ] = myhilb( n, m )
  if nargin == 1
    m = n;
  end;

  for i = 1:n, 
      for j = 1: m
          H(i, j) = 1/(i+j-1);
      end,
  end;
end
% call function
% H = myhilb(sym(4), 3)
% H = myhilb(3)

例2-25 阶乘的递归调用 n! = n(n-1)!

% n!
function k = my_fact(n)
  if nargin ~= 1
    error('Error: Only one input variable accepted');
  end;
  if abs(n - floor(n)) > eps || n < 0
    error('n should be a non-negative integer');
  end;
  if n > 1, k = n * my_fact(n-1);
  elseif  any([0, 1] == n)
    k =1;
  end;
end

阶乘的不同计算方法

   factorial(n)

   prod(1:n)

   gamma(1+n)

   gamma(1 + sym(n))

可变输入输出

   输入输出变量: varargin, varargout

   变量的提取 - 单元数组(cell)

   varargin{1}, varargin{2}, …, varargin{n}

例2-27 任意多输入变量

   conv()可计算两个多项式的积,试用varargin实现任意多个多项式的积

function a = convs(varargin), a = 1;
  for i = 1: nargin;
    a = conv(a, varargin{i});
  end;
end
% call convs
% P = [1 2 4 0 5]; Q = [1 2];
% F = [1 2 3]; D = convs(P, Q, F), E = conv(conv(P, Q), F)
% G = convs(P, Q, F, [1, 1], [1, 3], [1, 1])

匿名函数

   f = @(list of variables) function_contents

伪代码与代码保密

   伪代码:提高程序的执行速度，保密，把ascii的.m文件转换成二进制文件

   伪代码语句:

pcode mytest
pcode *.m
pcode mytest -inplace

一定要先保存好.m源文件,pcode没有对应恢复命令。没办法通过.p文件恢复成.m文件。

02.07 二维曲线的绘制

例2-28 函数曲线绘制与检验

   绘制函数: 𝑦=sin(tan𝑥)−tan(sin𝑥),𝑥∈[−𝜋,𝜋]  y=sin⁡(tan⁡x)−tan⁡(sin⁡x),x∈[−π,π]

x = [-pi: 0.05: pi];
y = sin(tan(x)) - tan(sin(x)); plot(x, y)

不同步距效果

x = [-pi: 0.0001: pi];
y = sin(tan(x)) - tan(sin(x)); plot(x, y)

例2-29 分段函数

   绘制饱和函数方程𝑦={1.1𝕤𝕚𝕘𝕟(𝑥),|𝑥|>1.1𝑥,|𝑥|≤1.1  y= {    1.1  s i g n  ( x ) ,  |  x  |  > 1.1     x ,    |  x  |  ≤ 1.1      

   互斥条件:

x = [-2: 0.02: 2];
y = 1.1 * sign(x) .* (abs(x) > 1.1) + ...
    x .* (abs(x) <= 1.1); plot(x, y)

更简单的命令–折线

plot([-2, -1.1, 1.1, 2], [-1.1, -1.1, 1.1, 1.1])

plot其它调用格式

   t仍为向量，而y为矩阵,
𝑦=𝑦11𝑦21⋮𝑦𝑛1𝑦12𝑦22⋮𝑦𝑛2⋯⋯⋱⋯𝑦1𝑛𝑦2𝑛⋮𝑦𝑛𝑛  y= [     y  11      y  12     ⋯    y  1 n        y  21      y  22     ⋯    y  2 n       ⋮   ⋮   ⋱   ⋮      y  n 1      y  n 2     ⋯    y  n n      ]  
plot(t, y)

   t和y为矩阵,且t和y的行和列数均相同

   假设有多对这样的向量或矩阵
(𝑡1,𝑦1),(𝑡2,𝑦2),⋯,(𝑡𝑚,𝑦𝑚),  ( t 1 , y 1 ),( t 2 , y 2 ),⋯,( t m , y m ),

   plot(𝑡1,𝑦1,𝑡2,𝑦2,⋯,𝑡𝑚,𝑦𝑚)  ( t 1 , y 1 , t 2 , y 2 ,⋯, t m , y m )

更一般的调用格式

   改变曲线性质

   h = plot(𝑡1,𝑦1   t 1 , y 1  , option 1, 𝑡2,𝑦2   t 2 , y 2  , option 2, ⋯,𝑡𝑚,𝑦𝑚  ⋯, t m , y m  , option m)

图形修饰与属性设置

   每个窗口、曲线、坐标轴都是一个对象

   对象的属性可以通过set()来设置

   可以通过函数get()来获取

   set(handle, ‘p_name1’, …, p_value1, ‘p_name2’, p_value2, …)

   v = get (object, ‘p_name’)

   图形对象的属性可以通过快捷菜单(鼠标右键)直接修改

多纵轴曲线的绘制

   例2-30 试绘制曲线 𝑦1=sin𝑥,𝑦2=0.01cos𝑥   y 1 =sin⁡x, y 2 =0.01cos⁡x

x = 0: 0.01: 2*pi;
y1 = sin(x); y2 = 0.01 * cos(x);
plot(x, y1, x, y2, '--')

x = 0: 0.01: 2*pi;
y1 = sin(x); y2 = 0.01 * cos(x);
plotyy(x, y1, x, y2)

三、四纵轴图形可以下载相应函数绘制

   plotyyy(), plot4y(), 从MathWorks File Exchange 下载，还可以使用 plotxx()函数

02.08 特殊二维图形

其它二维图形绘制语言

例2-31 极坐标曲线

   绘制极坐标函数:𝜌=5sin(4𝜃/3)𝜌=5sin(𝜃/3)  ρ=5sin⁡(4θ / 3)ρ=5sin⁡(θ / 3)

theta = 0: 0.01: 6*pi;
rho = 5*sin(4*theta/3);
polar(theta, rho)

rho = 5*sin(theta/3);
polar(theta, rho)

例2-32 不同的二维图形

t = 0: .2: 2*pi; y = sin(t);
subplot(2, 2, 1), stairs(t, y)
subplot(2, 2, 2), stem(t, y)
subplot(2, 2, 3), bar(t, y)
subplot(2, 2, 4), semilogx(t, y)

绘制区域设置

   分割不同区域: subplot(n, m, k)

隐函数绘制及应用

   隐函数: 𝑥2sin(𝑥+𝑦2)+𝑦2𝑒𝑥+𝑦+5cos(𝑥2+𝑦)=0  x2sin⁡(x+y2)+y2e x + y +5cos⁡(x2+y)=0

   隐函数绘图语句: ezplot(implicit function expression)

   用字符串表示隐函数

   默认区域是 [−2𝜋,2𝜋]  [−2π,2π]

   其它语法: ezplot(im_function, [𝑥𝑚,𝑥𝑀   x m , x M  ])

   ezplot(im_function, [𝑥𝑚,𝑥𝑀,𝑦𝑚,𝑦𝑀   x m , x M , y m , y M  ])

例2-33 隐函数曲线

   试绘制隐函数: 𝑥2sin(𝑥+𝑦2)+𝑦2𝑒𝑥+𝑦+5cos(𝑥2+𝑦)=0  x2sin⁡(x+y2)+y2e x + y +5cos⁡(x2+y)=0

% 默认范围 -2*pi, 2*pi
h = ezplot('x^2*sin(x+y^2) + y^2*exp(x+y)+5*cos(x^2+y)');
set(h, 'Color', 'b')

% 扩大范围
h = ezplot('x^2*sin(x+y^2) + y^2*exp(x+y) + 5*cos(x^2+y)', [-10, 10]);
set(h, 'Color', 'r');

图形修饰

   直接采用工具栏

   文字修饰

   特殊符号表 LaTex

   上、下标分别用 ^ 和 _ 表示

   a_2^2 + b_2^2 = c_2^2 ? 𝑎22+𝑏22=𝑐22   a 2 2 + b 2 2 = c 2 2  

数据文件的读取与存储

   可以采用save和load命令存储和读取数据

   save mydat A B C

   save /ascii mydat.dat A B C

   X = load(filename)

   matlab和excel交互数据

   X = xlsread(filename, range) ‘B6:C67’

   写文件: xlswrite()

例2-34 Excel文件读取

   已知Excel文件 census.xls 给出某省人口数

   第5-67行存储数据

   B列存储年份，C列存储人口数

   先读入matlab再绘图

X = xlsread('census.xls', 'B5:C67');
t = X(:1); p = X(:, 2);
plot(t, p)

更加单方法 – Copy & Paste

02.09 三维图形绘制

三维曲线绘制

   三维曲线绘制(空间质点的运动轨迹)

   plot3(x, y, z)

   plot3(𝑥1,𝑦1   x 1 , y 1  , option 1, 𝑥2,𝑦2   x 2 , y 2  , option 2, ⋯,𝑥𝑚,𝑦𝑚  ⋯, x m , y m   option m)

其它三维曲线绘制函数

   stem3, fill3, bar3等

   绘制三维轨迹图: comet3

例2-35 空间质点的位置

   试绘制参数方程:
𝑥(𝑡)=𝑡3sin(3𝑡)𝑒−𝑡𝑦(𝑡)=𝑡3cos(3𝑡)𝑒−𝑡𝑧=𝑡2   {    x ( t ) = t3 sin ⁡ ( 3 t ) e − t      y ( t ) = t3 cos ⁡ ( 3 t ) e − t      z = t2      

   质点的空间位置(随时间变化)

   其中, 𝑡∈[0,2𝜋]  t∈[0,2π]

t = 0: 0.01: 2*pi;
x = t .^ 3 .* sin(3*t) .* exp(-t);
y = t .^ 3 .* cos(3*t) .* exp(-t);
z = t .^ 2;
plot3(x, y, z), grid

其它曲线绘制

% 其它曲线绘制
stem3(x, y, z); hold on;
plot3(x, y, z); grid; hold off

% 动画效果
figure; comet3(x, y, z)

图形窗口的工具栏

   3D绘图和视角变换

   读取坐标值、局部放大

三维曲面绘制

   一般曲面绘制: z = f(x, y)

   [x, y] = meshgrid(𝑣1,𝑣2   v 1 , v 2  )

   z = …, for instance z = x .* y

   surf(x, y, z) or mesh(x, y, z)

   其它函数surfl(), surfc()

   等高线绘制

   contour(), contours()

例2-36 三维曲面

   给出二元函数如下, 绘制3D图形
𝑧=𝑓(𝑥,𝑦)=(𝑥2−2𝑥)𝑒−𝑥2−𝑦2−𝑥𝑦  z=f(x,y)=(x2−2x)e − x2 − y2 − x y  

[x, y] = meshgrid(-3: 0.1: 3, -2: 0.1: 2);
z = (x .^ 2 - 2*x) .* exp(-x .^ 2 -y .^ 2 - x .* y);
mesh(x, y, z)

% 表面图:
surf(x, y, z)

例2-37 试绘制出二元函数
𝑧=𝑓(𝑥,𝑦)=1(1−𝑥)2+𝑦2⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯√+1(1+𝑥)2+𝑦2⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯√  z=f(x,y)=  1  ( 1 − x )2 + y2   +  1  ( 1 + x )2 + y2    

   绘制3D图形

[x, y] = meshgrid(-2: .1: 2);
z = 1 ./ (sqrt((1-x).^2 + y .^2)) + 1 ./ (sqrt((1+x).^2 + y .^2));
surf(x, y, z), shading flat

例2-38 分段函数处理

   绘制如下分段函数的三维曲面
𝑝(𝑥1,𝑦1)=0.5457𝑒𝑥𝑝(−0.75𝑥22−3.75𝑥21−1.5𝑥1),𝑥1+𝑥2>10.7575𝑒𝑥𝑝(−𝑥22−6𝑥21),−1<𝑥1+𝑥2≤10.5457𝑒𝑥𝑝(−0.75𝑥22−3.75𝑥21+1.5𝑥1),𝑥1+𝑥2≤−1  p( x 1 , y 1 )= {    0.5457 e x p ( − 0.75  x 2 2  − 3.75  x 1 2  − 1.5  x 1  ) ,  x 1  +  x 2  > 1     0.7575 e x p ( −  x 2 2  − 6  x 1 2  ) , − 1 <  x 1  +  x 2  ≤ 1     0.5457 e x p ( − 0.75  x 2 2  − 3.75  x 1 2  + 1.5  x 1  ) ,  x 1  +  x 2  ≤ − 1      

   分段函数求值: 互斥的不等式、点运算

[x1, x2] = meshgrid(-1.5: .1: 1.5, -2: .1: 2);
p = 0.5457*exp(-0.75*x2.^2 - 3.75*x1.^2-1.5*x1).*(x1+x2>1) + ...
0.7575*exp(-x2.^2-6*x1.^2).*((x1+x2>-1)&(x1+x2<=1)) +...
0.5457*exp(-0.75*x2.^2-3.75*x1.^2+1.5*x1).*(x1+x2<=-1);
surf(x1, x2, p), xlim([-1.5 1.5]);

参数方程的表面图

   三维函数参数方程
𝑥=𝑓𝑥(𝑢,𝑣),𝑦=𝑓𝑦(𝑢,𝑣),𝑧=𝑓𝑧(𝑢,𝑣)  x= f x (u,v),y= f y (u,v),z= f z (u,v)

   参数区间: 𝑢𝑥≤𝑢≤𝑢𝑀𝑣𝑚≤𝑣≤𝑣𝑀   u x ≤u≤ u M  v m ≤v≤ v M  

   ezsurf(𝑓𝑥,𝑓𝑦,𝑓𝑧,[𝑢𝑚,𝑢𝑀,𝑣𝑚,𝑣𝑀]   f x , f y , f z ,[ u m , u M , v m , v M ] )

   默认区间:(−2𝜋,2𝜋)  (−2π,2π)

例2-40 Mobius带的绘制

   数学模型:
𝑥=cos𝑢+𝑣cos𝑢cos𝑢/2𝑦=sin𝑢+𝑣sin𝑢cos𝑢/2𝑧=𝑣sin𝑢/2   {    x = cos ⁡ u + v cos ⁡ u cos ⁡ u  /  2     y = sin ⁡ u + v sin ⁡ u cos ⁡ u  /  2     z = v sin ⁡ u  /  2      

% Mobius带
syms u v; x= cos(u) + v*cos(u)*cos(u/2);
y = sin(u) + v*sin(u)*cos(u/2);
z = v*sin(u/2);
ezsurf(x, y, z, [0, 2*pi, -0.5, 0.5])

球面绘制

   matlab函数生成数据[x, y, z] = shpere(n)

   生成矩阵(n+1) x (n+1)

例2-41 球面绘制

   绘制两个球, 一个圆心在原点, 半径为1，另一个圆心(0.9, -0.8, 0.6),半径为0.3

[x, y, z] = sphere(50); surf(x, y, z), hold on
x1 = 0.3*x + 0.9; y1 = 0.3*y - 0.8; z1 = 0.3*z + 0.6;
surf(x1, y1, z1), hold off

柱面绘制

   某曲线r沿z轴旋转一周得出柱面

   [x, y, z] = cylinder(r, n)

例2-42 柱面曲线方程
𝑟(𝑧)=𝑒−𝑧2/2sin𝑧,𝑧∈(−1,3)  r(z)=e − z2  /  2 sin⁡z,z∈(−1,3)

z0 = -1: 0.1: 3; r = exp(-z0.^2/2).*sin(z0);
[x, y, z] = cylinder(r); z = -1 + 4*z; surf(x, y, z)

02.10 特殊三维图形

等高线绘制

   各种等高线绘制方法

   直接绘制: contour(x, y, z, n)

   带有标注的等高线:

   [C, h] = countour(x, y, z, n)

   clabel(C, h)

   三维等高线countour3(x, y, z, n) coutourf(x, y, z, n)

例2-43 等高线

   例2-38的分段函数，绘制各种等高线

   生成数据，再绘制图形
𝑝(𝑥1,𝑦1)=0.5457𝑒𝑥𝑝(−0.75𝑥22−3.75𝑥21−1.5𝑥1),𝑥1+𝑥2>10.7575𝑒𝑥𝑝(−𝑥22−6𝑥21),−1<𝑥1+𝑥2≤10.5457𝑒𝑥𝑝(−0.75𝑥22−3.75𝑥21+1.5𝑥1),𝑥1+𝑥2≤−1  p( x 1 , y 1 )= {    0.5457 e x p ( − 0.75  x 2 2  − 3.75  x 1 2  − 1.5  x 1  ) ,  x 1  +  x 2  > 1     0.7575 e x p ( −  x 2 2  − 6  x 1 2  ) , − 1 <  x 1  +  x 2  ≤ 1     0.5457 e x p ( − 0.75  x 2 2  − 3.75  x 1 2  + 1.5  x 1  ) ,  x 1  +  x 2  ≤ − 1      

[x1, x2] = meshgrid(-1.5: .1: 1.5, -2: .1: 2);
p = 0.5457*exp(-0.75*x2.^2 - 3.75*x1.^2-1.5*x1).*(x1+x2>1) + ...
0.7575*exp(-x2.^2-6*x1.^2).*((x1+x2>-1)&(x1+x2<=1)) +...
0.5457*exp(-0.75*x2.^2-3.75*x1.^2+1.5*x1).*(x1+x2<=-1);
[C, h] = contour(x1, x2, p); clabel(C, h)

contourf(x1, x2, p);
figure; contour3(x1, x2, p, 30)

三维隐函数图绘制

   三元隐函数: 𝑓(𝑥,𝑦,𝑧)=0  f(x,y,z)=0

   下载 ezimplot3()函数

   ezimplot3(fun, [𝑥𝑚,𝑥𝑀,𝑦𝑚,𝑦𝑀,𝑧𝑚,𝑧𝑀   x m , x M , y m , y M , z m , z M  ])

   默认范围:[−2𝜋,2𝜋]  [−2π,2π]

   fun可以为[隐函数字符串、匿名函数、M-函数]

例2-44 三元隐函数绘制

   已知三元隐函数
𝑥(𝑥,𝑦,𝑧)=𝑥sin(𝑦+𝑧2)+𝑦2cos(𝑥+𝑧)+𝑧𝑥cos(𝑧+𝑦2)=0  x(x,y,z)=xsin⁡(y+z2)+y2cos⁡(x+z)+zxcos⁡(z+y2)=0

f = 'x*sin(y+z^2) + y^2*cos(x+z)+z*x*cos(z+y^2)';
f = @(x, y, z)x*sin(y+z^2) + y^2*cos(x+z)+z*x*cos(z+y^2);
ezimplot3(f, [-1, 1]);
f1 = 'x^2+y^2+z^2-1'; ezimplot3(f1, [-1 1]);

三维图形视角设置

   两种方法改变图形视角

   直接采用工具栏

   命令语句:view()
𝑣𝑖𝑒𝑤(𝛼,𝛽)  view(α,β)

   读角度: [𝛼,𝛽]=𝑣𝑖𝑒𝑤(3)  [α,β]=view(3)
𝛼  α 定义为方位角, 𝛽  β 定义为仰角

例 2-45 三视图的设置

   函数: 𝑧=𝑓(𝑥,𝑦)=(𝑥2−2𝑥)𝑒−𝑥2−𝑦2−𝑥𝑦  z=f(x,y)=(x2−2x)e − x2 − y2 − x y  

   默认视角的提取: [a, b] = view(3);

   三视图绘制:

[x, y] = meshgrid(-3: 0.1: 3, -2: 0.1: 2);
z = (x .^2 - 2*x) .* exp(-x.^2-y.^2-x.*y);
subplot(224), surf(x,y,z),
subplot(221), surf(x,y,z), view(0, 90);
subplot(222), surf(x,y,z), view(90, 0);
subplot(223), surf(x,y,z), view(0, 0);

三维曲面旋转

   三维曲面旋转函数: rotate(h, v, 𝛼  α )

   其中: 三维曲线句柄，旋转基线，旋转角度

   如: 绕x轴正向旋转, v = [1, 0, 0]

例2-46 三维图形旋转

   分段函数曲面旋转绘制

[x, y] = meshgrid(-1: 0.04: 1, -2: 0.04: 2);
z = 0.5457*exp(-0.75*y.^2 - 3.75*x.^2-1.5*x).*(x+y>1) + ...
0.7575*exp(-y.^2-6*x.^2).*((x+y>-1)&(x+y<=1)) +...
0.5457*exp(-0.75*y.^2-3.75*x.^2+1.5*x).*(x+y<=-1);
h = surf(x, y, z);
% 绕x轴正方向旋转,基线 v = [1, 0, 0]
rot_ax = [1, 0, 0]; rotate(h, rot_ax, 15)
% 绕基线 v = [1, 1, 1] 旋转
figure;
h = surf(x, y, z); rot_ax = [1, 1, 1];
rotate(h, rot_ax, 15)

旋转的动画演示

   绕x轴正向旋转360度的动画

   每步(0.01秒)旋转1度

   循环结构旋转

figure;h = surf(x,y,z);
r_ax = [1, 0, 0]; axis tight
for i = 0:360
  rotate(h, r_ax, 1); pause(0.01);
end;

02.11 四维图形绘制

四维图形

   时间驱动的三维图形的动画

   三维图形的切面 – 体视化(volume visualization)

   体视化举例: CT、固体内部温度、密度，需要用切面观察

   流体的流速、液体的浓度分布

   matlab函数:slice(x, y, z, V, 𝑥1,𝑦1,𝑧1   x 1 , y 1 , z 1  )

例2-47 三维实体的切片

   三元函数
𝑉(𝑥,𝑦,𝑧)=𝑥𝑥+𝑦(𝑥+𝑦)/2+𝑧(𝑥+𝑦+𝑧)/3⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯√  V(x,y,z)= xx + y ( x + y )  /  2  + z ( x + y + z )  /  3  

   普通切面

   生成网格数据

   计算网格各点的函数值(点运算)

   绘图

[x, y, z] = meshgrid(0: 0.1: 2);
V = sqrt(x .^ x + y .^((x+y)/2) + z .^((x+y+z)/3) );
slice(x, y, z, V, [1 2], [1 2], [0 1]);

特殊切面

   先构造 z = 1 平面

   将该平面沿x轴旋转45度

   用该切面给原体视化图形切片

[x0, y0] = meshgrid(0: 0.1: 2);
z0 = ones(size(x0));
h = surf(x0, y0, z0);
rotate(h, [1,0,0], 45);
x1 = get(h, 'XData');
y1 = get(h, 'YData');
z1 = get(h, 'ZData');
slice(x, y, z, V, x1, y1, z1),
hold on, slice(x, y, z, V, 2, 2, 0), hold off

shading interp

编写体视化程序界面

   程序界面 vol_visual4d()

   数学函数: V = f(x, y, z)

   vol_visual.fig, vol_visual4d.m

   调用方法

   生成体视化数据：x, y, z, V

   调用该函数: vol_visual4d(x, y, z, V)

   通过界面任意设置切面

例2-48 利用界面体视化

   三元函数
𝑉(𝑥,𝑦,𝑧)=𝑥𝑥+𝑦(𝑥+𝑦)/2+𝑧(𝑥+𝑦+𝑧)/3⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯√  V(x,y,z)= xx + y ( x + y )  /  2  + z ( x + y + z )  /  3  

   体视化切面绘图

   通过滚动杆调整切面位置

   通过捡取框选择是否显示某轴切面

[x,y,z] = meshgrid(0: 0.1: 2);
V = sqrt(x .^ x + y .^ ((x+y)/2) + z .^ ((x+y+z)/3) );
vol_visual4d(x, y, z, V);


免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删