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

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

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

第四章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+2

poly2sym（[1-432]）

ans=

x^3-4*x^2+3*x+2

b）由根矢量创建多项式

通过调用函数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-144

ppa=

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

pr=

1.00001.10000.55000.1250

ppr=

x^3+11/10*x^2+11/20*x+1/8

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

MATLAB约定，多项式系数矢量用行矢量表示，根矢量用列矢量表示。

>>

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’——是指定变量为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

1.其他常用的多项式运算命令（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）之比，其部分分式展开为：

式中：

p1,p2,…,pn称为极点（poles），r1,r2,…,rn称为留数（residues），k（s）称为直项（directterms），假如a（s）含有m重根pj，则相应部分应写成：

例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

1.方阵系统：

（Squarearray）

最常见的是系数矩阵为方阵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）

实验数据较多，寻求他们的曲线拟合。

如在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）

不定系统为线性相关系统，其解都不唯一，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.

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=

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

四．矩阵分解（Matrixdecomposition）

MATLAB求解线性方程的过程基于三种分解法则：

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

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

（３）正交化，针对一般矩阵（行数≠列数）

这三种分解运算分别由chol,lu和qr三个函数来分解.

1．Cholesky分解（CholeskyDecomposition）

仅适用于对称和上三角矩阵

例：

cholesky分解。

a=pascal（6）

b=chol（a）

a=

111111

123456

136101521

1410203556

15153570126

162156126252

b=

111111

012345

0013610

0001410

000015

000001

2．LU分解（LUfactorization）.

用lu函数完成LU分解，将矩阵分解为上、下两个三角阵，其调用格式为：

[l,u]=lu（a）l代表下三角阵，u代表上三角阵。

例：

LU分解。

a=[472422;11440;303841]

[l,u]=lu（a）

a=

472422

11440

303841

l=

1.000000

0.23401.00000

0.63830.59091.0000

u=

47.000024.000022.0000

038.3830-5.1489

0030.0000

3．QR分解（Orthogonal-triangulardecomposition）.

函数调用格式：

[q,r]=qr（a）,q代表正规正交矩阵，r代表三角形矩阵。

原始阵a不必一定是方阵。

如果矩阵a是m×n阶的，则矩阵q是m×m阶的，矩阵r是m×n阶的。

例：

QR分解.

A=[22462020;30364644;398452];

[q,r]=qr（A）

q=

-0.4082-0.7209-0.5601

-0.5566-0.28980.7786

-0.72360.6296-0.2829

r=

-53.8981-44.6027-66.3289-34.1014

0-38.55640.5823-25.9097

0011.880022.4896

4.特征值与特征矢量（Eigenvaluesandeigenvectors）.

MATLAB中使用函数eig计算特征值和特征矢量，有两种调用方法：

*e=eig（a）,其中e是包含特征值的矢量；

*[v,d]=eig（a）,其中v是一个与a相同的n×n阶矩阵，它的每一列是矩阵a的一个特征值所对应的特征矢量，d为对角阵，其对角元素即为矩阵a的特征值。

例：

计算特征值和特征矢量。

a=[342515;18359;41219]

e=eig（a）

[v,d]=eig（a）

a=

342515

18359

41219

e=

68.5066

15.5122

-6.0187

v=

-0.6227-0.4409-0.3105

-0.49690.6786-0.0717

-0.6044-0.58750.9479

d=

68.506600

015.51220

00-6.0187

EIGEigenvaluesandeigenvectors.

五．数据分析（DataAnalyaia）

MATLAB对数据分析有两条约定：

（1）若输入量X是矢量，则不论是行矢量还是列矢量，运算是对整个矢量进行的；

（2）若输入量X是数组，（或称矩阵），则命令运算是按列进行的。

即默认每个列是有一个变量的不同“观察“所得的数据组成。

1.基本统计命令（表4-1）

例：

做各种基本统计运算。

A=[5-10-60;263-3;-95-1011;-22170-19;-16-44]

Amax=max（A）%找A各列的最大元素

Amin=min（A）%找A各列的最小元素

Amed=median（A）%找A各列的中位元素

Amean=mean（A）%找A各列的平均值

Astd=std（A）%求A各列的标准差

Aprod=prod（A）%求A各列元素的积

Asum=sum（A）%求A各列元素的和

S=cumsum（A）%求A各列元素的累积和

P=cumprod（A）%求A各列元素的累积积

I=sort（A）%使A的各列元素按递增排列

A=

5-10-60

263-3

-95-1011

-22170-19

-16-44

Amax=

517311

Amin=

-22-10-10-19

Amed=

-16-40

Amean=

-5.00004.8000-3.4000-1.4000

Astd=

10.83979.62815.079411.1490

Aprod=

-1980-3060000

Asum=

-2524-17-7

S=

5-10-60第一行就是原数据，

7-4-3-3第二行是前两行之和，

-21-138以此类推。

-2418-13-11

-2524-17-7

P=

5-10-60

10-60-180

-90-3001800

1980-510000

-1980-3060000

I=

-22-10-10-19

-95-6-3

-16-40

2604

517311

>>

求矩阵元素的最大值、最小值可用：

Amax=max（max（A））或Amax=max（A（:

）），

Amin=min（min（A））或Amin=min（A（:

））

六．插值:

（Interpolation）

在已知数据之间计算估计值的过程。

1．一维插值（1DInterpolation）

由interp1实现，用多项式技术计算插值点。

Yi=interp1（x,y,xi,method）y—函数值矢量，x—自变量取值范围，xi—插值点的自变量矢量，

Method—插值方法选项。

MATLAB6.1的4种方法：

*临近点插值：

method=‘nearest’

*线性插值：

method=‘linear’

*三次样条插值：

method=‘spline’

*立方插值：

method=‘pchip’or‘cubic’

选择插值方法时主要考虑因素：

运算时间、占用计算机内存和插值的光滑程度。

比较：

运算时间、占用计算机内存光滑程度。

*临近点插值：

快少差

*线性插值：

稍长较多稍好

*三次样条插值：

最长较多最好

*立方插值：

较长多较好

例1：

一维插值函数插值方法的对比。

x=0:

10;

y=sin（x）;

xi=0:

0.25:

10;

strmod={'nearest','linear','spline','cubic'}%将插值方法定义为单元数组

str1b={'（a）method=nearest','（b）method=linear',...

'（c）method=spline','（d）method=cubic'}%将图标定义为单元数组

fori=1:

4

yi=interp1（x,y,xi,strmod{i}）;

subplot（2,2,i）

plot（x,y,'ro',xi,yi,'b'）,xlabel（str1b（i））

end

strmod=

'nearest''linear''spline''cubic'

例2：

三次样条插值

x0=0:

10;

y0=sin（x0）;

x=0:

.25:

10;

y=spline（x0,y0,x）;

plot（x0,y0,'or',x,y,'k'）

与interp1结果一样

2．二维插值（2DInterpolation）

用于图形图象处理和三维曲线拟合等领域，由interp2实现，一般格式为：

ZI=interp2（X,Y,Z,XI,YI,method）X,Y—自变量组成的数组，尺寸相同。

XI,YI—插值点的自变量数组

Method—插值方法选项，4种

*临近点插值：

method=‘nearest’

*线性插值：

method=‘linear’该方法是interp2函数的缺省方法

*三次样条插值：

method=‘spline’

*立方插值：

method=‘pchip’or‘cubic’

例：

二维插值4种方法的对比。

[x,y,z]=peaks（7）;figure