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

1、常微分方程的數(shù)值解法專業(yè)班級：信息軟件姓名：吳中原學(xué)號：120108010002、實驗?zāi)康?、熟悉各種初值問題的算法，編出算法程序；2、明確各種算法的精度與所選步長有密切關(guān)系；通過計算更加了解各種算法的優(yōu)越性。、實驗題目1、根據(jù)初值問題數(shù)值算法，分別選擇二個初值問題編程計算；2、試分別取不同步長，考察某節(jié)點Xj處數(shù)值解的誤差變化情況;3、試用不同算法求解某初值問題，結(jié)果有何異常；4、分析各個算法的優(yōu)缺點。三、實驗原理與理論基礎(chǔ)一)歐拉法算法設(shè)計對常微分方程初始問題d二f(x,y)(6-1)(6-2)<dx、y(%)二yo用數(shù)值方法求解時，我們總是認(rèn)為(6-1)、(6-2)的解存在且唯歐拉

2、法是解初值問題的最簡單的數(shù)值方法。從(6-2)式由于y(x0)=y0已給定,因而可以算出y'(x)二f(x,y)。000設(shè)x1=h充分小，則近似地有：y0)(6-3),ni記y=y(x)i=0,1,ii從而我們可以取y二y二hf(x,y)1000作為y(xi)的近似值。利用y1及fg,yj又可以算出y(x2)的近似值:y2=yi+hf(片,yi)一般地，在任意點x=(n+1)h處y(x)的近似值由下式給出n+1(6-4)y=y+hf(x,y)n+1nnn這就是歐拉法的計算公式，h稱為步長。二)四階龍格-庫塔法算法設(shè)計：歐拉公式可以改寫為：廣1甘(1)，匕每步計算f(x,y)的值次，截1

3、1ii斷誤差為o)。y=y+(k+k)i+1i22改進(jìn)的歐拉公式可以改寫為：<k=hf(x,y)，它每一步要計算k=hf(x+h,y+k)2ii1f(x,y)的值兩次，截斷誤差為o(h3)。改進(jìn)的歐拉方法之所以比歐拉方法具有更高的精度，是因為在每一步它都比歐拉方法多計算了一次f(x,y)的值。因此，要進(jìn)一步提高精度，可以考慮在每一步增加計算f(x,y)的次數(shù)。如果考慮在每一步計算f(x,y)的值四次，則可以推得如下公式：y=y+1(k+2k+2k(xi+1i61231iih1k=hfx+,y+k2Ij2i21丿(h1k=hfx+,y+k3Ji2i22丿k=hf(x,y)k4hfi+h,y

4、+k)i3此公式稱為標(biāo)準(zhǔn)四階龍格-庫塔(Runge-Kutta)公式，它的截斷誤差為o(h5)雖然用龍格-庫塔方法每一步需要四次調(diào)用f，計算量較改進(jìn)的歐拉方法大一倍，這里由于龍格-庫塔方法的步長增大了一倍，因而兩種方法總的計算量相同，但龍格-庫塔方法精確度更高。所以龍格-庫塔公式兼顧了精度和計算工作量的較為理想的公式，在實際計算中最為常用。四、實驗內(nèi)容(一)問題重述：科學(xué)計算中經(jīng)常遇到微分方程(組)初值問題，需要利用Euler法，改進(jìn)Euler法，Rung-Kutta方法求其數(shù)值解，諸如以下問題：1),4xy=xy<yyG)=02)0<x<1分別取h=0.1,0.2,0.4時

5、數(shù)值解。初值問題的精確解y7+5e-x2。用r=3的Adams顯式和預(yù)-校式求解y'=X2-y23)、y(-1)=01<x<0取步長h=0.1,用四階標(biāo)準(zhǔn)R-K方法求值。用改進(jìn)Euler法或四階標(biāo)準(zhǔn)R-K方法求解y1G)=-1y2G)=0y3G)=10<x<1(、0.01,計算y血)，血10)；血15)數(shù)值解，參y615)q0.9880787,y(0.15)q0.1493359,y(0.15人0.8613125123(4)利用四階標(biāo)準(zhǔn)R-K方法求二階方程初值問題的數(shù)值解3y，+2y=0=0,y，G)=10<x<1h=0.02y1=y2<y2=y

6、1、y3=-y3取步長h=0.9880787,y考結(jié)果I)(II)y-0.1(1-y2)y'+y=0、yG)=1,y(0)=00<x<1,h=0.1(III)y'二-7<ex+1y(0)=1,y心)=00<x<2,h=0.1y"+siny=0(IV)M)=i,y心)=0(二)實驗代碼：1、歐拉法程序functiony=Euler(a,b,M,y0)%a=1,b=2,M=10,f=t*yA(1/3),y0=1;h=(b-a)/M;t=zeros(1,M+1);t=a:h:b;y=zeros(1,M+1);yy=zeros(1,M+1);y(

7、1)=y0;fork=1:My(k+1)=y(k)+h*t(k)*y(k)人(1/3);endyb=y(M+1);0<x<4,h=0.2yy=(t.A2+2)./3)41.5;det=yy-y;plot(t,y'r-',t,yy,'b:',t,det);2、改進(jìn)歐拉法程序functionH=heeuler(a,b,M,ya,f)%a=0,b=1,M=10,f=t*t+t-y,y0=0;h=(b-a)/M;t=zeros(1,M+1);y=zeros(1,M+1);p=0;q=0;t=a:h:b;y(1)=ya;fork=1:Mp=feval(f,t(

8、k),y(k);q=feval(f,t(k+1),y(k)+h*p);y(k+1)=y(k)+0.5*h*(p+q);endyy=t.*t-t+1-exp(-t);det=yy-y;plot(t,y'r-',t,yy,'b:',t,det);H=t',y',yy',det'functionf=ff(t,y);f=t.入2+t-y;3、四階龍格-庫塔法程序functionH=r_k4(a,b,M,ya,f)%a=0,b=1,M=10,f=t*t+t-y,y0=0;h=(b-a)/M;t=zeros(1,M+1);t=a:h:b;y=

9、zeros(1,M+1);K1=0;K2=0;K3=0;K4=0;y(1)=ya;fork=1:MK1=feval(f,t(k),y(k);K2=feval(f,t(k)+0.5*h,y(k)+0.5*h*K1);K3=feval(f,t(k)+0.5*h,y(k)+0.5*h*K2);K4=feval(f,t(k)+h,y(k)+h*K3);y(k+1)=y(k)+1/6*(K1+2*K2+2*K3+K4);endyy=t.*t-t+1-exp(-t);det=yy-y;plot(t,y,t,yy,t,det);H=t',y',yy',det'function

10、f=ff(t,y);f=t.入2+t-y;五、實驗結(jié)果,4xy=-xy<y1)yG)=0分別取h=0.1,0.2,0.4時數(shù)值解。初值問題的精確解y='''4+5e-x2oEuler('han',0,0,2,0.1)ans=000.100000.20000.04000.30000.11920.40000.23560.50000.38620.60000.56690.70000.77290.80000.99880.90001.23891.00001.48741.10001.73861.20001.98741.30002.22891.40002.4591

11、1.50002.67491.60002.87361.70003.05391.80003.21471.90003.35612.00003.4784Euler('han',0,0,2,0.2)ans=000.200000.40000.16000.60000.46720.80000.89111.00001.38861.20001.91081.40002.41221.60002.85681.80003.22262.00003.5025Euler('han',0,0,2,0.4)ans=000.400000.80000.64001.20001.71521.60002.81192.00003.57232、四階龍格-庫塔法結(jié)果y'=X2-y2V、y(1L0-1<x<0取步長h=0.1,用四階標(biāo)準(zhǔn)R-K方法求值。RK4('han',-1,0,1,0.1)ans=-1.00001.0000-0.90000.9909-0.80000.9672-0.70000.9331-0.60000.8921-0.50000.8468-0.