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

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

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

3、步長為步長為0.1的計(jì)算結(jié)果。的計(jì)算結(jié)果。步長為步長為0.01的計(jì)算結(jié)果的計(jì)算結(jié)果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歐拉中點(diǎn)公式歐拉中點(diǎn)公式 歐拉公式中我們利用了近似公式歐拉公式中我們利用了近似公式1()()()nnny xy

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

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

7、關(guān)鍵在于節(jié)點(diǎn)的函數(shù)值，關(guān)鍵在于選擇適當(dāng)?shù)倪x擇適當(dāng)?shù)姆e分公式計(jì)算積分積分公式計(jì)算積分！( )( , ( )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、這是所謂的這是所謂的歐拉梯形公式歐拉梯形公式。直接利用已經(jīng)求得的已知節(jié)點(diǎn)上的值計(jì)算未知節(jié)點(diǎn)上的函數(shù)直接利用已經(jīng)求得的已知節(jié)點(diǎn)上的值計(jì)算未知節(jié)點(diǎn)上的函數(shù)值的算法稱為值的算法稱為顯式法顯式法。例如：歐拉公式、歐拉中點(diǎn)公式例如：歐拉公式、歐拉中點(diǎn)公式計(jì)算未知節(jié)點(diǎn)上的函數(shù)值時(shí)，用到了未知節(jié)點(diǎn)上的函數(shù)值，計(jì)算未知節(jié)點(diǎn)上的函數(shù)值時(shí)，用到了未知節(jié)點(diǎn)上的函數(shù)值，這種算法稱為這種算法稱為隱式法。隱式法。例如：后退歐拉法、歐拉梯形公式例如：后退歐拉法、歐拉梯形公式顯然，利用隱式法求微分方程的數(shù)值解是，顯然，利用隱式法求微分方程的數(shù)值解是，需要從表達(dá)式中需要從表達(dá)式中反解未知節(jié)點(diǎn)上的函數(shù)值反解未知節(jié)點(diǎn)上的函數(shù)值。1

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

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

11、h1121211()2(,)(,)nnnnnnyyKKKhf xyKhf xh yK1nnxxh例子：例子：取步長為取步長為0.1計(jì)算，結(jié)果如圖。計(jì)算，結(jié)果如圖。(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同理，對于后退歐拉公式同理，對于后退歐拉公式有預(yù)測校正公式有預(yù)測校正公式或改寫為：或改寫為：111(,)nnnnyyf xyh(0)1(0)111(,)(,)nnnnnnnnyyhf xyyyf xyh11(,)(,)nnnnnnKhf xyyyf xyK h用此法解前面的例子用此法解前面的例子步長步長0.1步長步長0.011.4 誤差估計(jì)誤差估計(jì)定義：定義：利用第利用第n個(gè)節(jié)點(diǎn)或之前更多節(jié)點(diǎn)的函數(shù)精確值，利用近個(gè)節(jié)點(diǎn)或之前更多

13、節(jié)點(diǎn)的函數(shù)精確值，利用近似公式數(shù)值計(jì)算第似公式數(shù)值計(jì)算第n+1個(gè)節(jié)點(diǎn)的近似值，所引起的誤差，稱個(gè)節(jié)點(diǎn)的近似值，所引起的誤差，稱為第為第n+1個(gè)節(jié)點(diǎn)上的個(gè)節(jié)點(diǎn)上的局部截?cái)嗾`差局部截?cái)嗾`差。我們記我們記 為第為第n+1個(gè)節(jié)點(diǎn)上解的精確值，個(gè)節(jié)點(diǎn)上解的精確值， 為假設(shè)為假設(shè) 等條件下計(jì)算所得的近似值，等條件下計(jì)算所得的近似值，則局部截?cái)嗾`差為：則局部截?cái)嗾`差為：如局部截?cái)嗾`差為如局部截?cái)嗾`差為 ，稱為具有，稱為具有 p 階局部截?cái)嗾`差。階局部截?cái)嗾`差。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而而完全類似的可以得到完全類似的可以得到后退歐拉公式的局部截?cái)嗾`差為：后退歐拉公式的局部截?cái)嗾`差為：歐拉中點(diǎn)公式的局部截?cái)嗾`差為：歐拉中點(diǎn)公式的局部截?cái)嗾`差為：歐拉梯形公式的局部截?cái)嗾`差為：歐拉梯形公式的局部截?cái)嗾`差為：2111()()nnny xyO h3111()()nnny xyO h3111()()n

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

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

17、b fb fb fA1(1) ( )kkkiR yB hy類似地，如果取類似地，如果取 作為作為節(jié)點(diǎn)，可得節(jié)點(diǎn)，可得亞當(dāng)斯亞當(dāng)斯(Adams)隱式公式隱式公式局部截?cái)嗾`差：局部截?cái)嗾`差：11,iiii kxx xx*10112()iiiiiki khyyb fb fb fb fA*2(2)1 ( )kkkiR yBhy進(jìn)一步，如果我們將恒等式中的積分區(qū)間改為進(jìn)一步，如果我們將恒等式中的積分區(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)定性一個(gè)常微分方程數(shù)值解法應(yīng)用于方程一個(gè)常微分方程數(shù)值解法應(yīng)用于方程 時(shí)對給定步長時(shí)對給定步長 h 稱為稱為絕對穩(wěn)定絕對穩(wěn)定，是指在某一步（如，是指在某一步（如第一步）產(chǎn)生的誤差（如計(jì)算機(jī)的存儲(chǔ)誤差），在第一步）產(chǎn)生的誤差（如計(jì)算機(jī)的存儲(chǔ)誤差），在計(jì)算中會(huì)逐步減小。計(jì)算中會(huì)逐步減小。稱方程稱方程 為為試驗(yàn)方程試驗(yàn)方程，設(shè)，設(shè)在計(jì)算開始時(shí)產(chǎn)生誤在計(jì)算開始時(shí)產(chǎn)生誤差（存儲(chǔ)誤差），此誤差在以后差（存儲(chǔ)誤差），此誤差在以后會(huì)逐步減小會(huì)逐步減小，我們，我們稱該

19、算法相對于稱該算法相對于 是絕對穩(wěn)定的，這樣的是絕對穩(wěn)定的，這樣的 的全體稱為該算法的的全體稱為該算法的絕對穩(wěn)定域絕對穩(wěn)定域。 通常在復(fù)平面上通常在復(fù)平面上考慮絕對穩(wěn)定域。穩(wěn)定域與實(shí)軸的交集為考慮絕對穩(wěn)定域。穩(wěn)定域與實(shí)軸的交集為穩(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 局部截?cái)嗾`差的實(shí)用估計(jì)局部截?cái)嗾`差的實(shí)用估計(jì)（1）用兩種階數(shù)相同的算法求解，計(jì)算出）用兩種階數(shù)相同的算法求解，計(jì)

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

21、函數(shù)如函數(shù) y( x )有有p+1階導(dǎo)數(shù)，容易得到階導(dǎo)數(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公式中的導(dǎo)數(shù)用下面公式計(jì)算：公式中的導(dǎo)數(shù)用下面公式計(jì)算：例子：例子：2(,)(,)(,)( ):(,)2(,)(,)(,)nnnnxnnynnnnnxxnnxynnnyynnynnnyf xyyfxyfxyyyfxyfxyyfxyyfxyy (

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

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

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

25、yK1(/2,(,)/2)nnnnnnyyhf xhyhf xy例子例子( )cos04(0)1y xxxy三階三級庫塔法三階三級庫塔法局部截?cái)嗾`差為局部截?cái)嗾`差為4階，整體截?cái)嗾`差為階，整體截?cái)嗾`差為3階階11231123121(4)6(,)(,)22(,2)nnnnnnnnyyKKKKhf xyKhKhf xyKhf xh yKK最常用的四級四階公式，稱為最常用的四級四階公式，稱為龍格庫塔公式龍格庫塔公式：局部截?cái)嗾`差為局部截?cái)嗾`差為5階，整體截?cái)嗾`差為階，整體截?cái)嗾`差為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ū)域?yàn)椋?|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，當(dāng)0k時(shí)為顯式法，當(dāng)0k時(shí)為隱式法。通常選取系數(shù),jj使得格式具有盡可能高階的局部截?cái)嗾`差，可以用待定系數(shù)法來確定系數(shù)，用廣義Peano定理來給出簡化的局部截?cái)嗾`差表示式。微分方程組微分方程組

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