第四章 MATLAB的数值计算功能 内容参考

第四章 MATLAB 的数值计算功能一 多项式1多项式的表达与创建Matlab用矢量表达多项式系数，元素按降幂排列：P(x)=a0xn+a1xn-1+a2xn-2an-1x+a0其系数矢量为：P=a0 a1 an-1 an如将根矢量表示为：ar= ar1 ar2 arn则根矢量与系数矢量之间关系为：（x-ar1）(x- ar2) (x- arn)= a0xn+a1xn-1+a2xn-2an-1x+a0多项式系数矢量可通过调用函数 p=poly(ar)产生例 1：由根矢量创建多项式。将多项式(x-6)(x-3)(x-8)表示为系数形式 a=6 3 8pa=poly(a) %求系数矢量ppa=poly2sym(pa) %以符号形式表示原多项式ezplot(ppa,-50,50)pa = 1 -17 90 -144ppa =x3-17*x2+90*x-144注：（1）根矢量元素为n ，则多项式系数矢量元素为n+1；2）函数poly2sym(pa) 把多项式系数矢量表达成符号形式的多项式，缺省情况下自变量符号为x，可以指定自变量。（3）使用简单绘图函数可以直接绘制 符号形式多项式的曲线。例 2： 求三阶方阵A的特征多项式系数，并转换为多项式形式。a=6 3 8;7 5 6; 1 3 5Pa=poly(a) %求矩阵的特征多项式系数矢量Ppa=poly2sym(pa)Pa = 1.0000 -16.0000 38.0000 -83.0000Ppa =x3-17*x2+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.0000 1.1000 0.5500 0.1250pr = 1.0000 1.1000 0.5500 0.1250ppr =x3+11/10*x2+11/20*x+1/8注：（1）要形成实系数多项式，根矢量中的复数根必须共轭成对； （2）含复数根的根矢量所创建的多项式系数矢量中，可能带有很小的虚部，此时可采用取实部的命令(real)把虚部滤掉。例 4： 将多项式的系数表示形式转换为根表现形式，poly和 roots 互为逆函数。求 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+2s+2)(s+4)(s+1) （多个多项式相乘）c=conv(1,2,2,conv(1,4,1,1)cs=poly2sym(c,s) （指定变量为s）c = 1 7 16 18 8cs =s4+7*s3+16*s2+18*s+8例2：求多项式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 = 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称为极点，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)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=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，“”较多用，两者的关系为：(b/a)=(ab)系数矩阵a可能是m行n列的，有三种情况：*方阵系统: m=n 可求出精确解（a必须是非奇异，即满秩）*超定系统 ：mn 可求出最小二乘解*欠定系统： m求矩阵元素的最大值、最小值可用： 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=corrcoef(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 = 1.0000 -0.2561 -0.2561 1.00002. 微分与梯度 (表43)例1：按列求微分。x=1,10,20;2,12,23;3,14,26;3,16,29d=diff(x) %求一阶微分x = 1 10 20 2 12 23 3 14 26 3 16 29d = 1 2 3 1 2 30 2 3例2： 对于（u=x2+y2和2=4）求5点差分。x,y=meshgrid(-4:4,-3:3);u=x.2+y.2v4=4*del2(u) %求mn阶矩阵U的五点差分矩阵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 16 17 10 5 2 1 2 5 10 17 20 13 8 5 4 5 8 13 20 25 18 13 10 9 10 13 18 25v4 = 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4例3：产生一个二元函数偏导数和梯度。x=-2:0.2:2;y=-2:0.2:2;xx,yy=meshgrid(x,y);z=xx.*exp(-xx.2-yy.2);Gx,Gy=gradient(z,0.2,0.2); %Gx,Gy分别是二元函数的偏导contour(x,y,z,k),hold on,quiver(xx,yy,Gx,Gy,r),hold off（FX,FY = GRADIENT(F,HX,HY), when F is 2-D, uses the spacing specified by HX and HY. HX and HY can either be scalars to specify the spacing between coordinates or vectors to specify the coordinates of the points. If HX and HY are vectors, their length must match the corresponding dimension of F.）（QUIVER(X,Y,U,V) plots velocity vectors as arrows with components (u,v) at the points (x,y). The matrices X,Y,U,V must all be the same size and contain corresponding position and velocity components (X and Y can also be vectors to specify a uniform grid). QUIVER automatically scales the arrows to fit within the grid.）六 插值: 在已知数据之间计算估计值的过程。1 一维插值 由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 %将图标定义为单元数组for i=1:4 yi=interp1(x,y,xi,strmodi);subplot(2,2,i)plot(x,y, ro ,xi,yi, b),xlabel(str1b(i)endstrmod = 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 二维插值用于图形图象处理和三维曲线拟合等领域，由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(1), mesh(x,y,z)xi,yi=meshgrid(-3:0.2:3,-3:0.2:3);z1=interp2(x,y,z,xi,yi,nearest);z2=interp2(x,y,z,xi,yi,linear);z3=interp2(x,y,z,xi,yi,spline);z4=interp2(x,y,z,xi,yi,cubic);figure(2), subplot(2,2,1) mesh(xi,yi,z1) title(nearest)subplot(2,2,2) mesh(xi,yi,z2) title(linear)subplot(2,2,3) mesh(xi,yi,z3) title(spiine)subplot(2,2,4) mesh(xi,yi,z4) title(cubic)3 多维插值: 包括三维插值函数interp3和多维插值函数interpn,函数调用格式与一、二维插值基本相同。VI=interp3(X,Y,Z,V,XI,YI,ZI,method)其中： X,Y,Z自变量组成的数组；V三维函数数组XI,YI,ZI插值点的自变量数组Method -插值方法选项。（FLOW A simple function of 3 variables. FLOW, a function of three variables, is the speed profile of a submerged jet within a infinite tank. FLOW is useful for demonstrating SLICE and INTERP3.）(SLICE(X,Y,Z,V,XI,YI,ZI) draws slices through the volume V along the surface defined by the arrays XI,YI,ZI.) (SHADING FLAT sets the shading of the current graph to flat. SHADING INTERP sets the shading to interpolated. SHADING FACETED sets the shading to faceted, which is the default.)例：三维插值实例。x,y,z,v=flow(10);figure(1),slice(x,y,z,v,6 9.5,2,-2 .2)xi,yi,zi=meshgrid(.1:.25:10,-3:.25:3,-3:.25:3);vi=interp3(x,y,z,v,xi,yi,zi);figure(2),slice(xi,yi,zi,vi,6 9.5,2,-2 .2)shading flat七 数字信号处理初步MATLAB主要包括信号处理（Signal Processing Toolbox）和滤波器设计(Filter Design Toolbox)两部分.基本平台中提供了一些常用的信号处理函数（表 4-4），可实现数字信号处理的基本功能。1 快速傅立叶变换(FFT): MATLAB6.1 采用了新的快速傅立叶变换计算方法，速度高，可以作到实时处理。 函数fft的调用格式： *Y=fft(X) 返回应用快速傅立叶方法计算得到的矢量X的离散傅立叶变换（DFT）， 如果 X为矩阵，fft返回矩阵每一列的傅立叶变换，如果X为多维数组，fft运算从第一个非独立维开始执行。*Y=fft(X,n) 返回n点的离散傅立叶变换，如果X的长度小于n，X中补0使其与n的长度相同，；如果X的长度大于n，则X的多出部分将被删除；如X为矩阵，用同样方法处理矩阵列的长度。*Y=fft(X,dim) 和Y=fft(X,n,dim)沿dim维进行FFT操作。注：快速傅立叶变换的结果为复数。例： 产生一个正弦衰减曲线，进行快速傅立叶变换，并画出幅值图，相位图、实部图和虚部图。tp=0:2048; %时域数据点数Nyt=sin(0.08*pi*tp).*exp(-tp/80); figure(1), plot(tp,yt), axis(0,400,-1,1),t=0:800/2048:800; %频域点数Nff=0:1.25:1000;yf=fft(yt)；ya=abs(yf(1:801); %幅值yp=angle(yf(1:801)*180/pi; %相位yr=real(yf(1:801); %实部yi=imag(yf(1:801); %虚部figure(2),subplot(2,2,1)plot(f,ya),axis(0,200,0,60)subplot(2,2,2)plot(f,yp),axis(0,200,-200,10)subplot(2,2,3)plot(f,yr),axis(0,200,-40,40)subplot(2,2,4)plot(f,yi),axis(0,200,-60,10)2. 快速傅立叶变换的长度与运算速度 使用fft函数时，可输入第二个参数n以指定变换点的数量：y=fft(x,n)fft的运算速度取决于变换的长度。*如n是2 的整数次幂，则运算速度最快；*如n是合数，fft采用质因数分解的算法，速度取决于质因数的大小。*如n是质数，计算速度最慢。例： 创建700001阶的随机矢量x，取快速傅立叶变换的长度分别为质数65539、 2的16次方， 两个合数66000和5535，分别计算使用这些长度fft所占用的cpu时间。x=rand(70000,1);isprime(65539) %质数ans = 1216ans = 65536（FACTOR(N) returns a vector containing the prime factors of N.）factor(66000) %合数ans = 2 2 2 2 3 5 5 5 11factor(65535) %合数ans = 3 5 17 257 t=cputime; %开始运行时间y=fft(x,65539);e=cputime-t %运行结束与开始运行时的时间差 e= 0. 9900t=cputime;y=fft(x,65536);e=cputime-te= 0. 1600 t=cputime;y=fft(x,65535);e=cputime-te =0. 3800由上例可见，变换长度的差别对运算速度的影响很大，这正是大多数的数字信号处理系统中FFT的长度取为2 的n次方的原因。第五章 符号数学基础一 符号对象的创建1 创建符号变量和表达式创建符号变量和表达式的两个基本函数：sym, syms*x=sym(x) 创建一个符号变量x，可以是字符、字符串、表达式或字符表达式。*syms用于方便地一次创建多个符号变量，调用格式为： syms a b c d .例1： 使用sym函数创建符号变量.a=sym(a)b=sym( hello)c=sym( (1+sqrt(5)/2)y=sym( x3+5*x2+12*x+20)a =ab = helloC = (1+sqrt(5)/2 Y =x3+5*x2+12*x+20例2：用syms函数创建符号变量。syms a b c d书写简洁意义清楚，建议使用。2 创建符号矩阵 使用sym 和syms函数可以创建符号矩阵，可以直接输入或从数值矩阵转换。例1：创建一个循环矩阵。syms a b c dn=a b c d;b c d a;c d a b;d a b cn = a, b, c, d b, c, d, a c, d, a, b d, a, b