5.4節(jié)數(shù)字積分和微分方程數(shù)值解

1、5.4節(jié) 數(shù)值積分和微分方程數(shù)值解一數(shù)值定積分求面積【例5-4-1】 用數(shù)值積分法求由 ，y=0, x=0與x=10圍成的圖形面積，并討論步長和積分方法對精度的影響。解: 原理 用矩形法和梯形法分別求數(shù)值積分并作比較，步長的變化用循環(huán)語句實現(xiàn)。MATLAB中的定積分有專門的函數(shù)QUAD,QUADL等實現(xiàn)。為了弄清原理，我們先用直接編程的方法來計算，然后再介紹定積分函數(shù)及其調(diào)用方法。設x向量的長度取為n，即將積分區(qū)間分為n-1段，各段長度為 。算出各點的，則矩形法數(shù)值積分公式為：矩形和梯形定積分公式梯形法的公式為：比較兩個公式，它們之間的差別只是 。在MATLAB中，把向量中各元素疊加的命令是s

2、um。把向量中各元素按梯形法疊加的命令是trapz。梯形法的幾何意義是把被積分的函數(shù)的各計算點以直線相聯(lián)，形成許多窄長梯形條，然后疊加，我們把兩種算法都編入同一個程序進行比較。求面積的數(shù)值積分程序exn541for dx=2,1,0.5,0.1 % 設不同步長x=0:.1:10;y=-x.*x+115; % 取較密的函數(shù)樣本 plot(x,y),hold on % 畫出被積曲線并保持x1=0:dx:10;y1=-x1.*x1+115; % 求取樣點上的y1% 用矩形（歐拉）法求積分，注意末尾去掉一個點n=length(x1);s=sum(y1(1:n-1)*dx; q=trapz(y1)*dx

3、; % 用梯形法求積分stairs(x1,y1), % 畫出歐拉法的積分區(qū)域plot(x1,y1) % 畫出梯形法的積分區(qū)域dx,s,q,pause(1),hold off，end 程序exn541運行結果程序運行的結果如下：步長dx 矩形法解s 梯形法解q2 9108101 865815 .5841.25816.25 .1821.65816.65用解析法求出的精確解為2450/3=816.6666.。dx=2時矩形法和梯形法的積分面積見圖5-4-1.。在曲線的切線斜率為負的情況下，矩形法的積分結果一定偏大，梯形法是由各采樣點聯(lián)線包圍的面積，在曲線曲率為負（上凸）時，其積分結果一定偏小，因此精

4、確解在這兩者之間。由這結果也能看出，在步長相同時，梯形法的精度比矩形法高。矩形法數(shù)字積分的演示程序rsumsMATLAB中有一個矩形法數(shù)字積分的演示程序rsums，可以作一個對比。鍵入rsums(115-x.2,0,10)就得到右圖。圖中表示了被積函數(shù)的曲線和被步長分割的小區(qū)間，并按各區(qū)間中點的函數(shù)值構成了各個窄矩形面積。用鼠標拖動圖下方的滑尺可以改變步長的值，圖的上方顯示的是這些矩形面積疊加的結果。MATLAB內(nèi)的數(shù)值定積分函數(shù)在實際工作中，用MATLAB中的定積分求面積的函數(shù)quad和quadl可以得到比自編程序更高的精度，因為quad函數(shù)用的是辛普生法，即把被積函數(shù)用二次曲線逼近的算法，

5、而quadl函數(shù)采用了更高階的逼近方法。它們的調(diào)用格式如下：Q = QUADL(FUN,A,B,TOL)其中，F(xiàn)UN是表示被積函數(shù)的字符串， A是積分下限，B是積分上限。TOL是規(guī)定計算的容差，其默認值為1e-6 例如，鍵入S = quad(-x.*x+115,0,10)得到S = 8.166666666666666e+002二求兩條曲線所圍圖形的面積【例5-4-2】。設計算區(qū)間0,4上兩曲線所圍面積。解：原理：先畫出圖形， dx=input(dx= ) ;x=0:dx:4; f=exp(-(x-2).2.*cos(pi*x); g=4*cos(x-2); plot(x,f,x,g,:r)得到

6、右圖。從圖上看到，其中既有f(x)g(x)的區(qū)域，也有f(x)g(x)的區(qū)域，求兩條曲線所圍圖形的面積(1)若要求兩曲線所圍總面積（不管正負），則可加一條語句 s=trapz(abs(f-g)*dx, 在dx=0.001時，得到s = 6.47743996919702若要求兩曲線所圍的f(x)g(x)的正面積，則需要一定的技巧.方法一。先求出交點x1 ，再規(guī)定積分上下限。 x1=fzero(exp(-(x-2).2.*cos(pi*x)-4*cos(x-2),1)%把積分限設定為0 x1，求出積分結果再乘以2： x=0:dx:x1; f=exp(-(x-2).2.*cos(pi*x); g=4

7、*cos(x-2); s1=2*trapz(abs(f-g)*dx在設定dx=0.001時，得到s1 = 2.30330486000857求兩條曲線所圍圖形的面積(2)方法二。調(diào)用MATLAB中求面積函數(shù)quad。這里的關鍵是建立一個函數(shù)文件，把e1=f(x)-g(x)0的部分取出來。利用邏輯算式(e10)，它在e10處取值為1，在e10)與e1作元素群乘法，正的e1將全部保留，而負的e1就全部為零。因此編出子程序exn542f.m如下：function e = exn542f(x)e1=exp(-(x-2).2.*cos(pi*x)- 4*cos(x-2);e = (e1=0).*e1;將它

8、存入工作目錄下。于是求此積分的主程序語句為： s2=quad(exn542f,0,4)得到的結果為：s2= 2.30330244952618三求曲線長度【例5-4-3】設曲線方程及定義域為：用計算機做如下工作：(a) 按給定區(qū)間畫出曲線，再按n=2,4,8份分割并畫出割線。(b) 求這些線段長度之和，作為弧長的近似值。(c) 用積分來估算弧長，并與用割線計算的結果比較。解：原理：先按分區(qū)間算割線長度的方法編程，然后令分段數(shù)不斷增加求得其精密的結果，最后可以與解析結果進行比較。因此編程應該具有普遍性，能由用戶設定段數(shù)，并在任何分段數(shù)下算出結果。求曲線長度的程序exn543n=input(分段數(shù)目

9、n= ), % 輸入分段數(shù)目x=linspace(-1,1,n+1); % 設定x向量y=sqrt(1-x.2); % 求y向量plot(x,y), hold on % 繪圖并保持Dx=diff(x); % 求各段割線的x方向長度% x向量長度為n+1，Dx是相鄰x元素的差，其元素數(shù)為nDy=diff(y); % Dy是相鄰兩個y元素的差Ln=sqrt(Dx.2+Dy.2); % 求各割線長度L=sum(Ln) % 求n段割線的總長度程序exn543的運行結果程序運行后得到圖5-32，在不同的n下，其數(shù)值結果為：n=2, L =2.82842712474619n= 4,L = 3.035276

10、18041008n= 8,L = 3.10449606777484n= 1000L = 3.14156635621648我們已經(jīng)可以大致猜測出它將趨向于，精確的極限值可用下列符號數(shù)學語句導出。 x=syms x,y= sqrt(1-x2),L=int(sqrt(1+diff(y)2),-1,1)這個程序其實有相當?shù)耐ㄓ眯?，不同的被積函數(shù)，只要改變其中的一條函數(shù)賦值語句，并相應地改變自變量的賦值范圍就行了。 四求旋轉(zhuǎn)體體積【例5-4-4】求曲線與x軸所圍成的圖形分別繞x軸和y軸旋轉(zhuǎn)所成的旋轉(zhuǎn)體的體積。解：原理：由于旋轉(zhuǎn)對稱性，在圓周方向的計算只要乘以圓周長度，不需要積分運算。因此旋轉(zhuǎn)體的體積計算

11、實際上就退化單變量求積分。程序如下： %先畫出平面圖形 dx=input(dx= ) ;x=0:dx:pi; g=x.*sin(x).2; plot(x,g) 求筒形旋轉(zhuǎn)體體積(a)。繞y軸體積用薄圓柱筒形體作為微分體積單元，其半徑為x，厚度為dx，高度為g(x)，其立體圖見圖5-34左，此筒形單元的截面積為g*dx，薄環(huán)的微體積為：dv=2*pi*x*dx*g，旋轉(zhuǎn)體的體積為微分體積單元沿x方向的和，鍵入： v=trapz(2*pi*x.*g*dx)得：v = 27.53489480036561求盤形旋轉(zhuǎn)體體積(b)。繞x軸體積它繞x軸旋轉(zhuǎn)形成一個薄圓盤，其厚度為dx，而半徑為g(x) 。所

12、以此薄盤體的微體積為： dv1=pi*g.2.*dx，旋轉(zhuǎn)體的體積為微分體積單元沿x方向的和： v1=trapz(pi*g.2*dx)得：v1 = 9.86294784774497用符號數(shù)學工具箱的程序精確的理論結果可用符號數(shù)學工具箱函數(shù)求得如下： syms x, g=x*sin(x)2; v1t= int(pi*g2,0,pi),double(v1t)v1t = 1/8*pi4-15/64*pi2ans = 9.86294784774499大多數(shù)的定積分并不會有理論的解析結果，所以這樣的驗證一般是不必要的。 五多重積分【例5-4-5】 計算二重積分積分區(qū)域為由x=1,y=x及y=0所圍成的閉

13、合區(qū)域.解: 原理 先畫出積分區(qū)域，在任意x處取出沿y向的一個單元條，其寬度為dx，而高度為y=x，所以y是一個數(shù)組。其上的被積函數(shù)f也是一個數(shù)組，沿y向的積分可用trapz(f)完成，得到s1(k)，它是隨x而變的。用for循環(huán)求出所有的s1(k)。 再沿x方向用trapz函數(shù)積分。MATLAB的數(shù)組運算可以代替一個for循環(huán)，所以二重積分只用了一組for語句。二重積分的MATLAB程序exn545clear,format compactfill(0,1,1,0,0,0,1,0,y),hold% 畫出積分區(qū)域fill(0.55,0.6,0.6,0.55,0.55,0,0,0.6,0.55,0

14、,r) %畫出單元條dx=input(步長dx= );dy=dx;x=0:dx:1;lx=length(x);for k=1:lx x1=(k-1)*dx; y1=0:dy:x1; f=x1.2+y1.2; s1(k)=trapz(f)*dy;ends=trapz(s1)*dx用MATLAB函數(shù)求二重積分(1)運行的數(shù)值結果在步長dx=0.01時為:s=0.3334另一種方法是利用MATLAB中現(xiàn)成的二重積分函數(shù)dblquad，其調(diào)用格式為：Q = DBLQUAD(FUN,XMIN,XMAX,YMIN,YMAX,TOL)其中FUN是x,y的函數(shù)，接下來的四個變元是四個積分限，其中前兩個對應于x

15、，后兩個對應于y，TOL為允許誤差（默認值為1.e-16）。這四個積分限只許用常數(shù)代入，可見dblquad函數(shù)只能用于積分區(qū)域為矩形的情況。解決的方法之一是仍用矩形區(qū)積分，但把不屬于積分區(qū)域內(nèi)的函數(shù)置成零，其方法與上題有些類似。 用MATLAB函數(shù)求二重積分(2)在圖示的積分區(qū)域中，對角線左上方的白色區(qū)域滿足y-x0，邏輯式(y-x Q=dblquad(x.2+y.2).*(y-x0),0,1,0,1)得到Q = 0.33332245532028三重積分的計算【例5-4-6】計算三重積分積分區(qū)域為由x=1,y=x,z=xy及三個坐標面所圍成的區(qū)域.解: 方法 先畫出積分區(qū)域圖5-36,這個區(qū)域

16、在xy平面上的投影與圖5-35相仿，只是增加了z方向的高度，從而構成了一個三維的實體。 先畫出積分區(qū)域。這個區(qū)域在xy平面上的投影上例相仿，只是增加了z方向的高度，從而構成了一個三維的實體。程序exn546a用來畫這個立體空間。x=1,y=x都是沿z向的柱面，本題還用了plot3命令以畫出與z軸平行的輔助平面。頂面則是由z=xy構成的二次曲面，用mesh函數(shù)容易畫出。難點在于此區(qū)域有效的xy底面是一個三角形,不易用自變量網(wǎng)格表示。為了解決這個問題，采取了上題中對因變量乘以邏輯式的處理方法。將z乘以(y-x0)就可使y=x直線左上方所有z值都變成0。畫出三維圖后可以靠鼠標來拖動三維圖形旋轉(zhuǎn),以得

17、到一個最好的視覺效果。 繪制積分區(qū)域的程序exn546a%本程序給出由x=1,y=x,z=xy三個曲面圍成的積分區(qū)域.x,y=meshgrid(0:.05:1); % 確定矩形定義域網(wǎng)格z1=x.*y.*(y-x syms x, y=x2*atan(x), Z=int(y)得到y(tǒng) = x2*atan(x)Z = 1/3*x3*atan(x)-1/6*x2+1/6*log(x2+1)符號數(shù)學求不定積分例 548(b)例5-4-8(b)。解下列積分，畫出解的曲線。解：程序為 syms x Y=int(cos(x)2+sin(x)系統(tǒng)運行的結果為：Y =1/2*cos(x)*sin(x)+1/2*x

18、-cos(x)+C代入邊界條件，得： C=1-(1/2*cos(pi)*sin(pi)+1/2*pi-cos(pi)結果是C =-pi/2符號數(shù)學求不定積分例 548(c)例5-4-8(c)。求下列積分。解：程序為： syms x,Y=int(exp(-x2) % 不定積分：Y1=int(exp(-x2),0,1) % 定積分：運行結果為：Y = 1/2*pi(1/2)*erf(x)Y1 = 1/2*erf(1)*pi(1/2)【應用篇中與本節(jié)相關的例題】【例6-5-1】和例6-5-2磁場計算由元電流生成的磁場積分，得到全部線圈產(chǎn)生的磁場?！纠?-1-3】導彈跟蹤目標的軌跡，根據(jù)x,y兩個方向的速度積分得出軌跡?！纠?-1-5】有空氣阻力下的拋射體軌跡，是微分方程的數(shù)字積分求解問題。【例7-2-2】懸臂梁的橈度計算，這是純粹的兩次積分問題。【例7-2-3】簡支梁的橈度計算，也是純粹的兩次積分問題。【例7-3-1】，【例7-3-2】，一自由度自由振動和強迫振動，都是常微分方程求解的典型問題。這兩道題用的是解析解，只是用數(shù)學軟件來求根和繪制波形?！緫闷信c本節(jié)相關的例題】【例7-3-3】兩自由度振動，是四階的矩陣微分方程的問題。本題對