第四章数值计算

第四章数值计算功能

数值计算能力是MATLAB称雄世界的根本柱石。本章将以简要明了的方式，分述线性方程组的解、特征值问题的解、多项式和卷积、数据分析、泛函指令的使用、信号处理和系统分析等内容。4.1线性方程组的解4.1.1LU分解、行列式、逆和恰定方程的解（1）LU分解矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。线性代数中已经证明，只要方阵A是非奇异的，LU分解总是可以进行的。MATLAB提供的lu函数用于对矩阵进行LU分解，其调用格式为：[L,U]=lu(X)产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，使之满足X=LU。注意，这里的矩阵X必须是方阵。[L,U,P]=lu(X)

产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PX=LU。当然矩阵X同样必须是方阵。实现LU分解后，线性方程组Ax=b的解x=U\(L\b)或x=U\(L\Pb)，这样可以大大提高运算速度。（2）行列式命令：det（A）求矩阵A的行列式（3）逆命令：inv（A）求矩阵A的逆（4）除法指令求方程的解命令：x=A\b对恰定方程，采用LU法执行

注：此命令中“\”是“左除”符号，当矩阵A条件数很大时，机器会给出警告信息，提醒用户。

【例】“求逆”法和“左除”法解恰定方程的性能对比（1）为对比这两种方法的性能，先用以下指令构造一个条件数很大的高阶恰定方程randn('state',0);A=gallery('randsvd',100,2e13,2); x=ones(100,1); b=A*x; cond(A)ans=1.9990e+013（2）“求逆”法恰定方程的误差、残差、运算次数和所有时间tic xi=inv(A)*b; ti=toc eri=norm(x-xi) rei=norm(A*xi-b)/norm(b)

ti=0.4400eri=0.0469rei=0.0047（3）“左除”法解恰定方程的误差、残差、运算次数和所有时间tic;xd=A\b; td=toc,erd=norm(x-xd),red=norm(A*xd-b)/norm(b)td=0.0600erd=0.0078red=2.6829e-015计算结果表明：除法求解不但速度快，而且精度高得多。所以，在解方程时，尽量不要使用指令inv（A）*b。4.1.2奇异值分解和矩阵结构对于任意矩阵A，存在酉阵U，V，使得UTAV等于一个对角阵，那么我们就把这种分解称为奇异值分解。相应的MATLAB指令如下：

[U,S,V]=svd(A)给出矩阵A的奇异值分解三对组阵，使A=USVT

r=rank(A,tol)大于阈值tol的奇异值为矩阵A的“数值”秩norm(A,flag)计算矩阵A的范数，范数类型由flag指定（可取2，1等）cond(A)根据cond（A）σ1/σp计算条件数Pinv(A,tol)在指定阈值tol下，求矩阵A的广义逆null(A)求取矩阵A零空间的酉矩阵orth(A)求取矩阵A值空间的酉矩阵subspace(A,B)计算两矩阵之间的夹角4.1.3线性二乘问题的解对于线性模型y=Ax+η，式中η为服从正态分布N(0,1)的白噪音，求该超定方程最小二乘解有三种途径：(1)正则方程法得解x=(ATA)-1ATb(2)广义逆法得解x=A+b(3)用矩阵除法得解x=A\b【例4.1-2】对于超定方程，进行三种解法比较。其中取MATLAB库中的特殊函数生成。（1）生成矩阵A及y，并用三种方法求解A=gallery(5);A(:,1)=[];y=[1.77.56.30.83-0.082]';x=inv(A'*A)*A'*y,xx=pinv(A)*y,xxx=A\yWarning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=7.087751e-018.x=3.48115.15950.9534-0.0466xx=3.47595.19480.7121-0.1101Warning:Rankdeficient,rank=3tol=1.0829e-010.xxx=3.46055.29870-0.2974（2）计算三个解的范数nx=norm(x),nxx=norm(xx),nxxx=norm(xxx)

nx=6.2968e+000nxx=6.2918e+000nxxx=6.3356e+000（3）比较三种解法的方程误差e=norm(y-A*x),ee=norm(y-A*xx),eee=norm(y-A*xxx)

e=0.6985ee=0.0474eee=0.0474

4.2特征值分解和矩阵函数4.2.1特征值分解和矩阵函数求解Ax=λx的精良数值计算方法是:先对矩阵A实施一系列的Householder变换,产生一个准上三角阵(即对角线下还有一个非零的次对角线),然后再运用著名的QR法迭代使准上三角阵对角化.涉及矩阵特征值分解的常用MATLAB指令如下:d=eig(A)仅计算矩阵A的特征值(以向量形式d存放)[V,D]=eig(A)计算A的特征向量阵V和特征值对角阵D，使A*V=V*D成立[VR,DR]=cdf2rdf(VC,DC)把复数对角形转换成实数块对角形[VC,DC]=rsf2csf(VR,DR)把实数块对角形转换成复数对角形[V,J]=jordan(A)Jordan分解使A*V=V*J，J是块对角阵C=condeig(A)向量c中包含矩阵A关于各特征值条件数[V,D,c]=condeig(A)相当于[V，D]=eig(A)和c=condeig(A)两条指令的组合【例1】简单实阵的特征值问题。A=[1,-3;2,2/3];[V,D]=eig(A)

V=0.77460.77460.0430-0.6310i0.0430+0.6310iD=0.8333+2.4438i000.8333-2.4438i【例2】把例1中的复数特征值对角阵D转换成实数块对角阵，使VR*DR/VR=A。[VR,DR]=cdf2rdf(V,D)

VR=0.774600.0430-0.6310DR=0.83332.4438-2.44380.8333【例3】对亏损矩阵进行Jordan分解。A=gallery(5) [VJ,DJ]=jordan(A); [V,D,c_eig]=condeig(A);c_equ=cond(A);DJ,D,c_eig,c_equA=-911-2163-25270-69141-4211684-575575-11493451-138013891-38917782-23345933651024-10242048-614424572DJ=0100000100000100000100000D=-0.040800000-0.0119+0.0386i00000-0.0119-0.0386i000000.0323+0.0230i000000.0323-0.0230ic_eig=1.0e+010*2.12932.07962.07962.00202.0020c_equ=2.0253e+0184.2.2矩阵的谱分解和矩阵函数

(特征根各异的)矩阵A的特征值分解可重写为：据此定义矩阵函数。

常见矩阵函数的数学定义和指令表数学表达式指令形式功用A^p据谱分解，求ApP^A据谱分解，求pAexpm(A)一般常用expm1(A)用Pade近似求eAexpm2(A)用Taylor近似求eA，适用任何方差，精度稍差expm3(A)用特征值分解求eA，仅对非亏损阵使用sqrtm(A)据Schur分解，求。多解时，仅给出一个解logm(A)据谱分解，求lnAfunm(A,’FN’)FN是函数名(如sin等)【例1】数组乘方与矩阵乘方的比较。clear,A=[123;456;789];A_Ap=A.^0.3 A_Mp=A^0.3 A_Ap=1.00001.23111.39041.51571.62071.71181.79281.86611.9332A_Mp=0.6962+0.6032i0.4358+0.1636i0.1755-0.2759i0.6325+0.0666i0.7309+0.0181i0.8292-0.0305i0.5688-0.4700i1.0259-0.1275i1.4830+0.2150i【例2】标量的数组乘方和矩阵乘方的比较。（A取自例1）pA_A=(0.3).^A pA_M=(0.3)^A

pA_A=0.30000.09000.02700.00810.00240.00070.00020.00010.0000pA_M=2.93420.4175-1.0993-0.02780.7495-0.4731-1.9898-0.91841.1531【例3】sin的数组运算和矩阵运算比较。（A取自例1）A_sinA=sin(A) A_sinM=funm(A,'sin')

A_sinA=0.84150.90930.1411-0.7568-0.9589-0.27940.65700.98940.4121A_sinM=-0.6928-0.23060.2316-0.1724-0.1434-0.11430.3479-0.0561-0.46024.3多项式和卷积4.3.1多项式一、多项式表达方式的约定MATLAB约定降幂多项式用以下系数行向量表示：即把多项式的各项系数依降幂次序排放在行向量的元素位置上。其中，多项式中缺某幂次项，则应认为该幂次项的系数为零。二、多项式运算函数指令含义P=conv(p1,p2)p是多项式p1和p2的乘积多项式[q,r]=deconv(p1,p2)多项式p1被p2除的商多项式为q，而余多项式为rP=poly(AR)求方阵AR的特征多项式p；或求向量AR指定根所对应的多项式P=polyfit(x,y,n)求x，y向量给定的数据的n阶拟合多项式ppA=polyval(p,S)按数组运算规则计算多项式的值。P为多项式，S为矩阵PM=polyvalm(p,S)按矩阵运算规则计算多项式的值。P为多项式，S为矩阵[r,p,k]=residue(b,a)部分分式展开：b，a分别式分子分母多项式系数向量；r，p，k分别是留数、极点、直项R=roots(p)求多项式P的根例求的“商”及“余”多项式。p1=conv([1,0,2],conv([1,4],[1,1]));%计算分子多项式p2=[1011]; %注意缺项补零 [q,r]=deconv(p1,p2);cq='商多项式为';cr='余多项式为';disp([cq,poly2str(q,'s')]),disp([cr,poly2str(r,'s')])商多项式为s+5余多项式为5s^2+4s+3例求3阶方阵A的特征多项式。A=[111213;141516;171819];PA=poly(A)%A的特征多项式PPA=poly2str(PA,‘s’)%比较习惯的方式显示多项

PA=1.0000-45.0000-18.0000-0.0000PPA=s^3-45s^2-18s-2.8387e-015

例由给定根向量求多项式系数向量。R=[-0.5,-0.3+0.4*i,-0.3-0.4*i];%根向量P=poly(R)%R的特征多项式 PR=real(P)%求P的实部 PPR=poly2str(PR,'x')P=1.00001.10000.55000.1250PR=1.00001.10000.55000.1250PPR=x^3+1.1x^2+0.55x+0.125注：含复数的根向量所生成的多项式（如P）的系数有可能带有截断误差数量级的虚部，此时可采用实部的指令“real”把那很小的虚部滤掉。

例两种多项式求值指令的差别。S=pascal(4) %生成一个4阶方阵 P=poly(S);PP=poly2str(P,'s')PA=polyval(P,S) %独立变量取数组S元素时的多项式PM=polyvalm(P,S)%独立变量取矩阵S时的多项式

S=1111123413610141020PP=s^4-29s^3+72s^2-29s+1PA=1.0e+004*0.00160.00160.00160.00160.00160.0015-0.0140-0.05630.0016-0.0140-0.2549-1.20890.0016-0.0563-1.2089-4.3779PM=1.0e-010*0.00160.00330.00900.02050.00450.01010.02860.06970.00950.02100.06530.15960.01630.03870.12260.3019从理论上讲，PM应该为0，这是著名的“Caylay-Hamilton”定理：任何一个矩阵满足他自己的特征多项式方程.本例中的PM很小，这是截断误差造成的。例部分分式展开。a=[1,3,4,2,7,2]; %分母多项式系数向量b=[3,2,5,4,6]; %分子多项式系数向量[r,s,k]=residue(b,a)r=1.1274+1.1513i1.1274-1.1513i-0.0232-0.0722i-0.0232+0.0722i0.7916s=-1.7680+1.2673i-1.7680-1.2673i0.4176+1.1130i0.4176-1.1130i-0.2991k=[]三拟合和插值例对于给定数据对x0,y0，求拟合三阶多项式，并图示拟合情况。%给定数据对x0=0:0.1:1;y0=[-.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];%求拟合多项式系数n=3;P=polyfit(x0,y0,n)%图示拟合情况xx=0:0.01:1;yy=polyval(P,xx);plot(xx,yy,'-b',x0,y0,'.r','MarkerSize',20),xlabel('x')P=56.6915-87.117440.0070-0.9043采用三次多项式所得的拟合曲线

例以上例所给数据，研究一维插值，并观察插值与拟合的区别。x0=0:0.1:1;%给出数据对y0=[.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];

xi=0:0.02:1;%采用三次多项式进行插值yi=interp1(x0,y0,xi,'cubic');

plot(xi,yi,‘-b’,x0,y0,‘.r’,‘MarkerSize’,20),xlabel(‘x’)%绘图通过三次多项式插值所得的曲线

拟合和插值的区别：曲线拟和研究如何寻找“平滑”曲线最好地表现带噪声的“测量数据”，但是不要求拟合曲线穿过这些“测量数据”点。插值是在认定所给出“基准数据”完全正确的情况下，研究如何“平滑”地估算出“基准数据”之间其它点的函数值，因此插值所的曲线一定穿过“基准数据”。4.3.2卷积一离散序列的数值卷积设有如下两个任意序列，那么该卷积为它可准确计算的非平凡区间是[Nc1,Nc2].在此，Nc1=N1+N3，Nc2=min{N1+N4，N2+N3}二MATLAB的“卷积”指令1对序列A、B进行的卷积运算C=conv（A，B）2解卷运算，使B=conv（A，Q）+R成立[Q，R]=dconv（B，A）【例4.3-8】有序列和。（A）求这两个完整序列的卷积，并图示。（B）假设序列中最后4个非零值未知，而成为截尾序列，求卷积并图示。a=ones(1,10);n1=3;n2=12; %完整A（n）序列的非平凡值和区间端点b=ones(1,8);n3=2;n4=9;%完整B（n）序列的非平凡值和区间端点c=conv(a,b);nc1=n1+n3;nc2=n2+n4;%计算卷积和确定卷积非平凡区间端点 kc=nc1:nc2;%构成非平凡区间的序号自变量 %截尾序列卷积 aa=a(1:6);nn1=3;nn2=8; %截尾A（n）序列的非平凡值和区间端点cc=conv(aa,b);ncc1=nn1+n3;nx=nn2+n4; %“非平凡”区间右端点 ncc2=min(nn1+n4,nn2+n3);%截尾序列卷积被正确计算区间的右端点kx=(ncc2+1):nx;kcc=ncc1:ncc2;N=length(kcc);stem(kcc,cc(1:N),'r','filled')axis([nc1-2,nc2+2,0,10]),grid,holdonstem(kc,c,'b'),stem(kx,cc(N+1:end),'g'),holdoff

“完整”序列卷积和“截尾”序列卷积

对于“完整”序列，指令conv所给出的结果在整个非平凡区间上都是准确的，如本例中“蓝空心杆图”所示；对于“非完整”的截尾序列，指令conv给出的结果，只有部分非平凡区间上的结果是准确的，即本例的“红实心杆图”，而“绿空心杆图”虽也是由“截尾”序列算得的“非平凡”值，但不是真正的卷积。4.4数据分析函数MATLAB进行数据分析时的约定：（1）若输入宗量X是向量，那么不管是行向量还是列向量，运算是对整个向量进行的。（2）若输入宗量X是二维数组（或称矩阵），那么指令运算是按列进行的。即默认每个列是由一个变量的不同“观察”所得的数据组成。4.4.1随机数发生器和统计分析指令rand(n,m)

%产生（n×m）维的“[0，1]区间”均匀分布随机数组randn(n,m)%产生(n×m)维的“均值为0标准差为1”的正态分布随机数组min(X)

%对X矩阵各列分别求最小值max(X)

%对X矩阵各列分别求最大值median(X)

%对X矩阵各列分别求中位数mean(X)

%对X矩阵各列分别求均值std(X)

%对X矩阵各列分别求标准差var(X)

%对X矩阵各列分别求方差C=cov(X)

%对X矩阵各列间的协方差阵P=corrcoef(X)

%给出矩阵X各列间的相关数例基本统计示例。randn(‘state’,0) %为随机仿真可重复而设置 A=randn(1000,4);AMAX=max(A),AMIN=min(A)%因分别在+3，-3附近AMED=median(A) %应接近0 AMEAN=mean(A) %应接近0 ASTD=std(A) %应接近1

AMAX=2.73163.20253.41283.0868AMIN=-2.6442-3.0737-3.5027-3.0461AMED=-0.01310.05960.01220.0459AMEAN=-0.04310.04550.01770.0263ASTD=0.94351.03131.02480.9913例cov和corrcoef的使用示例。rand('state',1)X=rand(10,3);Y=rand(10,3);mx=mean(X);Xmx=X-ones(size(X))*diag(mx);CCX=Xmx‘*Xmx/(size(Xmx,1)-1) %它于CX相同CX=cov(X),CY=cov(Y)Cxy=cov(X,Y) %相当于cov(X(:),Y(:)) PX=corrcoef(X)Pxy=corrcoef(X,Y) CCX=0.08710.0242-0.00730.02420.08460.0056-0.00730.00560.0607CX=0.08710.0242-0.00730.02420.08460.0056-0.00730.00560.0607CY=0.07210.00130.01650.00130.0714-0.05350.0165-0.05350.0720Cxy=0.0761-0.0012-0.00120.0708PX=1.00000.2819-0.10050.28191.00000.0785-0.10050.07851.0000Pxy=1.0000-0.0168-0.01681.0000

4.4.2差分和累计指令del2(U,hx,hy)五点Laplacian，近似给出diff(X,m,n)沿第n维求m差分[DZx,DZy]=gradient(Z,dx,dy)对Z求x,y方向的梯度prod(X,n)沿第n维求乘积sum(X,n)沿第n维求和trapz(x,Y,n)采用梯形法沿第n维求函数Y关于资变量x的积分cumprod(X,n)沿第n维求累计乘积cumsum(X,n)沿第n为求累计和cumtrapz(x,Y,n)采用梯形法沿第n维求函数Y关于自变能量x的累计积分[XS,KK]=sort(X,n)沿第n维对X元素按模增大排列。KK给出XS元素的位置

例用一个简单矩阵表现diff和gradient指令计算方式。F=[1,2,3;4,5,6;7,8,9]Dx=diff(F) ％相邻行差分 Dx_2=diff(F,1,2) ％相邻列差分，第三输入宗量2表示“列”差分 [FX,FY]=gradient(F) ％数据点步长默认为1 [FX_2,FY_2]=gradient(F,0.5) ％数据点步长为0.5F=123456789Dx=333333Dx_2=111111FX=111111111FY=333333333FX_2=222222222FY_2=666666666

例函数的梯度和，用数值计算验证，并图示。

clear,dd=0.2;x=-1:dd:1;y=x;[X,Y]=meshgrid(x,y);Z=(X.^2)+(Y.^2);[DZx,DZy]=gradient(Z,dd,dd);DZ2=4*del2(Z,dd,dd);％与理论计算比较DDZx0=DZx-2*X;DDZy0=DZy-2*Y;DDZ20=DZ2-4;subplot(1,3,1),stem3(X,Y,DDZx0)subplot(1,3,2),stem3(X,Y,DDZy0)subplot(1,3,3),stem3(X,Y,DDZ20)axis([-1,1,-1,1,-0.4,0.4])xlabel('x'),ylabel('y')

理论计算和数值计算的差别图示

可见，数值梯度与理论梯度的区别仅发生在区间的端部；而Laplacian的数值结果与理论值几乎没有差别。例求积分,其中,dt=0.1;t=(0:dt:10)';y=exp(-0.8*t.*abs(sin(t)));ss=dt*cumsum(y); %矩形法求积 ss10=dt*sum(y);ssend=ss(end);%sum所给结果就是ss(end)st=cumtrapz(t,y); %梯形法求积 st10=trapz(t,y);stend=st(end);%trapz所给结果就是st(end)disp([blanks(5),'sum',blanks(6),'cumsum',blanks(4),'trapz',blanks(5),'cumtrapz'])disp([ss10,ssend,st10,stend])plot(t,y,'b:',t,ss,'r-',t,st,'r.')legend('y(x)','cunsum','cumtrapz',0)

sumcumsumtrapzcumtrapz2.70822.70822.65762.6576矩形法和梯形法求积比较4.5MATLAB泛函指令在MATLAB中，凡以“函数”为输入宗量的指令，都被统称为MATLAB泛函指令。最常见的有：

z=fzero(fun,x0)求一元函数零点的指令的最简格式

x=fsolve(fun,x0)解非线性方程组的最简单格式

x=fminbnd(fun,x1,x2)求函数在区间（x1,x2）中极小值的指令最简格式

x=fminsearch(fun,x0)单纯形法求多元函数极值点指令的最简格式

x=fminunc(fun,x0)拟牛顿法求多元函数极值点指令的最简格式

a=lsqnonlin(fun,a0)解非线性最小二乘问题指令的最简格式

q=quad(fun,a,b)采用递推自适应Simpson法计算积分

q=quadl(fun,a,b)采用递推自适应Lobatto法求数值积分

ss=dblquad(fun,inmin,inmax,outmin,outmax)

二重（闭型）数值积分指令[t,YY]=ode45(fun,tspan,Y0)采用4、5阶Runge-Kutta方程解算ODE初值问题4.5.1对于任意函数f(x)=0来说，可能有零点，也可能没有；可能有一个零点，也可能有多个甚至无数多个零点。因此，很难说出一个通用解法。解题步骤如下：（1）利用MATLAB作图指令获取初步近似解具体做法：现确定一个零点可能存在的自变量区间；然后利用plot指令画出f(x)在该区间中的图形；用眼观察f(x)与横轴的交点坐标，或者更细致些，用zoom对交点处进行局部放大再读数。借助ginput指令获得更精确些的交点坐标值。（2）利用fzero指令求精确界z=fzero(fun,x0)求一元函数零点指令的简单格式[z,f_z,exitflag]=fzero(fun,x0,options,P1,P2,…)求一元函数零点指令的最完整格式注1由于fzero是根据函数是否穿越横轴来决定零点的，因此本程序无法确定函数曲线仅触及横轴和不穿越的零点。注2第二输入宗量x0是表示零点初时猜测。注3输入宗量options是优化迭代所采用的参数选项。注4输入宗量P1，P2是向函数fun传递的参数。注5指令中，输出宗量z是所求零点的自变量值。注6第二个输出宗量f_z是函数值；第三个宗量是表示程序终止计算的条件。若exitflag>0，则表明找到零点后退出。例通过求的零点，综合叙述相关指令的用法。（1）使用字符串表示被处理函数P1=0.1;P2=0.5; ％按泛函指令要求，这里参数必须用P1,P2,…表示y_C=‘sin(x).^2.*exp(-P1*x)-P2*abs(x)’;％这里自变量必须用x表示（2）作图法观察函数零点分布x=-10:0.01:10;％对自变量采样，采样步长不宜过大Y=eval(y_C);％在采样点上计算函数值clf,plot(x,Y,‘r’);holdon,plot(x,zeros(size(x)),‘k’);％画坐标横轴xlabel('t');ylabel('y(t)'),holdoff

函数零点分布观察图（3）利用zoom和ginput指令获得零点的初始近似值zoomon %在MATLAB指令窗中运行，获局部放大图[tt,yy]=ginput(5);%在MATLAB指令窗中运行，用鼠标获5个零点猜测值zoomoff

tt%显示所得零点初始猜测值tt=-2.0046-0.5300-0.01150.61061.6590（4）求近似tt(4)的精确零点[t4,y4,exitflag]=fzero(y_C,tt(4),[],P1,P2)

Zerofoundintheinterval:[0.59333,0.6106].t4=0.5993y4=0exitflag=1局部放大和利用鼠标取值图

图中的十字是ginput运行后产生的取值图符4.5.2求函数极值点

许多科学研究和工程计算问题都可以归结为一个极值问题，而且所有的求极大值问题都可以化成求极小值的问题，所以MATLAB中只有处理极小值的指令求函数极值的三条指令：[f,fval,exitflag,output]=fminbnd(fun,x1,x2,options,P1,P2,…)求一元函数在区间（x1,x2）中极小值的指令最完整格式[f,fval,exitflag,output]=fminsearch(fun,x0,options,P1,P2,…)单纯形法求多元函数极值点指令的最完整格式[f,fval,exitflag,output,grad,hessian]=fminunc(fun,x0,options,P1,P2,…)拟牛顿法求多元函数极值点指令的最完整格式注：x1,x2为被研究区间的左、右边界；x0为极值点位置的初始猜测例求的极小值点。它即是著名的Rosenbrock‘s“Banana”测试函数，它的理论极小值是。该测试函数有一片浅谷，许多算法难以越过此谷。（1）本例采用内联函数表示测试函数如下：ff=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2','x');（2）用单纯形法求极小值点x0=[-1.2,1];[sx,sfval,sexit,soutput]=fminsearch(ff,x0)sx=1.00001.0000sfval=8.1777e-010sexit=1soutput=iterations:85funcCount:159algorithm:'Nelder-Meadsimplexdirectsearch'4.5.3数值积分数值积分有闭型、开型算法之分。两者的区别在于：是否需要计算积分区间端点处的函数值。MATLAB中仅提供闭型数值积分指令：q=quad(fun,a,b,tol,trace,p1,p2,…)采用递推自适应Simpson法计算积分q=quadl(fun,a,b,tol,trace,p1,p2,…)采用递推自适应Lobatto法求数值积分ss=dblquad(fun,inmin,inmax,outmin,outmax,tol,method)二重（闭型）数值积分指令注：tol是一个标量，用来控制绝对误差，缺省时，积分的绝对精度为1×10－6

trace取非0值时，将随积分的进程逐点画出被积函数p1,p2是向被积函数输送的参数method是积分方法，如”quad”,”quadl”例求，其精确值为.（1）编写被积函数的M文件[FUN.M]functiony=fun(x)y=exp(-x.*x);（2）利用函数句柄运行MATLAB指令quad和quadl进行积分Hf=@fun;Isim=quad(Hf,0,1),IL=quadl(Hf,0,1)Isim=0.7468IL=0.7468

4.5.4解常微分方程以例题说明：求微分方程在初始条件情况下的解，并图示。

（1）把高阶微分方程改写成一阶微分方程组令，于是原二阶方程可改写成如下一阶方程组（2）根据上述一阶微分方程组编写M函数文件[DyDt.m]functionydot=DyDt(t,y)mu=2;ydot=[y(2);mu*(1-y(1)^2)*y(2)-y(1)];（3）解算微分方程tspan=[0,30];y0=[1;0];[tt,yy]=ode45(@DyDt,tspan,y0); plot(tt,yy(:,1)),title('x(t)')（4）画相平面图

plot(yy(:,1),yy(:,2))

微分方程解相平面轨迹4.6信号处理MATLAB有专门的工具包（如SignalProcessingToolbox,DSPToolbox等）解决信号处理中的各种计算和仿真问题。本节只介绍MATLAB基本组件中的三个信号处理重要指令fit,ifft,filter。4.6.1快速Fourier变换和逆变换

离散序列x(k)的离散Fourier变换和逆离散Fourier变换的“教科书”定义为：而MATLAB出于软件的考虑：序列或向量元素下标从1开始记述，而不是从0开始。因此，上述两式在MATLAB中相应的表达式为：MATLAB根据以上定义给出了如下指令：

计算N点Xt序列的N点离散Fourier变换Xf

Xf=fft(Xt,N,DIM)计算N点Xt序列的N点离散Fourier逆变换Xt

Xt=ifft(Xt,N,DIM)注：使用N的目的：（1）使被变换序列的长度恰为2，以使指令fft较高速的运作。（2）当N大于被变换序列长度时，原序列的尾部将用0填补，使其长度为N后，再做变换。（3）缺省时，指令算法将以被变换序列的长度为N。例利用fft和ifft指令重新计算两序列的卷积。所给已知序列为

（1）产生给定序列a(n),b(n)cleara=ones(1,13);a([1,2,3])=0; ％a时间序列应是0时刻算起b=ones(1,10);b([1,2])=0;％b时间序列应是0时刻算起（2）直接卷积c=conv(a,b);

（3）通过变换求卷积M=32;AF=fft(a,M);BF=fft(b,M);CF=AF.*BF; ％使用点乘 cc=real(ifft(CF));

％过滤掉由于截断误差引起的虚部（4）图式计算结果nn=0:(M-1); ％保证绘图时，时间从0时刻开始c(M)=0; ％使直接卷积所得序列通过尾部补0，与cc长度相同error=c-cc; ％直接法和变换法所得卷积之差subplot(2,1,1),stem(nn,c,'fill'),grid,axis([0,31,0,9])xlabel('nn'),ylabel('cc')subplot(2,1,2),stem(nn,error,'fill'),axis([0,31,-1,1])ylabel('error')

变换法和直接法求卷积结果比较

例fft在信号分析中的应用。试用频谱分析方法从受噪声污染的信号中鉴别出有用信号。在此，。(1)构遭受噪音污染的信号clear,randn('state',0)t=linspace(0,10,512);y=3*sin(5*t)-6*cos(9*t)+5*randn(size(t)); plot(t,y)受噪声污染的信号(2)画幅频曲线Y=fft(y);Ts=t(2)-t(1) ％时间信号的采样周期 Ws=2*pi/Ts; ％时间信号的采样频率Wn=Ws/2 ％Nyquest频率w=linspace(0,Wn,length(t)/2);％半采样频率中相应的刻度Ya=abs(Y(1:length(t)/2));plot(w,Ya)Ts=0.0196Wn=160.5354受噪声污染信号的幅频谱(3)对感兴趣频段局部放大ii=find(w<=20);plot(w(ii),Ya(ii))grid,xlabel('FrequencyRad/s')受噪声污染信号幅频谱的局部放大4.6.2数字滤波命令：y=filter[B,A,x,[],dim]运用数字滤波对信号x进行滤波注：1)x可以是一维数组，当x是二维数组时，默认设置下，滤波沿列进行。当x是高维数组时，要用dim指定滤波沿哪维进行。2)该指令的第三和四个输入宗量可以缺省例设计一个低通滤波器，从受噪声干扰的多频率混合信号中获取10Hz的信号。在此，而。clear,randn('state',1)ws=1000;％采样频率t=0:1/ws:0.4; x=sin(2*pi*10*t)+cos(2*pi*100*t)+0.2*randn(size(t));％生成带噪音的多频率信号wn=ws/2; ％Nyquest频率 [B,A]=butter(10,30/wn); ％截至频率为30/wn的10阶ButterWorth低通滤波器y=filter(B,A,x); ％进行（初值为0的）滤波处理 plot(t,x,'b-',t,y,'r.','MarkerSize',10)legend('Input','Output',0)

注：butter指令在SignalProcessingToolbox工具包中。10阶Butterworth滤波器的滤波效果4.7系统分析4.7.1线性时不变对象LTIS_ss=ss(A,B,C,D)利用状态方程四对组生成LTI对象S_zpk=zpk(Z,P,K)利用零极点增益三对组生成LTI对象S_tf=tf(Num,Den)利用传递函数二对组生成LTI对象[A,B,C,D]=ssdata(S_Lti)从LTI对象获取状态四对组[Z,P,K]=zpkdata(S_Lti)从LTI对象获取零极点增益三对组[Num,Den]=tfdata(S_Lti)从LTI对象获取传递函数二对组注:(1)MATLAB中规定：LTI对象三个子类中，“状态方程子类”的优先级别最高，“零极点增益子类”次之，“传递函数子类”级别最低。(2)使用LTI对象的好处：①可以直接运用“方块图”见例2；②可以把MIMO多输入多输出系统作为一个整体进行处理。例1已知一个两输入两输出系统的状态方程四对组，据此建立LTI对象“状态方程子类”模型。然后再从此模型获取传递函数二对组。（1）建立LTI对象“状态方程子类”模型A=[0.50.50.7071;-0.5-0.50.7071;-6.364-0.7071-8];B=[004;101]';C=[0.7071-0.70710;00.5-0.5];D=0;S_ss=ss(A,B,C,D) ％建立状态方程类模型a=x1x2x3x10.50.50.7071x2-0.5-0.50.7071x3-6.364-0.7071-8b=u1u2x10