計算方法第八章(常微分方程數(shù)值解)

1、第八章第八章 常微分方程初值常微分方程初值問題數(shù)值解法問題數(shù)值解法本章主要研究常微分方程初值問題的數(shù)值求解：本章主要研究常微分方程初值問題的數(shù)值求解：通常，假設函數(shù)通常，假設函數(shù) f 關于第二個變量滿足李普希茨條件（關于第二個變量滿足李普希茨條件（L條條件），即為存在常數(shù)件），即為存在常數(shù) L 0，使得使得0( )( , ( )( )y xf x y xaxby ay1212( ,)( ,)f x yf x yL yy第一節(jié)第一節(jié) 一般概念一般概念1.1 歐拉法及其簡單改進歐拉法及其簡單改進0( )( , ( )( )y xf x y xaxby ay方法：選擇適當?shù)墓?jié)點，用差分近似微分，將方

2、程離散化，方法：選擇適當?shù)墓?jié)點，用差分近似微分，將方程離散化，從而求在這些節(jié)點上的解的近似值。從而求在這些節(jié)點上的解的近似值。012111NNnnnnnaxxxxxbhxxxx稱為到的步長（通常取為常數(shù)h）1()()( )|nnnx xy xy xy xh歐拉歐拉方法方法100(,),()nnnnyyhf xyy xy，記yn為 y(xn）的近似計算值，有例子：例子：(01)(0)1yyxy 10( , )(1)1nnnnf x yyyyhyh yy xye精確解為下面我們分別取步長為下面我們分別取步長為0.1與與0.01進行計算，進行計算，計算結果顯示在下面的圖中。計算結果顯示在下面的圖中。

3、步長為步長為0.1的計算結果。的計算結果。步長為步長為0.01的計算結果的計算結果0.01 0.99005 0.99 0.1 0.90484 0.90438 0.2 0.81873 0.81791 0.3 0.74082 0.7397 0.4 0.67032 0.66897 0.41 0.66365 0.66228 0.59 0.55433 0.55268 0.6 0.54881 0.54716 0.9 0.40657 0.40473 0.91 0.40252 0.40068 0.99 0.37158 0.36973 1 0.36788 0.36603 DOUBLE PRECISION h,y

4、(0：100) OPEN(20,FILE=OUTPUT.DAT,STATUS=UNKNOWN)h=1.0/100y(0)=1.0do 10 i=1,100 y(i)=y(i-1)*(1.0-h)write(20,*) i*h,y(i)10 continueEND我們可得精度更高的歐拉公式：我們可得精度更高的歐拉公式：11()()()2nnny xy xy xh22()()( )()()212nnny xhy xhyy xhO hh 11002(,),()nnnnyyhf xyy xy歐拉中點公式歐拉中點公式 歐拉公式中我們利用了近似公式歐拉公式中我們利用了近似公式1()()()nnny xy

5、xy xh光這個近似產生的誤差為光這個近似產生的誤差為1()()1()( )( )2nnny xy xy xyhO hh 利用利用利用中點公式求解微分方程時，有一個問題，利用中點公式求解微分方程時，有一個問題，就是計算時需要兩個迭代初值！這樣的算法就是計算時需要兩個迭代初值！這樣的算法稱為二步法。前面的歐拉法稱為單步法。稱為二步法。前面的歐拉法稱為單步法。對于這個問題，我們可以先用歐拉公式，通過對于這個問題，我們可以先用歐拉公式，通過給定的初值計算出給定的初值計算出 的值，然后再利用這兩個的值，然后再利用這兩個值（值（ y0 和 y1 ）進行計算，直到計算出全部節(jié)）進行計算，直到計算出全部節(jié)點

6、上的值。點上的值。111(,)nnnnnyxxyy一般的單步法：一般的單步法：一般的一般的 k 步法：步法：1(,)n knn knn kn kyxxyyy 1.2 歐拉方法的其他改進歐拉方法的其他改進微分方程數(shù)值解的關鍵在于對導數(shù)的處理，可以用差分來近似微分方程數(shù)值解的關鍵在于對導數(shù)的處理，可以用差分來近似導數(shù)，也可以通過積分，將導數(shù)項化掉。導數(shù)，也可以通過積分，將導數(shù)項化掉。對于方程：對于方程：首先，作出劃分首先，作出劃分設已經求出第設已經求出第 n 個個節(jié)點的函數(shù)值節(jié)點的函數(shù)值 ，在區(qū)間，在區(qū)間 上對上對方程兩邊積分方程兩邊積分容易看出，要求第容易看出，要求第 n+1 個個節(jié)點的函數(shù)值，

7、關鍵在于節(jié)點的函數(shù)值，關鍵在于選擇適當?shù)倪x擇適當?shù)姆e分公式計算積分積分公式計算積分！( )( , ( )y xf x y x0121NNaxxxxxb1,nnxxny11()()( , ( )nnxnnxy xy xf x y x dx（1）如選擇下矩形公式，則得）如選擇下矩形公式，則得這正是前面的這正是前面的歐拉公式歐拉公式。1(,)nnnnyyf xy h111(,)nnnnyyf xyh111 (,)(,)/2nnnnnnyyf xyf xyh（2）如選擇上矩形公式，則）如選擇上矩形公式，則得得這是所謂的這是所謂的后退歐拉公式后退歐拉公式。（3）如選擇梯形公式，則）如選擇梯形公式，則得得

8、這是所謂的這是所謂的歐拉梯形公式歐拉梯形公式。直接利用已經求得的已知節(jié)點上的值計算未知節(jié)點上的函數(shù)直接利用已經求得的已知節(jié)點上的值計算未知節(jié)點上的函數(shù)值的算法稱為值的算法稱為顯式法顯式法。例如：歐拉公式、歐拉中點公式例如：歐拉公式、歐拉中點公式計算未知節(jié)點上的函數(shù)值時，用到了未知節(jié)點上的函數(shù)值，計算未知節(jié)點上的函數(shù)值時，用到了未知節(jié)點上的函數(shù)值，這種算法稱為這種算法稱為隱式法。隱式法。例如：后退歐拉法、歐拉梯形公式例如：后退歐拉法、歐拉梯形公式顯然，利用隱式法求微分方程的數(shù)值解是，顯然，利用隱式法求微分方程的數(shù)值解是，需要從表達式中需要從表達式中反解未知節(jié)點上的函數(shù)值反解未知節(jié)點上的函數(shù)值。1

9、.3 隱式法的具體計算：隱式法的具體計算：例如歐拉梯形公式例如歐拉梯形公式用迭代法計算用迭代法計算 yn+1 的值。的值。（1）簡單迭代）簡單迭代收斂的條件：收斂的條件：11111 (,)(,)/2(,)/2(,)/2nnnnnnnnnnnyyf xyf xyhyf xyhf xyh(0)1(,)nnnnyyhf xy(1)( )111(,)/2(,)/2kknnnnnnyyf xyhf xyh/21Lh（2）牛頓迭代）牛頓迭代若用簡單迭代，而且只迭代一步，這樣組成的一組計算公式稱若用簡單迭代，而且只迭代一步，這樣組成的一組計算公式稱為為預測校正公式預測校正公式。（迭代初值。（迭代初值 稱為預

10、測，迭代步稱為預測，迭代步稱為校正）稱為校正）( )( )(1)( )11111( )11(,)/2(,)/21(,)/2kkkknnnnnnnnkynnyyf xyhf xyhyyfxyh(0)1ny(0)1(0)111(,) (,)(,)/2nnnnnnnnnnyyhf xyyyf xyf xyh預測校正公式也稱為改進的歐拉法，將上面的組合公式改預測校正公式也稱為改進的歐拉法，將上面的組合公式改寫為：寫為：注意到注意到 ，將上式進一步改寫為：，將上式進一步改寫為：這是我們最終使用的計算格式。這是我們最終使用的計算格式。11 (,)(,(,)/2nnnnnnnnyyf xyf xyhf xy

11、h1121211()2(,)(,)nnnnnnyyKKKhf xyKhf xh yK1nnxxh例子：例子：取步長為取步長為0.1計算，結果如圖。計算，結果如圖。(01)(0)1yyxy 11212111()2(,)(,)()nnnnnnnnyyKKKhf xyhyKhf xh yKh yK 圖：圖： DOUBLE PRECISION h,y(0:10),ak1,ak2 OPEN(20,FILE=OUTPUT1.DAT,STATUS=UNKNOWN)h=1.0/10y(0)=1.0do 10 i=1,10 ak1=-h*y(i-1) ak2=-h*(y(i-1)+ak1) y(i)=y(i-1

12、)+(ak1+ak2)/2.010 continue do 20 i=0,10 write(20,*) i*h,y(i),exp(-i*h)20 continueEND同理，對于后退歐拉公式同理，對于后退歐拉公式有預測校正公式有預測校正公式或改寫為：或改寫為：111(,)nnnnyyf xyh(0)1(0)111(,)(,)nnnnnnnnyyhf xyyyf xyh11(,)(,)nnnnnnKhf xyyyf xyK h用此法解前面的例子用此法解前面的例子步長步長0.1步長步長0.011.4 誤差估計誤差估計定義：定義：利用第利用第n個節(jié)點或之前更多節(jié)點的函數(shù)精確值，利用近個節(jié)點或之前更多

13、節(jié)點的函數(shù)精確值，利用近似公式數(shù)值計算第似公式數(shù)值計算第n+1個節(jié)點的近似值，所引起的誤差，稱個節(jié)點的近似值，所引起的誤差，稱為第為第n+1個節(jié)點上的個節(jié)點上的局部截斷誤差局部截斷誤差。我們記我們記 為第為第n+1個節(jié)點上解的精確值，個節(jié)點上解的精確值， 為假設為假設 等條件下計算所得的近似值，等條件下計算所得的近似值，則局部截斷誤差為：則局部截斷誤差為：如局部截斷誤差為如局部截斷誤差為 ，稱為具有，稱為具有 p 階局部截斷誤差。階局部截斷誤差。1ny1()ny x111()nnny xy()pO h()nnyy x歐拉方法的誤差分析：歐拉方法的誤差分析：12()()()()()()( )/2

14、()()( ) /2nnnnnnnny xy xy xhy xhhy xy x hyhy xy xyhh21()(, ()()()(, ()()nnnnnnny xf xy xy xy xhf xy xO h1()(, ()nnnnyy xhf xy x2111()()nnny xyO h而而完全類似的可以得到完全類似的可以得到后退歐拉公式的局部截斷誤差為：后退歐拉公式的局部截斷誤差為：歐拉中點公式的局部截斷誤差為：歐拉中點公式的局部截斷誤差為：歐拉梯形公式的局部截斷誤差為：歐拉梯形公式的局部截斷誤差為：2111()()nnny xyO h3111()()nnny xyO h3111()()n

15、nny xyO h定義：由初值計算到定義：由初值計算到第第 n 個節(jié)點的近似值與其精確值個節(jié)點的近似值與其精確值之間的誤差稱為第之間的誤差稱為第 n 個節(jié)點整體誤差。個節(jié)點整體誤差。定理：設下面求解微分方程的數(shù)值計算方法定理：設下面求解微分方程的數(shù)值計算方法局部截斷誤差為局部截斷誤差為p+1階，且函數(shù)階，且函數(shù) 關于關于 y 滿足滿足利普希茨條件，利普希茨條件，同時初值是準確的，則同時初值是準確的，則整體截斷誤差為整體截斷誤差為p階階。歐拉公式、后退歐拉公式的整體誤差為歐拉公式、后退歐拉公式的整體誤差為 1 階。階。歐拉中點公式、歐拉梯形公式的整體誤差為歐拉中點公式、歐拉梯形公式的整體誤差為

16、2 階。階。1(, )nnnnyyhxy h( , , )x y h( , , )( , , )x y hx y hL yy微分方程數(shù)值解法的進一步改進。再回到恒等式微分方程數(shù)值解法的進一步改進。再回到恒等式如果取如果取 作為節(jié)點，將被積函數(shù)用插值多作為節(jié)點，將被積函數(shù)用插值多項式來近似，用插值多項式帶到積分中去求出積分，則可以得項式來近似，用插值多項式帶到積分中去求出積分，則可以得到所謂的到所謂的亞當斯亞當斯(Adams)顯式公式顯式公式局部截斷誤差：局部截斷誤差：11()( )( , ( )iixiixy xy xf x y x dx1,iii kx xx1011()iiiiki khyy

17、b fb fb fA1(1) ( )kkkiR yB hy類似地，如果取類似地，如果取 作為作為節(jié)點，可得節(jié)點，可得亞當斯亞當斯(Adams)隱式公式隱式公式局部截斷誤差：局部截斷誤差：11,iiii kxx xx*10112()iiiiiki khyyb fb fb fb fA*2(2)1 ( )kkkiR yBhy進一步，如果我們將恒等式中的積分區(qū)間改為進一步，如果我們將恒等式中的積分區(qū)間改為 ，并在此區(qū)間上用辛普森公式，并在此區(qū)間上用辛普森公式，1,i iixx111111111()()( , ( )() (, ()4 ( , ()(, ()3iixiixiiiiiiiy xy xf x

18、 y x dxhy xf xy xf x y xf xy x111143iiiiihyyfff可得可得辛普森公式辛普森公式1.5 絕對穩(wěn)定性絕對穩(wěn)定性一個常微分方程數(shù)值解法應用于方程一個常微分方程數(shù)值解法應用于方程 時對給定步長時對給定步長 h 稱為稱為絕對穩(wěn)定絕對穩(wěn)定，是指在某一步（如，是指在某一步（如第一步）產生的誤差（如計算機的存儲誤差），在第一步）產生的誤差（如計算機的存儲誤差），在計算中會逐步減小。計算中會逐步減小。稱方程稱方程 為為試驗方程試驗方程，設，設在計算開始時產生誤在計算開始時產生誤差（存儲誤差），此誤差在以后差（存儲誤差），此誤差在以后會逐步減小會逐步減小，我們，我們稱該

19、算法相對于稱該算法相對于 是絕對穩(wěn)定的，這樣的是絕對穩(wěn)定的，這樣的 的全體稱為該算法的的全體稱為該算法的絕對穩(wěn)定域絕對穩(wěn)定域。 通常在復平面上通常在復平面上考慮絕對穩(wěn)定域。穩(wěn)定域與實軸的交集為考慮絕對穩(wěn)定域。穩(wěn)定域與實軸的交集為穩(wěn)定區(qū)間穩(wěn)定區(qū)間yy hh( Re( )0 )yy h歐拉法的絕對穩(wěn)定區(qū)域歐拉法的絕對穩(wěn)定區(qū)域 后退歐拉法的絕對穩(wěn)定區(qū)域后退歐拉法的絕對穩(wěn)定區(qū)域11h11h絕對穩(wěn)定域越大，絕對穩(wěn)定域越大，h 的選取范圍就越大，算法就越好！的選取范圍就越大，算法就越好！1.6 局部截斷誤差的實用估計局部截斷誤差的實用估計（1）用兩種階數(shù)相同的算法求解，計算出）用兩種階數(shù)相同的算法求解，計

20、算出n+1步的近似值，步的近似值，從而得到局部截斷誤差估計。從而得到局部截斷誤差估計。（2）用同樣的公式，用不同步長計算出）用同樣的公式，用不同步長計算出n+1步的近似值，從步的近似值，從而得到局部截斷誤差估計。而得到局部截斷誤差估計。1.7 隱式法隱式法隱式法具有較好的絕對穩(wěn)定性隱式法具有較好的絕對穩(wěn)定性！只不過在使用隱式法的時候，需要進行迭代，或者使用預測只不過在使用隱式法的時候，需要進行迭代，或者使用預測校正計算格式。校正計算格式。第二節(jié)第二節(jié) 泰勒級數(shù)法與龍格庫塔法泰勒級數(shù)法與龍格庫塔法對于方程：對于方程：取計算步長取計算步長為為 h ，則則 ，將函數(shù)進行泰勒展開，將函數(shù)進行泰勒展開如

21、函數(shù)如函數(shù) y( x )有有p+1階導數(shù)，容易得到階導數(shù)，容易得到p階泰勒級數(shù)展開法：階泰勒級數(shù)展開法：0( )( , ( )( )y xf x y xaxby ay1nnxxh21()()()()()/2!nnnnny xy xhy xhy xh y x( )21(1)12!( )(, )(1)!ppnnnnnppnyyyyhyhhpyE x hhp公式中的導數(shù)用下面公式計算：公式中的導數(shù)用下面公式計算：例子：例子：2(,)(,)(,)( ):(,)2(,)(,)(,)nnnnxnnynnnnnxxnnxynnnyynnynnnyf xyyfxyfxyyyfxyfxyyfxyyfxyy (

22、)cos04(0)1y xxxycos,sin,cosnnnnnnyxyxyx 231cossin/2cos/6nnnnnyyhxhxhx步長步長0.1 步長步長0.01龍格庫塔法：龍格庫塔法：對于常微分方程的數(shù)值解法，一個關鍵在于選擇精度高的算法對于常微分方程的數(shù)值解法，一個關鍵在于選擇精度高的算法計算下面公式中的積分。計算下面公式中的積分。要高精度的計算積分，常用的方法是適當增加計算節(jié)點，不妨要高精度的計算積分，常用的方法是適當增加計算節(jié)點，不妨設用設用 m 個節(jié)點計算上面積分，節(jié)點為個節(jié)點計算上面積分，節(jié)點為則積分為則積分為11()()( , ( )nnxnnxy xy xf x y x

23、 dx*,1,2,inixxhim11( , ( )(, ()nnmxininixif x y x dxhf xh y xh將積分改寫為：將積分改寫為：則得公式：則得公式：取取這樣的公式稱為這樣的公式稱為顯式龍格庫塔顯式龍格庫塔公式，其中系數(shù)選取公式，其中系數(shù)選取為使其具有盡可能高階的截斷誤差，可通過為使其具有盡可能高階的截斷誤差，可通過taylor展開，再比較得到。展開，再比較得到。11( , ( )nnmxiixif x y x dxK11mnniiiyyK111(,)(,),2,3,nniininijjjKhf xyKhf xh yKim確定二確定二階階 RK 法法:將將 在在 點展開，

24、并利用前點展開，并利用前面的公式面的公式 ，比較后得系數(shù)應滿足：，比較后得系數(shù)應滿足：2211(,)nnf xh yK11122122211(,)(,)nnnnnnyyKKKhf xyKhf xh yK12222121110,0,022 21222111,2 (,)nnxy此為四個未知數(shù)的三個方程，任意取此為四個未知數(shù)的三個方程，任意取 ，得，得( )特別，取特別，取 ，得到，得到通常也稱為變形歐拉法，也常寫為通常也稱為變形歐拉法，也常寫為它具有二階精度，也稱為它具有二階精度，也稱為二階二二階二級級 RK 方法方法。1/212121(,)(/2,/2)nnnnnnyyKKhf xyKhf xh

25、yK1(/2,(,)/2)nnnnnnyyhf xhyhf xy例子例子( )cos04(0)1y xxxy三階三級庫塔法三階三級庫塔法局部截斷誤差為局部截斷誤差為4階，整體截斷誤差為階，整體截斷誤差為3階階11231123121(4)6(,)(,)22(,2)nnnnnnnnyyKKKKhf xyKhKhf xyKhf xh yKK最常用的四級四階公式，稱為最常用的四級四階公式，稱為龍格庫塔公式龍格庫塔公式：局部截斷誤差為局部截斷誤差為5階，整體截斷誤差為階，整體截斷誤差為4階。階。1123411223431(22)6(,)(,)22(,)22(,)nnnnnnnnnnyyKKKKKhf x

26、yKhKhf xyKhKhf xyKhf xh yK顯式顯式Lunge-Kutta法的法的絕對穩(wěn)定區(qū)域絕對穩(wěn)定區(qū)域S階的穩(wěn)定區(qū)域為：2|1| 12!sszzzsz 如右圖。穩(wěn)定區(qū)間為：1:( 2,0)2:( 2,0)3:( 2.51,0)4:( 2.78,0)ssss隱式龍格庫塔法：隱式龍格庫塔法：111221122(,)1,2,nnmmininiiimmyyKKKKhf xh yKKKim常用的有二階常用的有二階RK法：法：111111(,)22nnnnyyKKhf xh yK隱式隱式Lunge-Kutta法一般是無條件穩(wěn)定的，穩(wěn)定性法一般是無條件穩(wěn)定的，穩(wěn)定性比較好！比較好！例子：例子：用

27、隱式用隱式二階二階RK法：法：(01)(0)1yyxy 111110( , ),(/2)1/21/21nnnnnnhf x yyKh yKKyhhyyKyyhy 210(1/2)1nnyyhhy用顯式用顯式二階二階RK法：法： DOUBLE PRECISION h,y(0:100),yy(0:100)double precision ak OPEN(20,FILE=OUTPUT5.DAT,STATUS=UNKNOWN)h=1.0/100y(0)=1.0yy(0)=1.0do 10 i=1,100 y(i)=y(i-1)-(h/(1-h/2)*y(i-1) yy(i)=yy(i-1)*(1-h+

28、h*h/2)10 continue do 20 i=0,100 write(20,100) i*h,y(i),yy(i),exp(-i*h)100 format (1x,f4.2,f8.5,f8.5,f8.5)20 continueENDh=0.1h=0.01半隱式龍格庫塔法：半隱式龍格庫塔法：111221122,1111,11(,)(,)1,2,nnmmininiii iiiininii iiyyKKKKhf xh yKKKg K hfxh yKKim1122,11111,11(,)1(,)ininiii iiininii iiKhf xh yKKKhg fxh yKK最常用的二級三階半隱式

29、龍格庫塔公式：最常用的二級三階半隱式龍格庫塔公式：1121112110.413154321.4131543261 (1)(,)(,)661 (1)(,)6(,)0.17378667nnynnnnynnnnyyKKKhhfxyf xyKhhfxh yKf xh yK待定系數(shù)法待定系數(shù)法 線性多步法線性多步法線性 k 步法的一般形式為：100kkn kjnjjnjjjyyhf其中(,)njnjnjff xy，當0k時為顯式法，當0k時為隱式法。通常選取系數(shù),jj使得格式具有盡可能高階的局部截斷誤差，可以用待定系數(shù)法來確定系數(shù)，用廣義Peano定理來給出簡化的局部截斷誤差表示式。微分方程組微分方程組

30、一階方程組一階方程組11122212121010202000( )( ,( ),( ),( )( )( ,( ),( ),( )( )( ,( ),( ),( )(),(),()mmmmmmmy xf x y xyxyxyxfx y xyxyxyxfx y xyxyxy xyyxyyxy00( )( , ( )()y xf x y xy xy 1212010200(,) ,(,) ,(,)mmmyy yyffffyyyy龍格庫塔公式龍格庫塔公式1123411223431(22)6(,)(,)22(,)22(,)nnnnnnnnnnyyKKKKKhf xyKhKhf xyKhKhf xyKhf

31、xh yK寫成分量形式：寫成分量形式：,1,12341121112121221222312411322331(22)6(,)(,)2222(,)2222(,)1,2,i ni niiiiiinnnmnmiinnnmnmiinnnmniinnnmnmyyKKKKKhf xyyyKKKhKhf xyyyKKKhKhf xyyyKhf xh yKyKyKim例子：例子：1122212(01)(0)0,(0)1yyyyyxyy 1,11,1112131411121121121212221312141132231(262)() ()22 ()22 ()nnnnnnnnnnyyKKKKKhyyKKKhyy

32、KKKhyyKhyKyK2,12,212223242122122222232242231(262)()2()2()nnnnnnyyKKKKKhyKKh yKKh yKh yK DOUBLE PRECISION h,y1(0:100),y2(0:100)double precision ak1,ak2,ak3,ak4,bk1,bk2,bk3,bk4 OPEN(20,FILE=OUTPUT6.DAT,STATUS=UNKNOWN)h=1.0/10y1(0)=0.0y2(0)=1.0do 10 i=1,10 ak1=h*(-y1(i-1)+y2(i-1) bk1=-h*y2(i-1) ak2=h*(-(y1(i-1)+ak1/2.0)+y2(i-1)+bk1/2.0) bk2=-h*(y2(i-1)+bk1/2.0) ak3=h*(-(y1(i-1)+ak2/2.0)+y2(i-1)+bk2/2.0) bk3=-h*(y2(i-1)+bk2/2.0) ak4=h*(-(y1(i-1)+ak3)+y2(i-1)+bk3) bk4=-h*(y2(i-1)+bk3) y1(i)=y1(i-1)+(ak1+2*ak2+2*ak3+ak4)/6.0 y2(i)=y2