第四章 MATLAB的数值计算功能Word下载.docx

第四章 MATLAB的数值计算功能Word下载.docx

《第四章 MATLAB的数值计算功能Word下载.docx》由会员分享，可在线阅读，更多相关《第四章 MATLAB的数值计算功能Word下载.docx（42页珍藏版）》请在冰点文库上搜索。

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=

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=

注：

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

展开（s2+2s+2）（s+4）（s+1）（多个多项式相乘）

c=conv（[1,2,2],conv（[1,4],[1,1]））

cs=poly2sym（c,’s’）%（指定变量为s）

1716188

s^4+7*s^3+16*s^2+18*s+8

求多项式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（A）,otherwise

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）

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）m<

n可尝试找出含有最少m个基解或最小范数解

MATLAB对不同形式的参数矩阵，采用不同的运算法则来处理，它会自动检测参数矩阵，以区别下面几种形式：

*三角矩阵（TriangularMatrix）

*对称正定矩阵（symmetricalpositivedeterminedmatrix）

*非奇异方阵（Nonsingularmatrix）

*超定系统（Overdetermindsystem）

*欠定系统（Underdetermindsystem）

1.方阵系统：

（Squarearray）m=n

最常见的是系数矩阵为方阵a，常数项b为列矢量,其解x可写成x=a\b,x和b大小相同。

求方阵系统的根。

a=[1167;

5139;

1718]

b=[16134]’

x=a\b

a=

1167

5139

1718

b=

16

13

4

x=

3.9763

5.4455

-8.6303

假如a,b为两个大小相同的矩阵，求方阵系统的根。

a=[459;

18195;

1413]

b=[1512;

31519;

7610]

C=a*x

459

18195

1413

1512

31519

7610

-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之间距离的最小平方和。

求上述数据的最小二乘解。

将数据带入方程式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

0.4744

0.3434

带入方程得：

y（t）~0.4744+0.3434e-t

用此方程可绘制曲线：

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

111

11-1

10

6

8.0000

0

2.0000

q=

4.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

814

13

96

-0.0661-0.04020.1743

0.10450.0406-0.0974

1.0000-0.0000

-0.00001.0000

0.93340.24720.0317

0.24720.0817-0.1177

0.0317-0.11770.9849

d=

8.000014.0000

1.00003.0000

9.00006.0000

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求解线性方程的过程基于三种分解法则：

（１）Cholesky分解，针对对称正定矩阵；

（２）LU分解，高斯消元法，针对一般矩阵；

（３）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）

111111