第四章MATLAB的数值计算功能

第四章MATLAB的数值计算功能Chapter4:NumericalcomputationofMATLAB一、多项式(Polynomial)`.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函数直接输入多项式的系数矢量，就可方便的建立符号形式的多项式。例：创建多项式x3-4x2+3x+2poly2sym([1-432])ans=x^3-4*x^2+3*x+2b)由根矢量创建多项式通过调用函数p=poly(ar)产生多项式的系数矢量,再利用poly2sym函数就可方便的建立符号形式的多项式。注：（1）根矢量元素为n，则多项式系数矢量元素为n+1；（2）函数poly2sym(pa)把多项式系数矢量表达成符号形式的多项式，缺省情况下自变量符号为x，可以指定自变量。（3）使用简单绘图函数ezplot可以直接绘制符号形式多项式的曲线。例1：由根矢量创建多项式。将多项式(x-6)(x-3)(x-8)表示为系数形式a=[638]%根矢量pa=poly(a)%求系数矢量ppa=poly2sym(pa)%以符号形式表示原多项式ezplot(ppa,[-50,50])pa=1-1790-144ppa=x^3-17*x^2+90*x-144注：含复数根的根矢量所创建的多项式要注意：（1）要形成实系数多项式，根矢量中的复数根必须共轭成对；（2）含复数根的根矢量所创建的多项式系数矢量中，可能带有很小的虚部，此时可采用取实部的命令(real)把虚部滤掉。进行多项式的求根运算时，有两种方法，一是直接调用求根函数roots，poly和roots互为逆函数。另一种是先把多项式转化为伴随矩阵，然后再求其特征值，该特征值即是多项式的根。例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.1250pr=1.00001.10000.55000.1250ppr=x^3+11/10*x^2+11/20*x+1/8c)特征多项式输入法用poly函数可实现由矩阵的特征多项式系数创建多项式。条件：特征多项式系数矢量的第一个元素必须为一。例2：求三阶方阵A的特征多项式系数，并转换为多项式形式。a=[638;756;135]Pa=poly(a)%求矩阵的特征多项式系数矢量Ppa=poly2sym(pa)Pa=1.0000-16.000038.0000-83.0000Ppa=x^3-16*x^2+38*x-83注：n阶方阵的特征多项式系数矢量一定是n+1阶的。注：（1）要形成实系数多项式，根矢量中的复数根必须共轭成对；（2）含复数根的根矢量所创建的多项式系数矢量中，可能带有很小的虚部，此时可采用取实部的命令(real)把虚部滤掉。进行多项式的求根运算时，有两种方法，一是直接调用求根函数roots，poly和roots互为逆函数。另一种是先把多项式转化为伴随矩阵，然后再求其特征值，该特征值即是多项式的根。例4：将多项式的系数表示形式转换为根表现形式。求x3-6x2-72x-27的根a=[1-6-72-27]r=roots(a)r=12.1229-5.7345-0.3884MATLAB约定，多项式系数矢量用行矢量表示，根矢量用列矢量表示。1.多项式的乘除运算(Multiplicationanddivisionofpolynomial)多项式乘法用函数conv(a,b)实现，除法用函数deconv(a,b)实现。例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’——是指定变量为sc=413282718cs=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=1716188cs=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=1342r1=00000q2=1446r2=0000-10cc=17161818test=111111.其他常用的多项式运算命令(Othercomputationcommandofpolynomial)pa=polyval(p,s)按数组运算规则计算给定s时多项式p的值，p为多项式系数矢。pm=polyvalm(p,s)按矩阵运算规则计算给定s时多项式p的值。[r,p,k]=residue(b,a)部分分式展开，b,a分别是分子分母多项式系数矢量，r,p,k分别是留数、极点和直项矢量p=polyfit(x,y,n)用n阶多项式拟合x，y矢量给定的数据。polyder(p)多项式微分。注：对于多项式b(s)与不重根的n阶多项式a(s)之比，其部分分式展开为：)()()(2211skpsrLpsrpsrsasbnn式中：p1,p2,…,pn称为极点(poles)，r1,r2,…,rn称为留数(residues)，k(s)称为直项(directterms)，假如a(s)含有m重根pj，则相应部分应写成：112()()jjjmmjjjrrrspspsp例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=bMALAB中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)mn可求出最小二乘解*不定系统：(Underdetermindsystem)mn可尝试找出含有最少m个基解或最小范数解1.方阵系统：(Squarearray)最常见的是系数矩阵为方阵a，常数项b为列矢量,其解x可写成x=a\b,x和b大小相同。例1：求方阵系统的根。a=[1167;5139;1718]b=[16134]’x=a\ba=116751391718b=16134x=3.97635.4455-8.6303例2：假如a,b为两个大小相同的矩阵，求方阵系统的根。a=[459;18195;1413]b=[1512;31519;7610]x=a\bC=a*xa=459181951413b=1512315197610x=-3.6750-0.73332.97083.72501.4667-2.1292-0.32500.06671.1958C=1.00005.000012.00003.000015.000019.00007.00006.000010.0000若方阵a的各个行矢量线性相关(linearcorrelation)，则称方阵a为奇异矩阵。这时线性方程将有无穷多组解。若方阵是奇异矩阵，则反斜线运算因子将发出警告信息。2．超定系统(Overdetermindsystem)实验数据较多，寻求他们的曲线拟合。如在t内测得一组数据y：ty0.00.820.30.720.80.631.10.601.60.552.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.00001.00000.74081.00000.44931.00000.33291.00000.20191.00000.1108c=0.47440.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)]*cplot(t1,y1,’b’,t,y,’ro’)如果一个矩阵的行矢量是线性相关的，则它的最小二乘解并不唯一，因此，a\b运算将给出警告，并产生含有最少元素的基解。3.不定系统:(Underdetermindsystem)不定系统为线性相关系统，其解都不唯一，MATLAB会计算一组构成通解的基解，而方程的特解则用QR分解法决定。两种解法：最少元素解a\b，最小范数解pinv(a)*b.例：用两种方法求解欠定系统。对a和矢量b分别用a\b和pinv(a)*b求解:a=[111;11-1]b=[106]’p=a\bq=pinv(a)*ba=11111-1b=106p=8.000002.0000q=4.00004.00002.0000三．逆矩阵及行列式(Reversanddeterminantofmatrix)1．方阵的逆和行列式(Reversanddeterminantofsquarematrix)若a是方阵，且为非奇异阵，则方程ax=I和xa=I有相同的解X。X称为a的逆矩阵，记做a-1，在MATLAB中用inv函数来计算矩阵的逆。计算方阵的行列式则用det函数。DETDeterminant.INVMatrixinverse.例：计算方阵的行列式和逆矩阵。a=[3-31;-35-2;1-21];b=[14135;5112;6145];d1=det(a)x1=inv(a)d2=det(b)x2=inv(b)d1=1x1=1.00001.00001.00001.00002.00003.00001.00003.00006.0000d2=-1351x