1、第六章第六章 MATLAB 數(shù)值計(jì)算數(shù)值計(jì)算6.1 線性方程組的解線性方程組的解6.1.1 LU 分解、行列式、逆和恰定方程的解分解、行列式、逆和恰定方程的解【例 6.1-1】“求逆”法和“左除”法解恰定方程的性能對(duì)比（1）randn(state,0);A=gallery(randsvd,100,2e13,2);x=ones(100,1);b=A*x;cond(A) ans = 1.9990e+013 （2）ticxi=inv(A)*b;ti=toceri=norm(x-xi)rei=norm(A*xi-b)/norm(b) ti = 0.4400eri = 0.0469rei = 0.004

2、7 （3）tic;xd=Ab;td=toc,erd=norm(x-xd),red=norm(A*xd-b)/norm(b) td = 0.0600erd = 0.0078red = 2.6829e-015 6.1.2 奇異值分解和矩陣結(jié)構(gòu)奇異值分解和矩陣結(jié)構(gòu)6.1.3 線性二乘問題的解線性二乘問題的解【例 6.1-2】對(duì)于超定方程，進(jìn)行三種解法比較。其中取 MATLAB 庫(kù)中的特殊Axy A函數(shù)生成。A=gallery(5);A(:,1)=;y=1.7 7.5 6.3 0.83 -0.082;x=inv(A*A)*A*y,xx=pinv(A)*y,xxx=Ay Warning: Matrix

3、is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.087751e-018.x = 3.4811 5.1595 0.9534 -0.0466xx = 3.4759 5.1948 0.7121 -0.1101Warning: Rank deficient, rank = 3 tol = 1.0829e-010.xxx = 3.4605 5.2987 0 -0.2974 nx=norm(x),nxx=norm(xx),nxxx=norm(xxx) nx = 6.2968e+000nxx = 6.291

4、8e+000nxxx = 6.3356e+000 e=norm(y-A*x),ee=norm(y-A*xx),eee=norm(y-A*xxx) e = 1.1020e+000ee = 4.7424e-002eee = 4.7424e-002 6.2 特征值分解和矩陣函數(shù)特征值分解和矩陣函數(shù)6.2.2 特征值分解問題特征值分解問題【例 6.2-1】簡(jiǎn)單實(shí)陣的特征值問題。A=1,-3;2,2/3;V,D=eig(A) V = -0.7728 + 0.0527i -0.7728 - 0.0527i 0 + 0.6325i 0 - 0.6325iD = 0.8333 + 2.4438i 0 0 0.

5、8333 - 2.4438i 【例 6.2-2】把例 6.2-1 中的復(fù)數(shù)特征值對(duì)角陣 D 轉(zhuǎn)換成實(shí)數(shù)塊對(duì)角陣，使 VR*DR/VR=A。VR,DR=cdf2rdf(V,D) VR = -0.7728 0.0527 0 0.6325DR = 0.8333 2.4438 -2.4438 0.8333 【例 6.2-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_equ A = -9 11 -21 63 -252 70 -69 141 -421 1684

6、-575 575 -1149 3451 -13801 3891 -3891 7782 -23345 93365 1024 -1024 2048 -6144 24572DJ = 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0D = Columns 1 through 4 -0.0328 + 0.0243i 0 0 0 0 -0.0328 - 0.0243i 0 0 0 0 0.0130 + 0.0379i 0 0 0 0 0.0130 - 0.0379i 0 0 0 0 Column 5 0 0 0 0 0.0396 c_eig = 1.0e+

7、010 * 2.1016 2.1016 2.0251 2.0251 1.9796c_equ = 5.2129e+017 6.2.3 矩陣的譜分解和矩陣函數(shù)矩陣的譜分解和矩陣函數(shù)【例 6.2-4】數(shù)組乘方與矩陣乘方的比較。clear,A=1 2 3;4 5 6;7 8 9;A_Ap=A.0.3A_Mp=A0.3 A_Ap = 1.0000 1.2311 1.3904 1.5157 1.6207 1.7118 1.7928 1.8661 1.9332A_Mp = 0.6962 + 0.6032i 0.4358 + 0.1636i 0.1755 - 0.2759i 0.6325 + 0.0666i

8、0.7309 + 0.0181i 0.8292 - 0.0305i 0.5688 - 0.4700i 1.0259 - 0.1275i 1.4830 + 0.2150i 【例 6.2-5】標(biāo)量的數(shù)組乘方和矩陣乘方的比較。（A 取自例 6.2-4）pA_A=(0.3).ApA_M=(0.3)A pA_A = 0.3000 0.0900 0.0270 0.0081 0.0024 0.0007 0.0002 0.0001 0.0000pA_M = 2.9342 0.4175 -1.0993 -0.0278 0.7495 -0.4731 -1.9898 -0.9184 1.1531 【例 6.2-6】

9、sin 的數(shù)組運(yùn)算和矩陣運(yùn)算比較。（A 取自例 6.2-4）A_sinA=sin(A)A_sinM=funm(A,sin) A_sinA = 0.8415 0.9093 0.1411 -0.7568 -0.9589 -0.2794 0.6570 0.9894 0.4121A_sinM = -0.6928 -0.2306 0.2316 -0.1724 -0.1434 -0.1143 0.3479 -0.0561 -0.4602 6.3 多項(xiàng)式和卷積多項(xiàng)式和卷積6.3.2 多項(xiàng)式多項(xiàng)式6.3.2.1多項(xiàng)式表達(dá)方式的約定多項(xiàng)式表達(dá)方式的約定6.3.2.2多項(xiàng)式運(yùn)算多項(xiàng)式運(yùn)算函數(shù)函數(shù)【例 6.3-1】

10、求的“商”及“余”多項(xiàng)式。1) 1)(4)(2(32sssssp1=conv(1,0,2,conv(1,4,1,1);p2=1 0 1 1;q,r=deconv(p1,p2);cq=商多項(xiàng)式為商多項(xiàng)式為 ; cr=余多項(xiàng)式為余多項(xiàng)式為 ;disp(cq,poly2str(q,s),disp(cr,poly2str(r,s) 商多項(xiàng)式為 s + 5余多項(xiàng)式為 5 s2 + 4 s + 3 【例 6.3-2】求 3 階方陣 A 的特征多項(xiàng)式。A=11 12 13;14 15 16;17 18 19;PA=poly(A) PPA=poly2str(PA,s) PA = 1.0000 -45.0000

11、 -18.0000 -0.0000PPA = s3 - 45 s2 - 18 s - 2.8387e-015 【例 6.3.1.2-3】由給定根向量求多項(xiàng)式系數(shù)向量。R=-0.5,-0.3+0.4*i,-0.3-0.4*i;P=poly(R)PR=real(P) PPR=poly2str(PR,x) P = 1.0000 1.1000 0.5500 0.1250PR = 1.0000 1.1000 0.5500 0.1250PPR = x3 + 1.1 x2 + 0.55 x + 0.125 【例 6.3-4】?jī)煞N多項(xiàng)式求值指令的差別。S=pascal(4)P=poly(S);PP=poly2

12、str(P,s)PA=polyval(P,S)PM=polyvalm(P,S) S = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20PP = s4 - 29 s3 + 72 s2 - 29 s + 1PA = 1.0e+004 * 0.0016 0.0016 0.0016 0.0016 0.0016 0.0015 -0.0140 -0.0563 0.0016 -0.0140 -0.2549 -1.2089 0.0016 -0.0563 -1.2089 -4.3779PM = 1.0e-011 * -0.0077 0.0053 -0.0096 0.0430 -0.0068

13、 0.0481 -0.0110 0.1222 0.0075 0.1400 -0.0095 0.2608 0.0430 0.2920 -0.0007 0.4737 【例 6.3-5】部分分式展開。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.

14、2991 k = 6.3.2.3擬合和插值擬合和插值【例 6.3-6】 對(duì)于給定數(shù)據(jù)對(duì) x0 , y0 ，求擬合三階多項(xiàng)式，并圖示擬合情況。（見圖 6.3-1）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.1174 40.0070 -0.904302 4 8 1 2468

15、x圖6.3-1 采用三次多項(xiàng)式所得的擬合曲線【例 6.3-7】以上例所給數(shù)據(jù)，研究一維插值，并觀察插值與擬合的區(qū)別。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) x圖6.3-2 通過(guò)三次多項(xiàng)式插值所得的曲線6.3.3 卷積卷積6.3.3.1離散序列的數(shù)值卷積離散序列的數(shù)值卷積6.3.3.2MATLAB 的的“卷積卷積”指令指

16、令【例 6.3-8】有序列 和 。elsennA12, 4 , 301)(elsennB9 , 3 , 201)(（A）求這兩個(gè)完整序列的卷積，并圖示。（B）假設(shè)序列中最后 4 個(gè)非零值未知，而成A為截尾序列，求卷積并圖示。（見圖 6.3-3）a=ones(1,10);n1=3;n2=12;b=ones(1,8);n3=2;n4=9;c=conv(a,b);nc1=n1+n3;nc2=n2+n4;kc=nc1:nc2;aa=a(1:6);nn1=3;nn2=8;cc=conv(aa,b);ncc1=nn1+n3;nx=nn2+n4;ncc2=min(nn1+n4,nn2+n3);kx=(ncc

17、2+1):nx;kcc=ncc1:ncc2;N=length(kcc);stem(kcc,cc(1:N),r,filled)axis(nc1-2,nc2+2,0,10),grid,hold onstem(kc,c,b),stem(kx,cc(N+1:end),g),hold off 6 圖6.3-3 “完整”序列卷積和“截尾”序列卷積6.4 數(shù)據(jù)分析函數(shù)數(shù)據(jù)分析函數(shù)6.4.2 隨機(jī)數(shù)發(fā)生器和隨機(jī)數(shù)發(fā)生器和 統(tǒng)計(jì)分析指令統(tǒng)計(jì)分析指令【例 6.4-1】基本統(tǒng)計(jì)示例。randn(state,0)A=randn(1000,4);AMAX=max(A),AMIN=min(A)AMED=median(A)

18、AMEAN=mean(A)ASTD=std(A) AMAX = 2.7316 3.2025 3.4128 3.0868AMIN = -2.6442 -3.0737 -3.5027 -3.0461AMED = -0.0131 0.0596 0.0122 0.0459AMEAN = -0.0431 0.0455 0.0177 0.0263ASTD = 0.9435 1.0313 1.0248 0.9913 【例 6.4-2】cov 和 corrcoef 的使用示例。rand(state,1)X=rand(10,3);Y=rand(10,3);mx=mean(X);Xmx=X-ones(size(X

19、)*diag(mx);CCX=Xmx*Xmx/(size(Xmx,1)-1)CX=cov(X),CY=cov(Y)Cxy=cov(X,Y)PX=corrcoef(X)Pxy=corrcoef(X,Y) CCX = 0.0871 0.0242 -0.0073 0.0242 0.0846 0.0056 -0.0073 0.0056 0.0607CX = 0.0871 0.0242 -0.0073 0.0242 0.0846 0.0056 -0.0073 0.0056 0.0607CY = 0.0721 0.0013 0.0165 0.0013 0.0714 -0.0535 0.0165 -0.05

20、35 0.0720Cxy = 0.0761 -0.0012 -0.0012 0.0708PX = 1.0000 0.2819 -0.1005 0.2819 1.0000 0.0785 -0.1005 0.0785 1.0000Pxy = 1.0000 -0.0168 -0.0168 1.0000 6.4.3 差分和累計(jì)指令差分和累計(jì)指令【例 6.4-3】用一個(gè)簡(jiǎn)單矩陣表現(xiàn) diff 和 gradient 指令計(jì)算方式。F=1,2,3;4,5,6;7,8,9Dx=diff(F)Dx_2=diff(F,1,2)FX,FY=gradient(F)FX_2,FY_2=gradient(F,0.5) F

21、 = 1 2 3 4 5 6 7 8 9Dx = 3 3 3 3 3 3Dx_2 = 1 1 1 1 1 1FX = 1 1 1 1 1 1 1 1 1FY = 3 3 3 3 3 3 3 3 3FX_2 = 2 2 2 2 2 2 2 2 2FY_2 = 6 6 6 6 6 6 6 6 6 【例 6.4-4】函數(shù)的梯度和，用數(shù)值22),(yxyxz22),(yxyxz4),(2yxz計(jì)算驗(yàn)證，并圖示。（見圖 6.4-1）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);DZ

22、2=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) 01 01 0 2 4 01 01 0 2 4 01 0 2 4 x 圖6.4-1理論計(jì)算和數(shù)值計(jì)算的差別圖示【例 6.4-5】求積分，其中，。xdttyxs0)()(ttetysin8.0)(100 xdt=

23、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);st=cumtrapz(t,y);st10=trapz(t,y);stend=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) sum cumsum

24、 trapz cumtrapz 2.7082 2.7082 2.6576 2.65760246 0 5 1 2 5 圖6.4-2矩形法和梯形法求積比較6.5 MATLAB 泛函指令泛函指令6.5.2 求函數(shù)零點(diǎn)求函數(shù)零點(diǎn)【例 6.5-1】通過(guò)求的零點(diǎn)，綜合敘述相關(guān)指令的用法。tbettfat)(sin)(2P1=0.1;P2=0.5;y_C=sin(x).2.*exp(-P1*x)-P2*abs(x); x=-10:0.01:10;Y=eval(y_C);clf,plot(x,Y,r);hold on,plot(x,zeros(size(x),k);xlabel(t);ylabel(y(t),

25、hold off 圖6.5-1 函數(shù)零點(diǎn)分布觀察圖zoom ontt,yy=ginput(5);zoom off圖6.5-2 局部放大和利用鼠標(biāo)取值圖tt tt = -2.0046 -0.5300 -0.0115 0.6106 1.6590 t4,y4,exitflag=fzero(y_C,tt(4),P1,P2) Zero found in the interval: 0.59333, 0.6106.t4 = 0.5993y4 = 0exitflag = 1 6.5.3 求函數(shù)極值點(diǎn)求函數(shù)極值點(diǎn)【例 6.5-2】求的極小值點(diǎn)。它即是著名的 Rosenbrocks 222)1 ()(100),

26、(xxyyxfBanana 測(cè)試函數(shù)，它的理論極小值是。該測(cè)試函數(shù)有一片淺谷，許多算法難1, 1yx以越過(guò)此谷。ff=inline(100*(x(2)-x(1)2)2+(1-x(1)2,x); x0=-1.2,1;sx,sfval,sexit,soutput=fminsearch(ff,x0) Optimization terminated successfully: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 and F(X) satisfies the conver

27、gence criteria using OPTIONS.TolFun of 1.000000e-004 sx = 1.0000 1.0000sfval = 8.1777e-010sexit = 1soutput = iterations: 85 funcCount: 159 algorithm: Nelder-Mead simplex direct search 6.5.4 數(shù)值積分?jǐn)?shù)值積分【例 6.5-3】求，其精確值為 。dxeIx10274684204. 0（1）fun.mfunction y=fun(x)y=exp(-x.*x);（2）Hf=fun;Isim=quad(Hf,0,1)

28、,IL=quadl(Hf,0,1) Isim = 0.7468IL = 0.7468 6.5.5 解常微分方程解常微分方程【例 6.5-4】求微分方程在初始條件情況下0)1 (222xdtdxxdtxd0)0(, 1)0(dtdxx的解，并圖示。（見圖 6.5-3 和 6.5-4）（1）（2）DyDt.mfunction ydot=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) 5 5 5 圖6.5

29、-3 微分方程解（4）plot(yy(:,1),yy(:,2) 0 5 1 5 2 5 0123圖6.5-4 相平面軌跡6.6 信號(hào)處理信號(hào)處理6.6.2 快速快速 Fourier 變換和逆變換變換和逆變換【例 6.6-1】利用 fft 和 ifft 指令重新計(jì)算兩序列的卷積。所給已知序列為 ， 。elsenna12, 4 , 301)(elsennb9 , 3 , 201)(cleara=ones(1,13);a(1,2,3)=0;b=ones(1,10);b(1,2)=0; c=conv(a,b); M=32;AF=fft(a,M);BF=fft(b,M);CF=AF.*BF; cc=re

30、al(ifft(CF); nn=0:(M-1);c(M)=0;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) 8 5 0 圖6.6-1 變換法和直接法求卷積結(jié)果比較【例6.6-2】fft在信號(hào)分析中的應(yīng)用。試用頻譜分析方法從受噪聲污染的信號(hào)中鑒別出有用信號(hào)。在此，。)(9cos65sin3)(tNttty)5 , 0()(NtNclear,randn(s

31、tate,0)t=linspace(0,10,512);y=3*sin(5*t)-6*cos(9*t)+5*randn(size(t);plot(t,y) 8 0 0 5 圖6.6-2 受噪聲污染的信號(hào)Y=fft(y);Ts=t(2)-t(1)Ws=2*pi/Ts;Wn=Ws/2w=linspace(0,Wn,length(t)/2);Ya=abs(Y(1:length(t)/2);plot(w,Ya) Ts = 0.0196Wn = 160.5354 , 0 ? 0 0 00 0 0 0 圖6.6-3 受噪聲污染信號(hào)的幅頻譜ii=find(w=20);plot(w(ii),Ya(ii)gri

32、d,xlabel(Frequency Rad/s) 05 00 0 0 0 圖6.6-4 受噪聲污染信號(hào)幅頻譜的局部放大6.6.3 數(shù)字濾波數(shù)字濾波【例 6.6-3】設(shè)計(jì)一個(gè)低通濾波器，從受噪聲干擾的多頻率混合信號(hào)中獲取 10Hz 的信)(tx號(hào)，見圖 6.6-5。在此，)()1002cos()102sin()(tntttx而。)2 . 0 , 0()(Ntnclear,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;B,A=butter(10,30/wn

33、);y=filter(B,A,x);plot(t,x,b-,t,y,r.,MarkerSize,10)legend(Input,Output,0) 1 5 2 5 圖6.6-5 10 階 Butterworth 濾波器的濾波效果6.7 系統(tǒng)分析系統(tǒng)分析6.7.2 線性時(shí)不變對(duì)象線性時(shí)不變對(duì)象 LTI【例 6.7-1】已知一個(gè)兩輸入兩輸出系統(tǒng)的狀態(tài)方程四對(duì)組，據(jù)此建立 LTI 對(duì)象“狀態(tài)方程子類”模型。然后再?gòu)拇四Ｐ瞳@取傳遞函數(shù)二對(duì)組。A=0.5 0.5 0.7071;-0.5 -0.5 0.7071;-6.364 -0.7071 -8;B=0 0 4;1 0 1;C=0.7071 -0.7071 0;0 0.5 -0.5;D=0;S_ss=ss(A,B,C,D) a = x1 x2 x3 x1 0.5 0.5 0.7071 x2 -0.5 -0.5 0.7071 x3 -6.364 -0.7071 -8 b = u1 u2 x1 0 1 x2 0 0 x3 4 1 c = x1 x2 x3 y1 0.7071 -0.7071 0 y2 0 0.5 -0.5 d = u1 u2 y1 0 0 y2 0 0 Continuous-time model. S_tf=tf(minr