常微分方程數(shù)值解實驗報告

1、常微分方程數(shù)值解實驗報告 數(shù)學(xué)與信息科學(xué) 信息與計算科學(xué) 鄭思義 學(xué)號 課程 201216524 常微分方程數(shù)值解 實驗一：常微分方程的數(shù)值解法 1、分別用Euler法、改進的Euler法（預(yù)報校正格式）和S K法求 解初值問題。（h=0.1 ）并與真解作比較。 y 1.1實驗代碼: %歐拉法 fun cti onx,y =n aeuler(dyfu n,xspan,y 0,h) %dyfun是常微分方程，xspan是x的取值范圍, x=xs pan( 1):h:xs pan( 2); y(1)=y0; yo 是初值, 曰卜【/. 疋步長 歡迎下載13 forn=1:le ngth(x)-1

2、y(n+1)=y( n)+h*feval(dyfu n, x( n), y( n); end %改進的歐拉法 fun cti onx,m, y=n aeuler2(dyfu n,xspan,y 0,h) %dyfun是常微分方程，xspan是x的取值范圍, %返回值x為x取值，m為預(yù)報解，y為校正解 x=xs pan( 1):h:xs pan( 2); y(1)=y0; m=zeros(le ngth(x)-1,1); yo 是初值, 曰卜【/. 疋步長。 forn=1:le ngth(x)-1 k1=feval(dyfu n,x(n ),y (n); y( n+1)=y( n)+h*k1;

3、m( n)=y( n+1); k2=feval(dyfu n,x( n+1),y( n+1); y( n+1)=y( n)+h*(k1+k2)/2; end %四階SK法 fun cti onx,y=rk(dyfu n,xspan,y 0,h) %dyfun是常微分方程，xspan是x的取值范圍, x=xs pan( 1):h:xs pan( 2); y(1)=y0; for n=1:le ngth(x)-1 yo 是初值, 曰卜【/. 疋步長。 k1=feval(dyfu n,x(n ),y (n); k2=feval(dyfu n,x( n)+h/2,y( n)+(h*k1)/2); k3

4、=feval(dyfu n,x( n)+h/2,y( n)+(h*k2)/2); k4=feval(dyfu n,x( n)+h,y( n)+h*k3); y(n+1)=y( n)+(h/6)*(k1+2*k2+2*k3+k4); end %主程序 x=0:0.1:1;y=ex p(-x)+x; dyfu n=in li ne( -y+x+1); x1,y1=naeuler(dyfu n,0,1,1,0.1); x2,m,y2=naeuler2(dyfu n,0,1,1,0.1); x3,y3=rk(dyfu n, 0,1,1,0.1); P lot(x,y, xlabel( r,x1,y1,

5、 x );ylabel( y為真解，y1 *,x3,y3,o); lege nd( ,Location, NorthWest +,x2,y2, y); 為歐拉解,y2為改進歐拉解,y3為S K解 )； 1.2實驗結(jié)果: x 真解y 歐拉解y1 預(yù)報值m 校正值y2 S K 解 y3 0.0 1.0000 1.0000 1.0000 1.0000 0.1 1.0048 1.0000 1.0000 1.0050 1.0048 0.2 1.0187 1.0100 1.0145 1.0190 1.0187 0.3 1.0408 1.0290 1.0371 1.0412 1.0408 0.4 1.070

6、3 1.0561 1.0671 1.0708 1.0703 0.5 1.1065 1.0905 1.1037 1.1071 1.1065 0.6 1.1488 1.1314 1.1464 1.1494 1.1488 0.7 1.1966 1.1783 1.1945 1.1972 1.1966 0.8 1.2493 1.2305 1.2475 1.2500 1.2493 0.9 1.3066 1.2874 1.3050 1.3072 1.3066 1.0 1.3679 1.3487 1.3665 1.3685 1.3679 2、選取一種理論上收斂但是不穩(wěn)定的算法對問題1進行計算，并與 真解作比較

7、。（選改進的歐拉法） 2.1實驗思路：算法的穩(wěn)定性是與步長h密切相關(guān)的。而對于問題 一而言，取定步長 h=0.1不論是單步法或低階多步法都是穩(wěn)定的算 法。所以考慮改變 h取值范圍，借此分析不同步長會對結(jié)果造成什么 影響。故依次米用 h=2.0、2.2、2.4、2.6的改進歐拉法。 2.2實驗代碼: %主程序 x=0:3:30;y=ex p(-x)+x; dyfu n=in li ne(-y+x+1 )； x1,m1,y1=naeuler2(dyfu n,0,20,1,2); x2,m2,y2=naeuler2(dyfu n,0,22,1,2.2); x3,m3,y3=naeuler2(dyfu

8、 n,0,24,1,2.4); x4,m4,y4=naeuler2(dyfu n,0,26,1,2.6); sub plot(2,2,1) p lot(x,y, r,x1,y1,+ );xlabel( sub plot(2,2,2) h=2.0 ）； plot(x,y,r ,x2,y2, + );xlabel( h=2.2 sub plot(2,2,3) plot(x,y,r ,x3,y3, + );xlabel( h=2.4 sub plot(2,2,4) plot(x,y,r ,x4,y4. + );xlabel( h=2.6 ); ); ); 2.3實驗結(jié)果: x h=2.0 h=2.2

9、 h=2.4 h=2.6 0.0 1.0000 1.0000 1.0000 1.0000 0.1 3.0000 3.4200 3.8800 4.3800 0.2 5.0000 5.8884 6.9904 8.3684 0.3 7.0000 8.4158 10.4418 13.4398 0.4 9.0000 11.0153 14.3979 20.4388 0.5 11.0000 13.7027 19.1008 30.8690 0.6 13.0000 16.4973 24.9092 47.4068 0.7 15.0000 19.4227 32.3536 74.8161 0.8 17.0000 22

10、.5077 42.2194 121.5767 0.9 19.0000 25.7874 55.6687 202.7825 1.0 21.0000 29.3046 74.4217 345.3008 實驗結(jié)果分析： 從實驗 1 結(jié)果可以看出，在算法滿足收斂性和穩(wěn)定性的前提下， Eluer 法雖然計算并不復(fù)雜，凡是精度不足，反觀改進的 Eluer 法和 s K法雖然計算略微復(fù)雜但是結(jié)果很精確。 實驗 2 改變了步長，導(dǎo)致算法理論上收斂但是不滿足穩(wěn)定性。 結(jié) 果表示步長 h 越大，結(jié)果越失真。對于同一個問題，步長 h 的選取變 得尤為重要， 這三種單步法的絕對穩(wěn)定區(qū)間并不一樣， 所以并沒有一 種方法是萬

11、能的，我們應(yīng)該根據(jù)不同的步長來選取合適的方法。 實驗二： Ritz-Galerkin 方法與有限差分法 1、用中心差分格式求解邊值問題 y y x 0,0 x 1 y(0) 0, y(1) 0 取步長 h=0.1 ，并與真解作比較。 1.1 實驗代碼： %中心差分法 function U=fdm(xspan,y0,y1,h) %xspan 為 x 取值范圍， y0,y1 為邊界條件， h 為步長 N=1/h; d=zeros(1,N-1); for i=1:N x(i)=xspan(1)+i*h; q(i)=1; f(i)=x(i); end for i=1:N-1 d(i)=q(i)*h*h

12、+2; end a=diag(d); b=zeros(N-1); c=zeros(N-1); for i=1:N-2 b(i+1,i)=-1; end for i=1:N-2 c(i,i+1)=-1; end A=a+b+c; for i=2:N-2 B(i,1)=f(i)*h*h; end B(1,1)=f(1)*h*h+y0; B(N-1,1)=f(N-1)*h*h+y1; U= in v(A)*B; %主程序 x=0:0.1:1; y=x+(ex p(1)*ex p(-x)/(ex p(2)-1)-(ex p(1)*ex p(x)/(ex p( 2)-1); y1=fdm(0,1,0,0

13、,0.1); y1=0,y1,0; P lot(x,y. r,x,y1, +) xlabel(x );ylabel( y); lege nd( y 真解，y1 中心差分法 ,Location, NorthWest ) 1.2實驗結(jié)果： x y真解 y1中心差分法 0.0 0.0000 0.0000 0.1 0.0148 0.0148 0.2 0.0287 0.0287 0.3 0.0409 0.0408 0.4 0.0505 0.0504 0.5 0.0566 0.0565 0.6 0.0583 0.0582 0.7 0.0545 0.0545 0.8 0.0443 0.0443 0.9 0.

14、0265 0.0265 1.0 0.0000 0.0000 2、用Ritz-Galerkin方法求解上述問題，并與真值作比較，列表畫 圖。 2.1實驗代碼: %Ritz Galerkin 法 fun cti on vu=Ritz_Galerki n(x0,y0,x1,y1,h) %x0,x1為x取值范圍，y0,y1為邊界條件，h為步長 N=1/h; syms x; for i=1:N fai(i)=x*(1-x)*(xA(i-1); dfai(i)=diff(x*(1-x)*(xA(i-1); end for i=1:N for j=1:N fun=dfai(i)*dfaia)+fai(i)*

15、fai(j); A(i,j)=int(fun, x,0,1); end fun=x*fai(i)+dfai(i); f(i)=in t(fu n, x,0,1); end c=i nv(A)*f： p roduct=c.*fai： sum=0; for i=1:N sum=sum+product(i); end vu=; for y=0:h:1 v=subs(sum,x,y); vu=vu,v; end y=0:h:1; yy=0:0.1:1; u=si n(yy)/si n(1)-yy; u=v pa(u,5); vu=v pa(vu,5); %主程序 x=0:0.1:1; y=x+(ex

16、p(1)*ex p(-x)/(ex p(2)-1)-(ex p(1)*ex p(x)/(ex p( 2)-1); y1=Ritz_Galerki n(0,0,1,0,0.1); y1=double(y1); P lot(x,y. r ,x,y1,+) xlabel( lege nd( x );ylabel( y); y 為真解,y1 為 R G 法，Location ,NorthWest ); 2.2實驗結(jié)果: x y真解 y1R G 法 0.0 0.0000 0.0000 0.1 0.0148 0.0148 0.2 0.0287 0.0287 0.3 0.0409 0.0409 0.4 0.0505 0.0505 0.5 0.0566 0.0566 0.6 0.0583 0.0583 0.7 0.0545 0.0545 0.8 0.0443 0.0443 0.9 0.0265 0.0265 1.0 0.000