第四章MATLAB地数值计算功能.docx
- 文档编号:5357444
- 上传时间:2023-05-08
- 格式:DOCX
- 页数:49
- 大小:58.70KB
第四章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可求出最小二乘解
*欠定系统:
m MATLAB对不同形式的参数矩阵,采用不同的运算法则来处理,它会自动检测参数矩阵,以区别下面几种形式: *三角矩阵 *对称正定矩阵 *非奇异方阵 *超定系统 *欠定系统 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求解线性方程的过程基于三种分解法则: (1)Cholesky分解,针对对称正定矩阵; (2)高斯消元法,针对一般矩阵; (3)正交化,针对一般矩阵(行数≠列数) 这三种分解运算分别由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=
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第四 MATLAB 数值 计算 功能