1.5 MATLAB的数值运算
MATLAB具有强大的数值运算能力,它不仅能对矩阵和向量进行相应的运算,而且也可处理多项式的解、数据分析、函数的极值、线性方程组的解、函数的微积分和函数绘图等问题。
1.5.1 矩阵运算
MATLAB的基本数据单元是不需要指定维数的复数矩阵,它提供了各种矩阵的运算与操作。它既可以对矩阵进行整体的处理,也可以对矩阵的某个或某些元素进行单独的处理,所以在MATLAB环境下矩阵的操作同数的操作一样简单。
1.矩阵的实现
在MATLAB语言中不必描述矩阵的维数和类型,它们是由输入的格式和内容来确定的。例如,
当时,把A当做一个2 ×2维的矩阵;
当A=[12]时,把A当做一个2维行向量;
当时,把A当做一个2维列向量;
当A=5时,把A当做一个标量;
当A=1+2j时,把A当做一个复数。
(1)矩阵的赋值
在MATLAB中,矩阵可以用以下几种方式进行赋值。
●直接列出元素的形式;
●通过命令和函数产生;
●建立在文件中;
●从外部的数据文件中装入。
①简单矩阵的输入
对于比较小的简单矩阵可以使用直接排列的形式输入,把矩阵的元素直接排列到方括号中,每行内的元素间用空格或逗号分开,行与行的内容用分号隔开。例如,矩阵
在MATLAB下的输入方式为
>>A=[1,2,3;4,5,6;7,8,9]
或
>>A=[123;456;789]
都将得到相同的结果:
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矩阵将一直保存在工作空间中,除非被替代或清除,在MATLAB的命令窗口中可随时输入矩阵名查看其内容。
②利用语句或函数产生矩阵
在MATLAB中,向量可利用下面的语句来产生:
s1:s2:s3
其中,s1为起始值,s3为终止值,s2为步矩。
使用这样的命令就可以产生一个由s1开始,以步距s2自增,并终止于s3的行向量。例如
>>y=[0:pi/4:pi;0:10/4:10]
结果显示:
y= 0 0.7854 1.5708 2.3562 3.1416 0 2.5000 5.0000 7.5000 10.0000
如果s2省略,则认为自增步距为1。例如,利用以下命令可产生一个行向量:
>>x=1:5
结果显示:
x= 1 2 3 4 5
利用linspace()函数也可产生一个行向量,该函数的调用格式为
x=linspace(n,m,z)
其中,x为产生的等间隔行向量;参数n和m分别为行向量中的起始和终止元素的值;z为需要产生行向量的元素数。例如,对于以下命令:
>>x=linspace(0,2*pi,5)
结果显示:
x= 0 1.5708 3.1416 4.7124 6.2832
利用size()函数可测取一个矩阵的维数,该函数的调用格式为
[n,m] =size(A)
其中,A为要测试的矩阵名;而返回的两个参数n和m分别为A矩阵的行数和列数。
当要测试的变量是一个向量时,仍可由size()函数来得出其大小;更简洁地,用户可以使用length()函数来求出,该函数的调用格式为
n=length(x)
其中,x为要测试的向量名;而返回的n值为向量x的元素个数。
如果对一个矩阵A用length(A)函数测试,则返回该矩阵行列的最大值,即该函数等效于max(size(A))。
(2)矩阵的元素
MATLAB的矩阵元素可用任何表达式来描述,它既可以是实数,也可以是复数。例如
>>B=[ -1/3 1.3; sqrt(3)(1+2+3)*j]
结果显示:
B= -0.3333 1.3000 1.7321 0+6.0000i
MATLAB允许把矩阵作为元素来建立新的矩阵。例如
>>A=[1 2 3;4 5 6;7 8 9];C=[A;[10 11 12]]
结果显示:
C= 1 2 3 4 5 6 7 8 9 10 11 12
MATLAB还允许对一个矩阵的单个元素进行赋值和操作。例如,如果想将A矩阵的第2行第3列的元素赋值为100,则可通过下面的语句来完成:
>>A=[1 2 3;4 5 6;7 8 9];A(2,3)=100
这时将只改变此元素的值,而不影响其他元素的值。
如果给出的行数或列数大于原来矩阵的范围,则MATLAB将自动扩展原来的矩阵,并将扩展后未赋值的矩阵元素置为0。例如,想把以上矩阵A的第4行第5列元素的值定义为8,可以通过下面语句来完成:
>>A(4,5)=8
结果显示:
A= 1 2 3 0 0 4 5 100 0 0 7 8 9 0 0 0 0 0 0 8
利用上面的语句除了对单个矩阵元素进行定义之外,MATLAB还允许对子矩阵进行定义和处理。例如:
A(1:3,1:2:5) %表示取A矩阵的第1行到第3行内,且位于第1,3,5 列上的所有元素构成的子矩阵;
A(2:3,:) %表示取A矩阵的第2行和第3行的所有元素构成的子矩阵; A(:,j) %表示取A矩阵第j列的全部元素构成的子矩阵; B(:,[3,5,10])=A(:,1:3) %表示将A矩阵的前3列,赋值给B矩阵的第3,5和第10列; A(:,n:-1:1) %表示由A矩阵中取n至1反增长的列元素组成一个新的 矩阵。
特别地,当A(:)在赋值语句的右边时,表示将A的所有元素按列在一个长的列向量中展成串。例如:
>>A=[12;34],B=A(:)
结果显示:
B= 1 3 2 4
(3)特殊矩阵的实现
在MATLAB中特殊矩阵可以利用以下函数来建立。
①单位矩阵函数eye()
基本格式:A=eye(n) %产生一个n阶的单位矩阵A A=eye(size(B))%产生与B矩阵同阶的单位矩阵A A=eye(n,m) %产生一个主对角线的元素为1,其余元素全为0的n ×m矩阵
②零矩阵函数zeros()
基本格式:A=zeros(n,m) %产生一个n ×m零矩阵A A=zeros(n) %产生一个n ×n零矩阵A A=zeros(size(B))%产生一个与B矩阵同阶的零矩阵A
③1矩阵函数ones()
基本格式:A=ones(n,m) A=ones(n) A=ones(size(B))
④随机元素矩阵函数rand()
随机元素矩阵的各个元素是随机产生的。如果矩阵的随机元素满足[0,1]区间上的均匀分布,则可以由MATLAB函数rand()来生成。该函数的调用格式为
A=rand(n,m) A=rand(n) A=rand(size(B))
⑤对角矩阵函数diag()
如果用MATLAB提供的方法建立一个向量V=[a1,a2,…,an],则可利用函数diag(V)来建立一个对角矩阵。例如:
>>V=[1,2,3,4];A=diag(V)
如果矩阵A为一个方阵,则调用V=diag(A)将提取出A矩阵的对角元素来构成向量V,而不管矩阵的非对角元素是何值。
⑥伴随矩阵函数compan()
假设有一个多项式
p(s)=sn+a1 sn-1 +a2 sn-2 +…+an-1 s+an
则可写出一个伴随矩阵
生成伴随矩阵的函数的调用格式为
A=compan(p)
其中,p=[1,a1,a2,…,an]为一个多项式向量。
例如有一个向量p=[1 2 3 4 5],则可通过下面的命令构成一个伴随矩阵。
>>p=[1 2 3 4 5];A=compan(p)
⑦上三角矩阵函数triu()和下三角矩阵函数tril()
调用格式为
A=triu(B)和A=tril(B)
其中,B为一个矩阵。例如:
>>B=[1 2 3;4 5 6;7 8 9];A=tril(B)
2.矩阵的基本运算
矩阵运算是MATLAB的基础,MATLAB的矩阵运算功能十分强大,并且运算的形式和一般的数学表示十分相似。
(1)矩阵的转置
矩阵转置的运算符为“′”。例如:
>>A=[1 2 3;4 5 6];B=A′
如果A为复数矩阵,则A′为它的复数共轭转置,非共轭转置使用A.′,或者用conj(A)实现。
(2)矩阵的加和减
矩阵的加减法的运算符为“+”和“-”。只有同阶矩阵方可进行加减运算,标量可以和矩阵进行加减运算,但应对矩阵的每个元素施加运算。例如:
>>A=[1 2 3;4 5 6;7 8 9];B=A+1
(3)矩阵的乘法
矩阵乘法的运算符为“*”。当两个矩阵中前一矩阵的列数和后一矩阵的行数相同时,可以进行乘法运算,这与数学上的形式是一致的。例如:
>>A=[1 2 3;4 5 6;7 8 9];B=[1 1;1 1;1 1];C=A*B
在MATLAB中还可将矩阵和标量进行相乘,其结果为标量与矩阵中的每个元素分别相乘。
(4)矩阵的除法
矩阵的除法有两种运算符“\”和“/”,分别表示左除和右除。一般地讲,x=A\B是A*x=B的解,x=B/A是x*A=B的解。通常A\B≠B/A,而A\B=inv(A)*B,B/A=B*inv(A)。
(5)矩阵的乘方
矩阵乘方的运算符为“^”。一个方阵的乘方运算可以用A^P来表示。P为正整数,则A的P次幂即为A矩阵自乘P次。如果P为负整数,则可以将A自乘P次,然后对结果进行求逆运算,就可得出该乘方结果。如果P是一个分数,如P=m\n,其中n和m均为整数,则首先应该将A矩阵自乘n次,然后对结果再开m次方。例如:
>>A=[1 2 3;4 5 6;7 8 9];B=A^2,C=A^0.1
(6)矩阵的翻转
MATLAB还提供了一些矩阵翻转处理的特殊命令,对n ×m维矩阵A,如
B=fliplr(A) %命令将矩阵A进行左右翻转再赋给矩阵B,即bij =ai,m+1-j C=flipud(A) %命令将矩阵A进行上下翻转再赋给矩阵C,即cij =an+1-i,j D=rot90(A) %命令将矩阵A进行旋转90度后赋给矩阵D,即dij =aj,m+1-i
例如:
>>A=[1 2 3;4 5 6;7 8 9];B=fliplr(A),C=flipud(A)
(7)矩阵的超越函数
MATLAB中exp(),sqrt(),sin(),cos()等基本函数命令可以直接使用在矩阵上,这种运算只定义在矩阵的单个元素上,即分别对矩阵的每个元素进行运算。超越数学函数,可以在函数后加上m而成为矩阵的超越函数,例如expm(A),sqrtm(A),logm(A)分别为矩阵指数、矩阵开方和矩阵对数。矩阵的超越函数要求运算的矩阵必须为方阵。例如
>>A=[1 2 3;4 5 6;7 8 9];B=expm(A),C=sqrtm(A)
3.矩阵的特殊运算
矩阵的特殊运算包括以下内容。
(1)矩阵行列式
矩阵A={aij}的行列式定义为
|A |=det(A)∑(-1)ka1k1a2k2…ankn
式中,k1,k2,…,kn是将序列1,2,3,…,n的元素交换k次所得出的一个序列,每个这样的序列称为一个置换,而∑表示对k1,k2,…,kn取遍1,2,3,…,n所有的排列的求和。
MATLAB求矩阵行列式函数的调用格式为
det(A)
计算矩阵的行列式有多种算法,在MATLAB中采用的方法为LU分解法。
(2)矩阵求逆
对于一个已知的n ×n维非奇异方阵A来说,如果有一个同样大小的C矩阵满足
AC =CA=I
式中,I为单位阵,则称C矩阵为A矩阵的逆矩阵,并记为C=A-1。
矩阵求逆的算法是多种多样的,比较常用的有行(列)主元素Gauss消去法、LU分解法和基于奇异值分解的方法等,MATLAB提供了一个求取逆矩阵的函数inv(),其调用格式为
inv(A)
如果A矩阵为奇异的或接近奇异,则利用此函数有可能产生错误的结果。
(3)矩阵的迹
假设一个方阵为A={aij},(i,j=1,2,…,n),则矩阵A的迹定义为
亦即矩阵的迹为该矩阵对角线上各个元素之和。由代数理论可知矩阵的迹和该矩阵的特征值之和是相同的。在MATLAB中提供了求取矩阵迹的函数trace(),其调用格式为
trace(A)
(4)矩阵的秩
对于n ×m维的矩阵A,若矩阵所有的列向量中共有rc 个线性无关,则称矩阵的列秩为rc;如果rc =m,则称A为列满秩矩阵。相应地,若矩阵A的行向量中有rr个是线性无关的,则称矩阵A的行秩为rr;如果rr=n,则称A为行满秩矩阵。可以证明,矩阵的行秩和列秩是相等的,记
rank{A} =rc =rr
这时矩阵的秩为rank{A}。矩阵的秩也表示该矩阵中行列式不等于0的子式的最大阶次。所谓子式,即为从原矩阵中任取k行及k列所构成的子矩阵。
矩阵求秩的算法也是多种多样的,有的算法是稳定的,而有的算法可能因矩阵的条件数变化而不是很稳定。MATLAB中采用的算法是基于矩阵的奇异值分解的算法,首先对矩阵做奇异值分解,得出矩阵A的n个奇异值σi(i=1,2,…,n),在这n个奇异值中找出大于给定误差限ε的个数r,这时r就可以认为是矩阵A的秩。
MATLAB提供了一个内部函数rank()来用数值方法求取一个已知矩阵的秩,其调用格式为
k=rank(A,tol)
其中,A为要求秩的矩阵;k为所求矩阵A的秩;而tol=ε为判0用的误差限,一般可以取默认值eps。这样调用格式就可以简化为
k=rank(A)
这里的判0用误差限就取做机器中的默认值eps。如果eps取的不合适,则求出的数值秩可能和原矩阵的秩不同。所以在使用数值秩时应当引起注意。
(5)矩阵的三角分解
矩阵的三角分解又称为LU分解,它的目的是将一个矩阵A分解成一个下三角矩阵L和一个上三角矩阵U的乘积,亦即可以写成A=LU,其中L和U矩阵分别可以写成
这样产生的矩阵与原来的矩阵A的关系可以写成
a11 =u11,a12 =u12,…,a1n =u1n
a21 =l21 u11,a22 =l21 u12 +u22,…,a2n =l21 u1n +u2n
︙
其中
由上式可以立即得出求lij和uij的递推计算公式
该公式的递推初值为u1i =a1i,i=1,2,…,n。
注意,在上述的算法中并未对主元素进行任何选取,因此该算法并不一定数值稳定。
在MATLAB下也给出了矩阵的LU分解函数lu(),该函数的调用格式为
[L,U] =lu(A)
其中,L和U分别为变换后的下三角矩阵和上三角矩阵。在MATLAB的lu()函数中考虑了主元素选取的问题,所以该函数一般会给出可靠的结果。由该函数得出的下三角矩阵L并不一定是一个真正的下三角矩阵。因为选取它可能进行了一些元素行的交换,这样主对角线的元素可能不是1。例如:
>>A=[1 2 3;4 5 6;7 8 0];[L,U] =lu(A)
LU分解使用的算法是高斯变量消去法,这种分解是矩阵求逆、矩阵求行列式和线性方程求解的基础,也是方阵“/”或“\”两种矩阵除法的基础。
(6)矩阵的奇异值分解
矩阵的奇异值也可以看成是矩阵的一种测度,对任意的n ×m矩阵A来说,总有ATA≥0,AAT≥0,且有
rank{ATA} =rank{AAT} =rank(A)
进一步可以证明,ATA与AAT有相同的非零特征值λi,且相同的非零特征值总是为正数。在数学上把这些非零的特征值的平方根称做矩阵A的奇异值,记为
矩阵的奇异值大小通常决定矩阵的性态。如果矩阵的奇异值变化特别大,则矩阵中某个参数有一个微小的变化将严重影响到原矩阵的参数,如其特征值的大小,这样的矩阵又称为病态矩阵,有时也称为奇异矩阵。
假设矩阵A为n ×m矩阵,且rank{A} =r,则矩阵A可以分解为
其中,L和M为正交矩阵;Δ=diag{σ1,…,σr}为对角矩阵,且其对角元素均不为0。
MATLAB提供了直接求取矩阵奇异值分解的函数,其调用格式为
[U,S,V] =svd(A)
其中,A为原始矩阵;S为对角矩阵,其对角元素就是A的奇异值;而U和V均为正交矩阵,并满足
A=USVT
矩阵最大奇异值σmax和最小奇异值σmin的比值又称为该矩阵的条件数,记为cond{A},即cond{A} =σmax/σmin。矩阵的条件数越大,则对参数变化越敏感。在MATLAB下也提供了求取矩阵条件数的函数cond(),其调用格式为
cond(A)
(7)矩阵的范数
矩阵的范数也是对矩阵的一种测度。在介绍矩阵的范数之前,首先要介绍向量范数的基本概念。
如果对线性空间中的一个向量x存在一个函数ρ(x)满足下面三个条件:
① ρ(x)≤0且ρ(x)=0的充要条件是x=0;
② ρ(ax)=|a|ρ(x),a为任意标量;
③对向量x和y有ρ(x+y)≤ρ(x)+ρ(y)。
则称ρ(x)为x向量范数。
范数的形式是多种多样的。可以证明,下面给出的一族式子都满足上述的三个条件:
且
这里用到了向量范数的记号‖x‖p。
对于任意的非零向量x,矩阵A的范数为
和向量的范数一样,矩阵的范数定义如下:
其中,s{A}为矩阵A的特征值,而smax{ATA}即为矩阵A的最大奇异值的平方。换句话说,‖A‖2 为矩阵A的最大奇异值。
MATLAB提供了求取矩阵范数的函数norm(),它允许求各种意义下的矩阵范数,该函数的调用格式为
N=norm(A,选项)
其中,选项如表1-8所示。
表1-8 矩阵范数函数的选项定义
(8)矩阵的特征值与特征向量
对一个矩阵A来说,如果存在一个非零的向量x,且有一个标量λ满足
Ax =λx
则称λ为矩阵A的一个特征值,而x为对应于特征值λ的特征向量。严格说来,x应该称为A的右特征向量。如果矩阵A的特征值不包含重复的值,则对应的各个
特征向量为线性独立的,这样由各个特征向量可以构成一个非奇异的矩阵,如果用它对原始矩阵做相似变换,则可以得出一个对角矩阵。
矩阵的特征值与特征向量可以由MATLAB提供的函数eig()容易地求出,该函数的调用格式为
[V,D] =eig(A)
其中,A为要处理的矩阵;D为一个对角矩阵,其对角线上的元素为矩阵A的特征值,而每个特征值对应的V矩阵的列为该特征值的特征向量。该矩阵是一个满秩矩阵,它满足AV=VD,且每个特征向量各元素的平方和(即2范数)均为1。如果调用该函数时只给出一个返回变量,则将只返回矩阵A的特征值。即使A为复数矩阵,也同样可以由eig()函数得出其特征值与特征向量矩阵。
(9)矩阵的特征多项式、特征方程和特征根
对于给定的n ×n阶矩阵A,称多项式
f(s)=det(sI-A)=a0 sn+a1 sn-1 +…+an-1 s+an
为矩阵A的特征多项式,其中系数a0,a1,…,an称为矩阵的特征多项式系数。
MATLAB提供了求取矩阵特征多项式系数的函数poly(),其调用格式为
p=poly(A)
其中,A为给定的矩阵;返回值p为一个行向量,其各个分量为矩阵A的降幂排列的特征多项式系数。即p=[a0 a1 … an]。
令特征多项式等于零所构成的方程称为该矩阵的特征方程,而特征方程的根称为该矩阵的特征根。MATLAB中根据矩阵特征多项式求特征根的函数为roots(),其调用格式为
V=roots(p)
其中,p为特征多项式的系数向量,而V为特征多项式的解,即原始矩阵的特征根。
【例1-11】 求的特征多项式、特征根及其模值。
解 MATLAB命令如下
>>A=[1 2 3;4 5 6;7 8 9];p=poly(A),r=roots(p)′,abs(r)
结果显示:
p= 1.0000 -15.0000 -18.0000 -0.0000 r = 16.1168 -1.1168 -0.0000 ans = 16.1168 1.1168 0.0000
1.5.2 向量运算
虽然在MATLAB中向量和矩阵在形式上有很多的一致性,但它们实际上遵循着不同的运算规则。MATLAB向量运算符由矩阵运算符前面加一点“.”来表示,如“.*”、“./”和“.^”等。
1.向量的加减
向量的加减运算与矩阵的运算相同,所以“+”和“-”既可被向量接受又可被矩阵接受。
2.向量的乘法
向量乘法的操作符为“.*”。如果x和y两向量具有相同的维数,则x.*y表示向量x和y单个对应元素之间的相乘。例如:
>>x=[1 2 3];y=[4 5 6];z=x.*y
结果显示:
z= 4 10 18
可见向量的输入和输出与矩阵具有相同的格式,但它们的运算规则不同。例如,如果x是一个向量,则求取向量x平方时不能直接写成x*x,而必须写成x.*x,否则将给出错误信息。
但是对于矩阵可以使用向量运算符号,这时实际上就相当于把矩阵看成了向量来进行运算。例如,对于两个维数相同矩阵A和B,C=A.*B表示A和B矩阵的对应元素之间直接进行乘法运算,然后将结果赋给C矩阵,把这种运算称为矩阵的点积运算。两个矩阵之间的点积是它们对应元素的直接运算,它与矩阵的乘法是不同的。例如:
>>A=[1 2 3;4 5 6;7 8 9];B=[2 3 4;5 6 7;8 9 0];C=A.*B
结果显示:
C= 2 6 12 20 30 42 56 72 0
3.向量的除法
向量除法的操作符为“./”或“.\”。它们的运算结果一样。例如:
>>x=[123];y=[456];z=y./x
结果显示:
z= 4.0000 2.5000 2.0000
对于向量除法运算x.\y和y./x一样,将得到相同的结果,这与矩阵的左、右除是不一样的,因向量的运算是它们对应元素间的运算。
对于矩阵也可使用向量的除法操作符,这时就相当于把矩阵看成向量进行运算。
4.向量的乘方
向量乘方的运算符为“.^”。向量的乘方是对应元素的乘方,在这种底与指数均为向量的情况下,要求它们的维数必须相同。例如:
>>x=[1 2 3];y=[4 5 6];z=x.^y
结果显示:
z= 1 32 729
它相当于:z=[1 2 3].^[4 5 6] =[14 25 36]。
若指数为标量时,会得到如下结果。例如:
>>x=[1 2 3];z=x.^2
结果显示:
z= 1 4 9
以上运算相当于:z=[1 2 3].^2=[12 22 32]。
若底为标量时,则会得到如下结果。例如:
>>x=[1 2 3];y=[4 5 6];z=2.^[x y]
结果显示:
z= 2 4 8 16 32 64
以上运算相当于:z=2.^[1 2 3 4 5 6] =[21 22 23 24 25 26]。
同样对于矩阵也可以采用运算符“.^”。例如:
>>A=[1 2 3;4 5 6;7 8 0];B=A.^A
结果显示:
B= 1 4 27 256 3125 46656 823543 16777216 1
即矩阵B中的每个元素都是矩阵A元素的相应乘方。例如,3125=5^5。
可见如果对矩阵使用向量运算符号,实际上就相当于把矩阵看成了向量的运算,否则将作为矩阵运算。
1.5.3 关系和逻辑运算
1.关系运算
MATLAB常用的关系操作符如表1-9所示。
表1-9 关系操作符
MATLAB的关系操作符可以用来比较两个大小相同的矩阵,或者比较一个矩阵和一个标量。比较两个元素大小时,结果是1表明为真,结果是0表明为假。函数find()在关系运算中很有用,它可以在矩阵中找出一些满足一定关系的数据元素。例如:
>>A=1:9;B=A>4
结果显示:
B= 0 0 0 0 1 1 1 1 1 >>A=1:9;C=A(A>4)
或
>>A=1:9;C=find(A>4)
结果显示:
C= 5 6 7 8 9
2.逻辑运算
MATLAB的逻辑操作符有 &(与)、|(或)和~(非)。它们通常用于元素或0-1矩阵的逻辑运算。
与运算符和或运算符可比较两个标量或两个同阶矩阵。对于矩阵,逻辑运算符是作用于矩阵中的元素。逻辑运算结果信息也用“0”和“1”表示,逻辑操作符认定任何非零元素都表示为真。给出1为真,0为假。
非是一个元操作符,当A非零时,~A返回的信息为0;当A为零时,~A返回的信息为1。因而就有:P|(~P)返回值为1,P&(~P)返回值为0。例如:
>>A=1:9;C=~(A>4)
结果显示:
C= 1 1 1 1 0 0 0 0 0 >>A=1:9;C=(A>4)&(A<7)
结果显示:
C= 0 0 0 0 1 1 0 0 0
3.关系和逻辑运算函数
除了上面介绍的关系和逻辑运算符外,MATLAB中还提供了一些关系和逻辑运算函数,如表1-10所示。
表1-10 关系和逻辑运算函数
对于矩阵,any()和all()命令按列对其进行处理,并返回带有处理列所得结果的一个行向量。
1.5.4 多项式运算
多项式运算是数学中最基本的运算之一。在MATLAB中同样可对多项式进行相应的一系列运算。
1.多项式的表示
在高等数学中,多项式一般可表示成以下形式:
f(x)=a0 xn+a1 xn-1 +…+an-1 x+an
其中,a0,a1,…,an称为多项式的系数。
所以多项式很容易用其系数组成的行向量p=[a0 a1 … an]来表示,其中行向量是按其系数降幂排列组成的系数向量。在MATLAB中,构造多项式正是采用了把多项式的各项系数依降幂次序排放在行向量的对应元素位置,直接输入其系数向量的方法来实现的。对于缺项的系数一定要进行补零。例如,对于多项式
f(x)=x4 +5x3 +3x+2
可用以下MATLAB命令来表示:
>>p=[1 5 0 3 2]
在MATLAB中,利用函数poly2str()可将多项式的系数向量表示成相应多项式的习惯表示形式,该函数的调用格式为
f=poly2str(p,′s′)
其中,p为多项式的系数向量,s为多项式的变量名,f为相应的多项式。例如:
>>p=[1 5 0 3 2];f=poly2str(p,′x′)
结果显示:
f= x^4+5 x^3+3x+2
2.多项式的四则运算
多项式的四则运算是指多项式的加、减、乘和除运算。其中多项式的加、减运算要求两个相加、减的多项式的系数向量维数的大小必须相等。
(1)多项式的加减
在MATLAB中,当两个相加、减的多项式阶次不同时,低阶多项式的系数向量必须用首零填补,使其与高阶多项式的系数向量有相同维数。
【例1-12】 求以下两个多项式的和。
f1(x)=x4 +5x3 +3x+2,f2(x)=x2 +6x+5
解 MATLAB命令如下:
>>p1=[1 5 0 3 2];p2=[0 0 1 6 5];p=p1+p2
结果显示:
p= 1 5 1 9 7
(2)多项式的乘法
在MATLAB中,多项式的乘法运算,利用函数conv()来实现。函数conv()相当于执行两个数组的卷积,其调用格式为
p=conv(p1,p2)
其中,p1和p2为多项式的系数按降幂排列构成的系数向量;p为多项式p1和p2的乘积多项式,按其系数降幂排列构成的多项式积的系数向量。
【例1-13】 求以下两个多项式的乘积。
f1(x)=x4 +5x3 +3x+2,f2(x)=x2 +6x+5
解 MATLAB命令如下:
>>p1=[1 5 0 3 2];p2=[1 6 5];p=conv(p1,p2)
结果显示:
p= 1 11 35 28 20 27 10
需要说明的是,当对多个多项式执行乘法运算时,可重复使用conv()函数。
【例1-14】 求多项式f(x)=(x+1)2(x2 +6x+5)的展开式。
解 MATLAB命令如下:
>>p=conv([1 1],conv([1 1],[1 6 5]))
结果显示:
p = 1 8 18 16 5
(3)多项式的除法
在MATLAB中,多项式的除法运算是利用函数deconv()来实现的,其调用格式为
[p,r] =deconv(p1,p2)
其中,p1,p2为多项式的系数按降幂排列构成的系数向量;p为多项式p1被p2除的商多项式,按其系数降幂排列构成的多项式商的系数向量,而余多项式为r。函数deconv()相当于两个数组的解卷运算,使p1=conv(p2,p)+r成立。
3.多项式的值及其导数
如果多项式函数为
f(x)=a0 xn+a1 x(n-1)+…+an-1 x+an
则可以求出该函数的导数函数为
f′(x)=na0 xn-1 +(n-1)a1 xn-2 +…+an-1
在MATLAB中提供了多项式求值函数polyval()和多项式求导函数polyder(),它们的调用格式分别为
f0=polyval(p,x0)和dp=polyder(p)
其中,p为由多项式系数按降幂排列构成的系数向量;x0为求值点的x值。该函数将返回f(x)多项式在x=x0的值f0。函数polyder(p)返回多项式导数的系数向量,亦即向量
dp=[na0(n-1)a1 … an-1]
同样,MATLAB也提供了多项式矩阵的求值函数polyvalm(),其调用格式为
fA=polyvalm(p,A)
其中,p为矩阵多项式函数按降幂排列构成的向量,即
p=[a0 a1 … an]
而A为一个给定矩阵,返回值fA为下面的矩阵多项式的值:
f(A)=a0An+a1An-1 +…+an-1A+anI
4.多项式的求解
MATLAB中多项式的求解运算同样可利用前面介绍的roots()函数来实现,其调用格式为
r=roots(p)
其中,p为多项式的系数向量;r为多项式的解。
【例1-15】 求方程f(x)=x2 +5x+6=0的解。
解 MATLAB命令如下:
>>p=[1 5 6];x=roots(p)
结果显示:
x= -3.0000
-2.0000
5.多项式的拟合
在MATLAB中,多项式的拟合可用polyfit()函数来实现,其调用格式为
[p,s] =polyfit(x,y,n)
其中,x,y为利用最小二乘法进行拟合的数据;n为要拟合的多项式的阶次;p为要拟合的多项式的系数向量;s为使用该函数获得的错误预估计值。一般来说,多项式拟合的阶数越高,拟合的精度就越高。
【例1-16】 用6阶多项式对[0,2*pi]区间上的sint函数进行拟合。
解 MATLAB命令如下:
>>x=linspace(0,2*pi,80);y=sin(x);p=polyfit(x,y,6),y1=polyval(p,x); >>plot(x,y,′ro′,x,y1)
拟合结果如图1-10所示。
图1-10 正弦曲线及6阶拟合曲线
6.多项式的插值
在MATLAB中,多项式的插值方法包括一维线性插值、二维线性插值、三维线性插值和三次样条插值等。
(1)一维线性插值
对于一维线性数据,可以通过插值或查表来求得离散点之间的数据值。在MATLAB中,一维线性插值可用函数interp1()来实现,其调用格式为
yi=interp1(x,y,xi,′method′)
其中,返回值yi为在插值向量xi处的函数值向量,它是根据向量x与y插值而来的。如果y是一个矩阵,那么对y的每一列进行插值,返回的矩阵yi的大小为length(xi)×length(y,2)。如果xi中有元素不在x的范围内,则与之相对应的yi返回NaN。如果x省略,表示x=1:n,此处n为向量y的长度或为矩阵y的行数,即size(y,1)。参数′method′表示插值方法,它可采用以下方法:最近插值(′nearest′)、线性插值(′linear′)、三次插值(′cubic′)和三次样条插值(′spline′)。′method′的默认值为线性插值。
需要注意,上述函数interp1()的所有插值调用格式,都要求向量x为单调。当x为单调且等间距时,可以使用快速插值法,此时可将′method′参数的值设置为′nearest′,′linear′或′cubic′。例如:
>>x=linspace(0,2*pi,80);y=sin(x); >>x1=[0.5 1.4 2.6 4.2],y1=interp1(x,y,x1,′cubic′)
结果显示:
x1= 0.5000 1.4000 2.6000 4.2000 y1= 0.4794 0.9855 0.5155 -0.8716
(2)二维线性插值和三维线性插值
与一维线性插值一样,二维线性插值也可以通过插值或查表来求得离散点之间的数据值。在MATLAB中,二维线性插值可用函数interp2()来实现,其调用格式为
zi=interp2(x,y,z,xi,yi,′method′)
其中,返回值zi为在插值向量xi,yi处的函数值向量,它是根据向量x,y与z插值而来的。x,y和z也可以是矩阵。如果xi,yi中有元素不在x,y的范围内,则与之相对应的zi返回NaN。如果x,y省略,表示x=1:m,y=1:n,此处n与m为矩阵z的行数和列数,即[n,m] =size(z)。参数′method′的定义同上。
需要注意,上述的所有插值调用格式,都要求向量x与y为单调且具有相同的格式。当x与y为单调且等间距时,同样可以使用快速插值法。例如:
>>z=[3 5 7 9 11 10 4 15;1 3 2 6 11 5 7 13]; >>z1=interp2(z,[1,3],[1,2])
结果显示:
z1= 3 2
与一维线性插值和二维线性插值一样,三维线性插值也可以通过插值或查表来求得离散点之间的数据值。在MATLAB中,三维线性插值可用函数interp3()来实现,其使用方法与in-terp1()函数和interp2()函数基本类似,这里就不详细介绍了。
(3)三次样条插值
使用高阶多项式的插值常常会产生病变的结果,而三次样条插值就能消除这种病变。在三次样条插值中,要寻找三次多项式,以逼近每对数据点之间的曲线。在MATLAB中,三次样条插值利用spline()函数来实现,其调用格式为
yi=spline(x,y,xi)
其中,返回值yi为在插值向量xi处的函数值向量,它是根据向量x与y插值而来的。此函数的作用等同于interp1(x,y,xi,′spline′)。
1.5.5 数据分析
MATLAB强大的数组运算功能,决定了它很容易对一大批数据进行一般的数据分析,如求数组的极值、平均值、中值、和、积、标准差、方差、协方差和排序等。
1.随机数
在MATLAB中提供了两个用来产生随机数组的函数rand()和randn()。其调用格式为
A=rand(n,m)和A=randn(n,m)
其中,n和m分别为将要产生随机数组的行和列;A为n ×m维的随机数组。函数rand(n,m)用来产生一个n ×m维在[0,1]区间上均匀分布的随机数组;函数randn(n,m)用来产生一个n ×m维的均值为0、标准差为1的正态分布的随机数组。
2.最大值和最小值
如果给定一组数据{xi},i=1,2,…,n,则可利用MATLAB将这些数据用一个向量表示出来,即
x=[x1,x2,…,xn]
利用MATLAB中的函数max()和min()便可求出这组数据的最大值和最小值,调用格式为
[xmax,i] =max(x)和 [xmin,i] =min(x)
其中,返回的xmax及xmin分别为向量x的最大值和最小值,而i为最大值或最小值所在的位置。当然这两个函数均可以只返回一个参数,而不返回i。如果给出的x不是向量而是矩阵,则采用min()和max()函数得出的结果将不是数值,而是一个向量。它的含义是得出由每一列构成的向量的最大值或最小值所构成的行向量。当x为复数时,通过计算max(abs(x))返回结果。
对于函数max()和min()也可以采用以下格式:
z=max(x,y)和z=min(x,y)
其中,x与y为向量或矩阵,它们的大小一样;返回结果z是一个包含最大或最小元素,且与它们大小一样的向量或矩阵。
3.平均值
利用MATLAB中的mean()函数可求出向量或矩阵的平均值,调用格式为
y=mean(x)
其中,x为向量或矩阵;y为向量或矩阵x的平均值。若x为向量,则返回向量元素的平均值。若x为矩阵,则返回一行向量,包括矩阵每列元素的平均值。
4.中值
利用MATLAB中的median()函数可求出向量或矩阵的中值,调用格式为
y=median(x)
其中,x为向量或矩阵;y为向量或矩阵x的中值。若x为向量,则返回向量元素的中值;若x为矩阵,则返回一行向量,包括矩阵每列元素的中值。
5.求和
利用MATLAB中的sum()函数可求出向量或矩阵中元素的和,调用格式为
y=sum(x)
其中,x为向量或矩阵;y为向量或矩阵x中元素的和。若x为向量,则返回向量元素的总和;若x为矩阵,则返回一行向量,包括矩阵每列元素的和。
6.求积
利用MATLAB中的prod()函数可求出向量或矩阵中元素的积,调用格式为
y=prod(x)
其中,x为向量或矩阵;y为向量或矩阵x中元素的积。若x为向量,则返回向量元素的积;若x为矩阵,则返回一行向量,包括矩阵每列元素的积。
7.标准差
利用MATLAB中的std()函数可求出向量或矩阵中元素的标准差,调用格式为
y=std(x)
其中,x为向量或矩阵;y为向量或矩阵x中元素的标准差。若x为向量,则返回向量元素的标准差;若x为矩阵,则返回一行向量,包括矩阵每列元素的标准差。
8.方差
利用MATLAB中的var()函数可求出向量或矩阵中元素的方差,调用格式为
y=var(x)
其中,x为向量或矩阵;y为向量或矩阵x中元素的方差。若x为向量,则返回向量元素的方差;若x为矩阵,则返回一行向量,包括矩阵每列元素的方差。
9.协方差
利用MATLAB中的cov()函数可求出向量或矩阵中元素的协方差,调用格式为
y=cov(x)
其中,x为矩阵;y为矩阵x中各列间的协方差阵。函数cov(x,y)相当于
cov(x(:),y(:))
10.按实部或幅值对特征值进行排序
利用MATLAB中的esort()和dsort()函数,可对特征值按实部或幅值进行排序,函数的调用格式为
[s,ndx] =esort(p)和 [s,ndx] =dsort(p)
其中,esort(p)针对连续系统,根据实部按递减顺序对向量p中的复特征值进行排序;ndx为索引矢量,对于连续特征值,先列出不稳定特征值。dsort(p)针对离散系统,根据幅值按递减顺序对向量p中的复特征值进行排序;ndx为索引矢量。对于离散特征值,先列出不稳定特征值。
另外,sort()函数用于对元素按升序进行排序。
1.5.6 函数极值
在MATLAB中,提供了两个基于单纯形算法求解多元函数极小值的函数fmins()(仅适用于MATLAB 6.5及以前的版本)和fminsearch(),以及一个基于拟牛顿法求解多元函数极小值的函数fminunc(),其调用格式分别为
x=fmins(fun,x0,options) x=fminsearch(fun,x0,options) x=fminunc(fun,x0,options)
其中,fun表示函数名;x0表示函数的初值,其大小往往能决定最后解的精度和收敛速度;x为所求极小值点;options表示选项,它是由一些控制变量构成的向量,比如它的第一个分量不为0表示在求解时显示整个动态过程(其默认值为0),第二个分量表示求解的精度(默认值为1e-4)。可以指定这些参数来控制求解的条件。在调用以上函数时,首先应该写一个描述f(·)的函数,其格式为
y=fun(x)
【例1-17】 求函数f(x)=3x-x3的最小值。
解 首先根据方程编写一个函数ex1 -17.m。
%ex1 17.m function y=ex1 17(x) y=3*x-x^3;
然后利用下面的命令调用fminsearch()函数求方程的解。
>>x=fminsearch(′ex1 17′,0),ymin=ex1 17(x)
运行结果:
x= -1.0000 ymin= -2
结果表明,函数在x=-1时有最小值ymin=-2。
以上函数的极值问题,也可直接利用以下命令得到与上面相同的结果。
>>f=′3*x-x^3′,x=fminsearch(f,0),ymin=3*x-x^3
因为函数f(x)的极小值问题等价于函数-f(x)的极大值问题,所以函数fminsearch()也可用来求解函数 f(x)的极大值,这时只需在所求函数前加一负号即可。例如,求函数f(x)=3x-x3 的最大值,可利用以下命令:
>>f=′-(3*x-x^3)′,x=fminsearch(f,0),ymax=3*x-x^3
运行结果:
x = 1.0000 ymax= 2
结果表明,函数在x=1时有最大值ymax=2。
1.5.7 代数方程求解
利用MATLAB中求函数f(.)零点的函数fzero()和fsolve(),可以方便地求得非线性方程f(·)=0的解,它们的调用格式分别为
x=fzero(fun,x0)和x=fsolve(fun,x0)
其中,fun表示函数名,函数名定义同前;x0表示函数的初值,是为求解过程所假设的起始点,当x0为标量时,该命令将在它两侧寻找一个与之最靠近的解,当x0为二元向量[a,b]时,该命令将在区间[a,b]内寻找一个解;x为返回的解。fzero()用来对一元方程求解,fsolve()用来对多元方程求解。当然利用函数fzero()和fsolve()也可对线性方程进行求解。
【例1-18】 试求函数f=sinx的零点。
解 MATLAB命令如下:
>>y=fzero(′sin(x)′,[0,pi])
结果显示:
y= 0
【例1-19】 试求以下非线性方程组的解。
解 首先根据三元方程编写一个函数ex1-19.m。
%ex1 19.m function q=ex1 19(p) q(1)=sin(p(1))+p(2)^2+log(p(3))-7; q(2)=3*p(1)+2*p(2)-p(3)^3+1; q(3)=p(1)+p(2)+p(3)-5;
然后利用下面的命令在初值x0 =1,y0 =1,z0 =1下调用fsolve()函数直接求出方程的解。
>>x=fsolve(′ex1 19′,[1 1 1])
结果显示:
x= 0.5991 2.3959 2.0050
【例1-20】 试求以下线性方程组的解。
解 首先根据两元方程编写一个函数ex1-20.m。
%ex1 20.m function z=ex1 20(p) z(1)=p(1)-p(2); z(2)=p(1)+p(2)-6;
然后利用下面的命令在任意给的初值x0 =0,y0 =0下调用fsolve()函数直接求出方程的解。
>>x=fsolve(′ex1 20′,[0 0])
结果显示:
x= 3.0000 3.0000
特别指出,对线性方程Ax=B,在A的逆存在的条件下,有更简单的求解方法。在MAT-LAB下对该线性方程,可利用以下命令进行求解。
x=inv(A)*B
例如,对于例1-20也可利用以下命令,得到上述结果。
>>A=[1 -1;1 1];B=[0;6];x=inv(A)*B
1.5.8 微分方程求解
MATLAB中提供了求解常微分方程的函数ode45(),其调用格式为
x=ode45(fun,[t0,tf],x0,tol)
其中,fun为函数名,其定义同前;[t0,tf]为求解时间区间;x0为微分方程的初值;tol用来指定误差精度,其默认值为10 -3;x为返回的解。
【例1-21】 求下列微分方程在初始条件x(0)=1,下的解。
解 首先将微分方程写成一阶微分方程组。
令x1 =x,,则可得
然后根据以上微分方程组编写一个函数ex1-21.m。
%ex1 21.m function dx=ex1 21(t,x) dx=[x(2);(1-x(1)^2)*x(2)-x(1)];
最后利用以下的MATLAB命令,即可求出微分方程在时间区间[0,30]上的解曲线,如图1-11所示。
图1-11 解曲线
>>[t,x] =ode45(′ex1 21′,[0,30],[1;0]); >>plot(t,x(:,1),t,x(:,2));xlabel(′t′);ylabel(′x(t)′)
1.5.9 函数积分
对于函数f(x)的定积分
在被积函数f(x)相当复杂时,一般很难采用解析的方法求出定积分的值来,而往往要采用数值方法来求解。求解定积分的数值方法是多种多样的,如简单的梯形法、Simpson法和Rom-berg法等,它们的基本思想都是将整个积分空间[a,b]分割成若干个子空间[xi,xi+1],i=1,2,…,n。其中x1 =a,xn+1 =b。这样整个积分问题就分解为下面的求和形式:
而在每一个小的子空间上都可以近似地求解出来。例如可以采用下面给出的Simpson方法来求解出[xi,xi+1]上积分近似值:
式中,hi =xi+1 -xi。
MATLAB基于此算法,采用自适应变步长方法用quad()函数来求取定积分,该函数的调用格式为
y=quad(fun,a,b,tol)
其中,fun为函数名,其定义和其他函数一致;a,b分别为定积分的上下限;tol为变步长用的误差限,如不给出误差限,则将自动地假定tol=1 ×10 -3;y为返回的值。
【例1-22】 试求积分的解。
解 这一无穷定积分可以由有限积分来近似,一般情况下,选择积分的上下限为±15就能保证相当的精度。
首先根据给定的被积函数编写下面的函数ex1-22.m。
%ex1 22.m function f=ex1 22(x) f=1/sqrt(2*pi)*exp(-x.^2/2);
然后通过下面的MATLAB语句可求出所需函数的定积分。
>>format long;y=quad(′ex1 22′,-15,15)
结果显示:
y= 1.00000007247356
以上函数的积分问题,也可直接利用以下命令得到与上面相同的结果。
>>format long;y=quad(′1/sqrt(2*pi)*exp(-x.^2/2)′,-15,15)
或
>>format long;f=′1/sqrt(2*pi)*exp(-x.^2/2)′;y=quad(f,-15,15)
另外,MATLAB给出一种利用插值运算来更精确更快速地求出所需要的定积分的函数quadl(),该函数的调用格式为
y=quadl(fun,a,b,tol)
其中,tol的默认值为10 -6,其他参数的定义及使用方法和quad()几乎一致。该函数可以更准确地求出积分的值,且一般情况下函数调用的步数明显小于quad()。
对于上例,使用quadl()函数可得如下结果:
>>format long;f=′1/sqrt(2*pi)*exp(-x.^2/2)′;y=quadl(f,-15,15)
结果显示:
y= 1.00000000000338
在二维情况下求积分实质上是求函数曲线与坐标轴之间所夹的封闭图形的面积,故利用MATLAB中的trapz()函数,也可求取定积分。例如,对上例有
>>format long;x=-15:0.01:15;y=ex1 21(x);area=trapz(x,y)
结果显示:
area= 1.00000000000000