第四章MATLAB地数值计算功能.docx
《第四章MATLAB地数值计算功能.docx》由会员分享，可在线阅读，更多相关《第四章MATLAB地数值计算功能.docx（49页珍藏版）》请在冰点文库上搜索。
第四章MATLAB地数值计算功能
第四章MATLAB的数值计算功能
一．多项式`
1．多项式的表达与创建
Matlab用矢量表达多项式系数，元素按降幂排列：
P（x）=a0xn+a1xn-1+a2xn-2…an-1x+a0
其系数矢量为：
P=[a0a1…an-1an]
如将根矢量表示为：
ar=[ar1ar2…arn]
则根矢量与系数矢量之间关系为：
（x-ar1）（x-ar2）…（x-arn）=a0xn+a1xn-1+a2xn-2…an-1x+a0
多项式系数矢量可通过调用函数p=poly（ar）产生
例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）根矢量元素为n，则多项式系数矢量元素为n+1；
2）函数poly2sym（pa）把多项式系数矢量表达成符号形式的多项式，缺省情况下自变量符号为x，可以指定自变量。
（3）使用简单绘图函数可以直接绘制符号形式多项式的曲线。
例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阶的。
例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）把虚部滤掉。
例4：
将多项式的系数表示形式转换为根表现形式，poly和roots互为逆函数。
求x3-6x2-72x-27的根
a=[1-6-72-27]
r=roots（a）
r=
12.1229
-5.7345
-0.3884
MATLAB约定，多项式系数矢量用行矢量表示，根矢量用列矢量表示。
>>
1.多项式的乘除运算
多项式乘法用函数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’）
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.其他常用的多项式运算命令
pa=polyval（p,s）按数组运算规则计算给定s时多项式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称为极点，r1,r2,…,rn称为留数，k（s）称为直项，假如a（s）含有m重根pj，则相应部分应写成：
例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：
对一组实验数据进行多项式最小二乘拟合
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）
二．线性代数
解线性方程就是找出是否存在一个唯一的矩阵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列的，有三种情况：
*方阵系统:
m=n可求出精确解（a必须是非奇异，即满秩）
*超定系统：
m>n可求出最小二乘解
*欠定系统：
mMATLAB对不同形式的参数矩阵，采用不同的运算法则来处理，它会自动检测参数矩阵，以区别下面几种形式：
*三角矩阵
*对称正定矩阵
*非奇异方阵
*超定系统
*欠定系统
1.方阵系统：
最常见的是系数矩阵为方阵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的各个行矢量线性相关，则称方阵a为奇异矩阵。
这时线性方程将有无穷多组解。
若方阵是奇异矩阵，则反斜线运算因子将发出警告信息。
2．超定系统实验数据较多，寻求他们的曲线拟合。
如在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.欠定系统:
欠定系统为线性相关系统，其解都不唯一，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
三．逆矩阵及行列式
1．方阵的逆和行列式
若a是方阵，且为非奇异阵，则方程ax=I和xa=I有相同的解X。
X称为a的逆矩阵，记做a-1，在MATLAB中用inv函数来计算矩阵的逆。
计算方阵的行列式则用det函数。
例：
计算方阵的行列式和逆矩阵。
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．广义逆矩阵（伪逆）
一般非方阵无逆矩阵和行列式，方程ax=I和xa=I至少有一个无解，这种矩阵可以求得特殊的逆矩阵，成为广义逆矩阵（或伪逆）。
矩阵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
四．矩阵分解
MATLAB求解线性方程的过程基于三种分解法则：
（１）Cholesky分解，针对对称正定矩阵；
（２）高斯消元法，针对一般矩阵；
（３）正交化，针对一般矩阵（行数≠列数）
这三种分解运算分别由chol,lu和qr三个函数来分解.
1．Cholesky分解
例：
cholesky分解。
a=pascal（6）
b=chol（a）
a=
111111
123456
136101521
1410203556
15153570126
162156126252
b=
111111
012345
0013610
0001410
000015
000001
2．LU分解
用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分解
函数调用格式：
[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.特征值与特征矢量
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
5.奇异值分解.
如存在两个矢量u,v及一常数c,使得矩阵A满足：
Av=cu,A’u=cv
称c为奇异值，称u,v为奇异矢量。
将奇异值写成对角方阵∑，而相对应的奇异矢量作为列矢量则可写成两个正交矩阵U，V，使得：
AV=U∑，A‘U=V∑因为U，V正交，所以可得奇异值表达式：
A=U∑V’。
一个m行n列的矩阵A经奇异值分解，可求得m行m列的U,m行n列的矩阵∑和n行n列的矩阵V.。
奇异值分解用svd函数实现，调用格式为；
[u,s,v]=svd（a）
（SVDSingularvaluedecomposition.
[U,S,V]=SVD（X）producesadiagonalmatrixS,ofthesame
dimensionasXandwithnonnegativediagonalelementsin
decreasingorder,andunitarymatricesUandVsothat
X=U*S*V'.）
例：
奇异值分解。
a=[85;73;46];
[u,s,v]=svd（a）%s为奇异值对角方阵
u=
-0.6841-0.1826-0.7061
-0.5407-0.52280.6591
-0.48950.83270.2589
s=
13.76490
03.0865
00
v=
-0.8148-0.5797
-0.57970.8148
五．数据分析
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各列元素的累积j积
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（maxA））或Amax=max（A（:
）），
Amin=min（min（A））或Amin=min（A（:
））
2．协方差阵和相关阵（表4—2）
例：
计算协方差和相关阵。
x=rand（10,3）;
y=rand（10,3）;
cx=cov（x）%求协方差阵
cy=cov（y）
cxy=cov（x,y）%求两随机变量的协方差
px=corrcoef（x,y）%求相关阵
pxy=corrcoef（x,y）%求两随机变量的（2×2）相关系数
cx=
0.0893-0.0586-0.0320
-0.05860.07190.0298
-0.03200.02980.0617
cy=
0.0805-0.03080.0099
-0.03080.0761-0.0548
0.0099-0.05480.0667
cxy=
0.0978-0.0211
-0.02110.0696
px=
1.0000-0.2561
-0.25611.0000
pxy=
1.0000-0.2561
-0.25611.0000
2.微分与梯度（表4—3）
例1：
按列求微分。
x=[1,10,20;2,12,23;3,14,26;3,16,29]
d=diff（x）%求一阶微分
x=
11020
21223
31426
31629
d=
123
123
023
例2：
对于（u=x2+y2和Δ2=4）求5点差分。
[x,y]=meshgrid（-4:
4,-3:
3）;
u=x.^2+y.^2
v4=4*del2（u）%求m×n阶矩阵U的五点差分矩阵
u=
25181310910131825
2013854581320
1710521251017
16941014916
1710521251017
2013854581320
25181310910131825
v4=
444444444
444444444
444444444
444444444
444444444
444444444
444444444
例3：
产生一个二元函数偏导数和梯度。
x=-2:
0.2:
2;
y=-2:
0.2:
2;
[xx,yy]=meshgrid（x,y）;
z=
升级会员 