1.3 变量、数据与函数
1.3.1 变量
像任何其它计算机语言一样,MATLAB也有变量名规则。变量名必须是不含有空格的单个词。变量命名规则如下:
变量名区分字母大小写,如ltems, items, itEms及ITEMS都是不同的变量。
变量名最多不超过19个字符。第19个字符之后的字符将被忽略,如
howaboutthisvariablename只能表示为howaboutthisvariabl。
变量名必须以字母打头,之后可以是任意字母、数字或下划线,如x51488,a_b_c_d_e。标点符号在MATLAB中具有特殊含义,所以变量名中不允许使用标点符号。
除了这些命名规则,MATLAB还有几个特殊变量,见表1-3。
表1-3 MATLAB特殊变量表
特殊变量
取值
ans
pi
eps
flops
Inf
NnN
用于结果的缺省变量名
圆周率
计算机是最小数,当和1相加就产生一个比1大的数
浮点运算数
无穷大,如1/0
不定量,如0/0
i(和)j
nargin
nargout
realmin
realmax
所有函数的输入变量数目
所有函数的输出变量数目
最小可用正实数
最大可用正实数
表1-3中的特殊变量在启动MATLAB之后,自动赋予表中取值。如果定义了相同名字的变量,原始特殊取值将会丢失,直到清除所有变量或重新启动MATLAB。一般来讲,应当尽量避免重新定义特殊变量。
1.3.2 注释和标点
百分号后所有的文字为注释。注释语句不能执行。如键入
x=4.5 %在y=5 时,%注释该x取值是在y=5时有效
显示
多条命令可以放在同一行,中间用逗号或分号隔开。逗号要求显示结果,分号禁止显示结果。如键入
x=4.5; y=5, f=1.9
显示
在PC机上运行时,可以随时按下CTRL+C键中断MATLAB的运行。
1.3.3 复数表示
MATLAB对复数不需要特殊处理,用i、j和sqrt(-x)(X是任意整数、实数)表示。复数的数学运算可以写成与实数同样的形式。如键入
a=1-2i, b=sqrt(-2), c=a+b
显示
MATLAB还可以用real, imag, abs, angle命令来表示一个复数的实部、虚部、幅值和相角。如
1.3.4 数学函数
MATLAB所支持的常用数学函数见到表1-4。注意,MATLAB只对孤度操作。
表1-4 常用函数
命令
说明
abs(x)
acos(x)
acosh(x)
angle(x)
asin(x)
asinh(x)
atan(x)
atan2(x, y)
atanh(x)
ceil(x)
conj(x)
cos(x)
cosh(x)
exp(x)
fix(x)
floor(x)
gcd(x, y)
imag(x)
lcm(x, y)
log(x)
log10(x)
real(x)
rem(x, y)
round(x)
sign(x)
sin(x)
sinh(x)
sqrt(x)
tan(x)
tanh(x)
绝对值或复数的辐值
反余弦
反双曲余弦
四象限内取复数相角
反正弦
反双曲正弦
反正切
四象限内反正切
反双曲正切
对+方向取整数
复数共轭
余弦
双曲余弦
指数函数ex
对零方向取整数
对-方向取整数
整数x和y的最大公约数
复数虚部
整数x和y的最小公倍数
自然对数
常用对数
复数实部
除后余数;rem(x, y)给出x/y的余数
四舍五入到最接近的整数
符号函数;返回自变量的符号,例如sign(1, 5)=1, sign(-2.4)=-1, sigh(0)=0
正弦
双曲正弦
平方根
正切
双曲正切
1.4 数据的输入与输出
1.4.1 save与load命令
MATLAB可以通过计算机文件来保存或加载数据。File菜单的Save Workspace As…菜单命令打开一个标准的文件对象框来保存所有当前变量。保存变量并不会将其从MATLAB工作空间中删除。类似地,File菜单的Load Workspace As…菜单命令打开一个对话框,在工作空间中加载以前保存的变量。加载MATLAB工作空间中已有的同名变量将会把变量值改为从文件中加载的值。
如果File菜单方法不可行或不能满足用户需要,MATLAB还有save和load两个命令,可以提供更大的灵活性。Save命令允许自己选择文件格式保存一个或多个变量,而load命令将加载自己选择文件格式保存的文件,取出所需的变量。对于大多数用户,MATLAB命令load和save为装载和存储数据提供了足够的工具。利用以扩展名为.mat结尾的文件名,load和save假定数据是以与软件平台无关的二进制格式保存,或者用称之为flat的简单的ASCⅡ文件格式保存。例如:
(1)指令save以二进制格式将所有变量存入到matlab. mat文件中。
(2)指令save data以二进制格式将所有变量存入到date. mat文件中。
(3)指令save data x y z以二进制格式将变量x, y, z存入到data, mat文件中。
指令load data将data. Mat中的所有变量加载到工作空间中,供用户使用。
1.4.2 低级文件输入输出
当flat ASCⅡ或. mat这两种格式还不够时,MATLAB提供了基于C语言的低级文件I/O命令,可以读出你所知道的任意文件格式。表1-5列出了MATLAB中基本的低级文件I/O命令。
表1-5 MATLAB低级文件I/O命令
命令
说明
fclose
feof
ferror
fgetl
fgets
fopen
fprintf
fread
frewind
fscanf
fseek
ftell
fwrite
关闭文件
测试文件结束
查询文件I/O的错误状态
读文件的行,忽略回行符
读文件的行,包括回行符
打开文件
按照格式要求把数据写到文件或屏幕上
从文件中读二进制数据
返回到文件开始
按照格式要求从文件中读数据
设置文件位置指示符
获取文件位置指示符
把二进制数据写到文件里
表1-5的命令都可以通过help功能得到应用的详细解释。这里只对部分常用的命令进行解释。比如打开文件函数fopen( )的语句格式为
文件句柄=fopen(文件名,文件类型)
其中文件句柄为一个整数,供从该文件中读取数据时使用。文件名应该是用单引号括起来的字符串,而文件类型可以由一个字符串来描述,其意义和C语言的几乎一致。如它可以采用‘r’来表示一个只读型的文件。而‘a’表示一个可添加的文件。例如,如果想打开一个名为myfile, xdy的文件,但不想改变其中的内容,只想从中读出一些数据,则可以把它按一个只读型文件打开,这样就要使用下面的命令:
myf=fopen(‘myfile. xdy’,‘r’)
如果该文件存在,则返回一个整数myf句柄,可以调用fread()或fscanf()等命令从中读取数据。其中,fread()命令从该文件中按照二进制格式读取数据,fscanf命令按照用户给定的格式读取数据。如
S=fscanf(myf, ‘%s’) %按文字型读取变量,赋值于S
A=fscanf(myf, ‘%5d’) %读5个十进制数
fprintf命令的格式举例如下:
fprintf (fid, ‘%6.2f %12.8e\n’, y);
其中,fid是被写数据文件的句柄,单引号‘’中间的部分是写数据的格式要求,与C语言的数据格式相同。每一种不同格式要用%分开,小数点前、后的数字表示数据的整数和小数部分的位数。最后的变量是被写的数据变量名。
下面给出一个例题说明用低级文件读出数据的过程。设已有数据文件Test. tst,内容如下:
a_d_l
3.0
0.6
-2.00 .05372 .00707 -.02525
-1.00 .09308 .04892 .01294
运行文件inout.m内容如下:
% input and output a data file
fl=fopen(‘test. txt’,‘r’) %打开test. m文件,作为只读文件,句柄fl
p=fscanf(fl,‘%c’) %读取全部数据,按照数据原有格式,赋值于p
f2=fopen(‘name. m’,‘w’) %打开name.m文件,作为写文件,句柄f2
fprintf(f2,‘%s15\n, %9.5f\n, %9.5f\n, %4(9.5f)\n’,p) %写p,数据格式:
%第一行为文字型,15位;第2,3行为实数型,5位小数;第4行写4个实数型数据
fclose(‘all’) %关闭所有的文件
运行该文件时,键入
inout
屏幕显示结果
结果文件为name. m,内容为
a_d_l
3.0
0.6
-2.00 .05372 .00707 -.02525
-1.00 .09308 .04892 .01294
得到的结果文件内容与输入文件内容相同。
注意,fsconf命令读入的数据按文字型变量存储,里面的数据不能用来进行加法、乘法等运算。需要运算的数据可以用M文件形式输入,参看有关M文件的章节。
1.5 数组与矩阵运算
由于数据组可以定义为只有一行(或一列),因此,所有有关矩阵的运算都可以用于数组运算,除非有区别,今后不再区分数组与矩阵。
1.5.1 矩阵表示与块操作
1.矩阵表达式
数组和矩阵是MATLAB的基础。MATLAB的数组与矩阵用[]表示,程序可以自行解读矩阵的行、列标志和元素。定义矩阵的原则是:矩阵元素间用空格或逗号隔开,行用分号隔开。例如,键入
a=[1 2 3;4 5 6],b=[7 8 9]
显示
2.矩阵转置
矩阵的转置用A‘表示,例如,键入
aa=a‘ % 求矩阵转置
显示结果
若矩阵为复数矩阵,求转置时首先对矩阵元素进行转置,然后再逐项求取其共轭数值,这种转置方式又称为Hermit转置。例如,已知
x=
5.000+1.000i -2.000+1.0000i
4.0000 +3.000i
键入
xx=x’ % 求Hermit转置
结果为
3.矩阵大小
MATLAB具有矩阵大小的查询功能。查询矩阵的大小可以用表1-6中的命令。
表1-6 矩阵大小查询
命令
说明
Whos
size(A)
length(A)
find(A)
显示工作空间中存在的变量及其大小
返回A矩阵的行数和列数
返回A矩阵的最大行、列数
给出特殊要求的矩阵元素的行、列标记
·size命令
该函数的调用格式为[n,m]=size(A)。其中A为要测试的矩阵名,而返回的两个参数n和m分别为A矩阵的行数和列数。例如键入
[n,m]=size(A) %查询A的行、列数
可得
·length命令
如果要测试的变量是一个数组而不是矩阵时,仍可以由size()函数来得出大小。更简洁地,可以使用length()函数来求出。该函数的调用格式为n= length(A),其中A为要测试的数组名,返回的n为A数组的元素个数。如果A为矩阵,则将返回A的行、列数的最大值,即该函数等效于max(size(A))。例如键入
n=length(A) %查询A的最大维数
可得
·find命令
用find命令,MATLAB可以进行特殊要求的矩阵元素定位。如键入
[i,j]=find(A>3) % 指出矩阵元素中大于3的元素的行、列位置
显示结果
该结果表明,在A矩阵中,第2行的第1,2,3列元素均满足条件要求。又如,定义x1矩阵:
x1
2 -1 0 1 2
k=find(abs(x1)>1) %找出绝对值大于1的元素
k=
1 5
结果表明,x1矩阵中,第1,5个元素满足要求。
4.矩阵的块操作
MATLAB中提供了很多简便、智能的方式,可以对矩阵进行元素更改、插入子块、提取子块、重排子块、扩大维数等操作。这里,重要的是冒号“:”的应用。在MATLAB中,冒号“:”表示“全部”。如已知
A=
1 2 3
4 5 6
b=
7 8 9
如键入
a(1,:)=b %将a的第1行中所有元素用b的元素替代
如键入
a(:,:)=1 %矩阵所有元素设为1
如键入
a(2,3)=10 %第2行第3列元素等于10
MATLAB已定义的矩阵的维数可以扩大,但不能缩小,除非利用clear命令删除该矩阵。如果输入的同名矩阵的维数小于原矩阵维数,MATLAB认为是原矩阵修改了部分元素或子块。增加矩阵的维数时,可以只给出非零MATLAB自动将未定义元素改为0。如键入
a(5,5)=2 %定义矩阵A的第5行,第5列元素
结果显示
结果表示A矩阵扩展为5×5维矩阵。
另外,MATLAB可以进行3维数组操作,用一系列矩阵表示,所有矩阵维数必须相等。
如键入:
p(1,:,:)=[1 2; 3 4], p(2,:,:)=[5 6; 7 8]
结果显示了3维数组P的两个子矩阵的第1列、第2列元素。在多维数组的插值运算中,这种表示方式是有用的。
5.矩阵的翻转操作
MATLAB提供了几种函数可以进行矩阵的翻转操作。
·矩阵上下翻转
flipud命令可以将n行矩阵A作上下翻转,将A的行按照n…1的顺序重新排列,A的列保持不变。如有
a=
1 2 3
4 5 6
7 8 9
10 11 12
键入
b=flipud(a) %进行上下翻转
则有
·矩阵左右翻转
fliplr命令可以将m列矩阵A作左右翻转,将A的列按照m…1的顺序重新排列,A的列保持不变。如对上述a矩阵,键入
c=fliplr(a) %进行左右翻转
则有
·矩阵逆时针90°旋转
rot90命令将矩阵逆时针转90°,如对于上述A阵,键入
rot90(a) % 逆时针旋转90°
则有
另外,还有矩阵的对角化等操作。MATLAB命令见表1 - 7
表1 – 7 矩阵操作函数
命 令
说 明
flipud(A)
fliplr(A)
rot90(A)
diag(A)
diag(V)
tril(A)
triu(A)
矩阵作上下翻转
矩阵作左右翻转
矩阵逆时针翻转90°
提取矩阵A的对角元素,返回列向量
以列向量V作对角元素创建对角矩阵
提取A的下三角矩阵
提取A的上三角矩阵
1.5.2 矩阵运算
矩阵运算包括矩阵与标量、矩阵与矩阵的运算;矩阵函数和稀疏矩阵应用等。下面分别加以说明。
矩阵与标量的运算
运算包括+、—、×、÷和乘方等运算。矩阵和标量运算完成矩阵的每个元素对该标量的运算。如已知
a=
2 3
5 6
则
MATLAB用符号“^”表示乘方。求矩阵乘方时要求矩阵为方矩阵。已知矩阵
b=
4
5
若键入
b^2 % 其平方为b×b
则有
若键入
b^(-1) % 实际是求b的逆矩阵
则有
b^(0.2) %实际是将b矩阵开5次方
p=b^(0.2)
则有
矩阵与矩阵的运算
·矩阵加减法运算
矩阵A和B维数完全相同时,可以进行矩阵 加减法运算。它会自动地使得A和B矩阵的相应元素相加减。如果A和B的维数不相等,则MATLAB将自动地给出错误信息,提示两个矩阵的维数不相等。如已知
a=[1 2 3;4 5 6],b=[7;8;9],c=[10;11;12]
键入
a+b
则出现
由于A,B的维数不等,程序给出了错误信息。而如果键入
b+c
则有
·矩阵乘法运算
两个矩阵A,B的维数相容时(A的列数等于B的行数),可以进行C=A×B的运算。如对于上述定义的a,b阵,若键入
cc=a*b
则有
另外,MATLAB可以进行kronecker乘法运算。指令形式为C=kron(A,B),表示An×m和Bp×q矩阵的C=AB运算。结果为增广矩阵Cnp×mq,它表示A矩阵的每个元素依次与B矩阵的所有元素相乘,组成矩阵子块,n×m个子块共同组成新的矩阵C。如若
a=
1 2 3
4 5 6
b=
7
8
9
键入
cl=kron(a,b) % 求A,B的kronecker乘积
则有
·矩阵除法运算
矩阵的除法运算包括左除和右除两种运算。其中,
左除:A\B=A-1B,A为方矩阵;
右除:A/B=AB-1,B为方矩阵。
如已知
a=[1 2; 3 4],b=[1 3 5; 2 4 6],键入
cl=a\b %左除,结果为a-1×b
则有
若已知
c=
1 1 3
1 2 3
4 5 6
键入
b/c %右除,结果为b×c-1
则有
矩阵的除法运算实际是求AX=B的解的过程。当A为非奇异矩阵时,结果是最小二乘解,即矩阵除法可找到使||AX-B||误差绝对值最小的X。
·矩阵的点运算
MATLAB中定义了一种特殊的运算,即所谓的点运算。两个矩阵之间的点运算是该矩阵对应元素的直接运算,例如C=A.×B表示A和B矩阵的相应元素之间直接进行乘法运算,然后将结果赋给C矩阵。注意,点乘积运算要求A和B矩阵的维数相同。这种点乘积又称为Hadamard乘积。可以看出,这种运算和普通乘法运算是不同的,例如,已知A,B矩阵
a=
1 2
3 4
b=
2 2
1 2
若键入
c=a*b
则有
若键入
cc=a. *b
则有
可以看出,这两种乘积结果是不同的。前者是普通矩阵乘积,而后者是两个矩阵对应元素之间的乘积,形成了新的矩阵[aij * bij]。点运算在MATLAB中起着很重要的作用,例如,如果x是一个向量,则求取函数x的模值时不能直接写成x*x,而必须写成x. * x。在进行矩阵的点运算时,要求运算的两个矩阵的维数一致,其实一些特殊的矩阵函数,如sin()也是由点运算的形式来进行的,因为它要对矩阵的每个元素求取正弦值。
矩阵点运算不仅可以用于点乘积运算,还可以用于其它运算。比如对前面给出的α矩阵作a. * a运算,则将得出下面的结果:
键入
a. * a %结果形成新矩阵[aij2]
则有
·矩阵求幂
矩阵求幂的运算包括矩阵与常数和矩阵与矩阵的幂运算,用点运算的形式表示。具体解释如下:
a. ^3=[aij2],α矩阵的3次方——A矩阵的每个元素的3次方形成的新矩阵;
3. ^a=[3aij],3的a次方——新矩阵的每个矩阵元素都是以3为底,以A矩阵的对应元素为幂指数,形成的新矩阵。
a. ^b=[aij bij],a的b次方——新矩阵的每个元素都以A的元素为底,以B的对应元素为幂指数。例如有
a=
1 2
3 4
键入
a. ^3
则有
键入
3. ^a
则有
若有
b=
2 1
3 2
键入
a. ^b
则有
1.5.3 矩阵函数
MATLAB定义了一些特殊矩阵,不必一一赋值定义。特殊矩阵定义见表1-8。
表1-8 特殊矩阵
命令
说明
A=[ ]
A=eye(n)
A=ones(n,m)
A=rand(n,m)
A=randn(n,m)
A=zeros(n,m)
空矩阵
n维单位矩阵
全部元素都为1的矩阵
元素服从0和1之间均匀分布的随面矩阵
元素服从零均值单位方差正态分布的随机矩阵
全部元素都为0的矩阵
MATLAB还提供了很多用于求解线性代数数值问题的矩阵函数。表1-9给出了大部分矩阵函数的简短描述。
表1-9 矩阵函数
命令
说明
d=eig(A)
[v,d]=eig(A)
det(A)
expm(A)
inv(A)
logm(A)
norm(A)
norm(A,1)
norm(A,2)
norm(A, inf)
norm(A,p)
norm(A,’fro’)
null(A)
orth(A)
pinv(A)
poly(A)
schur(A)
sqrtm(A)
svd(A)
TRACE(A)
矩阵特征值
矩阵特征值与特征向量
行列式计算
矩阵求幂
矩阵的逆
矩阵的对数
矩阵和向量的范数
1——范数
2——范数(欧几里德范数)
无穷大范数
P——范数(只对向量)
F——范数
零空间
正交化
伪逆
特征多项式
分解
矩阵平方根
奇异值分解
对角元素之和
注意,上述矩阵函数,如矩阵求幂等运算是通过级数求出的,与矩阵元素的点运算结果不同。若希望求矩阵每个元素的相应函数,如sin,tan,log等运算,可以直接采用标量运算命令,如有
a=
1 2
3 4
如键入
sin(a)
则有
如键入
log(a)
则有
如键入
exp(a) %求[eaij]
即有
ans=
2.7183 7.3891
20.0855 54.5982
矩阵作为指数的运算:
如键入
expm(a) %求矩阵的幂级数
则有
给出了不同的结果。这是由于expm(a)是矩阵的幂级数,其计算公式为:
其它矩阵函数,如矩阵的对数也同样,使用时应注意。
1.5.4 稀疏矩阵
稀疏矩阵是一个阶级很高,只有极少部分元素非零的矩阵。用一般的矩阵表示时,由于要定义每一个元素,因此占用很大的存储空间。而如果用稀疏矩阵定义,可以只定义非零元素,因此可以节约大量的空间,如100×100维的单位矩阵,只有100个对角元素非零,用一般矩阵定义时需要10000个存储单元,而用稀疏矩阵定义,只有100个单位。考虑到高维矩阵表示不便,本文只有10维矩阵作为例题。
speye命令可以建立单位稀疏矩阵。例如,键入
a=speye(5) %建立5阶单位稀硫矩阵
则有
a=
(1,1) 1
(2,2) 1
(3,3) 1
(4,4) 1
(5,5) 1
如键入b=eye(5) %建立5阶普通矩阵
则有
若要查询a,b矩阵状态,键入
whos
则有
Grand total is 30 elements using 284 bytes
可以看出,采用稀疏矩阵定义,只有84个字节,而用一般矩阵定义,则需要200个字节。
Sparse(i,j,s)命令可以产生最大行为i,最大列为j,最后一个元素为s的稀疏矩阵。
如键入
aa=sparse(10,10,1) %建立10×10阶稀疏矩阵,最后一个元素为1
则有
如键入
ap=sparse(1:5,6,0.5) %建立5×6稀疏矩阵,第6列元素为0.5
则有
如键入
ap=sparse(1:5,1:5,0.5) %建立5×5稀疏矩阵,对角线元素为0.5
则有
注意,用最后一种指令形式时,行、列数要相等。
另外,B=sparse(A)命令可以将普通A矩阵转换为稀疏矩阵B。A=full(B)命令可以将稀疏矩阵B转换为普通矩阵A。其它稀疏矩阵的定义和运算命令在表1-10中列出,可以根据需要自行应用。
表1-10 稀疏矩阵函数
命令
说明
find
gplot
nnz
nonzeros
normest
nzmax
randperm
speye
spones
sprandn
sprandsym
sprank
spy
symmd
求非零项的下标
按图论画出稀疏矩阵
非零元素的数目
寻找非零元素
估计2—范数
为非零元素分配的内存数
随机排列向量
稀疏单位矩阵
用1代替非零元素
随机稀疏矩阵
随机稀疏对称矩阵
稀疏矩阵的轶
显示稀疏结构
最小对称度
1.6 M函数与M文件
MATLAB除了可以进行前面讲过的数学函数、矩阵函数运算之外,还提供了M函数、M函数文件和M文件功能,可以利用所有的已知函数编制自己的M函数或M文件,完成更为复杂的运算。实际上,MATLAB的许多复杂函数,如求特征值、特征向量、曲线拟合、插值运算等都是由建立在M函数基础上的M文件完成的。它使得自己的需求可以任意扩展和抽象化,给MATLAB增加了抽象思维能力。
1.6.1 M函数与M函数文件
1.M函数
MATLAB的M函数是function语句引导的,其基本格式如下:
Function[y1,y2,…]=ff(x1,x2,…)
其中,ff为函数名,xi和yi分别为输入和输出变量。它们可以是标量、数组、矩阵或字符串,例如,可用M函数完成的运算。给出M函数ff,内容如下:
function[p]=ff(x)
n=length(x);
for I=1:n
pp(i)=sqrt(x(i)?2+10);
end
p=pp*2-5;
然后在工作空间定义x:
x=1:5
x=
1 2 3 4 5
调用ff函数,可得:
在function命令中也可以完成没有输入输出变量,简单定义为:
function文件名
它执行该文件下指定的操作,例如,M函数的文件test1内容为:
function test1
a=’function test’
b=[1 2; 3 4]
它执行定义a为符号变量,b为矩阵的简单运算。在工作空间中输入文件名:
test1
结果显示:
应当注意,M函数中使用的变量,除输入和输出变量以外,所有变量都是局部变量,即在该函数返回之后,这些变量会自动在MATLAB的工作空间中清除掉。如果想使这些中间变量在工作空间中起作用,则应该把它们设置成全局变量。全局变量是由MATLAB提供的global命令来设置的,一般在M函数的开头定义。命令形式为:
global a b c
不同的全局变量名a b c用空格隔开。Global命令应当在工作空间和M函数中都出现。如果只在一方出现,则不被承认为全局变量。例如,上面的例题中,如果要增加全局变量z1,z2参与运算,则首先定义
global z1 z2 % 在工作空间中定义全局变量
键入
z1=1: -0.1: 0.6
则有
若键入
z2=0:0.5:2
则有
然后修改M函数:
function [p]=ff(x)
global z1 z2 % 增加全局变量
n=length(x)
for i=1:n
pp(i)=sqrt(x(i)?2+10)+z1(i)-z2(i); %z1,z2加入运算
end
p=pp*2-5;
运行M函数
y=ff(x)
y=
3.6332 3.2833 3.3178 3.5980 4.0322
结果为
2.M函数文件
除了MATLAB提供的M函数,所有自己定义的M函数都要通过一个M文件来产生,M文件是指扩展名为.m的文件。M函数文件必须按照一定的规则来编写。它的基本规则和属性如下:
·函数名和文件名必须相同,如上面定义的M函数test1,它是一个名为test1,m的文件。开头应当以function语句开始,第2条以后可以加入注释行(以%开始)和MATLAB运算语句。
·M函数文件有自己的工作空间,与MATAB的工作空间分开。M函数的变量都是内部变量,不会送工作空间去。M函数与工作空间的联系只有输入、输出变量。
·M函数中若有return命令,函数将中断运行,返回工作空间。
·MATLAB在执行一次M函数之后,将其编为机器码,下一次运行时直接调用,运算速度很快。
·M函数文件中可以调用其它一般M文件(或称为脚本文件)。这时,脚本文件中的变量都在M函数文件中有效,而在MATLAB的工作空间中无效。
·M函数文件可以重复调用自己,但应当尽量避免这种操作,因为容易形成死循环。
1.6.2 M文件
MATLAB提供的M文件(也叫脚本文件)是普通的ASCⅡ码构成的文件,只能由MATLAB语言所支持的语句组成。它类似于DOS下的批处理文件,可以直接在工作空间中执行。M文件以扩展名,m结尾,执行时只需键入文件名,MATLAB会自动执行该M文件中的各条语句。与M函数文件不同,它不以function语句开始,可以直接以任意MATLAB语句开始。如果要定义全局变量,则应以全局变量定义开始。
在工作空间中调用的M文件中的所有变量都可以在工作空间中直接引用,并可以在工作空间中的其它M文件或SIMLINK文件进行时调用。但如果要在某个M函数中引用,则必须定义为全局变量。M文件可以调用M函数文件,一般来讲,M文件应当在M函数的上层。
MATLAB的M文件的功能非常强大。它允许自由编制充分复杂的程序,调用各种已有的函数、其它M文件等。完成复杂的数值、逻辑和符号运算,与绘图功能、SIMULINK工具箱等联合使用,还可以仿真和绘制复杂的运动过程,是一个非常有用的工具。
例如上面的例子可以简单地用一个M文件完成。该文件名为f2.m,文件内容如下:
global z1 z2
x=1:5
z1=1:-0.1:0.6
z2=0:0.5:2
y=ff(x)
为了运行M函数,建立M函数文件ff.m,内容如下:
function[p]=ff(x)
global z1 z2
n=length(x)
for i=1:n
pp(i)=sqrt(x(i)?2+10)+z1(i)-z2(i);
end
p=pp/2-5;
运行f2.m文件,键入f2,结果显示:
有关M文件的操作命令,在表1-11中给出。
表1-11 有关M文件的函数
命令
说明
disp(ans)
echo
input
keyboard
pause
pause(n)
waitforbuttonpress
显示未指定变量名的结果
控制命令窗口对脚本文件命令的回响
提示用户输入
暂时把控制权交给键盘(按reture退出)
暂停,直至用户按任意键
暂停ns
暂停,直至用户按鼠标键或键盘键
1.7 多项式运算
多项式运算是线性代数、线性系统分析中的重要内容。MATLAB提供了多条命令,可以进行多项式运算。多项式用一个行向量表示,多项式的系数降幂排列。例如,多项式:
p(x)=x4+2x3-5x+6
在MATLAB工作空间输入为:
p=[1 2 0 -5 6]
则有
1.7.1 求根及其逆运算
roots命令可以求解多项式p的根,求出的根按列向量存储。如求上面给出的多项式的根,则键入
rr=roots(p)
结果为
poly命令可以由根的列向量表示求多项式系数,得到的多项式系数按行向量存储。如对于上面求出的根,键入
pp=poly(rr)
得到
当根有共轭虚部时,多项式系数中可能出现小的虚部,对于线性系统,这是不可能的。这时,可以利用pp=rela(pp)命令类消除虚假的虚部。
1.7.2 加法、减法与乘法
两个多项式的加、减法为多项式对应元素的加、减运算。两个多项式的阶数可以不同,但在多项式定义时应当补充0元素使其行向量元素数目相等,否则不能相加、减。如已知
p(x)=x4+2x3-5x+6
s(x)=x2+2x+3
键入
p=
1 2 0 -5 6
s=
0 0 1 2 3
则有
结果为
结果为
多项式系数定义为
a=
1 1 3
b=
1 2 4
键入
c=conv(a,b) % 求多项式的乘积
则有
1.7.3 微分与赋值运算
多项式的微分由polyder命令完成。如
p=
1 2 0 -5 6
键入
d=polyder(p) % 对多项式p(x)求微分
得
给出p(x)中自变量的范围,polyval命令可以计算多项式的值。如键入下列语句:
x=-1:0.1:2; % x由-1到2,步长0.1,共30个点
y=polyval (p,x);
可以得到一组对应于x的y的值(30个点),可以利用绘图命令plot(x,y),画出x和y的关系。
1.7.4 有理多项式
在线性系统的Fourier变换、Laplace变换和Z变换中,经常用到有理多项式。MATLAB提供了一些命令可以进行有理多项式的运算。MATLAB中的有理多项式是由分子多项式和分母多项式表示的,可以用residue命令进行部分分式展开。该命令形式为
[r,p,k]=residue(mum,den)
式中,mum和den分别表示分子和分母多项式的系数行向量。
分解的结果形式:
式中k(s)为常数项或纯微分项。
例如,对已知传递函数
进行部分分式分解,运算过程如下:
定义、输入数据为
d=
1 2 4 3 % 定义分母多项式
n=
10 30 % 定义分子多项式
键入
[r,p,k]=residue(n,d) % 进行部分分式展开
结果为:
其中,k=[][表示没有常数项。
结果表明,传递函数G(s)被分解为下面的部分分式形式。
注意,该命令得到的低阶分式都是一阶形式,分子上是常值。显然,由于原传递函数有复特征值,利用该命令得到的展开式不是常用的实数形式。当然,可以利用多项式的乘法命令conv进行通分得到实数形式如下。
若键入
r1=r(1),r2=r(2) % 数值定义成多项式r1,r2
则显示
键入
p1=[1—p(1)],p2=[1—p(2)] % 定义两个极点多项式p1(s)=s—p(1),p2(s)=s—p(2)
显示
键入den=vonv(p1,p2) % 求分母多项式den=p1(s)*p2(s)
显示
键入
mum=conv(r1,p2) + vonv(r2,p1) % 求分子多项式
显示
得到的部分分式形式为
根据给出的r,p,k的值,可以用一个命令求出传递函数的有理多项式。如利用上面求出的结果,键入命令:
[num,den]=residue(r,p,k)
可以得到
可以看到,求出的分子多项式和分母多项式与给定的传递函数G(s)形式相同。多项式运算的其它命令和功能,在表1—12中列出。
命令
说明
Conv(a,b)
[q,r]=deconv(a,b)
poly(r)
polyder(a)
polyfit(x,y,n)
polyval(p,x)
[r,p,k]=residue(a,b)
[a,b]=residue(r,p,k)
roots(a)
乘法
除法
用根构造多项式
对多项式或有理多项式求导
多项式数据拟合
计算x对应的多项式的值
部分分式展开式
部分分式组合
求多项式的根
1.8 控制语句与逻辑运算
控制语句包括循环与条件语句。它们决定了运算过程和路径,因此被称为控制语句。循环和条件语句包含在每一种可以用于进行科学计算的计算机语言中。它们更适合人的思维,扩展了计算功能,并节省了语句,使程序看来更为简洁、清晰。
MATLAB的循环和条件语句中,由于经常包括了饿大量的MATLAB命令,一般用于M文件中更为合适(可以查阅M文件说明),而不常采用直接在工作空间中键入方式。下面介绍的在工作空间中键入的程序,可以直接写入M文件调用。
1.8.1 for循环
for循环语句允许按照给出的范围或固定的次数重复完成一个(或一组)运算。它从for开始,用end结束,也叫做for-end结构。for语句的基本格式为:
for 循环变量=数组范围
命令串
end
执行for语句时,循环变量按照数组指定的范围逐步取值,每一步执行一次命令串,直至循环变量按照数组指定全部取值完毕。例如,运行下面的for循环子程序。
键入
for n=1:5 % 循环变量取值从1到5,每步按1递增
x(n)=n?2; % 运算命令
end % 结束循环运算
再键入
x
结果为
循环变量的范围可以是任意数组,如:
1 2 3 4
5 6 7 8
9 10 11 12
执行如下for循环运算:
for i=a
y=i(1)-i(2)+i(3)
end
结果为
由于A是3行4列矩阵,计算时将a按列分步赋值给变量i,所以,在每一步循环时,i取a的一列进行运算。很明显,如执行下面的循环运算:
for n=a
y=n
end
结果为
这种用法扩展了for语句的应用范围。另外,for语句可以嵌套,例如:
for i=1:3
for j=5:-1:1
a(i,j)=i?2+j?2
end
end
执行完上述for语句,若想知道结果,键入a
得到
1.8.2 While循环
for循环以固定的次数求一组命令的值。MATLAB提供了另一种循环语句while。它根据给出的条件,以不定的次数求一组命令的值。该循环语句的结构为:
while 条件表达式
命令串
end
其执行方式为,若条件表达式中的条件成立,则执行命令串;如果表达式不成立,则跳出循环,向下继续执行。例如执行下面的循环语句:
s=0; n=1;
while n<=10 % 循环条件
s=s+n; n=n+1 % 命令
end
若想知道结果,键入
s
得到
该循环只进行到n=10为止。
1.8.3 条件语句
除了前面介绍的循环语句结构之外,MATLAB还提供了条件转移语句,得到MATLAB更易于使用。条件语句的格式为
if 条件表达式
命令串
end
当给出的条件表达式成立时,执行命令语句,然后继续向下执行;若条件不成立,则跳出条件块而直接向下执行。循环语句和条件语句中的条件表达式用逻辑关系符号表示。如“大于等于”用“>=”表示,“等于”用“= =”表示,可参看1.8.4节“关系与逻辑运算”。
例如:
y=0;
for n=1:4
if n>2
end
end
结构为
很明显,当n=1,2时,不满足条件,程序转而执行下面的循环。
条件语句还有if-else-end结构,结构形式为:
if 条件表达式1
命令串1
elseif 条件表达式2
命令串2
┇
else
命令串3
end
执行上述语句时,如果条件表达式1的条件成立,那么就执行命令串1;如果条件1不成立,条件2成立,则执行命令串2;否则执行命令串3。MATLAB允许多层不相交的条件语句嵌套。
另外,在执行for和while循环语句时,可以利用if+break语句中止循环运算。如:
sum=0;
for m=1:100 % 循环变量从1到100
if (sum>100)
m % 当SUM大于100时,显示m
break; end % 中止运算
sum=sum+m;
end
运行结果:
结果表明,当运行到第15个循环时,满足给定条件,程序中止运行。
1.8.4 关系和逻辑运算
除了传统的数学运算,MATLAB还支持关系和逻辑运算,目的是提供求解真/假命题的答案。对于所有关系和逻辑表达式的输入,ATLAB把任何非零数值当做真,把零当做假。而对于所有关系和逻辑表达式的结果为真时输出,1假时输出0。MATLAB给出的关系操作符如表1-13。
表1-13 关系操作符
关系操作符
说明
逻辑操作符
说明
<
<=
>
>=
= =
≈=
小于
小于等于
大于
大于等于
等于
不等于
&
|
~
与
或
非
1.关系操作符
关系操作符用来比较两个同样大小的数组,或比较一个数组与一个标量。数组和标量比较时,数组的每一个元素与标量比较,结果数组与原数组大小一样,例如已有数组a,b分别为
a=
1 2 3 4 5
b=
3 4 5 6 7
键入
tl=(a>3) % 所有a>3的位置赋值1,否则赋值0
得到
键入
t2=(a>3)-b % 完成上面运算,结果数组减b的对应元素
得到
结果表明,由于关系运算的结果是0或1组成的数组,因此也可以用于一般数学运算中。
2.逻辑操作符
逻辑操作符提供了一种按照逻辑“与”。“或”,“非”形成的关系表达式,并可以用于运算,例如命令
t3=~(t1) % 取t1的“非”(0变为1,1变为0)
的结果为
键入
t4=(a>1)&(b<6) %当a>1和b<6同时成立到赋值1,否则赋值0
得到
用逻辑运算还可以截取函数的部分段,产生不连续信号。方法是首先由逻辑运算产生0或1,要保留的部分与1乘,不保留的部分与0乘。例如执行下面的程序
x=0:0.1:10; % 定义x范围
y=sin(x) % 计算y
z=(y>=0). * y; % 将y中的负值用0代替,正值保持不变
plot(x,y,’k’,x,z,’o’) % 绘图,y用一表示,x用O表示
绘图结果见图1-4。
3.NaN与空矩阵
MATLAB提供了函数NaN, Inf和空矩阵运算,表示一些特殊的概念。其中,NaN表示“Not-a-Number”,专门表示一类由0/0或Inf-Inf产生的数值;Inf表示“无穷大”,是由1/0或计算中的数值上溢产生的。
按照IEEE数字标准,对NaN的所有运算都产生NaN,例如健入
a=[1 2 nan inf] % 定义数组中有元素取值为NaN和Inf
得到
键入
b=2*a % 对NaN,Inf运算
得到
如取
c=sqrt(b)
则
如取
c(4)=c(4)-inf % Inf-Inf产生NaN
则
如取
d=(a= = NaN) % 寻找数组中的NaN
则
d中元素全部为0,由此可见,NaN与NaN相比较时,产生全部为0的结果。用户可以用isnan命令找取值为NaN的元素。键入
g=isnan(a)
可得到
其中,第3个元素值为1,表示“真”,则a矩阵中第3个元素是NaN。
空矩阵用符号[ ]表示,它不是元素取值为0的矩阵,而是一个行列数为0的矩阵。换句话说,空矩阵是一个标志,表示逻辑上的“无”或“不存在”。例如健入
size ([ ]) % 查询空矩阵的维数
可得到
可知,空矩阵的行、列数都是0。
在运算中若没有结果时,返回空矩阵,执行下面的程序:
x=1:5
x=
1 2 3 4 5
y=find(x>6)
y=
由于没有x>6的值,所以返回空矩阵。
1.8.5 字符运算
MATLAB可以给一串文字定义并进行有些字符串的处理和运算。字符串一般是ASCII码的数值数组,用单引号括起来,定义为字符表达式。字符表达式中每一个字符占用2个字节存储。例如:
键入
x=’matlab is a good software’
显示
键入查询命令
whos
结果为
经查询,x是文字型数组,共有25个元素,占用50个字节。
MATLAB定义了一些字符串转换函数,如表1-14。
表1-14 字符串转换函数
命令
说明
Abs
Fprintf
Int2str
Lower
Num2str
Setstr
Sprintf
Sscanf
Str2mat
Str2num
upper
字符串到ASCII码转换
按照给定格式把文本写到文件中或显示屏上
整数转换成字符串
字符串变为小写
数字转换成字符串
ASCⅡ码转换成字符串
按照给定格式,数字转换成字符串
按照给定格式交字符串转换成数字
字符串转换成一个文本矩阵
字符串转换成数字
字符串转换成大写
下面仅对表1-14中几个常用的命令加以介绍。
Num2str命令可以将数值变量转换成字符串。它可以与其它字符串组成新的文本。例名键入
a=2.7;
xx=[‘there are’, num2str(a), ‘kg eggs.’]
显示
xx=
there are 2.7kg eggs.
键入命令
upper(xx)
结果为
ans=
THERE ARE 2.7KG EGGS.
fprintf在前面数据的输入输出部分已有解释,这里可以用来写字符串。如
键入
fprintf(‘pi= %. 0e\n’, pi) % e型数,写pi,不写小数部分
显示
键入
fprintf(‘pi=%. 5e\n’, pi) % e型数,写pi,小数部分取5位
显示
键入
fprintf(‘pi= %. 0f\n’, pi) % f型数,写pi,不写小数部分
显示
键入
fprintf(‘pi= %. 5f\n’, pi) % f型数,写pi,小数部分取5位
显示
键入
fprintf(‘pi= %. 0g\n’, pi) % g型数,写pi,不写小数部分
显示
键入
fprintf(‘pi= %. 5g\n’, pi) % g型数,写pi,加小数部分共取5位
显示
键入
fprintf(‘pi= %. 10g\n’, pi) % g型数,写pi,加小数部分共取10位
显示
其中,g型数为固定小数点,不论取多少位,小数点位置不变。
MATLAB还提供了字符串运算函数,使得MATLAB的计算功能更加强大。如eval命令和feval命令,可以运行已有的M函数,计算并赋值给其它变量。如:
键入命令
a=eval(‘sqrt(2)’) % 计算2的平方根,赋于a
结果为
如果已建立了一个M函数,名为f1·m,内容如下:
function[y]=f1(x)
y=x?2+2
该程序输入x,计算并输出y。利用eval命令可以调用M函数f1并计算y值。如键入
a=eval(‘fl(2)’) % 调用M函数f1,其输入变量x=2,结果赋值给a
结果为
设另一个M函数,名为f2,内容如下:
function[p1,p2]=f2(x)
pl=x(1)?2+x(2)?2+x(3)?2;
p2=sqrt(p1);
该函数计算数组x各元素的平方和并开方,结果存入p1,p2。运行该程序,首先定义
a=
1 2 3
键入
[x y ]=eval(‘f2(a)’) % 调用M函数f2,以a为输入变量,x,y为输出变量
结果为
x=
14
y=
3.7417
feval用法与eval类似,约束更多一些。这两条命令应用很广泛,在很多Simulink程序中可以看到它们,读者可以自行研究并开发自己的程序。
1.9 曲线拟合与插值运算
1.9.1 曲线拟合
在许多应用领域中,人们经常需要从一系列已知离散点上的数据集[(x1,y1),(x2,y2),…(xn,yn)]得到一个解析函数y=f(x)。得到的解析函数f(x)应当在原离散点xi上尽可能接近给寂静的yi的值。这一过程称为曲线拟合。最常用的曲线拟合是最小二乘法曲线拟合,拟合结果可使误差的平方和最小,即使出使
最小的f(x)。
如果自己去编程序,求解最小二乘解,过程是比较复杂的。MATLAB提供的函数polyfit,根据给定的自变量数组x和函数数组y,按照拟合的阶数要求自动求解满足最小二乘意义的一阶或高阶解析函数f(x),使用很方便。为了说明这个问题,我们取以下函数为例。
Y=0.5-2*x?2
end
键入
y
得到
令y的每一项数据有一些偏差,构成待拟合的数据集:
y=[0.52 0.4 0.4 0.35 0.18 0.02 -0.25 -0.4 -0.81 -1.1 -1.5]
y=
Columns 1 through 7
0.5200 0.4500 0.4000 0.3500 0.1800 0.0200 -0.2500
Columns 8 through 11
-0.4000 0.8100 -1.1000 -1.5000
分别进行一维和二维曲线拟合:
键入命令
m=1;
fxy1=polyfit(x,y,m) % 一维拟合
fxy1=
-1.9873 0.7991
拟合出的多项式为fxy1=-1.9873x+0.7991,是一个线性方程。
键入命令:
m=2;
fxy2=polyfit(x,y,m) % 二维拟合
fxy2=
-2.0350 0.0477 0.4938
拟合出的多项式为fxy2=-2.035x?2+0.0477x+0.4938,是一个二阶方程。得明显,二维曲线拟合得到的结果与给定的多项式y=0.5-2*x?2更为接近,一维拟合效果较差。利用多项式估计命令分别计算一、二维拟合结果的数值,并与原函数值绘在一张图上。
y1=polyval(fxy1,x);
y2=polyval(fxy2,x);
plot(x,y,’o’, x,yl, ‘k:’ , x, y2, ‘k’)
由图1-5可以看出,二维拟合结果与原函数值基本重合。
为了使用polyfit命令,我们必须给出自变量数据组和期望的最佳拟合数据的多项式的阶数。选择不同的阶数,会得到不同的拟合效果,这一点从上面的举例中可以明显看出,如同选择阶数,需要在系统辨识方面有更多的知识,这里不再详细讨论。
1.9.2 插值运算
与曲线拟合不同,插值运算不是试图找出适合于所有自变量数组x的全局最优拟合函数y=f(x),而是要找到一个解析函数连接自变量相邻的两个点(xi,xi+1),由此还可以找到两点间的数值。根据自变量的维数不同,插值方法可以分为一维插值y=f(x)和二维插值z=f(x,y)等。在许多工程问题上,我们只能获得无规律的离散点上的数值,插值可以帮助我们得到近似的连续过程,便于用数学的解析方法对已有数据进行处理和运算,因此是一个很有用的工具。
应当注意的是,使用插值命令时,有两个约束:只能在自变量x的取值范围内进行插值运算;自变量x必须是单调的。
1.一维插值
一维插值是指对一个自变量的插值。MATLAB的一维插值命令可以进行线性插值、三次样条插值和三次插值等。
Interp1是一维插值命令。命令形式如下:
z=interpl(x,y,x0, ‘method’)
式中,x,y为给定数组,x0是要计算的x的位置,可以是一个或一组数据,全部数据应当在x的取值范围内。Method是所用的插值方法。线性插值使用下列方法:
‘nearest’——用最接近的相邻点插值;
‘linear’线性插值(缺省);
‘spline’——三次样条插值;
‘cubic’——三次插值。
·线性插值
线性插值的最简单的例子是MATLAB的作图。按缺省设置,MATLAB作图时用直线连接所有的数据点,它利用线性插值估计数据点之间的取值。例如:
x1=0:0.01:10; % 定义x1取值范围,共1000个点
y1=sin(x1);
x2=0:10; % 定义x2取值范围,共10个点
y2=sin(x2);
plot(x1, y1, ‘k’, x2, y2, ‘b*’) % 绘图,y1用‘—’,y2用“*”
plot(x2, y2, ‘k’) % 只绘 y2, 用‘—’
绘图结果如图1-6、图1-7。
图1-6表明,计算y2所用步长较大,y2数组只是y1的子集。图1-7表明,单独用直线画y2时,任意两点间用直线连接,程序是用线性插值来得到两点间的数值并绘出的。很明显,取不同的步长,显示的曲线的近似程度不同。对于一维插值,步长越小越好。
例如,x,y取值如下
x=
0 1 2 3 4 5 6 7 8 9 10
y=
Columns 1 thorugh 7
0.5200 0.4500 0.4000 0.3500 0.1800 0.0200 -0.2500
Columns 8 thorugh 11
-0.4000 -0.8100 -1.1000 -1.5000
进行下列线性插值运算:
键入
p=interpl(x,y,8.5) % 用线性插值求x=8.5处的值
可得
定义
x0=[3.4 4.7 7 6.1 7.6] % 定义数组x0
显示出
键入命令
p1=interpl(x,y,x0) % 用线性插值求数组x0对应点上的值
结果为
得出的一组数据为x0对应点上的y的数值。
·三次样条插值
三次样条插值用一个3阶多项式连接相邻的两个点。在每一个点上,相邻两个连续函数的1阶、2阶导数相同。这样可以使连接后的曲线光滑起来,不会出现线性插值的折线连接效果。使用三次样条插值时,在interp1命令中加入‘spline’标志。例如,对上例键入命令
pp=interpl(x,y,8.5, ‘spline’) % 用三次样条插值求x=8.5处的值
结果为
可以看出,用三次样条插值得到的结果与用线性插值得到的结果有所不同。一般来讲,用三次样条插值会得到更为精确的结果,但其计算较为复杂,因为必须找出3阶多项式。如果步长足够小,采用线性插值更为合理。
用‘cubic’标志的三次函数插值与三次样条插值不同,也是要找到3阶多项式,但使用的方法不同。
2.二维插值
二维插值基于与一维插值同样的思想,但是针对两个自变量的函数z=f(x,y)进行插值。二维插值简单地可以理解为连续三维空间函数的取值运算,如求解随平面位置变化的温度、湿度、气压等。
二维插值命令形式为
zi=interp2(x,y,z,x0,y0, ‘method’)
式中,x,y为自变量数组,z为测量数组,x0,y0为指定的自变量插值计算点数组,method是二维插值使用的方法,包括
‘nearest’——最近点插值
‘linear’——双线性插值(缺省)
‘cubic’——三次函数插值
二维插值要求所有自变量取值都是单调的。
例如,考虑一个平面范围内的温度变化,在多个离散点上测量出当地的温度,数据如下:
x=
1 2 3 4 5 6 % 自变量x的取值范围
y=
1 2 3 4 % y的取值范围
t= % 温度t的测量值
12 10 11 11 13 15
16 22 28 35 27 20
18 21 26 32 28 25
20 25 30 33 32 30
t以矩阵形式给出,t的列数与x的维数与x的维数相同,行数与y的维数相同。用三维绘图命令
mesh(x,y,t)
可以画出温度变化的三维空间网格图形,如图1-8用直线连接相邻两点。
如果只考虑一个方向上的温度变化,如取:
x0=1:0.1:6; % 细化x范围
y0= % 指定y为常值
2
t1=interp2(x,y,t,x0,y0); % 线性插值
t2=interp2(x,y,t,x0,y0, ‘cubic’) % 三次插值
plot(x0, t1, ‘k’ , x0, t2, ‘:’) % 线性插值用‘—’,三次插值用“:”
比较两个插值结果如图1-9。
如果考虑两个方向上插值,可以运行下列程序:
x0=1:0.1:6; % 将x0细化为51个点
y0=1:0.06:4; % 将y0细化为51个点
编制M函数f1:
function[t3]=f1(x,y,t,x0,y0)
n=length(y0) % 测量y0的维数
for i=1:n
tt=interp2(x,y,t,x0,y0(i)); % 对每一个y0(i),x0进行二维插值
for j=1: length(t)
t3(i,j)=tt(j); % 插值结果赋于t3矩阵第i行
end
end
运行该M函数,结果存在t3(x0,y0)矩阵,绘图:
t3=f1(x,y,t,x0,y0);
mesh(x0,y0,t3) % 绘x0,y0,t3的三维图形
绘出的图形如图1-10。与图1-8比较可以看出,在x,y范围细化以后,采用线性二维插值,结果比一维插值更为精确。
有关曲线拟合和插值的命令见表1-15。
表1-15 曲线拟合和插值函数
命令
说明
Polyfit(x,y,n)
Interpl(x,y,no)
Interpl(x,y,no, ′spline′)
Interpl(x,y,xo, ′cubic′)
Interp2(x,y,z,xi,yi)
Interp2(x,y,z,xi,yi, ′cubic′)
Interp2(x,y,z,xi,yi, ′nearest′)
对描述n阶多项式y=f(x)的数据进行最小二乘曲线拟合
一维线性插值
一维三次样条插值
一维三次插值
二维线性插值
二维三次插值
二维最近邻插值
1.10 符号运算
MATLAB具有符号数学工具箱,可以对符号表达式进行运算和处理,基本运算包括复合、化简、微分、积分以及求解代数方程式、微分方程式等。另外,还可以求解线性代数问题,如求解符号矩阵的逆、行列式、正则行的精确结果,找出符号矩阵的特征值而没有由数值计算引在运算最后将数字代入结果,因此避免了中间运算的误差,能够以指定的精度返回结果。
符号数学工具箱的功能建立在Maple软件的基础上。该软件最初是由加拿大的滑铁卢(Waterloo)大学开发的。当要求MATLAB进行符号运算时,它就转入Maple去计算并将结果返回到MATLAB命令窗口。因此,MATLAB中的符号运算是MATLAB处理数字功能的自然扩展。
1.10.1 符号表达式
符号表达式是包括数字、代数或有理运算和符号变量的MATLAB字符串。它不要求变量有预先确定的值。符号方程式是含有符号的符号表达式。符号矩阵是数组,其元素是符号表达式。符号运算是使用已知的数字规则和给定的符号恒等式求解这些符号方程,它与代数和微积分所学到的求解方法完全一样。
1.创建符号表达式
MATLAB用sym(‘’)命令建立符号表达式。符号表达式表示成字符串,用单引号‘’括起来,以与数字变量运算相区别。否则,这些符号表达式几乎完全像MATLAB命令。表1-16列有几则符号表达和MATLAB等效表达式的例子。
表1-16 符号表达式与等效的MATLAB表达式
符号表达式
MATLAB符号表达式
‘1/(2*x?n)’
‘1/sprt(2*x)’
‘cos(x?2)-sin(2*x)’
m=sym(‘[a,b;c,d]’)
f=int(‘x?3/sqrt(1-x)’,‘a’,‘b’)
MATLAB符号函数可用多种方法来操作这些表达式。例如,定义符号表达式:
Y=sym(‘cos(x)’)
结果为
对符号表达式进行运算,键入
c1=diff(y) % 求微分
结果为
某些简单的表达式运算可以不用sym进行定义,直接对符号表达式进行运算,如键入
c2=diff(‘cos(x)’)
结果为
符号运算的结果也产生符号变量,键入查询命令
whos
得到
表明y是1×6维文字型数组,用12个字节。c1,c2都是一维符号表达式,各占138个字节,良种方法得到的结果一样。
如果直接定义可能与MATLAB数值函数运算混淆,则必须用sym命令定义。例如,下列命令
m=[a,b;c,d]
定义数值矩阵m,其元素为a,b,c,d。由于矩阵元素没有事先赋值,程序给出错误信息:
???Undefined function or variable ‘a’。
而下列命令定义m为字符串,m=‘[a,b;c,d]’
结果为
只有下面的命令定义符号矩阵,
m=sym(‘[a,b;c,d]’)
结果为
另外,很多符号运算命令,如numden(提取分子分母多项式)只在sym定义下可以运行,在简单定义时出现错误信息。在使用时可能会遇到这个问题,这时需要重新用sym定义。
2.符号常量
不含变量的符号表达式叫做符号常量。例如键入命令
f=sym(‘2*4-6’) % 定义符号表达式,不含符号变量
结果为
若求f的数值(进行数值运算),键入
f1=numeric(f)
结果为
若进行符号运算,键入
f2=f+1
结果为
若对应符号为常量,其导数应为0。此时如果键入
diff(f)
则显示
ans=
0
检验上述变量类型,键入
whos
则显示
Name Size Bytes Class
ans 1×1 126 sym object
f 1×1 134 sym object
f1 1×1 8 double array
f2 1×1 126 sym object
Grand total is 11 elements using 394 bytes
f, f2和ans虽看来都是数字,但仍然是按照符号表达式存储的。只有f1是数值运算的结果,定义为8位的字节数。
结果表明,符号常量定义及运算结果仍属于符号变量范围,只有在数值运算之后,才变为数值变量。另外,符号运算占用较大的空间。
3.符号变量
当字符表达式中含有多于一个的变量时,只有一个变量是独立变量,其余的文字符号作为常量处理。如果不指定哪一个变量是独立变量,MATLAB将基于以下规则选择一个独立变量:
·除去i和j的小写字母,表达式中如果没有其它字母,选择x作为独立变量;
·如果有多个字符变量,选择在字母顺序中最接近x的字符变量;
·如果有相连的字母,选择在字母表中较后的那一个。
例如,键入:
diff(‘sin(x)+1’) % 只含有一个字符变量,该字符就是独立变量
则有
键入
diff(‘sin(a)+b’) % 只含有两个字符变量,字母表中靠后的是独立变量
则有
键入
diff(‘3*y+z’) % 含有两个字符变量,接近x的独立变量
则有
若键入
diff(‘sin(pi/2)-log(3/5)’) % 符号常量,取x为独立变量
由于表达式中不含x,所以得到该符号表达式的微分为
ans=
0
若键入
diff(‘sin(omiga)+a’) % omiga不是单个字符,也可以作为独立变量
则有
如果不熟悉或不愿意使用上述隐含规则,可以自己指定独立变量。例如键入:
diff(‘x?n’) % 未指定独立变量,将缺省变量x作为独立变量
则有
结果为对x的导数。
若键入
diff(‘x?n’, ‘n’) % 指定n为独立变量
结果为
结果为对n的导数。
1.10.2 符号表达式运算
一旦建立了一个符号表达式,符号运算功能可以完成,如提取表达式的一部分、合并两个表达式、求表达式的数值以及表达式的加、减、乘、除等运算。
1.提取分子和分母
如果表达式是一个有理分式(两个多项式之比),或是可以展开为有理分式(包括那些分母为1的分式)可利用numden命令来提取分子或分母,必要时还可以进行表达式合并。
例如,健入
m=sym(′x?2′) % 定义有理分式m
则有
键入
[nm,dm]=numden(m) % 求m的分子多项式nm(x)和分母多项式dm(x)
得到
由于m表达式中没有分母,结果dm=1。又如,键入
f=sym(′x?2/y?2) % 定义有理分式f
得到
键入
[n,d]=numden(f) % 求f的分子多项式n(x)和分母多项式d(x)
得到
numden命令还可以将有理分式进行通分,合并同类项并提取分子和分母。如,键入
g=sym(′3/2 * x?2+2/3 * x –3/5′) % 定义有理分式g
可得
键入
[ng, dg]=numden(g) % 求分子与分母表达式 ng(x), dg(x)
结果为
分母表达式为30,是通分的结果。
对于符号数组或符号矩阵,该命令对每一个元素进行上述计算。如键入。
a=sym(‘[4/5, 2 * x; 3/x, (x+3)/(2-x)]’) % 定义有理分式a
结果为
α为符号矩阵,键入
[na, da]=numden(a) % 求分子与分母表达式
得到
可以看出,分子和分母表达式分别对每一个元素求出。
2.代数运算
很多标准的代数运算可以在符号表达式上执行,命令symadd, symsub, symmul和symdiv可以进行两个表达式的加、减、乘、除运算,sympow将一个表达式表示为另一个表达式的幂次。例如,键入命令
f=‘2 * x-5’ % 定义符号表达式f
结果为
键入命令
g=‘x?2-x+7’ % 定义符号表达式g
结果为
键入命令
symadd(f,g) % 求f+g
结果为
键入命令
symsub (f,g) % 求f-g
结果为
键入命令
symsub (g,f) % 求g-f
结果为
键入
[nm, dm]=numden(m) % 求m的分子多项式nm(x)和分母多项式dm(x)
得到
由于m表达式中没有分母,结果dm=1。又如,键入
f=sym(′x?2/y?2′) % 定义有理分式f
得到
键入
[n,d]=numden (f) % 求f的分子多项式n(x)和分母多项式d(x)
得到
numden命令还可以将有理分式进行通分,合并同类项并提取分子和分母。如,键入
g=sym(?3/2* x?2+2/3 * x-3/5?) % 定义有理分式g
可得
键入
[ng,dg]=numden(g) % 求分子与分母表达式ng(x), dg(x)
结果为
分母表达式为30,是通分的结果。
对于符号数组或符号矩阵,该命令对每一个元素进行上述计算。如键入
a=sym(‘[4/5,2*x; 3/x, (x+3)/(2-x)]’) % 定义有理分式a
结果为
a为符号矩阵,键入
[na,da]=numden(a) % 求分子与分母表达式
得到
可以看出,分子和分母表达式分别对每一个元素求出。
2.代数运算
很多标准的代数运算可以在符号表达式上执行,命令symadd,symsub,symmul和symdiv可以进行两个表达式的加、减、乘、除运算,sympow将一个表达式表示为另一个表达式的冥次。例如,键入命令
f= ‘2*x-5’ % 定义符号表达式f
结果为
键入命令
g= ‘x^2-x+7’ % 定义符号表达式g
结果为
键入命令
symadd(f,g) % 求f+g
结果为
键入命令
symsub(f,g) %求f—g
结果为
键入命令
symsub(g,f) % 求g—f
结果为
键入命令
symdiv(g,f) % 求g/f
结果为
两个符号表达式相除,只能得出表达式。如果想得到表达式的结果,要用eval命令,如键入
y=symdiv(g,f) % 求y=g/f
得到
给定x=3
显示
x=
3
求y的值,键入
eval(y)
结果为
同样道理,sympow命令只产生结果表达式,若要求数值解,也要使用eval命令,例如键入
p=sympow(g,f) % 求p=g′
p=
(x?2-x+7)?(2 * x-5)
在上面给定的x取值下,求p值,键入
eval(p)
结果为
另外,symop命令可以将两个或多个表达式合并为一个,被合并的表达式由多个单引号括起来的运算符号连接起来,各单引号组之间用逗豆隔开。例如,键入
f=‘cos(x)’
显示
f=
cos(x)
键入
g=‘sin(3*x?2)’
显示
g=
sin(3*x?2)
键入 symop(f,‘/’,g,‘—’,‘5?x’) % 合并f和g,用给定的关系
得到
3.高级运算
符号表达式的高级运算包括表达式的复合、求逆函数、求前n-1项和等。其中,
·命令compose
把f(x)和g(x)复合成f(g(x)),例如键入
f=sym(‘cos(x)’)
显示
键入
g=sym( ‘3*x?2’)
显示
键入命令:compose(f,g) % 求f(g(x))
结果为
键入:compose(g,f) % 求g(f(x))
结果为
如果重新定义g为v的函数,键入
g=sym( ‘3*v?2)
显示
compose命令可以找到对v的复合函数:
键入:compose(f,g, ‘x’, ‘v’) % 求f(g(v))
得到
·命令finverse
求表达式的逆函数。表达式f(x)的逆函数g(x)满足g(f(x))=x。命令finverse给出表达式的逆函数。如果解不唯一,就给出警告。例如键入
p=sym( ‘e?x’)
显示
键入:finverse(p)
得到
结果相当于ln(x),因为ln(e?x)=x。又如已知
g=
3*v?2
键入:finverse(g)
则显示
由于g(v)的逆函数不唯一,程序给出了警告,并找到了另一个逆函数。
·命令symsum
求表达式的前n—1项和即。首先定义
f=sym(‘x?2’)
f=
x?2
键入
symsum(f) % 缺省,求前x—1项和y=12+22+……+(x—1)2
结果为
键入
symsum(f,‘s’) % 求对s的前s—1项和
结果为
结果表明由于函数f中不含有s,于是将x?2作为常值。结果为y=x?2(1+2……+1)从0到s—1共有s项,所以总和为x?2*s。
键入
symsum(f,0,10) % 求从0到10项的和,y=02+……+102
得到
键入
symsum(f,‘s’,0,10) % 求对s的从0到10项的和。
得到
如果定义y为x和v的函数:
y=sym( ‘x?2+3*v’)
显示
键入
symsum(y,‘v’,0,10) % 求对v的从0到10项的和。
该项和为sum=(x2+3×0)+(x2+3×1)+……+(x2+3×10),结果为
4.函数变换
MATLAB的符号运算功能,可以将符号表达式变换成数值或反之。有些符号函数可返回数值。
命令所以们可将MATLAB的一般(数值)表达式转换为符号表达式,前面已经作过介绍。函数numeric的功能正好相反,它把一个符号常数(无变量符号表达式)变换成一个数值。例如键入
numeric(‘1+2?2’)
结果为
命令eval可以计算符号表达式的值。例如,定义f和x为符号变量,如下:
f=
1+x?2
x=sym(‘2’)
显示
键入
p=eval(f) % 求符号表达式的值
结果为
结果p也是符号变量。
命令sym2poly可以将符号表达式转换成数值多项式的系数向量;命令poly2sym的功能相反。例如键入
f=sym(‘1+x?2’) % 定义符号表达式f
显示
键入
fp=sym2poly(f) % f转换成一般多项式fp
显示
键入
fs=poly2sym(fp) % 将fp转换为符号表达式fs
显示
转换后的符号表达式fs与f相同。
命令subs(f, ‘old’, ‘new’)可以在符号表达式中进行变量替换,用“new”字符串代替各个“old”字符串。例如给定
f=
1+x?2
键入
subs(f, ‘x’, ‘v’)
得到
1.10.3 微分和积分
微分和积分是微积分学研究和应用的核心,并广泛地用在许多工程学科。MATLAB符号运算功能可以解决许多这类问题。
1.微分
符号表达式的微分利用命令diff完成。有下面几种形式:
键入
p=‘a*x?2+b*x’ % 定义符号表达式p
显示
p=
a*x?2+b*x
键入
diff(p) % 对缺省变量x求微分
得到
键入
diff(p,‘a’) % 对指定变量a求微分
得到
键入
diff(p,2) % 对缺省变量x求二阶微分
得到
键入
diff(p,‘a’,2) % 对指定变量a求二阶微分
得到
另外,diff命令还可以对给定数组中的数值逐项求差值,例如,定义
a=[(1:8).?2]
显示
a=
1 4 9 16 25 36 49 64
键入命令
diff(a)
结果为
结果求出了a数组中每两项之间的差。
2.积分
符号运算的积分命令为int(f),其中f是一个符号表达式。积分的目的是求出另一个符号表达式p,使其微分满足diff(p)=f。积分命令有多种表达形式,例如键入
f=‘cos(2*x)+s?2’ % 建立符号表达式f
显示
f=
cos(2*x)+s?2
键入
int(f,‘s’) % 对s求积分
得到
键入
int(f, -pi/2, pi/2) %对x求从-pi/2到pi/2的定积分
得到
键入
int(f,‘s’,-2,2) % 对s求从-2到2的定积分
得到
键入
int(f,‘s’,‘m’,‘n’) % 对s求从m到n的定积分
得到
一般来讲,积分比微分复杂得多。一个函数的积分或逆求导的解不一定以封闭形式存在;或也许存在,但软件找不到;或虽然可以明显地求解,但超过内存或时间限制。出于上述种种原因,当MATLAB不能找到逆导数时,它将返回未经计算的命令。例如,
int(‘sin(x)/exp(x?2)’)
结果显示
它首先警告;期望的积分不能得到。然后指出出错的位置,最后返回原命令。
符号微分和积分命令都可以对符号矩阵进行微分和积分。它对矩阵的每一个元素进行运算和处理。
1.10.4 符号表达式画图
MATLAB提供了ezplot命令,将符号表达式可视化。对于一个自变量的函数,可视化实际上是求解自变量各点上的函数值并绘图的过程。例如,键入
y=‘2*x?2-3*x+10’
显示
y=
2*x?2-3*x+10
键入不同命令:
ezplot(y) % 程序自动选择自变量范围绘图
ezplot(y,[-10 10]) % 确定自变量范围绘图
绘图结果如图1-11,图1-12,表明了同一个函数,不同的自变量范围的图形。
1.10.5 符号表达式的简化
对于一些冗长繁复、难以理解的符号表达式,MATLAB提供了许多方法可以将其进行简化、约分、合并同类项等处理,使表达式变得更简洁易懂。包括:
1.命令pretty
以类似于数学课本上的形式(如有理分式)来显示符号表达式。例如键入
f=sym(‘cos(x)/(sin(3*x))?2-5?x’)
显示
f=
cos(x)/(sin(3*x)?2-5?x)
若键入
pretty(f)
则显示
cos(x) x
-5
2
sin(3x)
2.命令collect
可以合并同类项,给出降幂排列形式。若给定
f=(x?2-1)*(x-2) % 对于符号表达式f
键入
collect(f) % 合并同类项
得到
3.命令horner
可将降幂排列的多项式变成嵌套形式。对于上面的结果,键入
horner(ans)
得到
4.命令factor
将表达式分解成因式。对于给定的f,使用命令
foactor(f)
结果为
5.命令expand
展开表达式,给出降幂排列形式,使用命令
expand(f)
结果为
6.命令simplify
将表达式进行简化。Simplify是一个功能强大、通用的工具。它利用各种类型的代数恒等式,包括求和、分解、积分和幂、三角、指数和log函数等来简化表达式。例如定义
y=sym(‘sin(x)?2-2*x+cos(x)?2’) % 定义表达式y
显示
y=
sin(x)?2-2*x+cos(x)?2
使用命令
simplify(y) % 简化表达式y
得到
键入
y=sym(‘x?2-1)*(x+2)/(x+1)’) % 定义表达式y
显示
y=
(x?2-1)*(x+2)/(x+1)
使用命令
simplify(y) % 简化表达式y
结果为
7.命令simple
可在简化基础上进一步给出多种简化形式。符号表达式的简化是很有用的工具,但由于简化方式不同,有时会产生不同的答案。在符号工具箱内的simple命令可以试用几种方式进行简化,然后选择最简形式。它试用的简化命令除了已知的collect, factor, expand, simplify命令外,还包括如下功能:
combine(trig)——将cos(x)?2+sin(x)?2合并为1。
combing(trig)——将cos(x)?2-sin(x)?2合并为cos(2*x)。
radsimp——将cos(x)+(-sin(x)?2)?(1/2)表示为复数形式cos(x)+i*sin(x)。
convert(exp)——将cos(x)+i*sin(x)表示为指数形式exp(i*x)。
例如,定义函数
f=sym(‘(8/x?3)+12/x?2+6/x+1)?(1/3)’) % 定义函数
显示
f=
(8/x?3+12/x?2+6/x+1)?(1/3)
使用命令
simple(f) % 简化函数
程序会自动给出下列结果:
由于函数f中不含三角函数和指数函数,上述有关三角函数和指数函数的结果都相同。比较各种结果多项式,最终选择最简形式:
ans=
(x+2)/x
如果再作一次简化,键入
simple(ans)
显示
simplify:
(x+2)/x
factor:
(x+2)/x
expand:
1+2/x
collect(x):
1+2/x
在结果的显示中略去了有关三角函数和指数函数的结果。最终选择的最简形式为:
ans=
1+2/x
可以看出,再进行一次简化,结果形式还会有所不同。
1.10.6 可变精度算术运算
计算机内的数值计算精度受到每次计算结果所保留的位数(字长、字节数)的限制。如果保留位数是16位,则第17位以后的数据将被舍去,所以任何数值运算都会引入舍入误差。重复多次的数值运算还会造成累计误差。而MATLAB的符号运算是对符号表达式的运算,结果是非常准确的,因为它们不需要进行数值计算,所以无舍入误差。对符号运算结果用函数eval或numeric求其数值,仅在结果转换时会引入一次性的舍入误差。原理上,符号运算可以实现任何数位的运算,但当保留位数增加时,每次计算就需要增加时间和计算机内存。
Maple的缺省位数为16位精度。命令digits给出保留位数的当前值。命令digits(n)可以改变缺省位点,其中n是所期望精度的位数。用这种方法增加精度的副作用是,每个随后进行的Maple函数的计算都以新的精度为准,增加了计算时间。计算结果的显示不会改变,只有所用的Maple函数的缺省精度受到影响。
另外一个函数vpa可以用缺省的精度或任何指定的精度实行单个符号表达式的计算,以同样的精度来显示结果,而使全局的digits参数不齐。例如,已知π的值:
pi
ans=
3.1416
键入
format long % 用16位表示
显示
pi
ans=
3.14159265358979
键入
digits % 当前符号运算缺省位数
显示
Digits=32
键入
vpa(‘pi’) % 用缺省精度(32位)显示
显示
ans=
3.1415926535897932384626433832795
键入
vpa(‘pi’,20) % 用20位精度显示
显示
ans=
3.1415926535897932385
键入
vpa(‘pi’,50) % 用50位精度显示
显示
ans=
3.1415926535897932384626433832795028841971693993751
当然,还可以选择更多或更少的位数,以得到不同的精度。vpa命令也可以用于符号矩阵,对矩阵的每一个元素均进行上述处理。
1.10.7 符号方程求解
用MATLAB所具有的符号工具箱可以求解符号方程。
1.解单个代数方程
MATLAB用solve命令求解符号方程。如果表达式不是一个方程式(不含等号),则在求解之前自动将表达式置成0。它可以求解f(x)=0或y(x)=f(x)两种形式的代数方程。
如键入命令
solve(‘a*x?2+b*x+c’) % 求二次方程的根
则有
ans=
[1/2/a*(-b+(b?2-4*a*c)?(1/2))]
[1/2/a*(-b-(b?2-4*a*c)?(1/2))]
如键入
solve(‘a*x?2+b*x+c’,‘b’) % 对变量b求解
则有
如键入
solve(‘a*x?2=c+b’) % 解y(x)=f(x)形式的方程
则有
如键入
solve(‘a*x?2=c+b’,‘b’) % 对变量b求解
则有
注意,在求解带有周期函数的议程时,有无穷多个解。这时,solve命令在最接近于0的范围内(如±2π,±π,或±π/2等)搜索,返回一个在最小范围内得到的解。如果不能求得符号解,就返回可变精度数值解。例如键入
digits(5) % 设置5位精度
使用命令
solve(‘log(x)=tan(x)’)
则有
结果返回了数值解,精度为5位。
2.代数方程组求解
solve命令还可以同时求解若干代数方程,语句形式为:
solve(s1,s2,…sn)对缺省变量求解n个方程;
solve(s1,s2,…sn,‘v1,v2…vn’)对n个未知数v1,v2…vn求解n个方程。
例如,给出符号表达式f1,f2,f3如下:
键入
f1=‘x+y+z=2*b’ % 定义f1
显示
f1=
x+y+z=2*b
键入
f2=‘2* x+y+2*z=2*b’ % 定义f2
显示
f2=
2*x+y+2*z=2*b
键入
f3=‘2*x+2*y+z=5*b’ % 定义f3
显示
f3=
2*x+2*y+z=5*b
求解上述3阶代数方程,键入命令
[xx,yy,zz]=solve(f1,f2,f3)
结果为
3.解单个微分方程
用MATLAB符号工具箱的dsovle命令可以求解微分方程。由于用符号运算解微分方程可以得到解析解,比起一般的数值解更具有理论研究意义。
在微分方程的表达式中,包含微分符号。Dsovle命令中用大写字母D来表示求微分。D2,D3等等表示二次、三次重复求微分。例如方程dy2/d2x=0用符号表达式D2y=0来表示。独立变量缺省为t,也可以由用户指定。例如,给出一个一阶微分方程及初始条件:
求解该方程,键入
dsolve(‘Dy=1+y?2’) % 求一阶微分方程的通解
则有
式中,C1为积分常数。
键入
dsolve(‘Dy=1+y?2’,‘y(0) =1’) % 给定初始条件,求特解
则有
键入
dsolve(‘Dy=1+y?2’,‘y(0)=1’,‘t’) % 改变自变量为x
则有
如果求解二阶微分方程,需要给出两个初始条件Dy(0)和y(0),给出二阶导数方程
求解该方程,键入
doslve(‘D2y=cos(t)-y’,‘Dy(0)=0’,‘y(0)=1’) % 解二阶微分方程
则有
ans=
1/2*sin(t)*t+cos(t)
4.解微分方程组
命令dsolve也可同时处理若干个微分方程式。例如,两个线性一阶微分方程:
键入
[x,y]=dsolve(‘Dx=3*x+4*y’,‘Dy=-4*x+3*y’) % 求通解
则有
x=
exp(3*t)*cos(4*t)*Cl+exp(e*t)*sin(4*t)*C2
y=
-exp(3*t)*sin(4*t)*Cl+exp(3*t)*cos(4*t)*C2
结果中的C1,C2是积分常数。若加入初始条件,可求出特解,如键入
[x,y]=dsolve(‘Dx=3*x+4*y’,‘Dy=-4*x+3*y’,‘x(0)=0,y(0)=1’) %加入初始条件
则得
1.10.8 线性代数和矩阵
MATLAB可用线性代数求解符号矩阵问题
1.建立符号矩阵
符号矩阵和向量是数组,其元素为符号表达式,用sym命令产生。建立的符号矩阵,用行向量表示。例如,键入
s=sym(‘[cos(x),a-b;cd,2+x]’)
得到
如果要改变某个元素,可以直接赋值:
s(1,1)=‘sin(x)’ % 改变矩阵元素
得到
sym命令可以将数值矩阵转换为符号矩阵,便于进行符号运算。例如键入
x=[1.1, exp(1); 2.1, sqrt(2)] % 定义数值矩阵
则显示
x=
1.1000 2.7183
2.1000 1.4142
若键入
xx=sym(x) % 转换为符号矩阵
则显示
xx=
[11/10, 6121026514868073*2?(-51)]
[21/10, sqrt(2)]
得到了符号矩阵。由结果可以看出,对于一般常数,如1.1, 2.1,符号常数用11/10,21/10表示;对于无理数exp(1), sqrt(2)等,用符号浮点数或符号表达式表示。
2.矩阵代数运算
用symadd, symsub, symmul和symdiv命令,可以进行符号矩阵的加法、减法、乘法、除法等代数运算,用sympow可计算乘幂,这些都已经作过介绍。用transpose可以给出符号矩阵的转置。例如键入
y=sym(‘[cos(t), sin(t), -sin(t), cos(t)]’) % 定义符号矩阵y
则显示
y=
[cos(t), sin(t)]
[-sin(t), cos(t)]
键入
yt=transpose(y) % 求转置矩阵yt
结果为
键入
yy=symmul(y, yt) % 计算y*yt
结果为
yy=
[cos(t)?2+sin(t)?2 0]
[ 0,cos(t)?2+sin(t)?2]
键入
simplify(yy) % 简化yy
得到
由以上计算结果可知,该矩阵满足y×y′=1,是正交矩阵。
3.线性代数运算
在MATLAB5以上的版本中,符号矩阵运算的命令形式与数值矩阵运算命令相同,即同一命令既可以用于数值矩阵,也可以用于符号矩阵。但得出的解的形式不同,对于数值矩阵,得出数值解;而对于符号矩阵,得出符号表达式。例如,
·inv(a)可求矩阵的逆矩阵。对于前面给出的矩阵y,键入
inv(y)
结果为
·det(a)可计算符号矩阵的行列式。键入
det(y)
结果为
·poly可求解矩阵的特征多项式。键入
poly(y)
结果为
·[v,e]=eig(a)可以求得符号矩阵的特征根和特征向量,但要求矩阵元素是有理数。例如键入
a=sym(‘[2.1; 1.4]’)
显示
a=
[2,1]
[1,4]
键入
[p,e]=eig(a)
得到
解仍然是符号表达式。
·[v,j]=jordan(A)可以得到矩阵的约当(Jordan)标准型。得到的v是A的特征向量矩阵,j是特征值的对角矩阵。例如,对于上面给出的a矩阵,键入
[v,j]=jordan(a)
得到
·svd可求解矩阵奇异值。键入
svd(a)
可得
ans=
[3+2?(1/2)]
[3-2?(1/2)]
·jacobian(f,x)可求解函数f(x)对x的雅可比(Jacobia)矩阵。其结果的第(i,j)项是df(i)/dx(j)。缺省自变量是t。当f是标量时,f的雅可比值是f的梯度。例如键入
b=sym(‘[x1+x2,2-x2;3*x1,2*x1-x2]’) % 定义符号矩阵b
可得
键入
jacobian(b,‘[x1 x2]’) % 求b矩阵对x1, x2的Jacobia矩阵
结果为
结果的第一、二列分别是b对x1,x2的导数。
下面的表1-17~表1-23综合给出了符号工具箱的命令和应用。
表1-17 符号表达式的运算
命令
说明
numeric
pretty
subs
sym
symadd
symdiv
symmul
symop
sympow
symrat
symsub
符号到数值的转换
显示悦目的符号输出
替代子表达式
建立符号矩阵或表达式
符号加法
符号除法
符号乘法
符号运算
符号表达式的幂运算
有理近似
符号减法
表1-18 符号表达式的简化
命令
说明
colleet
xepand
factor
simple
simlify
symsum
合并同类项
展开
因式
求解最简形式
简化
和级数
表1-19 符号多项式
命令
说明
charpoly
horner
numden
poly2sym
sym2poly
特征多项式
嵌套多项式表示
分子或分母的提取
多项式向量到符号的转换
符号到多项式向量的转换
表1-20 符号微积分
命令
说明
diff
int
jordan
taylor
微分
积分
约当标准型
泰勒级数展开
表1-21 符号可变精度算术
命令
说明
digits
vpa
设置可变精度
可变精度计算
表1-22 求解符号方程
命令
说明
compose
dsolve
finverse
linsolve
solve
函数的复合
微分方程的求解
逆函数
齐次线性方程组的求解
代数方程的求解
表1-23 符号线性代数
命令
说明
ploy
det
eig
inv
jordan
linsolve
transose
特征多项式
矩阵行列式的值
特征值和特征向量
逆矩阵
约当标准型
齐次线性方程组的解
矩阵的转置
思考题及习题
1.用MATLAB语言可以求解线性代数中的哪些问题?
2.用符号工具箱能否解高等数学中的定积分?试举例说明。
3.M函数可以用来描述哪一类函数?可以描述动态环节吗?
4.用M函数或M文件描述下面的非线性函数:
5.已知矩阵
试求:(1)a+b, a-b, a+b*5, a-b+I (I为单位阵)
(2)a*b, a*b, 分析结果
(3)a/b, a\b, 分析结果
6.利用for和while循环语句求 当sum>1000时停止运算
7.求多项式x8+x6-2x5+x3+1=0的根,并求其对x的一阶、二阶微分。
8.已经y=x2,在x∈[-10, 10]上等距离取20个点,求出对应的y数组。给每个y的元素上加±0.05(或正或负),用曲线拟合方式求出y1=f(x),将结果与y=x2相比较。
9.给定下面矩阵a ,求a-1矩阵。将aij=i/j代入a和a-1矩阵,验证结果是否正确。
第2章 图形与可视化
2.1 概述
MATLAB具有很强的绘图功能,可以绘制多种二维、三维图形,也可以进行动画演示。本章将详细介绍各种图形的绘制命令。读者可以利用不同的色彩、线条和模块绘制出自己满意的图形。
2.2 二维绘图的plot命令
MATLAB最常用的二维绘图命令是plot命令。该命令将各个数据点用直线连接来绘制图形。MATLAB的其它二维绘图命令中的绝大多数是以plot为基础构造的。Plot命令打开一个默认的图形窗口,它还自动将数值标尺及单位标注加到两个坐标轴上。如果已经存在一个图形窗口,plot命令将刷新当前窗口的图形。
下面举一个简单的例子。
键入命令:
x=0:0.01:2; %图形的横坐标数据准备
y=sin(2*pi*x);%图形的纵坐标数据准备
plot(x, y); %绘制图形
grid
得到图形如图2-1所示,其中横坐标为x,纵坐标为y。
plot指令的调用格式:
plot(x1, y1,‘参数1’,x2, y2, ‘参数2’…)
plot可以用同一命令在同一坐标系中画多幅图形,x1, y1为第一条曲线x,y轴的坐标值,参数1为第一条曲线的先项参数;x2, y2为第二条曲线x, y轴的坐标值,参数2为第二条曲线的选项参数……。具体说明如下:
(1)x, y数据值
x, y可以是向量或矩阵
当x, y均为向量时,要求向量x与向量y的长度一致,则plot(x, y)绘制出以x为横坐标,y为纵坐标的二维图形。
当x为向量,y为矩阵时,plot(x, y)用不同颜色的图线绘制出y行或列对于x的图形。y矩阵的行或列的选择取决于x, y的维数。若y为方阵或y矩阵的列向量长度为x向量的长度一致,则绘制出y矩阵的各个列向量相对于x的一组二维图形;若y矩阵的行向量长度与x向量的长度一致,则绘制出y矩阵的各个行向量相对于x的一组二维图形。
若x为矩阵,y为向量,则按类似于(2)的规则。
若x, y是同维的距阵,则plot(x, y)绘制出y列向量相对于x的列向量之间的一组二维图形。
若x为向量,则plot(x)绘制出一个x元素排列序号之间关系的线性坐标图。
若x为矩阵,则plot(x)绘制出x的列向量相对于行号的二维图形。
(2)参数
参数选项为一个字符串,它决定了二维图形的颜色、线型及数据点的图标。表2-1~表2-3分别给出颜色、线型和标记的控制字符。
表2-1 颜色控制符
字符
颜色
字符
颜色
b
蓝色
m
紫红色
c
青色
r
红色
g
绿色
w
白色
k
黑色
y
黄色
表2-2 线型控制符
符号
线型
符号
线形
-
实线(默认)
:
点连线
-˙
点划线
--
虚线
表2-3 数据点标记字符
控制符
标记符
控制符
标记符
·
点
h
六角形
+
十字号
p
五角星
。
圆圈
^
下三角
*
星号
ˇ
上三角
×
叉号
>
右三角
s
正方形
<
左三角
d
菱形
如:plot(x, y, ‘r:o’)命令,字符串‘r:o’中,第一个字符“r”表示曲线颜色为红色;第二个字符“:”表示曲线线型采用点连线;“o”表示曲线上每一个数据点处用圆圈标出。当参数只指定数据点标记时,只按照标记字符画出孤立的数据点,不将数据点连接成线。
例如,键入命令
t=0:0.5:7;x=sin(t);plot(t, x, ’k:o’)得到的结果如图2-2所示。
2.3 图形修饰与控制
MATALB提供了一系列图形修饰函数,用于对plot命令绘制的图形进行修饰和控制。
2.3.1 坐标轴的调整
MATLAB用axis命令对绘制的图形的坐标轴进行调整。axis命令的功能非常丰富,可用它来控制轴的比例和特性。
1.坐标轴比例控制
命令:axis[xmin xmax ymin ymax]将图形的x轴范围限定在[xmin,xmax]之间,y轴的范围限定在[ymin,ymax]之间。MATLAB绘制图形时,按照给定的数据值确定坐标轴参数范围。对坐标轴范围参数的修改,也就相当于对原图形进行放大或缩小处理。
例如:在绘出图2-2之后,加一条命令:axis[0 3*pi -2 2],则图2-2变为图形2-3。
图2-3中,横轴为t,纵轴为y,整个图形的取值范围按照命令要求设置。
2.坐标轴特性控制
命令:axis(控制字符串)根据表2-4所示的功能控制图形。
表2-4 axis控制字符
控制字符
命令功能
auto
自动设置坐标系(默认):xmin=min(x); xmax=max(x)
ymin=min(x);ymax=max(y)
square
将图形设置为正方形图形
equal
将图形的x, y坐标轴的单位刻单设置为相等
mormal
关闭axis(square)和axis(equal)命令的作用
xy
使用笛卡儿坐标系
ij
使用“matrix”坐标系。即:坐标原点在左上方,x坐标从左向右增大,y坐标从上向下增大
on
打开所有轴标注、标记和背景
off
关闭所有轴标注、标记和背景
3.坐标刻度标示
命令:set (gca, ‘xtick’,标示向量)
set (gca, ‘ytick’,标示向量)
按照标示向量设置x, y轴的刻度标示。
命令:set(gca, ‘xticklabel’, ‘字符串∣字符串…’)
set(gca, ‘yticklabel’, ‘字符串∣字符串…’)
按照字符串设置x, y轴的刻度标注。
例如键入命令:
t=0:0.05:7;
plot(t, sin(t))
set(gca, ‘xtick’,[0 1.4 3.14 5 6.28]) %给定轴的标注点
得到结果如图2-4。图2-4中,横轴根据给定的数组标注刻度。
如果键入命令:
set (gca, ‘xticklabel’, ‘0∣1.4∣half∣5∣one’) %给定轴的标注点
得到结果如图2-4。图2-4中,横轴根据给定的数组标注刻度。
如果键入命令:
set(gca, ‘xticklabel’, ‘0∣1.4∣half∣5∣one’) %改变x轴的标注点
则横轴用字符串标注,结果如图2-5。
2.3.2 文字标示
有关图形的标题、轴线标注等指令有:
title(‘字符串’)——图形标题
xlable(‘字符串’)——x轴标注
ylable(‘字符串’)——y轴标注
text(x, y, ‘字符串’)——在坐标(x, y)处标注说明文字
gtext(‘字符串’)——用鼠标在特定处标注说明文字
输入特定的文字需要用反斜杠(\)开头,用法如表2-5所列。
表2-5 特殊字符
输入字符
表示的特殊字符
\pi
π
\alpha
α
\beta
β
\leftarrow
←
\rightarrow
→
\bullet
例如键入命令
t=0:0.05:2*pi;
plot(t, sin(t))
set(gca, ‘xtick’, [0 1.4 3.14 5 6.28])
xlabel(‘t(deg)’) %x轴单位
ylabel(‘magnitude(V)’) %y轴单位
title(‘This is a example 0\rightarrow 2\pi’) %图题
text(3.14, sin(3.14), ‘\leftarrow this is zero for\pi’) %在图中加文字
grid
结果如图2-6所示。
2.3.3 网格控制
设置或取消网格需要使用网格控制命令。网格是在坐标刻度标示上画出的格线。画出网格,便于对曲线进行观察和分析。
命令:grid on ——在所画的图形中添加网格线;
命令:grid off ——在所画的图形中去掉风格线。
也可以只输入命令grid添加网格线,再一次输入命令grid,则取消网格线。
图2-6中的网格就是执行grid命令后的结果。
2.3.4 图例注解
当在一个坐标系上画有多幅图形时,为区分各个图形,MATLAB提供了图例的注解说明命令。其格式如下:
legend(字符串1,字符串2,…,参数)
此命令在图形中开启一个注解视窗,依据绘图的先后顺序,依次输出字符串对各个图形进行注解说明。如字符串1表示第一个出现的线条,字符串2表示第二个出现的线条,参数字符串确定注解视察在图形中的位置,其含义见表2-6。同时,注解视窗可以用鼠标拖动,以便将其放置在一个合适的位置。
表2-6 参数字符串的含义
参数字符串
含义
0
尽量不与数据冲突,自动放置在最佳位置
1
放置在图形的右上角(默认)
2
放置在图形的左上角
3
放置在图形的左下角
4
放置在图形的右下角
-1
放置在图形视窗外右边
例如键入命令
x=0: 0.2: 12;
plot(x, sin(x), ‘-’, x, 1.5*cos(x), ‘:’);
legend(‘First’, ‘Second’);
结果如图2-7所示。
在图2-7的基础上,执行命令:
legend(‘First’, ‘Second’, -1) %强行将注解视察放置在图形视窗外的右方结果如图2-8所示。
2.3.5 图形的保持
hold命令用于保持当前图形。用plot命令绘图时,首先将当前图形窗口清屏,再绘制图形,所以我们只能见到最后一个plot命令绘制的图形。为了能利用多余plot命令绘制多幅图形,就需要保持窗口上的图形。
hold on——保持当前图形及轴系的所有特性。
Hold off——解除hold off命令。
例如键入命令:
x=0: .2: 12;
plot(x, sin(x), ‘-’)
hold on
plot(x, 1.5 * cos(x), ‘:’);
上例使用了两条plot命令绘制出两条曲线,如图2-9所示。
2.3.6 图形窗口的分割
MATLAB的绘图指令可以将绘图窗口分割成几个区域,在各个区域中分别绘图。
命令:subplot(m, n, p)将当前绘图窗口分割成m行n列区域,并指定第p个编号区别为以前绘图区域。区域的编号原则是“先上后下,先左后右”。MATLAB允许每个编号区域可以以不同的坐标系单独绘图。m, n和p前面的逗号可以省略。
例如键入命令:
x=0:0.05:7;
y1=sin(x);
y2=1.5*cos(x);
y3=sin(2*x);
y4=5*cos(2*x);
subplot(221)
plot(x, y1);
title(‘sin(x)’)
subplot(222)
plot(x, y2);
title(‘cos(x)’)
subplot(223)
plot(‘sin(2x)’)
xubplot(224)
plot(x, y4);
title(‘cos(2x)’)
2.3.7 图形的填充
fill命令用于填充二维封闭多边形。
命令:fill(x, y, ‘color’)在由数据x, y所构成的多边形内,用color所指定的颜色填充。如果该多边形不是封闭的,可以用初始点和终点的连线封闭。Color控制符见表2-1的说明。
例如键入命令
x=0:0.05:7;
y=sin(x);
fill(x, y, ‘k’);
结果如图2-11所示。可以看到,由于该图形不是封闭的,MATLAB用初试点(0,0)和终点(7,sin(7))的连线将其封闭,并填充黑色。
2.4 特殊坐标二维图形
MATLAB提供一些特殊坐标二维图形函数,如semilogx, semilogy, loglog和polar命令。这些命令与plot命令完全类似,不同的是将数据绘制到不同的图形坐标上。
2.4.1 对数坐标图形
semilogx(x, y, 参数)命令绘制半对数坐标图形,其x轴取以10为底的对数坐标,y轴为线性坐标。
Semilogy(x, y, 参数)绘制半对数坐标图形,其y轴取以10为低的对数坐标,x轴为线性坐标。
Loglog(x, y, 参数)绘制x, y轴都取以10为底的对数坐标图形。
例如键入命令
x=0:0.05:2;
y=10^x;
subplot(211), semilogx(x, y)
title(‘semilogx Graph’)
grid
subplot(212), semilogy(x, y)
title(‘semilogy Graph’)
grid
得到的结果如图2-12。
如果接着键入命令
subplot(211), olglog(x, y)
grid
title(‘loglog Graph’)
subplot(212), plot(x, y)
title(‘linear Graph’)
grid
得到的结果如图2-13。
2.4.2 极坐标图形
polar(theta, radius, 参数)命令绘制相角为theta,半径为radius的极坐标图形。
例如键入命令
t=0:0.01:2*pi;
r=2*cos(2*(t-pi/8));
polar(t, r)
结果如图2-14所示。
2.5 特殊二维图形
2.5.1 函数图形
fplot(‘函数运算式’,[xmin xmax])命令用来绘制给定函数在区间[xmin xmax]内的变化图形。
例如,要画出y=sin(3x),其中x在0到4之间变化的图形,只要执行如下命令:
fplot(‘sin(3*x)’, [0 4])
grid
即可得到如图2-5的结果。
2.5.2 饼图
饼图在统计中常用来表示各因素所占的百分比。
命令:pie(x)或pie(x,explode)根据矩阵或向量x绘制饼图,以表示各数据占sum(x)的百分比。若x为向量,该命令绘制x的每一元素占全部向量元素总和值的百分比的饼图;若x为矩阵,该命令绘制x的每一元素占全部矩阵元素总和值的百分比的饼图;参数explode表示某元素对应的扇形图是否从整个饼图中分离出来,若非零,则分离出来,它的维数应与x相同。
例如,键人命令
x=[15 35 10 15 25];
pie(x,[1 0 1 0 0])
运行结果如图2—16所示。该图形分为5块,按照命令,第l块和第3块分离出来。
2.5.3 条形图
条形图用来表示一些数据的对比情况。MATLAB提供了两类条形图的命令:一类是垂直方向的条形图;另一类是水平方向的条形图。
1. 垂直方向的条形图
命令:bar(x,width)或bar(x,‘参数’)根据矩阵或向量x绘制条形图。若x为向量,则以其各元素的序号为各个数据点的横坐标,以J向量的各个元素为纵坐标,绘出一个垂直方向的条形图。若x为矩阵,同时参数字符串为group或缺省,则以其各列序号为横坐标,每一列在其列序号坐标上分别以列的各元素为纵坐标,绘出一组垂直方向的条形图;若参数字符串为stack,则以其各列序号为横坐标,每一列在其列序号坐标上以列向量的累加值为纵坐标,绘出一个垂直方向的分组式条形图。Width给定条形的宽度,缺省值为0.8,若width大于1,则条形图重叠。
例如键人
x=[10,20,30;15,35,10;5,20,25];
subplot(121),bar(x,‘group’)
subplot(122),bar(x,‘stack’)
结果如图2—17。
2. 水平方向的条形图
水平方向的条形图命令形式与垂直方向条形图命令相同。命令为barh(x,width)或barh (x,‘参数,)。这类命令与上述bar命令的使用方法相同,只不过这类命令绘制的是水平方向的条形图。
2.5.4 梯形图
梯形图可以用来表示系统中的采样数据。
命令形式为
stairs(x)或stairs(x,y)
其中,x,y均为向量。stairs(x)命令绘制以 x的向量序号为横坐标,以x向量的各个对应元素为纵坐标的梯形图。stairs(x,y)命令绘制以向量x的各个对应元素为横坐标,以y向量的各个对应元素为纵坐标的梯形图。
例如键人命令:
x=0:0.1:7;
y=sin(x);
stairs(x,y)
结果如图2—18所示。
2.5.5 概率分布图
研究随机系统时,常常要用到概率分布图。MATLAB提供hist命令来绘制概率分布图。
命令:Hist(y,x)绘制y在以x为中心的区间中分布的个数条形图。
例如,键人命令:
y=randn(1,1000);
x=一2:0.1:2;
hist(y,x)
结果如图2—19所示。
2.5.6 向量图
1.原点向量图
命令compass(x)绘制相对原点的向量图,若工为复数,compass(x)命令相当于compass
(real(x),imag(x))。
命令compass(x,y)以复数坐标系的原点为起点,绘制出带箭头的一组复数向量,其中向量x表示复数的实部,向量y表示复数的虚部。
例如键人命令:
x=[一2+3*j3+4*j1-5*j-22];
subplot(121),compass(x)
y=[1,—2,3,—4,5,6];
z=[—2 3 —6 —5 —5 0];
subplot(122),compass(y,z)
结果如图2—20所示。
2.水平线向量图
命令形式为
feather(x)或feather(x,y)
feather命令与compass命令的功能相似,两者的区别是起点不同。compass命令起始于坐标原点,feather命令起始于x向量各元素的序号。
例如键人命令
x=[一2+3*j 3+4*i 1—5*j 一2—3*j 5—2*j];
subplot(121),feather(x)
grid
y=[l,—2,3,—4,5,6];
z=[—2 3 —6 —5 —5 1];
subplot(122),feather(y,z)
grid
结果如图2—21所示。
2.6 三维图形
MATLAB提供了各种各样的显示三维图形的函数。有些函数可在三维空间中画线,有些是画曲面或线格框架。下面首先讨论基本三维图形。
2.6.1 基本三维图形
MATLAB提供plot3命令绘制三维曲线。该命令将绘制二维图形的函数plot的特性扩展到三维空间。命令形式为:
plot3(x1,y1,z1,参数1,x2,y2,z2,参数2,…)
其功能和使用方法类似于绘制二维图形的plot命令。
例如键人命令:
t;0:pi/50;10*pi;
plot3(sin(t),cos(t),t);
grid
结果如图2—22所示。
2.6.2 三维图形的修饰与控制
三维图形的修饰与控制方式与二维图形类似。三维图形的控制命令与二维的相同;图形
的修饰命令中,增加了对z轴的修饰与标记。如zlabel(‘字符串’)命令。
例如键人命令:
t:0:pi/50;10‘pi;
plot3(sin(t),cos(t),t);
grid
axis([——2 2——2 2 0 40]);
t=0:pi/50:10‘pi;
plot3(sin(t)+1,cos(t)一1,t);
hold On
plot3(sin(2x t)一1,cos(2x t)+1,t);
grid
hold off
结果如图2—23所示。
2.6.3 特殊三维图形
1.三维网格曲面
三维网格曲面是连接三维空间的一些四边形所构成的曲面。在介绍三维网格曲面命令之 前,首先介绍产生三维网格数据点的函数meshgrid。
命令[X,Y]=meshgrid(x,y)将向量x(1×m),y(1×n)转换为三维网格数据矩阵X(n×m),Y(n×m)。
例如键入命令:
[X,Y]=meshgrid([1 2 3 4],[5 6 7])
执行后的结果为:
X=
1 2 3 4
1 2 3 4
1 2 3 4
Y=
5 5 5 5
6 6 6 6
7 7 7 7
三维网格曲面命令为mesh(x,y,z,c),mesh(x,y,z),mesh(z,c)或mesh(z)。
这四种命令格式都可以绘制三维网络曲线。当x(n×m),y(n×m)为矩阵时,且x矩阵的所有行向量相同,y矩阵的所有列向量相同时,mesh命令将自动执行meshgrid(x,y),将x,y转换为三维网格数据矩阵。z和c分别为(m×n)维矩阵,c表示网络曲面的颜色分布,若省略,则网络曲面的颜色与z方向上的高度值成正比。若x,y均省略,则三维网格数据矩阵取值x=1:n,y=1:m。
例如键入命令:
x=-10:0.5:10;
y=-8:0.5:8;
[X,Y]=meshgrid(x,y);
Z=sin(sqrt(X.?2+Y.?2))./sqrt(X.?2+Y.?2);
Mesh(X,Y,Z)
结果如图2-24所示。
2.带等高线的三维网络曲面
命令meshc(x,y,z,c),meshc(x,y,z),meshc(z,c)或meshc(z)绘制带有等高线(XY平面)的三维网格曲面。这些命令类似mesh命令,不同的是该命令还在XY平面上绘制曲面Z轴方向的等高线。
对于上述例题,如果接着键入命令:
meshc(X,Y,Z)
结果如图2-25所示。
3.带底座的三维网络曲面
命令meshz(x,y,z,c),meshz(x,y,z),meshz(z,c)或meshz(z)绘制带有底座的三维网格曲面。这些命令类似mesh命令,不同的是该命令还在XY平面上绘制曲面的底座。
例如上面例题,接着键入命令
meshz(X,Y,Z)
结果如图2-26所示。
4.填充颜色的三维网格曲面
命令mesh绘制的是连接三维空间的一些四边形所构成的曲面。该曲面只有四边形的边用某种颜色绘出,四边形的内部是透明的。Surf命令绘制的曲面也由一些四边形所构成,不同的是四边形的边是黑色的,其内部用不同的颜色填充。
命令形式为surf(x,y,z,c),surf(x,y,z),surf(z,c)或surf(z)。
这些命令使用方法类似mesh命令。
例如上例给出数据,键入
surf(X,Y,Z)
结果如图2-27。
5.三维直方图
与二维图形类似,MATLAB提供了两类绘制三维直方图的命令:
bar3(y,z,‘参数’)——绘制垂直的三维直方图;
bar3h(y,z,‘参数’)——绘制水平的三维直方图。
这些命令绘制以x=1:n为x坐标;以y(1×m)向量的各个元素为y坐标(若y省略,则y=1:m),且要求y是递增或递减的;以矩阵z(n×m)为z坐标,绘出垂直的或水平的三维直方图。参数可以选择grouped(分组式的)、detached(分离式的)或stacked(累加式的),与二维绘图命令的参数选择相同。
例如键入命令:
x=[10,20,30;15,35,10;5,20,25];
subplot(121), bar3(x,’grouped’)
subplot(122), bar3h(x)