第四章MATLAB地数值计算功能Word文件下载.docx

1、求 x3-6x2-72x-27的根a=1 -6 -72 -27r=roots（a）r = 12.1229 -5.7345 -0.3884MATLAB约定，多项式系数矢量用行矢量表示，根矢量用列矢量表示。1. 多项式的乘除运算多项式乘法用函数conv（a,b）实现， 除法用函数deconv（a,b）实现。例1：a（s）=s2+2s+3, b（s）=4s2+5s+6,计算 a（s）与 b（s）的乘积。a=1 2 3; b=4 5 6;c=conv（a,b）cs=poly2sym（c,s）c = 4 13 28 27 18cs =4*s4+13*s3+28*s2+27*s+18例2： 展开（s2+2

2、s+2）（s+4）（s+1） （多个多项式相乘）c=conv（1,2,2,conv（1,4,1,1）cs=poly2sym（c,s） （指定变量为s） 1 7 16 18 8s4+7*s3+16*s2+18*s+8求多项式s4+7*s3+16*s2+18*s+8分别被（s+4）,（s+3）除后的结果。c=1 7 16 18 8;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 = 1 3 4 2r1 = 0 0 0 0 0q2 = 1 4 4 6r2 =

3、0 0 0 0 -10cc = 1 7 16 18 18test = 1 1 1 1 11. 其他常用的多项式运算命令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称为极点，

4、r1,r2,rn 称为留数，k（s）称为直项，假如a（s）含有m重根pj，则相应部分应写成：例3：对 （3x4+2x3+5x2+4x+6）/（x5+3x4+4x3+2x2+7x+2） 做部分分式展开a=1 3 4 2 7 2;b=3 2 5 4 6;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.29

5、91 k = （分母阶数高于分子阶数时，k将是空矩阵，表示无此项）例5：对一组实验数据进行多项式最小二乘拟合x=1 2 3 4 5; % 实验数据y=5.5 43.1 128 290.7 498.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=bMALAB中x=ab 是方程 ax=b 的解， x=b/a是方程式xa=b的解。通常线性方程多写成ax=b，“”较多用，两者的关系为

6、：（b/a）=（ab）系数矩阵a可能是m行n列的，有三种情况：*方阵系统: m=n 可求出精确解（a必须是非奇异，即满秩）*超定系统 ：mn 可求出最小二乘解*欠定系统： mn 可尝试找出含有最少m个基解或最小范数解MATLAB对不同形式的参数矩阵，采用不同的运算法则来处理，它会自动检测参数矩阵，以区别下面几种形式：*三角矩阵*对称正定矩阵*非奇异方阵*超定系统*欠定系统1. 方阵系统：最常见的是系数矩阵为方阵a，常数项b为列矢量, 其解x可写成x=ab, x和b大小相同。 求方阵系统的根。a=11 6 7; 5 13 9; 17 1 8b=16 13 4x=aba = 11 6 7 5 13

7、 9 17 1 8b = 16 13 4x = 3.9763 5.4455 -8.6303假如a,b为两个大小相同的矩阵，求方阵系统的根。a=4 5 9; 18 19 5; 1 4 13b=1 5 12; 3 15 19; 7 6 10C=a*x 4 5 9 18 19 5 1 4 13 1 5 12 3 15 19 7 6 10 -3.6750 -0.7333 2.9708 3.7250 1.4667 -2.1292 -0.3250 0.0667 1.1958C = 1.0000 5.0000 12.0000 3.0000 15.0000 19.00007.0000 6.0000 10.00

8、00若方阵a的各个行矢量线性相关，则称方阵a为奇异矩阵。这时线性方程将有无穷多组解。若方阵是奇异矩阵，则反斜线运算因子将发出警告信息。2超定系统 实验数据较多，寻求他们的曲线拟合。 如在t内测得一组数据y：t y0.0 0.820.3 0.720.8 0.631.1 0.601.6 0.552.2 0.50这些数据显然有衰减指数趋势： y（t）c1+c2e-t此方程意为y矢量可以由两个矢量逐步逼近而得，一个是单行的常数矢量，一个是由指数e-t项构成，两个参数c1和c2可用最小二乘法求得，它们表示实验数据与方程y（t）c1+c2e-t之间距离的最小平方和。 求上述数据的最小二乘解。将数据带入方程

9、式y（t）c1+c2e-t中，可得到含有两个未知数的6个等式，可写成6行2 列的矩阵e.t=0 0.3 0.8 1.1 1.6 2.2;y=0.82 0.72 0.63 0.60 0.55 0.50;e=ones（size（t） exp（-t） %求6个y（t）方程的系数矩阵c=ey % 求方程的解e = 1.0000 1.0000 1.0000 0.7408 1.0000 0.4493 1.0000 0.3329 1.0000 0.2019 1.0000 0.1108 0.47440.3434带入方程得：y（t）0.4744+0.3434e-t用此方程可绘制曲线：t1=0:0.1:2.5;

10、y1=ones（size（t1）,exp（-t1）*cplot（t1,y1,b,t,y,ro）如果一个矩阵的行矢量是线性相关的，则它的最小二乘解并不唯一，因此，ab运算将给出警告，并产生含有最少元素的基解。3 .欠定系统: 欠定系统为线性相关系统，其解都不唯一，MATLAB会计算一组构成通解的基解，而方程的特解则用QR分解法决定。两种解法：最少元素解ab，最小范数解pinv（a）*b.例： 用两种方法求解欠定系统。对a和矢量b分别用ab和pinv（a）*b求解:a=1 1 1; 1 1 -1b=10 6p=abq=pinv（a）*b 1 1 1 1 1 -1 10 6 8.0000 0 2.0

11、000q = 4.0000三 逆矩阵及行列式1 方阵的逆和行列式若a是方阵，且为非奇异阵，则方程ax=I和 xa=I有相同的解X。X称为a的逆矩阵，记做a-1，在MATLAB中 用inv 函数来计算矩阵的逆。计算方阵的行列式则用det函数。计算方阵的行列式和逆矩阵。a=3 -3 1;-3 5 -2;1 -2 1;b=14 13 5; 5 1 12;6 14 5;d1=det（a）x1=inv（a）d2=det（b）x2=inv（b）d1 = 1x1 = 1.0000 1.0000 1.0000 1.0000 2.0000 3.0000 1.0000 3.0000 6.0000d2 = -135

12、1x2 = 0.1207 -0.0037 -0.1118 -0.0348 -0.0296 0.1058 -0.0474 0.0873 0.03772 广义逆矩阵（伪逆）一般非方阵无逆矩阵和行列式，方程ax=I 和xa=I至少有一个无解，这种矩阵可以求得特殊的逆矩阵，成为广义逆矩阵（或伪逆）。矩阵amn存在广义逆矩阵xnm，使得 ax=Imn， MATLAB用pinv函数来计算广义逆矩阵。计算广义逆矩阵。a=8 14; 1 3; 9 6x=pinv（a）b=x*ac=a*xd=c*a %d=a*x*a=ae=x*c %e=x*a*x=x 8 14 1 3 9 6 -0.0661 -0.0402

13、0.1743 0.1045 0.0406 -0.0974 1.0000 -0.0000 -0.0000 1.0000 0.9334 0.2472 0.0317 0.2472 0.0817 -0.1177 0.0317 -0.1177 0.9849d = 8.0000 14.0000 1.0000 3.0000 9.0000 6.00000.1045 0.0406 -0.0974四 矩阵分解MATLAB求解线性方程的过程基于三种分解法则：（）Cholesky分解，针对对称正定矩阵；（）高斯消元法， 针对一般矩阵；（）正交化， 针对一般矩阵（行数列数）这三种分解运算分别由chol, lu和 qr三

14、个函数来分解.1 Cholesky分解cholesky分解。a=pascal（6）b=chol（a） 1 1 1 1 1 1 1 2 3 4 5 6 1 3 6 10 15 21 1 4 10 20 35 56 1 5 15 35 70 126 1 6 21 56 126 252 0 1 2 3 4 5 0 0 1 3 6 10 0 0 0 1 4 10 0 0 0 0 1 50 0 0 0 0 12 LU分解用lu函数完成LU分解，其调用格式为：l,u=lu（a） l代表下三角阵，u代表上三角阵。LU分解。a=47 24 22; 11 44 0;30 38 41l,u=lu（a） 47 24

15、 22 11 44 0 30 38 41l = 1.0000 0 0 0.2340 1.0000 0 0.6383 0.5909 1.0000u = 47.0000 24.0000 22.0000 0 38.3830 -5.1489 0 0 30.00003 QR分解函数调用格式：q,r=qr（a）, q代表正规正交矩阵，r代表三角形矩阵。原始阵a不必一定是方阵。如果矩阵a是mn阶的，则矩阵q是mm阶的，矩阵r是mn阶的。QR分解.A=22 46 20 20; 30 36 46 44;39 8 45 2;q,r=qr（A） -0.4082 -0.7209 -0.5601 -0.5566 -0.

16、2898 0.7786 -0.7236 0.6296 -0.2829 -53.8981 -44.6027 -66.3289 -34.1014 0 -38.5564 0.5823 -25.9097 0 0 11.8800 22.48964. 特征值与特征矢量MATLAB中使用函数eig计算特征值和 特征矢量，有两种调用方法：*e=eig（a）, 其中e是包含特征值的矢量；*v,d=eig（a）, 其中v是一个与a相同的nn阶矩阵，它的每一列是矩阵a的一个特征值所对应的特征矢量，d为对角阵，其对角元素即为矩阵a的特征值。计算特征值和特征矢量。a=34 25 15; 18 35 9; 41 21 9

17、e=eig（a）v,d=eig（a） 34 25 15 18 35 9 41 21 9 68.5066 15.5122 -6.0187v = -0.6227 -0.4409 -0.3105 -0.4969 0.6786 -0.0717 -0.6044 -0.5875 0.9479 68.5066 0 0 0 15.5122 0 0 0 -6.01875. 奇异值分解.如存在两个矢量u,v及一常数c,使得矩阵A满足：Av=cu, Au=cv称c为奇异值，称u,v为奇异矢量。 将奇异值写成对角方阵，而相对应的奇异矢量作为列矢量则可写成两个正交矩阵U，V， 使得： AV=U， AU=V 因为U，V正

18、交，所以可得奇异值表达式： A=UV。一个m行n列的矩阵A经奇异值分解，可求得m行m列的U, m行n列的矩阵和n行n列的矩阵V.。奇异值分解用svd函数实现，调用格式为；u,s,v=svd（a） （ SVD Singular value decomposition. U,S,V = SVD（X） produces a diagonal matrix S, of the same dimension as X and with nonnegative diagonal elements in decreasing order, and unitary matrices U and V so th

19、at X = U*S*V.） 奇异值分解。a=8 5; 7 3;4 6;u,s,v=svd（a） % s为奇异值对角方阵 -0.6841 -0.1826 -0.7061 -0.5407 -0.5228 0.6591 -0.4895 0.8327 0.2589 13.7649 0 0 3.0865 0 0 -0.8148 -0.5797 -0.5797 0.8148 五 数据分析MATLAB对数据分析有两条约定：（1） 若输入量X是矢量，则不论是行矢量还是列矢量，运算是对整个矢量进行的； （2）若输入量X是数组，（或称矩阵），则命令运算是按列进行的。即默认每个列是有一个变量的不同“观察“所得的数

20、据组成。1. 基本统计命令 （表4-1） 做各种基本统计运算。A=5 -10 -6 0;2 6 3 -3;-9 5 -10 11;-22 17 0 -19;-1 6 -4 4Amax=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

21、） %使A的各列元素按递增排列A = 5 -10 -6 0 2 6 3 -3 -9 5 -10 11 -22 17 0 -19 -1 6 -4 4Amax = 5 17 3 11Amin = -22 -10 -10 -19Amed = -1 6 -4 0Amean = -5.0000 4.8000 -3.4000 -1.4000Astd = 10.8397 9.6281 5.0794 11.1490Aprod = -1980 -30600 0 0Asum = -25 24 -17 -7S = 7 -4 -3 -3 -2 1 -13 8 -24 18 -13 -11 -25 24 -17 -7P

22、 = 5 -10 -6 0 10 -60 -18 0 -90 -300 180 0 1980 -5100 0 0I = -9 5 -6 -3 -1 6 -4 0 2 6 0 4求矩阵元素的最大值、最小值可用： Amax=max（maxA） 或 Amax=max（A（:）， Amin=min（min（A） 或 Amin=min（A（:） 2 协方差 阵和相关阵（表 42） 计算协方差和相关阵。x=rand（10,3）;y=rand（10,3）;cx=cov（x） %求协方差阵cy=cov（y）cxy=cov（x,y） %求两随机变量的协方差px=corrcoef（x,y） %求相关阵pxy=co

23、rrcoef（x,y） %求两随机变量的（22）相关系数cx = 0.0893 -0.0586 -0.0320 -0.0586 0.0719 0.0298 -0.0320 0.0298 0.0617cy = 0.0805 -0.0308 0.0099 -0.0308 0.0761 -0.0548 0.0099 -0.0548 0.0667cxy = 0.0978 -0.0211 -0.0211 0.0696px = 1.0000 -0.2561 -0.2561 1.0000pxy =2. 微分与梯度 （表43）按列求微分。x=1,10,20;2,12,23;3,14,26;3,16,29d=d

24、iff（x） %求一阶微分 1 10 20 2 12 23 3 14 26 3 16 29 1 2 30 2 3 对于（u=x2+y2和2=4）求5点差分。x,y=meshgrid（-4:4,-3:3）;u=x.2+y.2v4=4*del2（u） %求mn阶矩阵U的五点差分矩阵 25 18 13 10 9 10 13 18 25 20 13 8 5 4 5 8 13 20 17 10 5 2 1 2 5 10 17 16 9 4 1 0 1 4 9 16v4 = 4 4 4 4 4 4 4 4 4产生一个二元函数偏导数和梯度。x=-2:0.2:2;y=-2:xx,yy=meshgrid（x,y）;z=