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

1、12022-5-1第三章第三章 MATLAB數(shù)值計(jì)算數(shù)值計(jì)算22022-5-1何時(shí)需采用數(shù)值計(jì)算？n解析解不存在；n解析解存在，但非常復(fù)雜；如：由Cramer法則，可將n元一次方程組化簡(jiǎn)為n個(gè)n-1元一次方程； n-1元一次方程組可化簡(jiǎn)為n-1個(gè)n-2元一次方程；結(jié)論： n元一次方程組的解析解是可以求出的。但求解需進(jìn)行方程(n-1)(n+1)!+n次基本運(yùn)算。一個(gè)20元一次方程組需進(jìn)行方程9.70731020次基本運(yùn)算。常見(jiàn)的問(wèn)題：力學(xué)中的有限元法；差分方程及微分方程的數(shù)值解法；FFT方法等。32022-5-1n3.1 線性方程組線性方程組n3.2 矩陣的非線性運(yùn)算矩陣的非線性運(yùn)算n3.3 數(shù)

2、值數(shù)值積分積分n3.4 數(shù)值微分?jǐn)?shù)值微分n3.5 常微分方程常微分方程n3.6 MEX技術(shù)簡(jiǎn)介技術(shù)簡(jiǎn)介42022-5-1 Ax = b 有x = A-1b，但實(shí)際上并不顯式求A-13.1 線性方程組例子： 7x = 21 x = 21/7=3如果求逆 x = 7-1 21 = .142857 21 = 2.99997這就需要一次除和一次乘，且精度更低52022-5-1運(yùn)算符 AX = BX = AB 左除 XA = BX = B/A 右除62022-5-13x3的例子6475156230710321xxx65 5 46237710 32132121xxxxxxxx5 . 21 . 6755 .

3、 2061 . 000710321xxx1 . 65 . 2761 . 0055 . 200710321xxx2 . 65 . 272 . 60055 . 200710321xxx72022-5-1算法矩陣表示2 . 60055 . 200710U104. 03 . 0015 . 0001L010100001PPALU 單位下三角陣上三角陣置換陣82022-5-1LU分解例子5156230710A1000100011P104. 000100012M0101000012P105 . 0013 . 00011M103 . 0015 . 00011L104. 000100012L92022-5-1線

4、性方程組求解Ax = b LU = PA Ly = Pb Ux = y x = Ab例： clearA=1 2 3 4 5 6 7 8 1;L,U,P = lu(A)Err=P*A - L*U102022-5-1矩陣的QR分解：nQR分解也稱(chēng)為長(zhǎng)方陣的正交分解，把長(zhǎng)方陣分解為正交矩陣和上三角矩陣的積。n格式：Q,R=qr(A)例：X=1 2 3 4 5 6 7 8 2 1 6 3;Q,R =qr(X)Q1,R1 =qr(X)112022-5-1矩陣的Cholesky分解：nCholesky分解把對(duì)稱(chēng)正定矩陣分解為上三角矩陣和其轉(zhuǎn)置的積，即X=CC，C為上三角陣。n格式：C=chol(A)例：X

5、=6 1 1 2 1 8 3 1 1 3 7 0 2 1 0 4;C=chol(X)C*C %驗(yàn)證X=CC122022-5-1分析在三位精度的算法中，x2是由兩個(gè)與舍入誤差相同階的數(shù)相除得到的，因此x2可以是任意的數(shù)。然后完全任意的x2代入第一個(gè)方程得到x1 ，所以第一方程的殘量是小的。我們也可以期望第二個(gè)方程的殘量也是小的，因?yàn)榫仃囀墙咏娈惖?。部分主元的高斯消去法可以確保產(chǎn)生小的殘量部分主元的高斯消去法可以確保產(chǎn)生小的殘量“小的”舍入誤差的階與三個(gè)數(shù)量有關(guān)：原來(lái)矩陣悉數(shù)矩陣單元的大小，消去過(guò)程中系數(shù)矩陣單元的大小和計(jì)算解的元素的大小。如果其中一個(gè)是“大的”，殘量在絕對(duì)意義上不必是小的。即使

6、殘量很小，但不能保證解的誤差很小。它們之間的關(guān)系可以部分由矩陣的條件數(shù)條件數(shù)來(lái)刻畫(huà)132022-5-1范數(shù)如果Ax = b, 如何測(cè)量x關(guān)于A和b的敏感性？pxxpnipip1 ,/11iiniiniixxxxxxmax , ,2/121211 x =(1:4)/5 n1 = norm(x,1) n2 = norm(x) ninf = norm(x,inf)142022-5-1擾動(dòng)分析bAx bbxxAxMb xmbbbAxx)(條件數(shù)是相對(duì)誤差放大因子，也是刻畫(huà)矩陣奇異性的測(cè)度。152022-5-1矩陣范數(shù)向量有三種常用范數(shù)，相對(duì)應(yīng)的矩陣范數(shù)的三種形式為： （A的行范數(shù)） （A的列范數(shù)） （

7、 是ATA 的最大特征值） （ A的2范數(shù)）11maxniji niAa 111maxnijjniAa 12A1162022-5-1A=1 2 3 4 5 6 7 8 9;k2=norm(A) k2_1=norm(A,2)k1=norm(A,1) k_inf=norm(A,inf)A=1 2 3 4 5 6 7 8 1;k2=norm (A) k2_1=norm (A,2)k1=norm (A,1) k_inf=norm (A,inf)例：172022-5-1矩陣條件數(shù)xAxmxAxMmin ,max條件數(shù)：1minmax)(AAxAxxAxa182022-5-1矩陣條件數(shù) cond(A) 或

8、cond(A,2):使用svd(A) cond(A,1)計(jì)算k1(A)：使用inv(A)，工作量比 cond(A,2)小 cond(A,inf): 同cond(A,1) condest(A) 估計(jì)k1(A)：用lu(A)用Higham和Tisseur的算法，適合大的稀疏矩陣192022-5-1A=1 2 3 4 5 6 7 8 9;k2=cond(A) k2_1=cond(A,2)k1=cond(A,1) k_inf=cond(A,inf)A=1 2 3 4 5 6 7 8 1;k2=cond(A) k2_1=cond(A,2)k1=cond(A,1) k_inf=cond(A,inf)例：2

9、02022-5-1非齊次線性方程組的求解n對(duì)于非齊次線性方程組Ax=b而言，要根據(jù)系數(shù)矩陣A的秩和增廣矩陣B=A b的秩和未知數(shù)個(gè)數(shù)n的關(guān)系，才能判斷方程組Ax=b的解的情況；1如果系數(shù)矩陣的秩=增廣矩陣的秩=n，則方程組有唯一解；2如果系數(shù)矩陣的秩=增廣矩陣的秩n，則方程組有無(wú)窮多解；3如果系數(shù)矩陣的秩增廣矩陣的秩，則方程組無(wú)解。212022-5-1n對(duì)于非齊次線性方程組Ax=b而言，首先應(yīng)判斷方程組解的情況，其次，若有解，先求出方程組的特解，再次，求出對(duì)應(yīng)齊次方程組的通解，最后，寫(xiě)出非齊次方程組的通解，即特解+對(duì)應(yīng)齊次方程組的通解。n求Ax=b對(duì)應(yīng)的齊次方程組Ax=0的通解，可以利用函數(shù)n

10、ull或?qū)ο禂?shù)矩陣A施行行變換；n求Ax=b的特解，方法較多，根據(jù)方程組中方程的個(gè)數(shù)m和未知數(shù)的個(gè)數(shù)n，可以把方程組Ax=b分為： 超定方程組（mn）； 恰定方程組（m=n）； 欠定方程組（mn）?？梢愿鶕?jù)系數(shù)矩陣A的性質(zhì)選用適當(dāng)?shù)挠?jì)算方法，如可利用左除左除 “” 、系數(shù)矩陣A的逆或偽逆、或利用lu分解、或cholesky分解等。222022-5-1clearA=4 2 -1;3 -1 2;11 3 0; b=2 10 8; B=A b; %增廣矩陣 n=3; RA=rank(A) RB=rank(B)if (RA=RB & RA=n) X=Ab else if (RA=RB &

11、; RA *Inner matrix dimensions must agree.372022-5-12牛頓柯特斯法牛頓柯特斯法基于牛頓柯特斯法，基于牛頓柯特斯法，MATLAB給出了給出了quadl函數(shù)來(lái)求定積函數(shù)來(lái)求定積分。分。（在（在6.5以前的版本中為：以前的版本中為： quad8 ） 該函數(shù)的調(diào)用格式為：該函數(shù)的調(diào)用格式為：I,n=quadl(fname,a,b,tol,trace)其中參數(shù)的含義和其中參數(shù)的含義和quad函數(shù)相似，只是函數(shù)相似，只是tol的缺省值取的缺省值取10-6。該函數(shù)可以更精確地求出定積分的值，且一般情況下函數(shù)該函數(shù)可以更精確地求出定積分的值，且一般情況下函數(shù)調(diào)

12、用的步數(shù)明顯小于調(diào)用的步數(shù)明顯小于quad函數(shù)，從而保證能以更高的效函數(shù)，從而保證能以更高的效率求出所需的定積分值。率求出所需的定積分值。382022-5-1例例3-2 求定積分。求定積分。(1) 被積函數(shù)文件被積函數(shù)文件fx.m。 function f=fx(x) f=x.*sin(x)./(1+cos(x).*cos(x);(2) 調(diào)用函數(shù)調(diào)用函數(shù)quadl求定積分。求定積分。 I=quadl(fx,0,pi) I = 2.4674注：在6.5以上的版本中將提示：Warning: QUAD8 is obsolete. We use QUADL instead.392022-5-1例例3-3

13、 分別用分別用quad函數(shù)和函數(shù)和quadl函數(shù)求定積分的近似函數(shù)求定積分的近似值，并在相同的積分精度下，比較函數(shù)的調(diào)用次數(shù)。值，并在相同的積分精度下，比較函數(shù)的調(diào)用次數(shù)。調(diào)用函數(shù)調(diào)用函數(shù)quad求定積分：求定積分：clearfx=inline(exp(-x) ); I,n=quad(fx,1,2.5,1e-10)I = 0.28579444254766n = 65構(gòu)造構(gòu)造inline函數(shù)，函數(shù)，可省略不必要的函可省略不必要的函數(shù)文件數(shù)文件402022-5-1 調(diào)用函數(shù)調(diào)用函數(shù)quadl求定積分：求定積分：clear fx=inline(exp(-x);I,n=quadl(fx,1,2.5,1

14、e-10)I=vpa(I,10)輸出：輸出：I = 0.2858n = 18I =.2857944425412022-5-13被積函數(shù)由一個(gè)表格定義被積函數(shù)由一個(gè)表格定義在在MATLAB中，對(duì)由表格形式定義的函數(shù)關(guān)系的求定積分中，對(duì)由表格形式定義的函數(shù)關(guān)系的求定積分問(wèn)題用問(wèn)題用trapz(X,Y)函數(shù)。其中向量函數(shù)。其中向量X,Y定義函數(shù)關(guān)系定義函數(shù)關(guān)系Y=f(X)。例：用例：用trapz函數(shù)（梯形法函數(shù)（梯形法trapezoidal ）計(jì)算定積分。）計(jì)算定積分。命令如下：命令如下：clearX=1:0.01:2.5;Y=exp(-X); %生成函數(shù)關(guān)系數(shù)據(jù)向量生成函數(shù)關(guān)系數(shù)據(jù)向量trapz(

15、X,Y)vpa(ans,11)ans = 0.2858ans = .28579682416422022-5-13.3.3 二重定積分的數(shù)值求解二重定積分的數(shù)值求解使用使用MATLAB提供的提供的dblquad函數(shù)就可以直函數(shù)就可以直接求出上述二重定積分的數(shù)值解。該函數(shù)接求出上述二重定積分的數(shù)值解。該函數(shù)的調(diào)用格式為：的調(diào)用格式為：I=dblquad(f,a,b,c,d,tol,trace)該函數(shù)求該函數(shù)求f(x,y)在在a,bc,d區(qū)域上的二重定區(qū)域上的二重定積分。參數(shù)積分。參數(shù)tol，trace的用法與函數(shù)的用法與函數(shù)quad完完全相同。全相同。432022-5-1例例3-5 計(jì)算二重定積分

16、計(jì)算二重定積分(1) 建立一個(gè)函數(shù)文件建立一個(gè)函數(shù)文件fxy.m：function f=fxy(x,y)global ki;ki=ki+1; %ki用于統(tǒng)計(jì)被積函數(shù)的調(diào)用次數(shù)用于統(tǒng)計(jì)被積函數(shù)的調(diào)用次數(shù)f=exp(-x.2/2).*sin(x.2+y);(2) 調(diào)用調(diào)用dblquad函數(shù)求解。函數(shù)求解。global ki; ki=0;I=dblquad(fxy,-2,2,-1,1) %單引號(hào)可有可無(wú)單引號(hào)可有可無(wú)kiI = 1.57449318974494ki = 1038442022-5-1clearfxy=inline(exp(-x.2/2).*sin(x.2+y);I=dblquad(fx

17、y,-2,2,-1,1)運(yùn)用運(yùn)用inline函數(shù)，更簡(jiǎn)潔函數(shù)，更簡(jiǎn)潔452022-5-13.4 數(shù)值微分?jǐn)?shù)值微分3.4.1 數(shù)值微分的實(shí)現(xiàn)數(shù)值微分的實(shí)現(xiàn)在在MATLAB中，沒(méi)有直接提供求數(shù)值導(dǎo)數(shù)的函數(shù)，中，沒(méi)有直接提供求數(shù)值導(dǎo)數(shù)的函數(shù)，只有計(jì)算向前差分的函數(shù)只有計(jì)算向前差分的函數(shù)diff，其調(diào)用格式為：，其調(diào)用格式為：DX=diff(X)：計(jì)算向量計(jì)算向量X的向前差分，的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,n-1。DX=diff(X,n)：計(jì)算計(jì)算X的的n階向前差分。例如，階向前差分。例如， diff(X,2)=diff(diff(X)。DX=diff(A,n,dim)：

18、計(jì)算矩陣計(jì)算矩陣A的的n階差分，階差分，dim=1時(shí)時(shí)(缺缺省狀態(tài)省狀態(tài))，按列計(jì)算差分；，按列計(jì)算差分；dim=2，按行計(jì)算差分。，按行計(jì)算差分。462022-5-1例例3-6 用不同的方法求函數(shù)用不同的方法求函數(shù)f(x)在在-3，3的數(shù)值導(dǎo)數(shù)，并在的數(shù)值導(dǎo)數(shù)，并在同一個(gè)坐標(biāo)系中做出同一個(gè)坐標(biāo)系中做出f(x)的圖像。的圖像。f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);g=inline(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5);132625326212(5)523411

19、522126(5)fxxxxxdfxxgdxxxxx472022-5-1例例3-6 程序如下：程序如下：clearf=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);g=inline(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5);x=-3:0.01:3;p=polyfit(x,f(x),5); %用用5次多項(xiàng)式次多項(xiàng)式p擬合擬合f(x)dp=polyder(p); %對(duì)擬合多項(xiàng)式對(duì)擬合多項(xiàng)式p求導(dǎo)數(shù)求導(dǎo)數(shù)dpdpx=polyval(dp,x); %求求dp在假設(shè)點(diǎn)的函數(shù)值在假設(shè)點(diǎn)的

20、函數(shù)值dx=diff(f(x,3.01)/0.01; %直接對(duì)直接對(duì)f(x)求數(shù)值導(dǎo)數(shù)求數(shù)值導(dǎo)數(shù)gx=g(x); %求函數(shù)求函數(shù)f的導(dǎo)函數(shù)的導(dǎo)函數(shù)g在假設(shè)點(diǎn)的導(dǎo)數(shù)在假設(shè)點(diǎn)的導(dǎo)數(shù)plot(x,dpx, k.,x,dx,.,x,gx, r-); %作圖作圖482022-5-13.5 常微分方程初值問(wèn)題的數(shù)值解法常微分方程初值問(wèn)題的數(shù)值解法3.5.1 龍格庫(kù)塔法簡(jiǎn)介龍格庫(kù)塔法簡(jiǎn)介Runge-Kutta方法：1905（德國(guó)）httsssshyyhsyhtfsshyhtfsshyhtfsytfsnnnnnnnnnnnn1432113423121)22(6),()2,2()2,2(),(若f不依賴(lài)y,則為

21、simpson公式步長(zhǎng) hRe(max)492022-5-13.5.2 龍格庫(kù)塔法的實(shí)現(xiàn)龍格庫(kù)塔法的實(shí)現(xiàn) 基于龍格庫(kù)塔法，基于龍格庫(kù)塔法，MATLAB提供了求提供了求非剛性非剛性常微分方程數(shù)值解的函常微分方程數(shù)值解的函數(shù)，一般調(diào)用格式為：數(shù)，一般調(diào)用格式為： t,y=ode23(fname,tspan,y0) t,y=ode45(fname,tspan,y0)其中其中fname是定義是定義f(t,y)的函數(shù)文件名，該函數(shù)文件必須返回一個(gè)列向量。的函數(shù)文件名，該函數(shù)文件必須返回一個(gè)列向量。tspan形式為形式為t0,tf,表示求解區(qū)間。表示求解區(qū)間。y0是初始狀態(tài)列向量。是初始狀態(tài)列向量。t和和

22、y分別給分別給出時(shí)間向量和相應(yīng)的狀態(tài)向量。出時(shí)間向量和相應(yīng)的狀態(tài)向量。剛性剛性常微分方程數(shù)值解的函數(shù)，一般調(diào)用格式為：常微分方程數(shù)值解的函數(shù)，一般調(diào)用格式為：ode15s及及ode23s502022-5-1例例3-10 試求初值問(wèn)題的數(shù)值解，并與精確解相比較。試求初值問(wèn)題的數(shù)值解，并與精確解相比較。 (1) 建立函數(shù)文件建立函數(shù)文件funt.m。function yp=funt(t,y)yp=(y2-t-2)/4/(t+1);(2) 求解微分方程。求解微分方程。cleart0=0; tf=10; y0=2;t,y=ode45(funt,t0,tf,y0); %求數(shù)值解求數(shù)值解y1=sqrt(t

23、+1)+1; %求精確解求精確解diary e:L.txt %生成日記文件生成日記文件t %輸出行向量輸出行向量yy1 diary off % y為數(shù)值解，為數(shù)值解，y1為精確值，顯然兩者近似。為精確值，顯然兩者近似。224(1)ytyt 512022-5-1nload-加載數(shù)據(jù)文件a.dat或b.txt中的數(shù)據(jù) 例：振動(dòng)測(cè)試的功率譜分析，powerp.mntype-在workspace中顯示a.dat或b.txt中的數(shù)據(jù) 例：type data.dat522022-5-1例例3-11 求解著名的求解著名的Van der Pol方程。方程。2(1)0yyyy變形：令變形：令12,xyxy122

24、2121(1)xxxxxx 532022-5-1function y=vdp_eq(t,x) mu=1;y=x(2); -mu*(x(1).2-1).*x(2)-x(1)clearx0=-0.2; -0.7; t_final=20;t,y=ode45(vdp_eq,0,t_final,x0);plot(t,y)figure; plot(y (:,1),y(:,2), -)1222121(1)xxxxxx 注意自變量必須寫(xiě)成注意自變量必須寫(xiě)成x(1)和和 x(2)542022-5-1例例3-12 Duffing方程，求時(shí)間相應(yīng)，試?yán)L制系方程，求時(shí)間相應(yīng)，試?yán)L制系統(tǒng)相平面圖。統(tǒng)相平面圖。30 xx

25、xx1232121xxxxxx 552022-5-1function y=duff_eq(t,x) mu=3; alpha=.3; k=-1;y=x(2); -k*x(1)-alpha*x(2)-mu*x(1).3;clear; closex0=-1;1; t_final=50;t,y=ode45(duff_eq,0,t_final,x0);plot(t,y);gridfigure; plot(y (:,1),y(:,2),:)1232121xxxxxx 562022-5-1n例例3-13 Lorenz系統(tǒng)，求時(shí)間相應(yīng)，試?yán)L制系統(tǒng)系統(tǒng)，求時(shí)間相應(yīng)，試?yán)L制系統(tǒng)相平面圖。相平面圖。()xxyyxz

26、rxyzxybz 1083b572022-5-1Lorenz系統(tǒng)方程：系統(tǒng)方程：function y=lorenzeq(t,x) sigma=10; b=8/3; r=28;y=-sigma*(x(1)-x(2); -x(1)*x(3)+r*x(1)-x(2); x(1)*x(2)-b*x(3);clear; close %混沌現(xiàn)象t_final=100; x0=0.01;0.01;0;t,x=ode45(lorenzeq,0,t_final,x0);plot(t,x),figure; plot3(x(:,1),x(:,2),x(:,3); grid()xxyyxzrxyzxybz 582022-5-13.6 M