第四章數(shù)值計(jì)算

第四章數(shù)值計(jì)算功能

數(shù)值計(jì)算能力是MATLAB稱(chēng)雄世界的根本柱石。本章將以簡(jiǎn)要明了的方式，分述線(xiàn)性方程組的解、特征值問(wèn)題的解、多項(xiàng)式和卷積、數(shù)據(jù)分析、泛函指令的使用、信號(hào)處理和系統(tǒng)分析等內(nèi)容。4.1線(xiàn)性方程組的解4.1.1LU分解、行列式、逆和恰定方程的解（1）LU分解矩陣的LU分解就是將一個(gè)矩陣表示為一個(gè)交換下三角矩陣和一個(gè)上三角矩陣的乘積形式。線(xiàn)性代數(shù)中已經(jīng)證明，只要方陣A是非奇異的，LU分解總是可以進(jìn)行的。MATLAB提供的lu函數(shù)用于對(duì)矩陣進(jìn)行LU分解，其調(diào)用格式為：[L,U]=lu(X)產(chǎn)生一個(gè)上三角陣U和一個(gè)變換形式的下三角陣L(行交換)，使之滿(mǎn)足X=LU。注意，這里的矩陣X必須是方陣。[L,U,P]=lu(X)

產(chǎn)生一個(gè)上三角陣U和一個(gè)下三角陣L以及一個(gè)置換矩陣P，使之滿(mǎn)足PX=LU。當(dāng)然矩陣X同樣必須是方陣。實(shí)現(xiàn)LU分解后，線(xiàn)性方程組Ax=b的解x=U\(L\b)或x=U\(L\Pb)，這樣可以大大提高運(yùn)算速度。（2）行列式命令：det（A）求矩陣A的行列式（3）逆命令：inv（A）求矩陣A的逆（4）除法指令求方程的解命令：x=A\b對(duì)恰定方程，采用LU法執(zhí)行

注：此命令中“\”是“左除”符號(hào)，當(dāng)矩陣A條件數(shù)很大時(shí)，機(jī)器會(huì)給出警告信息，提醒用戶(hù)。

【例】“求逆”法和“左除”法解恰定方程的性能對(duì)比（1）為對(duì)比這兩種方法的性能，先用以下指令構(gòu)造一個(gè)條件數(shù)很大的高階恰定方程randn('state',0);A=gallery('randsvd',100,2e13,2); x=ones(100,1); b=A*x; cond(A)ans=1.9990e+013（2）“求逆”法恰定方程的誤差、殘差、運(yùn)算次數(shù)和所有時(shí)間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）“左除”法解恰定方程的誤差、殘差、運(yùn)算次數(shù)和所有時(shí)間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計(jì)算結(jié)果表明：除法求解不但速度快，而且精度高得多。所以，在解方程時(shí)，盡量不要使用指令inv（A）*b。4.1.2奇異值分解和矩陣結(jié)構(gòu)對(duì)于任意矩陣A，存在酉陣U，V，使得UTAV等于一個(gè)對(duì)角陣，那么我們就把這種分解稱(chēng)為奇異值分解。相應(yīng)的MATLAB指令如下：

[U,S,V]=svd(A)給出矩陣A的奇異值分解三對(duì)組陣，使A=USVT

r=rank(A,tol)大于閾值tol的奇異值為矩陣A的“數(shù)值”秩norm(A,flag)計(jì)算矩陣A的范數(shù)，范數(shù)類(lèi)型由flag指定（可取2，1等）cond(A)根據(jù)cond（A）σ1/σp計(jì)算條件數(shù)Pinv(A,tol)在指定閾值tol下，求矩陣A的廣義逆null(A)求取矩陣A零空間的酉矩陣orth(A)求取矩陣A值空間的酉矩陣subspace(A,B)計(jì)算兩矩陣之間的夾角4.1.3線(xiàn)性二乘問(wèn)題的解對(duì)于線(xiàn)性模型y=Ax+η，式中η為服從正態(tài)分布N(0,1)的白噪音，求該超定方程最小二乘解有三種途徑：(1)正則方程法得解x=(ATA)-1ATb(2)廣義逆法得解x=A+b(3)用矩陣除法得解x=A\b【例4.1-2】對(duì)于超定方程，進(jìn)行三種解法比較。其中取MATLAB庫(kù)中的特殊函數(shù)生成。（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）計(jì)算三個(gè)解的范數(shù)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特征值分解和矩陣函數(shù)4.2.1特征值分解和矩陣函數(shù)求解Ax=λx的精良數(shù)值計(jì)算方法是:先對(duì)矩陣A實(shí)施一系列的Householder變換,產(chǎn)生一個(gè)準(zhǔn)上三角陣(即對(duì)角線(xiàn)下還有一個(gè)非零的次對(duì)角線(xiàn)),然后再運(yùn)用著名的QR法迭代使準(zhǔn)上三角陣對(duì)角化.涉及矩陣特征值分解的常用MATLAB指令如下:d=eig(A)僅計(jì)算矩陣A的特征值(以向量形式d存放)[V,D]=eig(A)計(jì)算A的特征向量陣V和特征值對(duì)角陣D，使A*V=V*D成立[VR,DR]=cdf2rdf(VC,DC)把復(fù)數(shù)對(duì)角形轉(zhuǎn)換成實(shí)數(shù)塊對(duì)角形[VC,DC]=rsf2csf(VR,DR)把實(shí)數(shù)塊對(duì)角形轉(zhuǎn)換成復(fù)數(shù)對(duì)角形[V,J]=jordan(A)Jordan分解使A*V=V*J，J是塊對(duì)角陣C=condeig(A)向量c中包含矩陣A關(guān)于各特征值條件數(shù)[V,D,c]=condeig(A)相當(dāng)于[V，D]=eig(A)和c=condeig(A)兩條指令的組合【例1】簡(jiǎn)單實(shí)陣的特征值問(wèn)題。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中的復(fù)數(shù)特征值對(duì)角陣D轉(zhuǎn)換成實(shí)數(shù)塊對(duì)角陣，使VR*DR/VR=A。[VR,DR]=cdf2rdf(V,D)

VR=0.774600.0430-0.6310DR=0.83332.4438-2.44380.8333【例3】對(duì)虧損矩陣進(jìn)行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矩陣的譜分解和矩陣函數(shù)

(特征根各異的)矩陣A的特征值分解可重寫(xiě)為：據(jù)此定義矩陣函數(shù)。

常見(jiàn)矩陣函數(shù)的數(shù)學(xué)定義和指令表數(shù)學(xué)表達(dá)式指令形式功用A^p據(jù)譜分解，求ApP^A據(jù)譜分解，求pAexpm(A)一般常用expm1(A)用Pade近似求eAexpm2(A)用Taylor近似求eA，適用任何方差，精度稍差expm3(A)用特征值分解求eA，僅對(duì)非虧損陣使用sqrtm(A)據(jù)Schur分解，求。多解時(shí)，僅給出一個(gè)解logm(A)據(jù)譜分解，求lnAfunm(A,’FN’)FN是函數(shù)名(如sin等)【例1】數(shù)組乘方與矩陣乘方的比較。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】標(biāo)量的數(shù)組乘方和矩陣乘方的比較。（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的數(shù)組運(yùn)算和矩陣運(yùn)算比較。（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多項(xiàng)式和卷積4.3.1多項(xiàng)式一、多項(xiàng)式表達(dá)方式的約定MATLAB約定降冪多項(xiàng)式用以下系數(shù)行向量表示：即把多項(xiàng)式的各項(xiàng)系數(shù)依降冪次序排放在行向量的元素位置上。其中，多項(xiàng)式中缺某冪次項(xiàng)，則應(yīng)認(rèn)為該冪次項(xiàng)的系數(shù)為零。二、多項(xiàng)式運(yùn)算函數(shù)指令含義P=conv(p1,p2)p是多項(xiàng)式p1和p2的乘積多項(xiàng)式[q,r]=deconv(p1,p2)多項(xiàng)式p1被p2除的商多項(xiàng)式為q，而余多項(xiàng)式為rP=poly(AR)求方陣AR的特征多項(xiàng)式p；或求向量AR指定根所對(duì)應(yīng)的多項(xiàng)式P=polyfit(x,y,n)求x，y向量給定的數(shù)據(jù)的n階擬合多項(xiàng)式ppA=polyval(p,S)按數(shù)組運(yùn)算規(guī)則計(jì)算多項(xiàng)式的值。P為多項(xiàng)式，S為矩陣PM=polyvalm(p,S)按矩陣運(yùn)算規(guī)則計(jì)算多項(xiàng)式的值。P為多項(xiàng)式，S為矩陣[r,p,k]=residue(b,a)部分分式展開(kāi)：b，a分別式分子分母多項(xiàng)式系數(shù)向量；r，p，k分別是留數(shù)、極點(diǎn)、直項(xiàng)R=roots(p)求多項(xiàng)式P的根例求的“商”及“余”多項(xiàng)式。p1=conv([1,0,2],conv([1,4],[1,1]));%計(jì)算分子多項(xiàng)式p2=[1011]; %注意缺項(xiàng)補(bǔ)零 [q,r]=deconv(p1,p2);cq='商多項(xiàng)式為';cr='余多項(xiàng)式為';disp([cq,poly2str(q,'s')]),disp([cr,poly2str(r,'s')])商多項(xiàng)式為s+5余多項(xiàng)式為5s^2+4s+3例求3階方陣A的特征多項(xiàng)式。A=[111213;141516;171819];PA=poly(A)%A的特征多項(xiàng)式PPA=poly2str(PA,‘s’)%比較習(xí)慣的方式顯示多項(xiàng)

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

例由給定根向量求多項(xiàng)式系數(shù)向量。R=[-0.5,-0.3+0.4*i,-0.3-0.4*i];%根向量P=poly(R)%R的特征多項(xiàng)式 PR=real(P)%求P的實(shí)部 PPR=poly2str(PR,'x')P=1.00001.10000.55000.1250PR=1.00001.10000.55000.1250PPR=x^3+1.1x^2+0.55x+0.125注：含復(fù)數(shù)的根向量所生成的多項(xiàng)式（如P）的系數(shù)有可能帶有截?cái)嗾`差數(shù)量級(jí)的虛部，此時(shí)可采用實(shí)部的指令“real”把那很小的虛部濾掉。

例兩種多項(xiàng)式求值指令的差別。S=pascal(4) %生成一個(gè)4階方陣 P=poly(S);PP=poly2str(P,'s')PA=polyval(P,S) %獨(dú)立變量取數(shù)組S元素時(shí)的多項(xiàng)式PM=polyvalm(P,S)%獨(dú)立變量取矩陣S時(shí)的多項(xiàng)式

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應(yīng)該為0，這是著名的“Caylay-Hamilton”定理：任何一個(gè)矩陣滿(mǎn)足他自己的特征多項(xiàng)式方程.本例中的PM很小，這是截?cái)嗾`差造成的。例部分分式展開(kāi)。a=[1,3,4,2,7,2]; %分母多項(xiàng)式系數(shù)向量b=[3,2,5,4,6]; %分子多項(xiàng)式系數(shù)向量[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=[]三擬合和插值例對(duì)于給定數(shù)據(jù)對(duì)x0,y0，求擬合三階多項(xiàng)式，并圖示擬合情況。%給定數(shù)據(jù)對(duì)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àng)式系數(shù)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采用三次多項(xiàng)式所得的擬合曲線(xiàn)

例以上例所給數(shù)據(jù)，研究一維插值，并觀(guān)察插值與擬合的區(qū)別。x0=0:0.1:1;%給出數(shù)據(jù)對(duì)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;%采用三次多項(xiàng)式進(jìn)行插值yi=interp1(x0,y0,xi,'cubic');

plot(xi,yi,‘-b’,x0,y0,‘.r’,‘MarkerSize’,20),xlabel(‘x’)%繪圖通過(guò)三次多項(xiàng)式插值所得的曲線(xiàn)

擬合和插值的區(qū)別：曲線(xiàn)擬和研究如何尋找“平滑”曲線(xiàn)最好地表現(xiàn)帶噪聲的“測(cè)量數(shù)據(jù)”，但是不要求擬合曲線(xiàn)穿過(guò)這些“測(cè)量數(shù)據(jù)”點(diǎn)。插值是在認(rèn)定所給出“基準(zhǔn)數(shù)據(jù)”完全正確的情況下，研究如何“平滑”地估算出“基準(zhǔn)數(shù)據(jù)”之間其它點(diǎn)的函數(shù)值，因此插值所的曲線(xiàn)一定穿過(guò)“基準(zhǔn)數(shù)據(jù)”。4.3.2卷積一離散序列的數(shù)值卷積設(shè)有如下兩個(gè)任意序列，那么該卷積為它可準(zhǔn)確計(jì)算的非平凡區(qū)間是[Nc1,Nc2].在此，Nc1=N1+N3，Nc2=min{N1+N4，N2+N3}二MATLAB的“卷積”指令1對(duì)序列A、B進(jìn)行的卷積運(yùn)算C=conv（A，B）2解卷運(yùn)算，使B=conv（A，Q）+R成立[Q，R]=dconv（B，A）【例4.3-8】有序列和。（A）求這兩個(gè)完整序列的卷積，并圖示。（B）假設(shè)序列中最后4個(gè)非零值未知，而成為截尾序列，求卷積并圖示。a=ones(1,10);n1=3;n2=12; %完整A（n）序列的非平凡值和區(qū)間端點(diǎn)b=ones(1,8);n3=2;n4=9;%完整B（n）序列的非平凡值和區(qū)間端點(diǎn)c=conv(a,b);nc1=n1+n3;nc2=n2+n4;%計(jì)算卷積和確定卷積非平凡區(qū)間端點(diǎn) kc=nc1:nc2;%構(gòu)成非平凡區(qū)間的序號(hào)自變量 %截尾序列卷積 aa=a(1:6);nn1=3;nn2=8; %截尾A（n）序列的非平凡值和區(qū)間端點(diǎn)cc=conv(aa,b);ncc1=nn1+n3;nx=nn2+n4; %“非平凡”區(qū)間右端點(diǎn) ncc2=min(nn1+n4,nn2+n3);%截尾序列卷積被正確計(jì)算區(qū)間的右端點(diǎn)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

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

對(duì)于“完整”序列，指令conv所給出的結(jié)果在整個(gè)非平凡區(qū)間上都是準(zhǔn)確的，如本例中“藍(lán)空心桿圖”所示；對(duì)于“非完整”的截尾序列，指令conv給出的結(jié)果，只有部分非平凡區(qū)間上的結(jié)果是準(zhǔn)確的，即本例的“紅實(shí)心桿圖”，而“綠空心桿圖”雖也是由“截尾”序列算得的“非平凡”值，但不是真正的卷積。4.4數(shù)據(jù)分析函數(shù)MATLAB進(jìn)行數(shù)據(jù)分析時(shí)的約定：（1）若輸入宗量X是向量，那么不管是行向量還是列向量，運(yùn)算是對(duì)整個(gè)向量進(jìn)行的。（2）若輸入宗量X是二維數(shù)組（或稱(chēng)矩陣），那么指令運(yùn)算是按列進(jìn)行的。即默認(rèn)每個(gè)列是由一個(gè)變量的不同“觀(guān)察”所得的數(shù)據(jù)組成。4.4.1隨機(jī)數(shù)發(fā)生器和統(tǒng)計(jì)分析指令rand(n,m)

%產(chǎn)生（n×m）維的“[0，1]區(qū)間”均勻分布隨機(jī)數(shù)組randn(n,m)%產(chǎn)生(n×m)維的“均值為0標(biāo)準(zhǔn)差為1”的正態(tài)分布隨機(jī)數(shù)組min(X)

%對(duì)X矩陣各列分別求最小值max(X)

%對(duì)X矩陣各列分別求最大值median(X)

%對(duì)X矩陣各列分別求中位數(shù)mean(X)

%對(duì)X矩陣各列分別求均值std(X)

%對(duì)X矩陣各列分別求標(biāo)準(zhǔn)差var(X)

%對(duì)X矩陣各列分別求方差C=cov(X)

%對(duì)X矩陣各列間的協(xié)方差陣P=corrcoef(X)

%給出矩陣X各列間的相關(guān)數(shù)例基本統(tǒng)計(jì)示例。randn(‘state’,0) %為隨機(jī)仿真可重復(fù)而設(shè)置 A=randn(1000,4);AMAX=max(A),AMIN=min(A)%因分別在+3，-3附近AMED=median(A) %應(yīng)接近0 AMEAN=mean(A) %應(yīng)接近0 ASTD=std(A) %應(yīng)接近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) %相當(dāng)于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差分和累計(jì)指令del2(U,hx,hy)五點(diǎn)Laplacian，近似給出diff(X,m,n)沿第n維求m差分[DZx,DZy]=gradient(Z,dx,dy)對(duì)Z求x,y方向的梯度prod(X,n)沿第n維求乘積sum(X,n)沿第n維求和trapz(x,Y,n)采用梯形法沿第n維求函數(shù)Y關(guān)于資變量x的積分cumprod(X,n)沿第n維求累計(jì)乘積cumsum(X,n)沿第n為求累計(jì)和cumtrapz(x,Y,n)采用梯形法沿第n維求函數(shù)Y關(guān)于自變能量x的累計(jì)積分[XS,KK]=sort(X,n)沿第n維對(duì)X元素按模增大排列。KK給出XS元素的位置

例用一個(gè)簡(jiǎn)單矩陣表現(xiàn)diff和gradient指令計(jì)算方式。F=[1,2,3;4,5,6;7,8,9]Dx=diff(F) ％相鄰行差分 Dx_2=diff(F,1,2) ％相鄰列差分，第三輸入宗量2表示“列”差分 [FX,FY]=gradient(F) ％數(shù)據(jù)點(diǎn)步長(zhǎng)默認(rèn)為1 [FX_2,FY_2]=gradient(F,0.5) ％數(shù)據(jù)點(diǎn)步長(zhǎng)為0.5F=123456789Dx=333333Dx_2=111111FX=111111111FY=333333333FX_2=222222222FY_2=666666666

例函數(shù)的梯度和，用數(shù)值計(jì)算驗(yàn)證，并圖示。

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);％與理論計(jì)算比較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')

理論計(jì)算和數(shù)值計(jì)算的差別圖示

可見(jiàn)，數(shù)值梯度與理論梯度的區(qū)別僅發(fā)生在區(qū)間的端部；而Laplacian的數(shù)值結(jié)果與理論值幾乎沒(méi)有差別。例求積分,其中,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所給結(jié)果就是ss(end)st=cumtrapz(t,y); %梯形法求積 st10=trapz(t,y);stend=st(end);%trapz所給結(jié)果就是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中，凡以“函數(shù)”為輸入宗量的指令，都被統(tǒng)稱(chēng)為MATLAB泛函指令。最常見(jiàn)的有：

z=fzero(fun,x0)求一元函數(shù)零點(diǎn)的指令的最簡(jiǎn)格式

x=fsolve(fun,x0)解非線(xiàn)性方程組的最簡(jiǎn)單格式

x=fminbnd(fun,x1,x2)求函數(shù)在區(qū)間（x1,x2）中極小值的指令最簡(jiǎn)格式

x=fminsearch(fun,x0)單純形法求多元函數(shù)極值點(diǎn)指令的最簡(jiǎn)格式

x=fminunc(fun,x0)擬牛頓法求多元函數(shù)極值點(diǎn)指令的最簡(jiǎn)格式

a=lsqnonlin(fun,a0)解非線(xiàn)性最小二乘問(wèn)題指令的最簡(jiǎn)格式

q=quad(fun,a,b)采用遞推自適應(yīng)Simpson法計(jì)算積分

q=quadl(fun,a,b)采用遞推自適應(yīng)Lobatto法求數(shù)值積分

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

二重（閉型）數(shù)值積分指令[t,YY]=ode45(fun,tspan,Y0)采用4、5階Runge-Kutta方程解算ODE初值問(wèn)題4.5.1對(duì)于任意函數(shù)f(x)=0來(lái)說(shuō)，可能有零點(diǎn)，也可能沒(méi)有；可能有一個(gè)零點(diǎn)，也可能有多個(gè)甚至無(wú)數(shù)多個(gè)零點(diǎn)。因此，很難說(shuō)出一個(gè)通用解法。解題步驟如下：（1）利用MATLAB作圖指令獲取初步近似解具體做法：現(xiàn)確定一個(gè)零點(diǎn)可能存在的自變量區(qū)間；然后利用plot指令畫(huà)出f(x)在該區(qū)間中的圖形；用眼觀(guān)察f(x)與橫軸的交點(diǎn)坐標(biāo)，或者更細(xì)致些，用zoom對(duì)交點(diǎn)處進(jìn)行局部放大再讀數(shù)。借助ginput指令獲得更精確些的交點(diǎn)坐標(biāo)值。（2）利用fzero指令求精確界z=fzero(fun,x0)求一元函數(shù)零點(diǎn)指令的簡(jiǎn)單格式[z,f_z,exitflag]=fzero(fun,x0,options,P1,P2,…)求一元函數(shù)零點(diǎn)指令的最完整格式注1由于fzero是根據(jù)函數(shù)是否穿越橫軸來(lái)決定零點(diǎn)的，因此本程序無(wú)法確定函數(shù)曲線(xiàn)僅觸及橫軸和不穿越的零點(diǎn)。注2第二輸入宗量x0是表示零點(diǎn)初時(shí)猜測(cè)。注3輸入宗量options是優(yōu)化迭代所采用的參數(shù)選項(xiàng)。注4輸入宗量P1，P2是向函數(shù)fun傳遞的參數(shù)。注5指令中，輸出宗量z是所求零點(diǎn)的自變量值。注6第二個(gè)輸出宗量f_z是函數(shù)值；第三個(gè)宗量是表示程序終止計(jì)算的條件。若exitflag>0，則表明找到零點(diǎn)后退出。例通過(guò)求的零點(diǎn)，綜合敘述相關(guān)指令的用法。（1）使用字符串表示被處理函數(shù)P1=0.1;P2=0.5; ％按泛函指令要求，這里參數(shù)必須用P1,P2,…表示y_C=‘sin(x).^2.*exp(-P1*x)-P2*abs(x)’;％這里自變量必須用x表示（2）作圖法觀(guān)察函數(shù)零點(diǎn)分布x=-10:0.01:10;％對(duì)自變量采樣，采樣步長(zhǎng)不宜過(guò)大Y=eval(y_C);％在采樣點(diǎn)上計(jì)算函數(shù)值clf,plot(x,Y,‘r’);holdon,plot(x,zeros(size(x)),‘k’);％畫(huà)坐標(biāo)橫軸xlabel('t');ylabel('y(t)'),holdoff

函數(shù)零點(diǎn)分布觀(guān)察圖（3）利用zoom和ginput指令獲得零點(diǎn)的初始近似值z(mì)oomon %在MATLAB指令窗中運(yùn)行，獲局部放大圖[tt,yy]=ginput(5);%在MATLAB指令窗中運(yùn)行，用鼠標(biāo)獲5個(gè)零點(diǎn)猜測(cè)值z(mì)oomoff

tt%顯示所得零點(diǎn)初始猜測(cè)值tt=-2.0046-0.5300-0.01150.61061.6590（4）求近似tt(4)的精確零點(diǎn)[t4,y4,exitflag]=fzero(y_C,tt(4),[],P1,P2)

Zerofoundintheinterval:[0.59333,0.6106].t4=0.5993y4=0exitflag=1局部放大和利用鼠標(biāo)取值圖

圖中的十字是ginput運(yùn)行后產(chǎn)生的取值圖符4.5.2求函數(shù)極值點(diǎn)

許多科學(xué)研究和工程計(jì)算問(wèn)題都可以歸結(jié)為一個(gè)極值問(wèn)題，而且所有的求極大值問(wèn)題都可以化成求極小值的問(wèn)題，所以MATLAB中只有處理極小值的指令求函數(shù)極值的三條指令：[f,fval,exitflag,output]=fminbnd(fun,x1,x2,options,P1,P2,…)求一元函數(shù)在區(qū)間（x1,x2）中極小值的指令最完整格式[f,fval,exitflag,output]=fminsearch(fun,x0,options,P1,P2,…)單純形法求多元函數(shù)極值點(diǎn)指令的最完整格式[f,fval,exitflag,output,grad,hessian]=fminunc(fun,x0,options,P1,P2,…)擬牛頓法求多元函數(shù)極值點(diǎn)指令的最完整格式注：x1,x2為被研究區(qū)間的左、右邊界；x0為極值點(diǎn)位置的初始猜測(cè)例求的極小值點(diǎn)。它即是著名的Rosenbrock‘s“Banana”測(cè)試函數(shù)，它的理論極小值是。該測(cè)試函數(shù)有一片淺谷，許多算法難以越過(guò)此谷。（1）本例采用內(nèi)聯(lián)函數(shù)表示測(cè)試函數(shù)如下：ff=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2','x');（2）用單純形法求極小值點(diǎn)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數(shù)值積分?jǐn)?shù)值積分有閉型、開(kāi)型算法之分。兩者的區(qū)別在于：是否需要計(jì)算積分區(qū)間端點(diǎn)處的函數(shù)值。MATLAB中僅提供閉型數(shù)值積分指令：q=quad(fun,a,b,tol,trace,p1,p2,…)采用遞推自適應(yīng)Simpson法計(jì)算積分q=quadl(fun,a,b,tol,trace,p1,p2,…)采用遞推自適應(yīng)Lobatto法求數(shù)值積分ss=dblquad(fun,inmin,inmax,outmin,outmax,tol,method)二重（閉型）數(shù)值積分指令注：tol是一個(gè)標(biāo)量，用來(lái)控制絕對(duì)誤差，缺省時(shí)，積分的絕對(duì)精度為1×10－6

trace取非0值時(shí)，將隨積分的進(jìn)程逐點(diǎn)畫(huà)出被積函數(shù)p1,p2是向被積函數(shù)輸送的參數(shù)method是積分方法，如”quad”,”quadl”例求，其精確值為.（1）編寫(xiě)被積函數(shù)的M文件[FUN.M]functiony=fun(x)y=exp(-x.*x);（2）利用函數(shù)句柄運(yùn)行MATLAB指令quad和quadl進(jìn)行積分Hf=@fun;Isim=quad(Hf,0,1),IL=quadl(Hf,0,1)Isim=0.7468IL=0.7468

4.5.4解常微分方程以例題說(shuō)明：求微分方程在初始條件情況下的解，并圖示。

（1）把高階微分方程改寫(xiě)成一階微分方程組令，于是原二階方程可改寫(xiě)成如下一階方程組（2）根據(jù)上述一階微分方程組編寫(xiě)M函數(shù)文件[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）畫(huà)相平面圖

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

微分方程解相平面軌跡4.6信號(hào)處理MATLAB有專(zhuān)門(mén)的工具包（如SignalProcessingToolbox,DSPToolbox等）解決信號(hào)處理中的各種計(jì)算和仿真問(wèn)題。本節(jié)只介紹MATLAB基本組件中的三個(gè)信號(hào)處理重要指令fit,ifft,filter。4.6.1快速Fourier變換和逆變換

離散序列x(k)的離散Fourier變換和逆離散Fourier變換的“教科書(shū)”定義為：而MATLAB出于軟件的考慮：序列或向量元素下標(biāo)從1開(kāi)始記述，而不是從0開(kāi)始。因此，上述兩式在MATLAB中相應(yīng)的表達(dá)式為：MATLAB根據(jù)以上定義給出了如下指令：

計(jì)算N點(diǎn)Xt序列的N點(diǎn)離散Fourier變換Xf

Xf=fft(Xt,N,DIM)計(jì)算N點(diǎn)Xt序列的N點(diǎn)離散Fourier逆變換Xt

Xt=ifft(Xt,N,DIM)注：使用N的目的：（1）使被變換序列的長(zhǎng)度恰為2，以使指令fft較高速的運(yùn)作。（2）當(dāng)N大于被變換序列長(zhǎng)度時(shí)，原序列的尾部將用0填補(bǔ)，使其長(zhǎng)度為N后，再做變換。（3）缺省時(shí)，指令算法將以被變換序列的長(zhǎng)度為N。例利用fft和ifft指令重新計(jì)算兩序列的卷積。所給已知序列為

（1）產(chǎn)生給定序列a(n),b(n)cleara=ones(1,13);a([1,2,3])=0; ％a時(shí)間序列應(yīng)是0時(shí)刻算起b=ones(1,10);b([1,2])=0;％b時(shí)間序列應(yīng)是0時(shí)刻算起（2）直接卷積c=conv(a,b);

（3）通過(guò)變換求卷積M=32;AF=fft(a,M);BF=fft(b,M);CF=AF.*BF; ％使用點(diǎn)乘 cc=real(ifft(CF));

％過(guò)濾掉由于截?cái)嗾`差引起的虛部（4）圖式計(jì)算結(jié)果nn=0:(M-1); ％保證繪圖時(shí)，時(shí)間從0時(shí)刻開(kāi)始c(M)=0; ％使直接卷積所得序列通過(guò)尾部補(bǔ)0，與cc長(zhǎng)度相同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')

變換法和直接法求卷積結(jié)果比較

例fft在信號(hào)分析中的應(yīng)用。試用頻譜分析方法從受噪聲污染的信號(hào)中鑒別出有用信號(hào)。在此，。(1)構(gòu)遭受噪音污染的信號(hào)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)受噪聲污染的信號(hào)(2)畫(huà)幅頻曲線(xiàn)Y=fft(y);Ts=t(2)-t(1) ％時(shí)間信號(hào)的采樣周期 Ws=2*pi/Ts; ％時(shí)間信號(hào)的采樣頻率Wn=Ws/2 ％Nyquest頻率w=linspace(0,Wn,length(t)/2);％半采樣頻率中相應(yīng)的刻度Ya=abs(Y(1:length(t)/2));plot(w,Ya)Ts=0.0196Wn=160.5354受噪聲污染信號(hào)的幅頻譜(3)對(duì)感興趣頻段局部放大ii=find(w<=20);plot(w(ii),Ya(ii))grid,xlabel('FrequencyRad/s')受噪聲污染信號(hào)幅頻譜的局部放大4.6.2數(shù)字濾波命令：y=filter[B,A,x,[],dim]運(yùn)用數(shù)字濾波對(duì)信號(hào)x進(jìn)行濾波注：1)x可以是一維數(shù)組，當(dāng)x是二維數(shù)組時(shí)，默認(rèn)設(shè)置下，濾波沿列進(jìn)行。當(dāng)x是高維數(shù)組時(shí)，要用dim指定濾波沿哪維進(jìn)行。2)該指令的第三和四個(gè)輸入宗量可以缺省例設(shè)計(jì)一個(gè)低通濾波器，從受噪聲干擾的多頻率混合信號(hào)中獲取10Hz的信號(hào)。在此，而。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));％生成帶噪音的多頻率信號(hào)wn=ws/2; ％Nyquest頻率 [B,A]=butter(10,30/wn); ％截至頻率為30/wn的10階ButterWorth低通濾波器y=filter(B,A,x); ％進(jìn)行（初值為0的）濾波處理 plot(t,x,'b-',t,y,'r.','MarkerSize',10)legend('Input','Output',0)

注：butter指令在SignalProcessingToolbox工具包中。10階Butterworth濾波器的濾波效果4.7系統(tǒng)分析4.7.1線(xiàn)性時(shí)不變對(duì)象LTIS_ss=ss(A,B,C,D)利用狀態(tài)方程四對(duì)組生成LTI對(duì)象S_zpk=zpk(Z,P,K)利用零極點(diǎn)增益三對(duì)組生成LTI對(duì)象S_tf=tf(Num,Den)利用傳遞函數(shù)二對(duì)組生成LTI對(duì)象[A,B,C,D]=ssdata(S_Lti)從LTI對(duì)象獲取狀態(tài)四對(duì)組[Z,P,K]=zpkdata(S_Lti)從LTI對(duì)象獲取零極點(diǎn)增益三對(duì)組[Num,Den]=tfdata(S_Lti)從LTI對(duì)象獲取傳遞函數(shù)二對(duì)組注:(1)MATLAB中規(guī)定：LTI對(duì)象三個(gè)子類(lèi)中，“狀態(tài)方程子類(lèi)”的優(yōu)先級(jí)別最高，“零極點(diǎn)增益子類(lèi)”次之，“傳遞函數(shù)子類(lèi)”級(jí)別最低。(2)使用LTI對(duì)象的好處：①可以直接運(yùn)用“方塊圖”見(jiàn)例2；②可以把MIMO多輸入多輸出系統(tǒng)作為一個(gè)整體進(jìn)行處理。例1已知一個(gè)兩輸入兩輸出系統(tǒng)的狀態(tài)方程四對(duì)組，據(jù)此建立LTI對(duì)象“狀態(tài)方程子類(lèi)”模型。然后再?gòu)拇四Ｐ瞳@取傳遞函數(shù)二對(duì)組。（1）建立LTI對(duì)象“狀態(tài)方程子類(lèi)”模型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) ％建立狀態(tài)方程類(lèi)模型a=x1x2x3x10.50.50.7071x2-0.5-0.50.7071x3-6.364-0.7071-8b=u1u2x10