第四章 MATLAB的数值计算功能.docx
《第四章 MATLAB的数值计算功能.docx》由会员分享,可在线阅读,更多相关《第四章 MATLAB的数值计算功能.docx(42页珍藏版)》请在冰豆网上搜索。
第四章MATLAB的数值计算功能
第四章MATLAB的数值计算功能
Chapter4:
NumericalcomputationofMATLAB
数值计算是MATLAB最基本、最重要的功能,是MATLAB最具代表性的特点。
MATLAB在数值计算过程中以数组和矩阵为基础。
数组是MATLAB运算中的重要数据组织形式。
前面章节对数组、矩阵的特征及其创建与基本运算规则等相关知识已作了较详尽的介绍,本章重点介绍常用的数值计算方法。
一、多项式(Polynomial)`
多项式在众多学科的计算中具有重要的作用,许多方程和定理都是多项式的形式。
MATLAB提供了标准多项式运算的函数,如多项式的求根、求值和微分,还提供了一些用于更高级运算的函数,如曲线拟合和多项式展开等。
1.多项式的表达与创建(ExpressionandCreatingofpolynomial)
(1)多项式的表达(expressionofpolynomial)_
Matlab用行矢量表达多项式系数(Coefficient),各元素按变量的降幂顺序排列,如多项式为:
P(x)=a0xn+a1xn-1+a2xn-2…an-1x+an
则其系数矢量(Vectorofcoefficient)为:
P=[a0a1…an-1an]
如将根矢量(Vectorofroot)表示为:
ar=[ar1ar2…arn]
则根矢量与系数矢量之间关系为:
(x-ar1)(x-ar2)…(x-arn)=a0xn+a1xn-1+a2xn-2…an-1x+an
(2)多项式的创建(polynomialcreating)
a)系数矢量的直接输入法
利用poly2sym函数直接输入多项式的系数矢量,就可方便的建立符号形式的多项式。
例1:
创建多项式x3-4x2+3x+2
poly2sym([1-432])
ans=
x^3-4*x^2+3*x+2
POLYConvertrootstopolynomial.
POLY(A),whenAisanNbyNmatrix,isarowvectorwith
N+1elementswhicharethecoefficientsofthecharacteristicpolynomial,DET(lambda*EYE(SIZE(A))-A).
POLY(V),whenVisavector,isavectorwhoseelementsare
thecoefficientsofthepolynomialwhoserootsarethe
elementsofV.Forvectors,ROOTSandPOLYareinverse
functionsofeachother,uptoordering,scaling,and
roundofferror.
b)由根矢量创建多项式
已知根矢量ar,通过调用函数p=poly(ar)产生多项式的系数矢量,再利用poly2sym函数就可方便的建立符号形式的多项式。
例2:
由根矢量创建多项式。
将多项式(x-6)(x-3)(x-8)表示为系数形式的多项式。
a=[638]%根矢量
pa=poly(a)%求系数矢量
ppa=poly2sym(pa)%以符号形式表示原多项式
ezplot(ppa,[-50,50])
pa=
1-1790-144
ppa=
x^3-17*x^2+90*x-144
说明:
(1)根矢量元素为n,则多项式系数矢量元素为n+1;
(2)函数poly2sym(pa)把多项式系数矢量表达成符号形式的多项式,缺省情况下自变量符号为x,可以指定自变量。
(3)使用简单绘图函数ezplot可以直接绘制符号形式多项式的曲线。
例3:
由给定复数根矢量求多项式系数矢量。
r=[-0.5-0.3+0.4i-0.3-0.4i];
p=poly(r)
pr=real(p)
ppr=poly2sym(pr)
p=
1.00001.10000.55000.1250
pr=
1.00001.10000.55000.1250
ppr=
x^3+11/10*x^2+11/20*x+1/8
说明:
含复数根的根矢量所创建的多项式要注意:
(1)要形成实系数多项式,根矢量中的复数根必须共轭成对;
(2)含复数根的根矢量所创建的多项式系数矢量中,可能带有很小的虚部,此时可采用取实部的命令(real)把虚部滤掉。
如果需要进行系数表示形式的多项式的求根运算,有两种方法可以实现,一是直接调用求根函数roots,poly和roots互为逆函数。
另一种是先把多项式转化为伴随矩阵,然后再求其特征值,该特征值即是多项式的根。
c)特征多项式输入法
用poly函数还可实现由矩阵的特征多项式系数创建多项式。
条件:
特征多项式系数矢量的第一个元素必须为一。
例2:
求三阶方阵A的特征多项式系数,并转换为多项式形式。
a=[638;756;135]
Pa=poly(a)%求矩阵的特征多项式系数矢量
Ppa=poly2sym(pa)
Pa=
1.0000-16.000038.0000-83.0000
Ppa=
x^3-17*x^2+90*x-144
注:
n阶方阵的特征多项式系数矢量一定是n+1阶的。
例4:
将多项式的系数表示形式转换为根表现形式。
求x3-6x2-72x-27的根
a=[1-6-72-27]
r=roots(a)
r=
12.1229
-5.7345
-0.3884
MATLAB约定,多项式系数矢量用行矢量表示,根矢量用列矢量表示。
>>
2.多项式的乘除运算(Multiplicationanddivisionofpolynomial)
向量的卷积与解卷积对应着多项式的乘除法,多项式乘法(卷积)用函数conv(a,b)实现,除法(解卷积)用函数deconv(a,b)实现。
长度为m的向量a和长度为n的向量b的卷积定义为:
C(k)=
C向量的长度为:
m+n-1
解卷积是卷积的逆运算,向量a对向量c进行解卷积将得到商向量q和余量r,并且满足:
例1:
a(s)=s2+2s+3,b(s)=4s2+5s+6,计算a(s)与b(s)的乘积。
a=[123];b=[456];%建立系数矢量
c=conv(a,b)
cs=poly2sym(c,’s’)%建立指定变量为s的符号形式多项式
c=
413282718
cs=
4*s^4+13*s^3+28*s^2+27*s+18
例2:
展开(s2+2s+2)(s+4)(s+1)(多个多项式相乘)
c=conv([1,2,2],conv([1,4],[1,1]))
cs=poly2sym(c,’s’)%(指定变量为s)
c=
1716188
cs=
s^4+7*s^3+16*s^2+18*s+8
例2:
求多项式s^4+7*s^3+16*s^2+18*s+8分别被(s+4),(s+3)除后的结果。
c=[1716188];
[q1,r1]=deconv(c,[1,4])%q—商矢量,r—余数矢量
[q2,r2]=deconv(c,[1,3])
cc=conv(q2,[1,3])%对除(s+3)结果检验
test=((c-r2)==cc)
q1=
1342
r1=
00000
q2=
1446
r2=
0000-10
cc=
17161818
test=
11111
3.其他常用的多项式运算命令(Othercomputationcommandofpolynomial)
pa=polyval(p,s)按数组运算规则计算给定s时多项式p的值。
pm=polyvalm(p,s)按矩阵运算规则计算给定s时多项式p的值。
[r,p,k]=residue(b,a)部分分式展开,b,a分别是分子(numerator)分母(denominator)多项式系数矢量,r,p,k分别是留数、极点和直项矢量
p=polyfit(x,y,n)用n阶多项式拟合x,y矢量给定的数据。
polyder(p)多项式微分。
注:
对于多项式b(s)与不重根的n阶多项式a(s)之比,其部分分式展开为:
式中:
p1,p2,…,pn称为极点(poles),r1,r2,…,rn称为留数(residues),k(s)称为直项(directterms),假如a(s)含有m重根pj,则相应部分应写成:
RESIDUEPartial-fractionexpansion(residues).
[R,P,K]=RESIDUE(B,A)findstheresidues,polesanddirecttermofapartialfractionexpansionoftheratiooftwopolynomialsB(s)/A(s).Iftherearenomultipleroots,
B(s)R
(1)R
(2)R(n)
----=--------+--------+...+--------+K(s)
A(s)s-P
(1)s-P
(2)s-P(n)
VectorsBandAspecifythecoefficientsofthenumeratoranddenominatorpolynomialsindescendingpowersofs.Theresidues
arereturnedinthecolumnvectorR,thepolelocationsincolumnvectorP,andthedirecttermsinrowvectorK.Thenumberofpolesisn=length(A)-1=length(R)=length(P).Thedirecttermcoefficientvectorisemptyiflength(B)length(K)=length(B)-length(A)+1.
IfP(j)=...=P(j+m-1)isapoleofmultplicitym,thentheexpansionincludestermsoftheform
R(j)R(j+1)R(j+m-1)
--------+------------+...+------------
s-P(j)(s-P(j))^2(s-P(j))^m
[B,A]=RESIDUE(R,P,K),with3inputargumentsand2outputarguments,convertsthepartialfractionexpansionbacktothepolynomialswithcoefficientsinBandA.
例3:
对(3x4+2x3+5x2+4x+6)/(x5+3x4+4x3+2x2+7x+2)做部分分式展开
a=[134272];
b=[32546];
[r,s,k]=residue(b,a)
r=
1.1274+1.1513i
1.1274-1.1513i
-0.0232-0.0722i
-0.0232+0.0722i
0.7916
s=
-1.7680+1.2673i
-1.7680-1.2673i
0.4176+1.1130i
0.4176-1.1130i
-0.2991
k=
[](分母阶数高于分子阶数时,k将是空矩阵,表示无此项)
例5:
对一组实验数据进行多项式最小二乘拟合(leastsquarefit)
x=[12345];%实验数据
y=[5.543.1128290.7498.4];
p=polyfit(x,y,3)%做三阶多项式拟合
x2=1:
.1:
5;
y2=polyval(p,x2);%根据给定值计算多项式结果
plot(x,y,’o’,x2,y2)
一.线性代数(LinearAlgebra)
解线性方程(Linearequation)就是找出是否存在一个唯一的矩阵x,使得a,b满足关系:
ax=b或xa=b
MALAB中x=a\b是方程ax=b的解,x=b/a是方程式xa=b的解。
通常线性方程多写成ax=b,“\”较多用,两者的关系为:
(b/a)’=(a’\b’)
系数矩阵a可能是m行n列的,有三种情况:
*方阵系统:
(Squarematrix)m=n可求出精确解(a必须是非奇异(nonsingular),即满秩(fullrank)。
*超定系统:
(Overdetermindsystem)m>n可求出最小二乘解
*欠定系统:
(Underdetermindsystem)mMATLAB对不同形式的参数矩阵,采用不同的运算法则来处理,它会自动检测参数矩阵,以区别下面几种形式:
*三角矩阵(TriangularMatrix)
*对称正定矩阵(symmetricalpositivedeterminedmatrix)
*非奇异方阵(Nonsingularmatrix)
*超定系统(Overdetermindsystem)
*欠定系统(Underdetermindsystem)
1.方阵系统:
(Squarearray)m=n
最常见的是系数矩阵为方阵a,常数项b为列矢量,其解x可写成x=a\b,x和b大小相同。
例1:
求方阵系统的根。
a=[1167;5139;1718]
b=[16134]’
x=a\b
a=
1167
5139
1718
b=
16
13
4
x=
3.9763
5.4455
-8.6303
例2:
假如a,b为两个大小相同的矩阵,求方阵系统的根。
a=[459;18195;1413]
b=[1512;31519;7610]
x=a\b
C=a*x
a=
459
18195
1413
b=
1512
31519
7610
x=
-3.6750-0.73332.9708
3.72501.4667-2.1292
-0.32500.06671.1958
C=
1.00005.000012.0000
3.000015.000019.0000
7.00006.000010.0000
若方阵a的各个行矢量线性相关(linearcorrelation),则称方阵a为奇异矩阵。
这时线性方程将有无穷多组解。
若方阵是奇异矩阵,则反斜线运算因子将发出警告信息。
2.超定系统(Overdetermindsystem)m﹥n
实验数据较多,寻求他们的曲线拟合。
如在t内测得一组数据y:
ty
0.00.82
0.30.72
0.80.63
1.10.60
1.60.55
2.20.50
这些数据显然有衰减指数趋势:
y(t)~c1+c2e-t
此方程意为y矢量可以由两个矢量逐步逼近而得,一个是单行的常数矢量,一个是由指数e-t项构成,两个参数c1和c2可用最小二乘法求得,它们表示实验数据与方程y(t)~c1+c2e-t之间距离的最小平方和。
例1:
求上述数据的最小二乘解。
将数据带入方程式y(t)~c1+c2e-t中,可得到含有两个未知数的6个等式,可写成6行2列的矩阵e.
t=[00.30.81.11.62.2]’;
y=[0.820.720.630.600.550.50]’;
e=[ones(size(t))exp(-t)]%求6个y(t)方程的系数矩阵
c=e\y%求方程的解
e=
1.00001.0000
1.00000.7408
1.00000.4493
1.00000.3329
1.00000.2019
1.00000.1108
c=
0.4744
0.3434
带入方程得:
y(t)~0.4744+0.3434e-t
用此方程可绘制曲线:
t=[00.30.81.11.62.2]’;
y=[0.820.720.630.600.550.50]’;
t1=[0:
0.1:
2.5]’;y1=[ones(size(t1)),exp(-t1)]*c
plot(t1,y1,’b’,t,y,’ro’)
如果一个矩阵的行矢量是线性相关的,则它的最小二乘解并不唯一,因此,a\b运算将给出警告,并产生含有最少元素的基解。
3.欠定系统:
(Underdetermindsystem)m﹤n
欠定系统为线性相关系统,其解都不唯一,MATLAB会计算一组构成通解的基解,而方程的特解则用QR分解法决定。
两种解法:
最少元素解a\b,最小范数解pinv(a)*b.
例:
用两种方法求解欠定系统。
对a和矢量b分别用a\b和pinv(a)*b求解:
a=[111;11-1]
b=[106]’
p=a\b
q=pinv(a)*b
a=
111
11-1
b=
10
6
p=
8.0000
0
2.0000
q=
4.0000
4.0000
2.0000
三.逆矩阵及行列式(Reversanddeterminantofmatrix)
1.方阵的逆和行列式(Reversanddeterminantofsquarematrix)
若a是方阵,且为非奇异阵,则方程ax=I和xa=I有相同的解X。
X称为a的逆矩阵,记做a-1,在MATLAB中用inv函数来计算矩阵的逆。
计算方阵的行列式则用det函数。
DETDeterminant.
DET(X)isthedeterminantofthesquarematrixX.
UseCONDinsteadofDETtotestformatrixsingularity.
INVMatrixinverse.
INV(X)istheinverseofthesquarematrixX.AwarningmessageisprintedifXisbadlyscaledornearlysingular.
例:
计算方阵的行列式和逆矩阵。
a=[3-31;-35-2;1-21];
b=[14135;5112;6145];
d1=det(a)
x1=inv(a)
d2=det(b)
x2=inv(b)
d1=
1
x1=
1.00001.00001.0000
1.00002.00003.0000
1.00003.00006.0000
d2=
-1351
x2=
0.1207-0.0037-0.1118
-0.0348-0.02960.1058
-0.04740.08730.0377
2.广义逆矩阵(伪逆)(Generalizedinversematrix)
一般非方阵无逆矩阵和行列式,方程ax=I和xa=I至少有一个无解,这种矩阵可以求得特殊的逆矩阵,成为广义逆矩阵(generalizedinversematrix)(或伪逆pseudoinverse)。
矩阵amn存在广义逆矩阵xnm,使得ax=Imn,MATLAB用pinv函数来计算广义逆矩阵。
例:
计算广义逆矩阵。
a=[814;13;96]
x=pinv(a)
b=x*a
c=a*x
d=c*a%d=a*x*a=a
e=x*c%e=x*a*x=x
a=
814
13
96
x=
-0.0661-0.04020.1743
0.10450.0406-0.0974
b=
1.0000-0.0000
-0.00001.0000
c=
0.93340.24720.0317
0.24720.0817-0.1177
0.0317-0.11770.9849
d=
8.000014.0000
1.00003.0000
9.00006.0000
e=
-0.0661-0.04020.1743
0.10450.0406-0.0974
PINVPseudoinverse.
X=PINV(A)producesamatrixXofthesamedimensionsasA'sothatA*X*A=A,X*A*X=XandA*XandX*AareHermitian.ThecomputationisbasedonSVD(A)andanysingularvalueslessthanatolerancearetreatedaszero.
ThedefaulttoleranceisMAX(SIZE(A))*NORM(A)*EPS.
PINV(A,TOL)usesthetoleranceTOLinsteadofthedefault.
四.矩阵分解(Matrixdecomposition)
通过矩阵分解的方法求解大型方程组非常有效的,这种方法可以使运算速度加快,节省磁盘空间,节省内存。
MATLAB求解线性方程的过程基于三种分解法则:
(1)Cholesky分解,针对对称正定矩阵;
(2)LU分解,高斯消元法,针对一般矩阵;
(3)QR分解,正交化,针对一般长方形矩阵(行数≠列数)
这三种分解运算分别由chol,lu和qr三个函数来分解.
1.Cholesky分解(CholeskyDecomposition)
仅适用于对称和上三角矩阵,如果A为对称正定矩阵,则Cholesky分解可将矩阵A分解为上三角矩阵和其转置的乘积,即:
A=R’*R,R为上三角阵。
方程A*X=b变成R’*R*X=b,所以有X=R\(R’\b)
例:
cholesky分解。
a=pascal(6)
b=chol(a)
a=
111111