常微分方程初值問題數值解法的比較

1、常微分方程初值問題數值解法的比較數值計算實踐一課程設計報告課題名稱常微分方程初值問題數值解法的比較完成時間2013-1-17姓名班級學號成績一.實驗目的及內容1實驗目的：(1)了解常微分方程初值問題的理論背景以及初值問題穩(wěn)定性、收斂性的研究；熟練掌握歐拉法、改進歐拉法、龍格-庫塔法以及截斷誤差分析；比較歐拉法、改進歐拉法及龍格-庫塔法，能夠選擇合適的方法進行問題的研究計算； ,2實驗內容：求微分方程yy (0 X 1)(歐拉法求解)y(0) = 1_ 2X求微分方程yy y (0 x 0,使1得| f (x, y1) | -f (x, y2) L | y1 - y2 |, yyy y2 w R

2、,則稱f關于y滿足利普希茨條件，L稱為利普希茨常數。對于常微分方程初值問題y = f (t, y), t _ 10 , 一、,00 ,考慮初值y0的擾動是問題的解y(t)發(fā)生偏差的情y(t0) = v。形。若t t 笛時y(t)的偏差被控制在有界范圍內，則稱該初值問題是穩(wěn)定的，否則該初值問題 不穩(wěn)定的。特別地，若t tm時y(t )的偏差收斂于零，則稱該初值問題是漸進穩(wěn)定的。對于初值問題Jy - y ,t _ 10 _ ，一(t 3),、,0穩(wěn)定性的研究，易知其準確解為 y(t) = y留*0),假定初值經過y(t0) = v。擾動后變?yōu)閥0十Ay0 ,對于擾動后的解為y(t) = (y 0

3、+ Ay0)e*,0)因此帶來的擾動誤差為y(t) = Ay0e其上0),因此考慮t t 餡時Ay(t )的值，它取決于 九。易知，若 九w 0 ,則原問題是穩(wěn)定的；若 九a 0,則原問題是不穩(wěn)定的；若九 0,則原問題是漸進穩(wěn)定的。實際遇到的大多數常微分方程初值問題都是穩(wěn)定的，因此在后面的討論數值解法時這常常是默認條件。依據：積分曲線上一點(x, y)的切線斜率等于函數值。方法：推進法，初始點 p0(x0, y0)出發(fā)，依照方向場在改點的方向推進到xi,x2，xn向前歐拉法的得到：y(x)y(xn i) = y(xn h) = y(xn)y (xn)hy,(xn)h2y (xn)h2!十取h的

4、線性部分，得Yn 1 二 Yn 八(。*)2 ) 將初值問題中得導數用差商來代替y(xn 1) - y(xn)y(Xn 1)- y(xn)xn 1 - xn=f(xn,yn)二ynhf(xn,yn)公dy將丁 =f(x,y)兩邊同 dx時對 x 的區(qū)間xn,xn+i上 積xn 1 y dxxnxn 1-x n _1=p(xn, yn)dx,即y(xn書)-y(xn) = Jf (xn, y(x)dx 對右端x nxnxn 1ff (xn, y(x)dx 用左xn矩形公式得yn書=yn +hf(xn,yn),此方法稱向前歐拉法，也叫 顯示歐拉法。Xn 1(4)對右端 jf (xn, y(x)dx

5、用右矩形公式得yn卡=yn + hf (xn#, y5)，也 xn叫隱式歐拉法。誤差分析：1.稱y(xnQ _yn4為計算yn邛時的局部截斷誤差； TOC o 1-5 h z 2.如果數值方法的局部截斷誤差為 qhp*),那么稱這種數值方法的階數是p,其實p為非負整數。通常情況下，步長 h越小，p越高，則局部截斷誤差越小，計算精,，, y (xn)hy (xn)h HYPERLINK l bookmark34 o Current Document .、生叱由療y(xn1) =y(xnh)二y(xn)y(xn)h-nn .y(x)在xn初泰勒展開有八 八n /八 y ”2!32.=ynhf(xn

6、, yn) qh )則有y(xn卡) -yn4i = Qh2)可見歐拉方法是一階方法，精度不是很高。2.改進歐拉方法：xn .1梯形公式：對右端Jf (xn, y(x)dx用梯形公式得+顯然梯形公式是隱式公式。 xn改進歐拉公式： 先用歐拉公式求的一個初步的近似值yn書，成為預測值，預測值 yn卡的精度可能達不到要求，在用梯形公式將他校正一次，記為yn41,這個結果成為校正值。預測：Yn .1 = Yn , hf(xn,yn)h校正：yn l = yn 21f(xn, yn) f(xn i, yn l)誤差分析：記Tn中為改進歐拉公式在xn 4i處的截斷誤差，Tn+1 = y( xn+)- y

7、(xn )丁一，、 h記 Tn 書=y(xn/ - y(xn) - - f ( x n , y( x n ) + f (x n 41 , y( x n+1 )因此有h 一Tn4=y(xn書)-y(xn)=，書+萬f ( x n如y ( x n由)f ( x n由，y( x n+),八中表不在*門書出的局到相應的計算公式。歐拉公式可以寫為Ki = f (Xn，yn)h 、y n 4i =yn ( Ki + K2 )改進歐拉公式可以寫成 Ki = f ( Xn, y n)因此推出一般的推出廣式K2 = f(Xn+i,yn + hKi)yn 平=yn +h(CiKi +C2K2 + . +CPKP)

8、Ki = f (Xn，yn)jK2 =f(xn +a2h,yn + hb21Kl)稱為p階龍格-庫塔方法，簡稱p階R-K方法。.pKP = f(Xn +aph, yn + h bp K )1i金r(Xn，yn，h) = Ci Ki， i i因為 Ki =f(Xn,yn),這里的6,兒，也均為常數。i JKi = f(Xn -h,ynhij Kj )j i因為給定的系數不唯一，因此這里的常數有無窮多個解，下面是特殊情況下和一般情況下得結果(一階龍格-庫塔)當r=i時，這就是歐拉法。,i,(二階龍格-庫塔)當r=2時，Ci = , % = = i ,這就是改進歐拉法。2j三階和四階龍格-庫塔也只是

9、在一般情況下得結果。hyn 書=yn + -(Ki +K2 + K3) 6Ki =f(Xn，yn)三階K2 =f(Xn + , +2Ki)K3 =f(Xn +h,y0- hKi +2h9)yn書=yn十h(K +26 +2K3 + K4) 6Ki =f(Xn，yn)rhh四階TK2 =f(Xn + 2，yn + 2 Ki)K2 =f(Xn + -2，yn + hKz) 2K3 =f(Xn + h,yn 十 h()Tn由=y(Xn中)-y(Xn) -h(4-iki + %kz)其局部截斷誤差是:三.程序代碼歐拉法代碼如下：(1)向前歐拉法：function y=Euler1(fun,x0,y0,

10、xN,N)%fun為一階微分方程，x0,y0為初始條件，xN為取值范圍的一個端點，N為區(qū)間個數x=zeros(1,N+1);x=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N; % h為區(qū)間步長for (n=1:N)x(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n);%艮據向前歐拉公式計算y 值endT=x,y(2)向后歐拉法：function y=Euler2(fun,x0,y0,xN,N)% fun為一階微分方程，x0 , y0為初始條件，xN為取值范圍的一個端點，N為區(qū)間個數x=zeros(1,N+1);x=z

11、eros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;% h為區(qū)間步長for (n=1:N)x(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n);for (k=1:3)z1=y(n)+h*feval(fun,x(n+1),z0);if (abs(z1-z0)1e-3) break ;endz0=z1;endy(n+1)=z1;%艮據向后歐拉公式計算y值endT=x,y改進歐拉代碼如下：function x,y=Gaijineuler(f,x0,y0,xZ,h)%f為一階微分方程的函數；x0,y0為初始條件；xZ為取值范圍的一個端點，h

12、為區(qū)間步長n=fix(xZ-x0)/h);%計算分點數y(1)=y0;x(1)=x0;for i=1:nx(i+1)=x0+i*h;yp=y(i)+h*feval(f,x(i),y(i);yc=y(i)+h*feval(f,x(i+1),yp);y(i+1)=(yp+yc)/2;endx=x;y=y;1.3.龍格-庫塔代碼如下： (三階龍格-庫塔)function R=Longgekuta3(f,a,b,aZ,h)%a, b為端點，h為步長，aZ為初值 n=(b-a)/h;T=zeros(1,n+1);Y=zeros(1,n+1);T=a:h:b;Y(1)=aZ;for j=1:nk1=fev

13、al(f,T(j),Y(j);k2=feval(f,TQ)+h/2,Y(j)+k1*h/2);k3=feval(f,T(j)+h,Y(j)-h*k1+k2*2*h);Y(j+1)=Y(j)+(k1+4*k2+k3)*h/6;endR=T Y;(四階龍格-庫塔)function R=Longgekuta4(f,a,b,aZ,h)% a, b為端點，h為步長，aZ為初值 n=(b-a)/h;T=zeros(1,n+1);Y=zeros(1,n+1);T=a:h:b;Y(1)=aZ;for j=1:nk1=feval(f,T(j),Y(j);k2=feval(f,T(j)+h/2,Y(j)+k1*h

14、/2);k3=feval(f,T(j)+h/2,Y(j)+k2*h/2);k4=feval(f,T(j)+h,Y(j)+k3*h);Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)*h/6;endR=T Y;四.數值結果：% 艮據改進歐拉公式計算結果%定義向量%計算各個分點%初值賦予%根據公式計算%定義向量%計算各個分點%初值賦予%根據公式計算輸入：定義M文件并保存ffx.mEuler1(ffx,0,1,1,10)結果：T =01.00000.10001.10000.20001.21000.30001.33100.40001.46410.50001.61050.60001.77160

15、.70001.94870.80002.14360.90002.35791.00002.5937ans =Columns 1 through 71.77161.00001.10001.21001.33101.46411.6105Columns 8 through 111.94872.14362.35792.5937向前歐拉公式結果：輸入：Euler2(ffx,0,1,1,10)結果T =01.00000.10001.11100.20001.23440.30001.37160.40001.52400.50001.69330.60001.88140.70002.09040.80002.32270.9

16、0002.58071.00002.8674ans =Columns 1 through 71.00001.11101.23441.37161.52401.69331.8814Columns 8 through 112.09042.32272.58072.8674改進歐拉公式結果：輸入：x,y=Gaijineuler(f,0,1,1,0.1)結果：x =00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000y =1.00001.09591.18411.26621.34341.41641.48601.55251.61651.678

17、21.7379龍格-庫塔計算結果：(三階)輸入：Longgekuta3(ff,0,1,1,0.1)結果：ans =Longgekuta3(ff,0,1,1,0.1)8ans =00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000(四階)輸入：結果:00.10000.20000.30000.40000.50000.60001.00001.10481.21881.34111.47111.60821.75201.90212.05802.21952.3863Longgekuta4(ff,0,1,1,0.1) ans =1.00001

18、.10481.21881.34111.47111.60821.75200.70001.90210.80002.05800.90002.21951.00002.3863五.計算結果分析方法顯示歐拉簡單；精度低隱式歐拉穩(wěn)定性最好精度低，計算量大梯形公式精度提高計算量大中點公式精度提高，顯式:多一個初值，可能影響精度1.收斂性：若某算法對于任意固定的 x = Xi = x0 + i h,當20 (同時i t8)時有yiT y( xi),則稱該算法是收斂的。以下討論的都是單步法(指在計算 yn 卡時只用到它前一步的信息yn)歐拉法，龍格-庫塔法都 是單步法的例子； V =,t 之 t C例：就初值問題

19、 y 0考察歐拉顯式格式的收斂性。了(t ) = y解：該問題的精確解為y(t) = y0e *工)歐拉公式為yi中=yi+h，y =(1十h九)yixi1對于固定的 X = Xi = X0+ i h,有 v =y0(1+Kh)h = y0( 1 + Ih)沖NNoe ，= y( Xi )設單步法具有中階精度，滿足利普希茨條件，其整體的截斷誤差：y(Xn) - yn = qhp)判斷單步法的收斂性，歸結為驗證增量函數邛能否滿足利普希茨條件。對于歐拉法，由于其增量函數 中就是f (x, y),因此當f (x, y)關于y滿足利普希茨條件 時候它就是收斂的。對于改進歐拉法和龍格-庫塔方法，找到增量

20、函數，滿足利普希茨條件的時候，找到利普希 茨常數，這些方法都是收斂的。2.穩(wěn)定性:r 例：考察初值問題 ”=一 y,在區(qū)間0,0.5上的解，分別用歐拉顯，隱式格式，改進歐拉格式計 7(0) = 1算數值解：10節(jié)點歐拉顯式改進歐拉法精 確 解30 xy = e-01.00001.00001.00001.00000.1-2.00002.5*10A(-1)2.50004.9787*10A(-2)0.24.00006.2500*10A(-2)6.25002.4788*10A(-3)0.3-8.00001.5625*10A(-2)1.5625*10A(1)1.2341*10A(-4)0.4.6000*

21、10A(1)3.9063*10A(-3)3.9063*10A(1)16.1442*10A(-6)0.513.2000*10A(1)9.7656*10A(-1)9.7656*10A(1)3.0590*10A(-7)則稱該算法是絕對穩(wěn)若某算法在計算過程中任一步產生的誤差在以后的計算中都逐步衰減, 定的.i 1yi 1 = V h yi =(1h)y, ； = yyi 1=(1+ h)i卡斯,由此可見,要保證初始誤差a以后逐步衰減h =Kh必須滿足:| 1 h | ： 1考察后退歐拉法：yi卅=yih yi 1玩由此可見,1梯形公式：模型y=y, yn .1yn，h,、+ -(y n + yn +)

22、,穩(wěn)定區(qū)域為| 2 + Ah |42 九h |相對穩(wěn)定區(qū)間為-二：h改進歐拉公式:yn相對穩(wěn)定區(qū)域為| 1四階龍格-庫塔:相對穩(wěn)定區(qū)域:1相對穩(wěn)定區(qū)間:yn 1yn 1ynhf(Xn,yn) = (1h)yn,hf(Xn，yn) f (Xn 1, yn 1 2,1, 2)=1 h 2( h)2yn1 ,. 、2 ,-, 八+ -(油)2 |M 1穩(wěn)定區(qū)間為2 Kh 021c 1a 1“H1 h 2( h)26( h)324( -h)4yn121314h 2(h)6(h)24(h)L- 2. 78 :二 h ：二 01 - h J據對絕對穩(wěn)定絕對穩(wěn)定 絕對穩(wěn)定的區(qū)域為| 1 - h | . 1節(jié)

23、點改進歐拉法四節(jié)龍格-庫塔法準確解y(0) = 1-12x,、，“在區(qū)間0,1上的解考察初值問題y1101110.21.1866671.1832291.1832160.41.3483121.3416671.3416410.61.4937041.4832811.4832400.81.6278611.6125141.61245211.7542051.7321421.732051變步長的龍格-庫塔法：實際計算中，步長h的選擇是一個十分重要的問題，若步長h取的太小，將增加許多不必要的的計算量，若步長取的太大，其計算競速很難保證，因此，與數值積分一樣，微分方程的數值解也有一個選擇步長的問題。由四階龍格-

24、庫塔法的性質，其局部截斷誤差為y(xn+1) y：?定ch5這h()一個近似值yn.,n里的c與y5(x)在xn,xn中內的值有關系。然后將步長折半，即取為步長，從xn跨兩步到xn中，在求h 5、一22)5 ,比較之后h 5(；)每跨一步的誤差為ci-J ,因此y(xn中)_yn1一誤差的大約減少16116h (夕) y(xn .1)- yn2.1y( x n .1) - y n 六.計算時出現的問題，解決方法及體會：對于初值問題，在用不同方法計算離散點的函數值的時候，有精確值給定的情況下，不同方法跟精確值 的差別很大，也就是不是一種好方法，這與實際問題的穩(wěn)定性有很大的關系，只要所求的初值問題

25、，能夠落 在使用方法的穩(wěn)定區(qū)間內，所得到的結果和精確值的離散點的函數值相差不會很大，并且在穩(wěn)定區(qū)間內，改 進歐拉法優(yōu)于歐拉法，三階之后的龍格-庫塔方法優(yōu)于改進歐拉法，離散點的函數值更加接近精確值。 所謂的歐拉法就是一階的龍格-庫塔方法，改進歐拉法也是在顯示歐拉法預測和梯形公式校正的基礎上作出的改進， 也是二階龍格-庫塔的特殊情況。所謂的三階，四階甚至n階龍格-庫塔方法，也是在待定系數有無數解得情況下的一般情況。一階常微分初值問題的收斂性，只要其增量函數關于y能滿足利普希茨條件，這就收斂。而起在計算過程過有避免不了的誤差， 歐拉法精度相對來說比較低，龍格-庫塔的階數越高，精度就越高，所得的結果和

26、精確值越接近，截斷誤差越小。同一個問題，想要減少截斷誤差，我們可以采用變步長的的方法，將步長變?yōu)閔h,這樣得出的結果，根據不同的方法，截斷誤差會有相應的改變。2而在每種方法的程序設計方面，主要思想還是根據每種方法的公式，找到步長，分點數，用循環(huán)語句， 得出相應的結果，只是在龍格 -庫塔方法里面的三階，四階的公式里面的待定系數有無數個解，我們給出的也 是一般情況的三階，四階龍格 -庫塔方法，精度，誤差相對于歐拉法，改進歐拉法，梯形公式都優(yōu)于他們。數值計算實踐一課程設計原始數據記錄實驗名稱：常微分方程初值問題數值解法研究 實驗時間： 2013年 1月 14 日12輸入：定義M文件并保存ffx.mE

27、uler1(ffx,0,1,1,10)結果：T =01.00000.10001.10000.20001.21000.30001.33100.40001.46410.50001.61050.60001.77160.70001.94870.80002.14360.90002.35791.00002.5937 ans =Columns 1 through 71.00001.10001.21001.33101.46411.61051.7716Columns 8 through 111.94872.14362.35792.5937向前歐拉公式結果：輸入：Euler2(ffx,0,1,1,10)結果T =

28、01.00000.10001.11100.20001.23440.30001.37160.40001.52400.50001.69330.60001.88140.70002.09040.80002.32270.90002.58071.00002.8674 ans =13Columns 1 through 71.00001.11101.23441.37161.52401.69331.8814Columns 8 through 112.09042.32272.58072.8674改進歐拉公式結果：輸入：x,y=Gaijineuler(f,0,1,1,0.1)結果：x =00.10000.20000.30000.40000.50000.600