數(shù)值積分和微分方程數(shù)值解.ppt

5.3節(jié) 數(shù)值積分和微分方程 數(shù)值解,一數(shù)值定積分求面積,【例5-3-1】 用數(shù)值積分法求由 ，y=0, x=0與x=10圍成的圖形面積，并討論步長和積分方法對精度的影響。 解: 原理 用矩形法和梯形法分別求數(shù)值積分并作比較，步長的變化用循環(huán)語句實(shí)現(xiàn)。MATLAB中的定積分有專門的函數(shù)QUAD,QUADL等實(shí)現(xiàn)。為了弄清原理，我們先用直接編程的方法來計(jì)算，然后再介紹定積分函數(shù)及其調(diào)用方法。設(shè)x向量的長度取為n，即將積分區(qū)間分為n-1段， 各段長度為 。算出各點(diǎn)的 ，則矩形法數(shù)值積分公式為：,矩形和梯形定積分公式,梯形法的公式為： 比較兩個(gè)公式，它們之間的差別只是 。 在MATLAB中，把向量中各元素疊加的命令是sum。把向量中各元素按梯形法疊加的命令是trapz。梯形法的幾何意義是把被積分的函數(shù)的各計(jì)算點(diǎn)以直線相聯(lián)，形成許多窄長梯形條，然后疊加，我們把兩種算法都編入同一個(gè)程序進(jìn)行比較。,求面積的數(shù)值積分程序exn531,for dx=2,1,0.5,0.1 % 設(shè)不同步長 x=0:.1:10;y=-x.*x+115; % 取較密的函數(shù)樣本 plot(x,y),hold on % 畫出被積曲線并保持 x1=0:dx:10;y1=-x1.*x1+115; % 求取樣點(diǎn)上的y1 % 用矩形（歐拉）法求積分，注意末尾去掉一個(gè)點(diǎn) n=length(x1);s=sum(y1(1:n-1)*dx; q=trapz(y1)*dx; % 用梯形法求積分 stairs(x1,y1), % 畫出歐拉法的積分區(qū)域 plot(x1,y1) % 畫出梯形法的積分區(qū)域 dx,s,q,pause(1), hold off，end,程序exn531運(yùn)行結(jié)果,程序運(yùn)行的結(jié)果如下： 步長dx 矩形法解s 梯形法解q 2 910 810 1 865 815 .5 841.25 816.25 .1 821.65 816.65 用解析法求出的精確解為2450/3=816.6666.。 dx=2時(shí)矩形法和梯形法的積分面積見圖5-4-1.。在曲線的切線斜率為負(fù)的情況下，矩形法的積分結(jié)果一定偏大，梯形法是由各采樣點(diǎn)聯(lián)線包圍的面積，在曲線曲率為負(fù)（上凸）時(shí)，其積分結(jié)果一定偏小，因此精確解在這兩者之間。由這結(jié)果也能看出，在步長相同時(shí)，梯形法的精度比矩形法高。,矩形法數(shù)字積分的演示程序rsums,MATLAB中有一個(gè)矩形法數(shù)字積分的演示程序rsums，可以作一個(gè)對比。鍵入 rsums(115-x.2,0,10) 就得到右圖。圖中表示了被積函數(shù)的曲線和被步長分割的小區(qū)間，并按各區(qū)間中點(diǎn)的函數(shù)值構(gòu)成了各個(gè)窄矩形面積。用鼠標(biāo)拖動(dòng)圖下方的滑尺可以改變步長的值，圖的上方顯示的是這些矩形面積疊加的結(jié)果。,MATLAB內(nèi)的數(shù)值定積分函數(shù),在實(shí)際工作中，用MATLAB中的定積分求面積的函數(shù)quad和quadl可以得到比自編程序更高的精度，因?yàn)閝uad函數(shù)用的是辛普生法，即把被積函數(shù)用二次曲線逼近的算法，而quadl函數(shù)采用了更高階的逼近方法。它們的調(diào)用格式如下： Q = QUADL(FUN,A,B,TOL) 其中，F(xiàn)UN是表示被積函數(shù)的字符串， A是積分下限，B是積分上限。TOL是規(guī)定計(jì)算的容差，其默認(rèn)值為1e-6 例如，鍵入 S = quad(-x.*x+115,0,10) 得到 S = 8.166666666666666e+002,二求兩條曲線所圍圖形的面積,【例5-3-2】。設(shè) 計(jì)算區(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),得到右圖。從圖上看到，其中既有f(x)g(x)的區(qū)域，也有f(x)g(x)的區(qū)域，,求兩條曲線所圍圖形的面積(1),若要求兩曲線所圍總面積（不管正負(fù)），則可加一條語句 s=trapz(abs(f-g)*dx, 在dx=0.001時(shí)，得到s = 6.47743996919702 若要求兩曲線所圍的f(x)g(x)的正面積，則需要一定的技巧. 方法一。先求出交點(diǎn)x1 ，再規(guī)定積分上下限。 x1=fzero(exp(-(x-2).2.*cos(pi*x)-4*cos(x-2),1) %把積分限設(shè)定為0x1，求出積分結(jié)果再乘以2： x=0:dx:x1; f=exp(-(x-2).2.*cos(pi*x); g=4*cos(x-2); s1=2*trapz(abs(f-g)*dx 在設(shè)定dx=0.001時(shí)，得到s1 = 2.30330486000857,求兩條曲線所圍圖形的面積(2),方法二。調(diào)用MATLAB中求面積函數(shù)quad。這里的關(guān)鍵是建立一個(gè)函數(shù)文件，把e1=f(x)-g(x)0的部分取出來。 利用邏輯算式(e10)，它在e10處取值為1，在e10)與e1作元素群乘法，正的e1將全部保留，而負(fù)的e1就全部為零。因此編出子程序exn542f.m如下： function e = exn542f(x) e1=exp(-(x-2).2.*cos(pi*x)- 4*cos(x-2); e = (e1=0).*e1; 將它存入工作目錄下。于是求此積分的主程序語句為： s2=quad(exn542f,0,4) 得到的結(jié)果為： s2= 2.30330244952618,三求曲線長度,【例5-3-3】設(shè)曲線方程及定義域?yàn)椋?用計(jì)算機(jī)做如下工作： (a) 按給定區(qū)間畫出曲線，再按n=2,4,8份分割并畫出割線。 (b) 求這些線段長度之和，作為弧長的近似值。 (c) 用積分來估算弧長，并與用割線計(jì)算的結(jié)果比較。 解：原理：先按分區(qū)間算割線長度的方法編程，然后令分段數(shù)不斷增加求得其精密的結(jié)果，最后可以與解析結(jié)果進(jìn)行比較。因此編程應(yīng)該具有普遍性，能由用戶設(shè)定段數(shù)，并在任何分段數(shù)下算出結(jié)果。,求曲線長度的程序exn533,n=input(分段數(shù)目n= ), % 輸入分段數(shù)目 x=linspace(-1,1,n+1); % 設(shè)定x向量 y=sqrt(1-x.2); % 求y向量 plot(x,y), hold on % 繪圖并保持 Dx=diff(x); % 求各段割線的x方向長度 % x向量長度為n+1，Dx是相鄰x元素的差，其元素?cái)?shù)為n Dy=diff(y); % Dy是相鄰兩個(gè)y元素的差 Ln=sqrt(Dx.2+Dy.2); % 求各割線長度 L=sum(Ln) % 求n段割線的總長度,程序exn533的運(yùn)行結(jié)果,程序運(yùn)行后得到圖5-32，在不同的n下，其數(shù)值結(jié)果為： n=2, L = 2.82842712474619 n= 4, L = 3.03527618041008 n= 8, L = 3.10449606777484 n= 1000 L = 3.14156635621648 我們已經(jīng)可以大致猜測出它將趨向于，精確的極限值可用下列符號數(shù)學(xué)語句導(dǎo)出。 syms x,y= sqrt(1-x2),L=int(sqrt(1+diff(y)2),-1,1) 這個(gè)程序其實(shí)有相當(dāng)?shù)耐ㄓ眯?，不同的被積函數(shù)，只要改變其中的一條函數(shù)賦值語句，并相應(yīng)地改變自變量的賦值范圍就行了。,四求旋轉(zhuǎn)體體積,【例5-3-4】求曲線與x軸所圍成的圖形分別繞x軸和y軸旋轉(zhuǎn)所成的旋轉(zhuǎn)體的體積。,解：原理：由于旋轉(zhuǎn)對稱性，在圓周方向的計(jì)算只要乘以圓周長度，不需要積分運(yùn)算。因此旋轉(zhuǎn)體的體積計(jì)算實(shí)際上就退化單變量求積分。程序如下： %先畫出平面圖形 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)形成一個(gè)薄圓盤，其厚度為dx，而半徑為g(x) 。所以此薄盤體的微體積為： dv1=pi*g.2.*dx, %旋轉(zhuǎn)體的體積為微分體積單元沿x方向的和： v1=trapz(pi*g.2*dx) 得：v1 = 9.86294784774497,用符號數(shù)學(xué)工具箱的程序,精確的理論結(jié)果可用符號數(shù)學(xué)工具箱函數(shù)求得如下： syms x, g=x*sin(x)2; v1t= int(pi*g2,0,pi),double(v1t) v1t = 1/8*pi4-15/64*pi2 ans = 9.86294784774499 大多數(shù)的定積分并不會(huì)有理論的解析結(jié)果，所以這樣的驗(yàn)證一般是不必要的。,五多重積分,【例5-3-5】 計(jì)算二重積分 積分區(qū)域?yàn)橛蓌=1,y=x及y=0所圍成的閉合區(qū)域.,解: 原理 先畫出積分區(qū)域，在任意x處取出沿y向的一個(gè)單元條，其寬度為dx，而高度為y=x，所以y是一個(gè)數(shù)組。其上的被積函數(shù)f也是一個(gè)數(shù)組，沿y向的積分可用trapz(f)完成，得到s1(k)，它是隨x而變的。用for循環(huán)求出所有的s1(k)。 再沿x方向用trapz函數(shù)積分。MATLAB的數(shù)組運(yùn)算可以代替一個(gè)for循環(huán)，所以二重積分只用了一組for語句。,二重積分的MATLAB程序exn535,clear,format compact fill(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,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; end s=trapz(s1)*dx,用MATLAB函數(shù)求二重積分(1),運(yùn)行的數(shù)值結(jié)果在步長dx=0.01時(shí)為: s =0.33337500000000 另一種方法是利用MATLAB中現(xiàn)成的二重積分函數(shù) dblquad，其調(diào)用格式為： Q = DBLQUAD(FUN,XMIN,XMAX,YMIN,YMAX,TOL) 其中FUN是x,y的函數(shù)，接下來的四個(gè)變元是四個(gè)積分限，其中前兩個(gè)對應(yīng)于x，后兩個(gè)對應(yīng)于y，TOL為允許誤差（默認(rèn)值為1.e-16）。這四個(gè)積分限只許用常數(shù)代入，可見dblquad函數(shù)只能用于積分區(qū)域?yàn)榫匦蔚那闆r。 解決的方法之一是仍用矩形區(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,三重積分的計(jì)算,【例5-3-6】計(jì)算三重積分 積分區(qū)域?yàn)橛蓌=1,y=x,z=xy及三個(gè)坐標(biāo)面所圍成的區(qū)域. 解: 方法 先畫出積分區(qū)域圖5-36,這個(gè)區(qū)域在xy平面上的投影與圖5-35相仿，只是增加了z方向的高度，從而構(gòu)成了一個(gè)三維的實(shí)體。,先畫出積分區(qū)域。這個(gè)區(qū)域在xy平面上的投影上例相仿，只是增加了z方向的高度，從而構(gòu)成了一個(gè)三維的實(shí)體。程序exn546a用來畫這個(gè)立體空間。,x=1,y=x都是沿z向的柱面，本題還用了plot3命令以畫出與z軸平行的輔助平面。頂面則是由z=xy構(gòu)成的二次曲面，用mesh函數(shù)容易畫出。難點(diǎn)在于此區(qū)域有效的xy底面是一個(gè)三角形,不易用自變量網(wǎng)格表示。為了解決這個(gè)問題，采取了上題中對因變量乘以邏輯式的處理方法。將z乘以(y-x0)就可使y=x直線左上方所有z值都變成0。畫出三維圖后可以靠鼠標(biāo)來拖動(dòng)三維圖形旋轉(zhuǎn),以得到一個(gè)最好的視覺效果。,繪制積分區(qū)域的程序exn536a,%本程序給出由x=1,y=x,z=xy三個(gè)曲面圍成的積分區(qū)域. x,y=meshgrid(0:.05:1); % 確定矩形定義域網(wǎng)格 z1=x.*y.*(y-x0); % 求z1=xy并構(gòu)成三角形定義域 mesh(x,y,z1);hold on; % 畫出積分區(qū)頂部 % 以下畫出積分區(qū)域的幾個(gè)側(cè)柱面 x1=0:0.02:1;y1=x1;sx1=length(x1); zd=zeros(1,sx1);x1.*y1; plot3(x1;x1,y1;y1,zd,*) line(ones(2,sx1),y1;y1,zeros(1,sx1);y1) plot3(ones(2,sx1),y1;y1,zeros(1,sx1);y1,o),編寫三重積分程序的思路,(1) 在任意點(diǎn)（x,y）處取出沿z向的一個(gè)單元條,其底面積為dx*dy,而高度為z=x*y,這一個(gè)細(xì)柱體上從z=0到z=x*y間的所有各點(diǎn)度屬于積分的區(qū)域，把它表為z向的一個(gè)數(shù)組。因?yàn)榇颂?x,y)固定，其上的被積函數(shù)f=f(z)是隨z而變的一個(gè)數(shù)組，沿z向的積分可用trapz(f(z)完成，得到s1。 (2) s1是隨x,y而變的,先固定x,用for循環(huán)求出沿y向所有的s1(j),用trapz函數(shù)求其和s2=trapz(s1); (3) s2(k)又是隨x而變的,再沿x方向用trapz函數(shù)積分.由于MATLAB的數(shù)組運(yùn)算可以代替一個(gè)for循環(huán),所以三重積分只用了兩組for語句.使本題的程序比較簡明。,三重積分程序exn536,dx=0.01;dy=dx;dz=dx;x=0:dx:1; for k=1:length(x) x1=(k-1)*dx; y=0:dy:x1; for j=1:length(y) y1=(j-1)*dy; z1=0:dz:x1*y1; %z1數(shù)組 f=x1.*y1.2.*z1.3;%f(z1) s1(j)=trapz(f)*dz; %沿z1積分 end s2(k)=trapz(s1)*dy; % 沿y1積分 end, s=trapz(s2)*dx % 沿x1積分,六微分方程的數(shù)值積分,MATLAB中用來進(jìn)行常微分方程數(shù)值積分的函數(shù)有好多種，例如ode23,ode45,等，ode是常微分方程(ordinary differential equation)的縮寫。它們都用來解形如 的一階微分方程組在給定初始值y0時(shí)的解。對入門者而言，會(huì)一種最簡單的ode23就行。它的最簡單的調(diào)用格式為： T,Y = ODE23(ODEFUN,TSPAN,Y0) 其中，輸入變元TSPAN= t0,tf 是自變量的初值和終值數(shù)組，Y0是輸出變量向量的初值，ODEFUN則是描述導(dǎo)數(shù)的函數(shù)f(y,t)。很大一類微分方程都可以用這種一階微分方程組（或向量形式的微分方程）描述，關(guān)鍵就是會(huì)列寫導(dǎo)數(shù)函數(shù)f(y,t)，下面舉例說明。,微分方程數(shù)值積分【例5-3-7】,用數(shù)值積分法求解微分方程 設(shè)初始時(shí)間t0=0;終止時(shí)間tf=3; 初始條件y(0)=1,y(0)=0. 解：先將方程化為兩個(gè)一階微分方程的方程組，其左端為兩維變量的一階導(dǎo)數(shù)。,微分方程化為標(biāo)準(zhǔn)形式,寫成矩陣形式為 其中 為取代變量y的變量向量， 為x的 導(dǎo)數(shù)，在程序中用xdot表示。x的初始條件為 這就是待積分的微分方程組的標(biāo)準(zhǔn)形式。 用MATLAB語句表述為： xdot=0, 1;-t, 0*x + 0; 1*(1-t2/pi2);,【例5-3-7】數(shù)值解的程序,將微分方程的右端寫成一個(gè)exn547f.m函數(shù)程序,內(nèi)容如下： function xdot=exn547f(t,x) u=1-(t.2)/(pi2); xdot=0, 1;-t, 0*x + 0; 1*u; % 向量導(dǎo)數(shù)方程 主程序exn547如下,它調(diào)用MATLAB中的現(xiàn)成的數(shù)值積分函數(shù)ode23進(jìn)行積分。 clf, t0=0; tf=3*pi; x0=1; 0; % 給出初始值 t,x=ode23(exn547f, t0,tf, x0) % 此處顯示結(jié)果 y=x(:,1); % y為x的第一列 plot(t,y) ,grid % 繪曲線 xlabel(t),ylabel(y(t),數(shù)值解程序exn537的運(yùn)行結(jié)果,程序運(yùn)行的結(jié)果見圖5-37。這個(gè)數(shù)值積分函數(shù)是按精度要求自動(dòng)選擇步長的。它的默認(rèn)精度為1.e-3，因此圖中的積分結(jié)果是可靠的。 若要改變精度要求，可在調(diào)用命令中增加備選變元，具體做法可鍵入help ode23查找。,exn537的運(yùn)行結(jié)果的討論,從物理意義看，這個(gè)方程表示了一個(gè)變系數(shù)的無阻尼振動(dòng)方程。如果這是一個(gè)機(jī)械振動(dòng)，則彈簧剛度隨時(shí)間成正比地增強(qiáng)，振動(dòng)頻率隨之逐漸提高。為了看得更為清楚，設(shè)彈簧剛度隨時(shí)間按三次方增強(qiáng)，即方程的第二項(xiàng)系數(shù)為y3，則只要把子程序exn537f中的核心語句改為 xdot=0, 1; -t.3, 0*x + 0; 1*u; 重新運(yùn)行程序exn537，就得到頻率迅速提高的波形,如圖5-37。如果我們在原來的方程中加進(jìn)y的一階導(dǎo)數(shù)項(xiàng)（阻尼項(xiàng)），也只要在函數(shù)子程序中把矩陣的系數(shù)作些改動(dòng)，馬上就可以得出新的結(jié)果。由此可見，用計(jì)算機(jī)解題的極大優(yōu)越性。,七常系數(shù)線性微分方程的數(shù)值解,第四章4.3.5節(jié)介紹了常系數(shù)線性微分方程用MATLAB求解的問題。其實(shí)這類方程是有解析解的，這個(gè)解析解取決于微分方程系數(shù)多項(xiàng)式的根。而四次以上多項(xiàng)式的求根卻沒有解析解，這就要依靠MATLAB用數(shù)值方法解決代數(shù)問題，這個(gè)函數(shù)稱為residue，根據(jù)微分方程左端系數(shù)多項(xiàng)式a和右端系數(shù)多項(xiàng)式b，就可求根p和求留數(shù)r。 r,p=residue(b,a) 讀者可以參閱4.3.5節(jié)的算例，并可參閱第七章中的機(jī)械振動(dòng)和第九章中求系統(tǒng)響應(yīng)的例題,什么是剛性問題？ 在用微分方程描述的一個(gè)變化過程中，若往往又包含著多個(gè)相互 作用但變化速度相差十分懸殊的子過程，這樣一類過程就認(rèn)為具有 “剛性”。描述這類過程的微分方程初值問題稱為“剛性問題”。例 如，宇航飛行器自動(dòng)控制系統(tǒng)一般包含兩個(gè)相互作用但效應(yīng)速度相差 十分懸殊的子系統(tǒng)，一個(gè)是控制飛行器質(zhì)心運(yùn)動(dòng)的系統(tǒng)，當(dāng)飛行器速 度較大時(shí)，質(zhì)心運(yùn)動(dòng)慣性較大，因而相對來說變化緩慢；另一個(gè)是控 制飛行器運(yùn)動(dòng)姿態(tài)的系統(tǒng)，由于慣性小，相對來說變化很快，因而整 個(gè)系統(tǒng)就是一個(gè)剛性系統(tǒng)。,八. 符號數(shù)學(xué)求不定積分,不定積分問題要用符號數(shù)學(xué)的公式推理功能來解決問題。而根據(jù)本書的指導(dǎo)思想和風(fēng)格，主要強(qiáng)調(diào)數(shù)值計(jì)算和計(jì)算的道理，不把公式推理放在主要的地位。但是工作中如果遇到這種需要，還是應(yīng)該利用符號數(shù)學(xué)的功能來解決，就算是查積分表也是應(yīng)該會(huì)查的。所以也大致地介紹一下例題的類型和解法。,符號數(shù)學(xué)解不定積分例 538(a),例5-3-8(a)。求不定積分 解：因?yàn)槭遣欢ǚe分，不能用數(shù)值方法計(jì)算，只能用符號數(shù)學(xué)工具箱。程序?yàn)椋?syms x, y=x2*atan(x), Z=int(y) 得到 y = x2*atan(x) Z = 1/3*x3*atan(x)-1/6*x2+1/6*log(x2+1),符號數(shù)學(xué)求不定積分例 538(b),例5-3-8(b)。解下列積分，畫出解的曲線。 解：程序?yàn)?syms x Y=int(cos(x)2