1、7: ODE1第第7章章 常微分方程（組）的數(shù)常微分方程（組）的數(shù)值解法值解法劉東毅劉東毅天津大學(xué)理學(xué)院數(shù)學(xué)系天津大學(xué)理學(xué)院數(shù)學(xué)系7: ODE2第第7章章 常微分方程（組）的數(shù)值解法常微分方程（組）的數(shù)值解法主要目標(biāo)：主要目標(biāo)： r掌握常微分方程初值問題數(shù)值解法的基本理論r掌握計算機(jī)上的常用算法主要內(nèi)容：主要內(nèi)容：r初值問題計算格式的建立rRunge-Kutta 方法r一階常微分方程組與高階方程的數(shù)值解法7: ODE3第第7章章 常微分方程（組）的數(shù)值解法常微分方程（組）的數(shù)值解法 在科學(xué)研究和工程實(shí)踐中會遇到很多微在科學(xué)研究和工程實(shí)踐中會遇到很多微分方程，雖然從理論上可以證明其解的存在分方程

2、，雖然從理論上可以證明其解的存在性，但其解的解析表達(dá)式往往是很難求解的性，但其解的解析表達(dá)式往往是很難求解的，或者即使可以寫出來，但也難于計算，此，或者即使可以寫出來，但也難于計算，此時，只能借助數(shù)值解來解決問題時，只能借助數(shù)值解來解決問題. 常微分方程（組）定解問題是自然科學(xué)常微分方程（組）定解問題是自然科學(xué)和工程技術(shù)領(lǐng)域中常見的數(shù)學(xué)模型和工程技術(shù)領(lǐng)域中常見的數(shù)學(xué)模型. 本章介紹本章介紹求解此類問題的基本理論和數(shù)值解法。求解此類問題的基本理論和數(shù)值解法。 7: ODE4yR定義定義7.0.1 若存在常數(shù)若存在常數(shù) L 0, 使得對一切的使得對一切的 xa , b 及及 y , , 均有均有

3、|( , )( , )|,f x yf x yL yy則稱則稱 f (x, y) 在在 D 上關(guān)于上關(guān)于 y 滿足滿足 Lipschitz 條件條件, 其中其中 L 稱為稱為 Lipschitz 常數(shù)常數(shù). 我們首先考慮一階常微分方程我們首先考慮一階常微分方程初值問題初值問題0( , ),( ).yf x yaxby ay 其中其中 f (x , y) 是區(qū)域是區(qū)域( , )|,Dx yaxbyR上的實(shí)值函數(shù)上的實(shí)值函數(shù).(7.0.1)7: ODE5 下面下面，我們給出常微分方程初值問題我們給出常微分方程初值問題解的存在惟一性定理。解的存在惟一性定理。定理定理 假設(shè)假設(shè) f (x, y)C(D

4、), ,且關(guān)于且關(guān)于 y 滿足滿足 Lipschitz 條件條件存在唯一解存在唯一解. 接下來，接下來，在此前提下在此前提下, 我們討論的我們討論的數(shù)值數(shù)值解法。解法。0( , ),( ).yf x yaxby ay 7: ODE6012,Naxxxxb然后在節(jié)點(diǎn)上建立逼近于原初值問題的計算格然后在節(jié)點(diǎn)上建立逼近于原初值問題的計算格式式 (或差分格式或差分格式), 由此計算出由此計算出原問題的解原問題的解 y( x ) 在節(jié)點(diǎn)在節(jié)點(diǎn) x1 , x2 , . . . , xN 處的近似值處的近似值: y1 , y2 , . . ., yN，稱它們?yōu)榉Q它們?yōu)槌Ｎ⒎址匠坛踔祮栴}的數(shù)值解常微分方程初值

5、問題的數(shù)值解. 本章，本章，y(xk)表示在節(jié)點(diǎn)表示在節(jié)點(diǎn)xk的真值，的真值，yk 表示相應(yīng)的近似表示相應(yīng)的近似值。值。相鄰兩個節(jié)點(diǎn)的距離相鄰兩個節(jié)點(diǎn)的距離 hn= xn+1 - xn 稱為稱為步長步長, 通常取定步長通常取定步長, 即節(jié)點(diǎn)即節(jié)點(diǎn) xn = x0 + nh , h 0, n = 0,1, , N. 其其基本思想是基本思想是在區(qū)間在區(qū)間 a, b 上引入一系列節(jié)點(diǎn)上引入一系列節(jié)點(diǎn)7: ODE77.1 初值問題計算格式的建立初值問題計算格式的建立1. 數(shù)值微分方法數(shù)值微分方法在等距節(jié)點(diǎn)下討論問題在等距節(jié)點(diǎn)下討論問題. 利用兩點(diǎn)數(shù)值微分公式利用兩點(diǎn)數(shù)值微分公式7.1.1 計算格式的建

6、立計算格式的建立將上式代入初值問題將上式代入初值問題 (7.0.1),)()(0yayx,yfy1()()(, ()() ,2nnnnny xy xhf xy xyh有有11()()()(),()2nnnnnnny xy xhy xyxxh(7.1.1) 00()(, () ,(),0,1,.nnny xf xy xy xy nN7: ODE8略去余項略去余項, 并以數(shù)值解并以數(shù)值解 yn , yn+1 代替代替 y (xn) 及及 y (xn+1), 則得則得差分方程差分方程上式稱為上式稱為 Euler 公式公式. 利用此式可由初值利用此式可由初值 y0 出發(fā)按出發(fā)按“步進(jìn)式步進(jìn)式” 方法方

7、法, 逐步求得數(shù)值逐步求得數(shù)值解解y1 , y2 , . . . , yN .)2 . 1 . 7(),(1nnnnyxhfyy由于計算由于計算yn+1時時, 只只用到它前一步的結(jié)用到它前一步的結(jié)果果yn , 這類公式稱為這類公式稱為單步法單步法. 又因?yàn)槠潢P(guān)又因?yàn)槠潢P(guān)于于yn+1是顯式形式是顯式形式, 故稱該故稱該Euler公式為公式為顯格式顯格式.7: ODE9如果如果利用下列數(shù)值微分公式利用下列數(shù)值微分公式111()()()(),()2nnnnnnny xy xhy xyxxh由由 類似的可導(dǎo)出類似的可導(dǎo)出)()(xx,yfxy),(111nnnnyxhfyy上述公式稱為上述公式稱為后退

8、的后退的 Euler 公式公式, 此公式為此公式為單步法公式單步法公式. 又因?yàn)樗P(guān)于又因?yàn)樗P(guān)于 yn+1 成隱式形式成隱式形式, 所以該公式為所以該公式為隱式公式隱式公式，簡稱，簡稱隱格式隱格式.7: ODE102111()()()(),()26nnnnnnny xy xhy xyxxh類似地，可導(dǎo)出類似地，可導(dǎo)出上述公式稱為上述公式稱為Euler兩步法兩步法公式公式. 這因?yàn)椋?dāng)這因?yàn)?，?dāng)計算計算 yn+1 時時, 要用到要用到 yn -1 與與 yn . 顯然它也是顯顯然它也是顯格式格式.如果如果利用下列三點(diǎn)數(shù)值微分公式利用下列三點(diǎn)數(shù)值微分公式112).(,nnnnyyhf xy(7.

9、1.3) 00()(, () ,(),0,1,.nnny xf xy xy xy nN7: ODE11設(shè)設(shè) y (x)C2a , b, 由由 Taylor 公式有公式有由于由于 故上式即為故上式即為 ( )( , ( ) ,y xf x y x略去余項略去余項, 并以并以 yn , yn+1 代替代替 y (xn) 及及 y (xn+1), 得到的差分方程正是得到的差分方程正是Euler 公式公式.211()()()()().2nnnnnnnhy xy xhy xyxx(7.1.4) )5 . 1 . 7()(! 2)(,()()(21nnnnnyhxyxhfxyxy 2. Taylor 展開

10、展開法法7: ODE123. 數(shù)值積分方法數(shù)值積分方法對對 ，在區(qū)間，在區(qū)間 xn , xn+1 上積分，得上積分，得( )( , )y xf x y11( )( , ( ),nnnnxxxxy x dxf x y x dx則有則有11()()( , ( ).nnxnnxy xy xf x y x dx對上式中的積分采用不同的數(shù)值積分公式可對上式中的積分采用不同的數(shù)值積分公式可得到不同的差分方程得到不同的差分方程. 例如例如, 對上式的積分采對上式的積分采用左矩形公式用左矩形公式, 可得到可得到 Euler 公式公式. 7: ODE1311131()() (, ()(, ()2(, ()().

11、12nnnnnnnnnnnhy xy xf xy xf xy xhfyxx若對此式的積分采用梯形公式若對此式的積分采用梯形公式, 11()()( , ( ).nnxnnxy xy xf x y x dx則有則有若略去余項若略去余項, 以以 yn , yn+1 代替代替 y (xn) 及及 y (xn+1), 得到的差分方程得到的差分方程111 (,)(,2.)nnnnnnhyyf xyf xy7: ODE14上式稱為上式稱為梯形公式梯形公式. 由于它關(guān)于由于它關(guān)于 yn+1 成隱式形成隱式形式式, 故其為故其為隱格式隱格式. 隱格式求解比較困難隱格式求解比較困難, 當(dāng)當(dāng) yn 已知時已知時,

12、要求要求 yn+1 ，需解關(guān)于需解關(guān)于 yn+1 的非線性方程的非線性方程. 在實(shí)際應(yīng)在實(shí)際應(yīng)用時用時, 上式常與上式常與 Euler 公式聯(lián)合使用公式聯(lián)合使用, 構(gòu)成如構(gòu)成如下計算格式下計算格式:111 (,)(,2.)nnnnnnhyyf xyf xy(7.1.6)(0)1(1)( )111(,), (,)(,)(0,1,)2nnnnkknnnnnnyyh f xyhyyf xyf xyk7: ODE15隱式梯形公式的迭代格式隱式梯形公式的迭代格式(0)1(1)( )111(,), (,)(,)(0,1,)2nnnnkknnnnnnyyh f xyhyyf xyf xyk(7.1.7)由由

13、上式可以得到一個序列上式可以得到一個序列: , k = 0,1, , 關(guān)于此序列的收斂性關(guān)于此序列的收斂性, 有如下的定理有如下的定理.( )1kny7: ODE16 設(shè)設(shè) f (x , y) 在區(qū)域在區(qū)域 D 上關(guān)于上關(guān)于 y 滿足滿足 Lipschitz 條件條件, 即即|( , )( , )|()|.f x yf x yLyy其中其中 L 為為 Lipschitz 常數(shù)常數(shù), 當(dāng)步長當(dāng)步長 時時, 對任意的初值對任意的初值 按格式（按格式（7.1.7）2hL(0)1ny(0)1(1)( )111(,), (,)(,)(0,1,)2nnnnkknnnnnnyyh f xyhyyf xyf

14、xyk生成的序列生成的序列 收斂于梯形公式（收斂于梯形公式（7.1.6）( )1kny111 (,)(,)2nnnnnnhyyf xyf xy1ny的解的解 .7: ODE17為了減少計算量為了減少計算量, 可采用可采用預(yù)測預(yù)測- -校正格式校正格式. 方法是先用方法是先用Euler公式求得一個初始近似值公式求得一個初始近似值 稱為稱為預(yù)測值預(yù)測值, 再把再把 帶入梯形公式右端計算一帶入梯形公式右端計算一次求得次求得 yn+1 稱之為稱之為校正值校正值, 即即1ny1ny預(yù)測預(yù)測:校正校正:1(,),nnnnyyhf xy111 (,)(,).2nnnnnnhyyf xyf xy上式稱為上式稱

15、為預(yù)測預(yù)測 - - 校正公式校正公式或或改進(jìn)的改進(jìn)的 EulerEuler 公式公式. 上式也可寫成如下形式上式也可寫成如下形式:11 (,)(,(,)2.nnnnnnnnhyyf xyf xyhf xy7: ODE18例例：利用利用Euler公式與改進(jìn)的公式與改進(jìn)的Euler公式求公式求解初值問題（步長解初值問題（步長h = 0.1）. 1)0(, 10,2yxyxyy解：由步長解：由步長h=0.1，知節(jié)點(diǎn)知節(jié)點(diǎn) 設(shè)數(shù)值解為設(shè)數(shù)值解為 利用利用Euler公公式得式得,1 . 00nnhxxn.10, 2 , 1 , 0n,nynnnnnnnnnnnyxyyxyyyxhfyy2 . 01 .

16、121 . 0),(17: ODE19計算結(jié)果見下表計算結(jié)果見下表(見書見書P227表表7.1) 此初值問題的解析解為此初值問題的解析解為 , 從從上表可以看出上表可以看出, 數(shù)值解數(shù)值解 yn與解析解與解析解 y(xn) 比較比較, yn精度較差精度較差. xy217: ODE20解此問題的改進(jìn)的解此問題的改進(jìn)的Euler公式為公式為 ),(1nnnnyxhfyy),(),(2111nnnnnnyxfyxfhyy.1 . 011. 0105. 12221 . 0,2 . 01 . 121 . 01111111nnnnnnnnnnnnnnnnnnnnnyxyxyyxyyxyyyyxyyxyyy

17、同同Euler公式比較公式比較, 改進(jìn)的改進(jìn)的Euler法顯然精度提高了法顯然精度提高了.由于誤差大小是評價計算格式優(yōu)劣的重要依據(jù)由于誤差大小是評價計算格式優(yōu)劣的重要依據(jù), 故需要給出有關(guān)誤差的概念故需要給出有關(guān)誤差的概念. 計算結(jié)計算結(jié)果見下果見下表表(見書見書P228表表7.2)7: ODE21 截斷誤差與方法的精度截斷誤差與方法的精度定義定義 7.1.1 稱誤差稱誤差 en+1 = y ( xn+1 ) - - yn+1為為數(shù)值方法在點(diǎn)數(shù)值方法在點(diǎn) xn+1 的截斷誤差的截斷誤差, 又稱又稱整體截斷整體截斷誤差誤差. 設(shè)設(shè) yk= y ( xk ) (k = 0,1,. . . , n)

18、，則則 為數(shù)值方法在點(diǎn)為數(shù)值方法在點(diǎn) xn+1 的的局部截斷誤差局部截斷誤差.111()nnny xy7: ODE22整體截斷誤差整體截斷誤差 en+1 是在沒有是在沒有引進(jìn)舍入誤差引進(jìn)舍入誤差的的情況下情況下, 純粹因?yàn)椴粶?zhǔn)確的計算格式造成的純粹因?yàn)椴粶?zhǔn)確的計算格式造成的, 故又稱為故又稱為方法誤差方法誤差.它不僅與它不僅與 x = xn+1 這一步這一步的計算有關(guān)的計算有關(guān), 而且和而且和 xn , xn-1 ,. . . , x1 這幾步的這幾步的計算都有關(guān)系計算都有關(guān)系. 局部截斷誤差是假設(shè)局部截斷誤差是假設(shè) xn 之前各數(shù)值解沒有誤之前各數(shù)值解沒有誤差差, 僅由僅由 xn 到到 xn

19、+1 這一步這一步計算由計算格式引計算由計算格式引起的誤差起的誤差.7: ODE23如如Euler公式公式1(,)nnnnyyhf xy在點(diǎn)在點(diǎn) xn+1 的整體截斷誤差的整體截斷誤差 en+1 = y (xn+1)- - yn+122()(, ()()(,)2(, ()(,)()2nnnnnnnnnnnnnhy xhf xy xyyhf xyheh f xy xf xyy局部截斷誤差局部截斷誤差2111()()2nnnnhy xyy7: ODE24定義定義7.1.2 若某數(shù)值方法的局部截若某數(shù)值方法的局部截斷斷誤差為誤差為 則稱該方法具有則稱該方法具有 P 階精度階精度, 或稱其為或稱其為

20、P 階方法階方法.11(),pnO h可以證明可以證明:Euler 方法的局部截斷誤差方法的局部截斷誤差 其具有其具有一階一階精度精度. 梯形方法的局部截斷誤差梯形方法的局部截斷誤差 其具有其具有二階二階精度精度.改進(jìn)的改進(jìn)的 Euler 方法方法的局部截斷誤的局部截斷誤差差 具有具有二階二階精度精度.21(),nO h31(),nO h31(),nO h7: ODE257.2 Runge-Kutta 方法方法 繼續(xù)討論前面的繼續(xù)討論前面的 Taylor 展開展開法。法。設(shè)設(shè) y (x)C2a , b, 由由 Taylor 公式有公式有211()()()()().2nnnnnnnhy xy x

21、hy xyxx 由由 故上式即為故上式即為 ( )( , ( ),y xf x y x)(! 2)(,()()(21nnnnnyhxyxhfxyxy 略去余項略去余項, 并以并以 yn , yn+1 代替代替 y (xn) 及及 y (xn+1), 得到得到Euler 公式公式.7: ODE26進(jìn)一步假設(shè)設(shè)進(jìn)一步假設(shè)設(shè) y (x)Cp+1a , b,由由 Taylor 公式有公式有) 1 . 2 . 7(,)(!)(! 2)()()()(21nnppnnnnRxyphxyhxyhxyxy )2 . 2 . 7.(),()()!1(11)1(1nnnpnppnxxhOyphR其中其中 由由 故式

22、故式(7.2.1)即為即為 ( )( , ( ),y xf x y x7: ODE277.2 Runge-Kutta 方法略去余項略去余項, 并以數(shù)值解并以數(shù)值解 yn , yn+1 y (xn) 及及 y (xn+1), 可得到一個差分方程，即可得到一個差分方程，即nnppnnnnRxyphxyhxyhxyxy )(!)(! 2)()()()(21) 3 . 2 . 7(,)(,(!)(,(! 2)(,()()()1(21nnnppnnnnnnRxyxfphxyxfhxyxhfxyxy., )()(,()!1(11)(1nnnpnnppnxxhOyfphR其中余項可寫成其中余項可寫成)(,(

23、),()()(xyxfdxdyxfkkk注：這里注：這里7: ODE28略去余項略去余項,用用yn , yn+1 代替代替 y (xn) 及及 y (xn+1)4 . 2 . 7(.! 2! 2)1(1)1(21pnpnnnpnpnnnnfphfhfhyfphfhfhyy, ),(nnnyxff . ) 1, 2 , 1(),()()(piyxffnniin其中其中Taylor的格式的格式 .21(1)()()(,()(,()2!(,(),(7.2.3)!nnnnnnppnnnhy xy xhf xy xfxy xhfxy xRp7: ODE29由于局部截斷誤差由于局部截斷誤差 ),(O)(1

24、111pnnnhyxy可知它是一個可知它是一個 p 階方法。當(dāng)階方法。當(dāng)p=1時，上式正是時，上式正是Euler 公式。但當(dāng)公式。但當(dāng) p 2 時，需要計算時，需要計算f (x, y(x) ) 的高階導(dǎo)數(shù)，特別是對于復(fù)雜函數(shù)的高階導(dǎo)數(shù)，特別是對于復(fù)雜函數(shù) f (x, y(x) 的的求導(dǎo)，這無疑是大大增加計算量，這是它最大求導(dǎo)，這無疑是大大增加計算量，這是它最大的缺點(diǎn)。因此高階的的缺點(diǎn)。因此高階的Taylor方法是不實(shí)用的。方法是不實(shí)用的。 德國數(shù)學(xué)家德國數(shù)學(xué)家C.Runge 及提出了一種改進(jìn)策及提出了一種改進(jìn)策略，得到了至今還被作為高精度的單步法廣泛略，得到了至今還被作為高精度的單步法廣泛使用

25、龍格使用龍格-庫塔法庫塔法 ( Runge- Kutta method )。 7: ODE30方法的基本思想方法的基本思想.! 2)1(11pnpnnnnfphfhfhyyRunge-Kutta方法是方法是利用利用 f 在在h 后面括號中后面括號中的因子來構(gòu)造差分方程的因子來構(gòu)造差分方程, 從而避免了高階導(dǎo)數(shù)的從而避免了高階導(dǎo)數(shù)的計算計算, 這就是這就是 Runge-Kutta 方法的基本思想方法的基本思想. )1(1! 2pnpnnfphfhf用用f 在在某些點(diǎn)處函數(shù)值的線性組合替代這一部分某些點(diǎn)處函數(shù)值的線性組合替代這一部分7: ODE31其一般形式為其一般形式為:11111,(,),(,

26、)(2,3, )rnniiinniininijjjyyhkkf xykf xh yhkir其中其中 r 是上式中調(diào)用是上式中調(diào)用 f 的個數(shù)的個數(shù), r 稱為級數(shù)，稱為級數(shù)， 為待定參數(shù)為待定參數(shù), 適當(dāng)確定這些參數(shù)適當(dāng)確定這些參數(shù), 可使可使之具有盡可能高的精度之具有盡可能高的精度. 如局部截斷誤差滿足如局部截斷誤差滿足,iii j . )(O)(1111rnnnhyxy7: ODE32 二階二階 Runge-Kutta 方法方法考慮考慮 r = 2 的情況的情況, 此時有此時有11 122121(),(,),(,).nnnnnnyyhkkkf xykf xh yh k利用二元函數(shù)的利用二元

27、函數(shù)的 一階一階Taylor 公式，即全微分公式公式，即全微分公式 希望適當(dāng)選擇參數(shù)希望適當(dāng)選擇參數(shù) 使上式的局部截斷誤差為使上式的局部截斷誤差為12, , 3111()()nnny xyO h即為二階方法即為二階方法.,)()()(,()(,(),(),(22nnnnnynnnxnnyyxxOyyyxfxxyxfyxfyxf下面將下面將yn+1與與y(xn+1)作比較作比較7: ODE332222( ,)(,)(,)()(,)()()()() ()() ()()()nnxnnnynnnnndefnxnnynnnnf x yf xyfxyxxfxyyyOxxyyffxxfyyOxxyy. )O

28、()()()(212hkfhfhfknynxn從而有, )(22111kkhyynn將上式代入),(12khyhxfknn再由7: ODE34得到. )()()()()O()()()()(322212121122111hOfffhfhyhkfhfhfhkhykkhyynynnxnnnynxnnnn在下面要將yn+1與y(xn+1)作比較，使它們的局部截斷誤差滿足. )O()()()(212hkfhfhfknynxn, )(O)(3111hyxynnn為此考慮y(xn+1)。7: ODE35再根據(jù) y(xn+1) 在點(diǎn) xn的一元 3 階 Taylor 展開式. )()()(! 2)(),(!

29、2)()(! 2)()()(323),(2321hOfffhhfyhOyxfhhfyhOxyhxyhxyxynynnxnnyxnnnnnnnn )()()()(322211hOfffhfhyynynnxnnn由剛才已得到的, )(O)(3111hyxynnn讓它們滿足7: ODE36即由122211212左式含有四個未知元三個左式含有四個未知元三個方程方程, 因此因此解不唯一解不唯一. 參數(shù)參數(shù)滿足左式的滿足左式的一族公式一族公式統(tǒng)稱統(tǒng)稱二階二階 Runge-Kutta 公式公式. 可得參數(shù)應(yīng)滿足下列方程組:. )()()(! 2)(321hOfffhhfyxynynnxnnn)()()()(

30、322211hOfffhfhyynynnxnnn7: ODE37 取上式稱為上式稱為中點(diǎn)公式中點(diǎn)公式 .1210,1,212,nnyyhk1(,),nnkf xy21(,).22nnhhkf xyk 取上式稱為上式稱為 Heun 公式公式 .112(3),4nnhyykk1(,),nnkf xy2122(,).33nnkf xh yhk12132,443 可見可見, 二階二階 Runge-Kutta 公式公式, 每計算一步需要每計算一步需要 兩次調(diào)用兩次調(diào)用 f 的函數(shù)值的函數(shù)值. 取121,1,2112121(),2(,),(,).nnnnnnhyykkkf xykf xh yhk得這正是這

31、正是改進(jìn)的改進(jìn)的 Euler 公式公式. 7: ODE38 四階四階Runge-Kutta方法方法當(dāng)當(dāng) r = 4 時時, 類似地可導(dǎo)出四階類似地可導(dǎo)出四階 Runge-Kutta 公公式式, 這種公式也有一族這種公式也有一族, 其中常用地有其中常用地有:標(biāo)準(zhǔn)標(biāo)準(zhǔn) (經(jīng)典經(jīng)典) 的的 Runge-Kutta 方法方法11234(22),6nnhyykkkk1(,),nnkf xy21(,),22nnhhkf xyk32(,),22nnhhkf xyk43(,).nnkf xh yk7: ODE3911234(22)(22),6nnhyykkkk1(,),nnkf xy21(,),22nnhhk

32、f xyk3122122(,),222nnhkf xyhkhk423222(,).22nnkf xh yhkhk Gill 公式公式 Gill 公式是標(biāo)準(zhǔn)的公式是標(biāo)準(zhǔn)的 Runge-Kutta 公式的改進(jìn)形公式的改進(jìn)形式式, 這種算法可節(jié)省存儲單元這種算法可節(jié)省存儲單元, 并能控制舍入誤并能控制舍入誤差的增長差的增長.四階四階 Runge-Kutta 公式公式, 每一步計算需四次調(diào)用每一步計算需四次調(diào)用 f 的函數(shù)值的函數(shù)值, 計算量較大計算量較大, 但其局部截斷誤差可達(dá)但其局部截斷誤差可達(dá) O(h5), 精度較高精度較高.7: ODE40例例7.2.1 用標(biāo)準(zhǔn)的用標(biāo)準(zhǔn)的 四階四階Rung-K

33、utta 法解初值問題法解初值問題,取取步長步長h=0.2.2, (0)1,(01).xyyyxy 解解: 解此問題的計算公式為解此問題的計算公式為112340.2(22),6nnyykkkk12,nnnxkyy21122,22()nnnhxhkykhyk32222,22()nnnhxhkykhyk4332,()nnnxhkyhkyhkxn0.20.40.60.81.0ynyn- y (xn)1.183 21.341 71.483 3 1.612 5 1.732 10.000 00.000 00.000 0 0.000 10.000 1計算結(jié)果如下計算結(jié)果如下:顯然在計算量大致相同的情顯然在計

34、算量大致相同的情況下況下, 標(biāo)準(zhǔn)的標(biāo)準(zhǔn)的 Runge-Kutta方方法比改進(jìn)的法比改進(jìn)的 Euler 方法精確方法精確度更高度更高. （參見（參見p227和和p228的表）的表）7: ODE417.5一階常微分方程組與高階方程一階常微分方程組與高階方程初值問題的數(shù)值解法初值問題的數(shù)值解法 一階常微分方程組一階常微分方程組初值問題初值問題 )y,y,y(x,fdxdy)y,y,y(x,fdxdy)y,y,y(x,fdxdym21mmm2122m21111122( ),( ),( ), , mmy a = s y a = s ya = s xa b7: ODE42寫成向量形式：寫成向量形式： 0(

35、 ,)( ) (7.5.2)dYF x YdxY a =Y,)()()(Yn21n21xyxyxyyyy.s ,s ,sYTn21011221212( ,)( ,)(),( ,)mmmmf x yyyfx yyyF x,Yfx yyy其中其中7: ODE43利用向量值函數(shù)的微積分理論利用向量值函數(shù)的微積分理論, 很容易推導(dǎo)很容易推導(dǎo)出出一階常微分方程組一階常微分方程組初值問題的數(shù)值解法初值問題的數(shù)值解法. )3 . 5 . 7(),(1nnnnYxhFYY如如Euler公式公式,),(T,2, 1nmnnnyyyY其中其中 .),(),(),(),(, 2, 1, 2, 12, 2, 11nm

36、nnnmnmnnnnmnnnnnyyyxfyyyxfyyyxfYxF7: ODE44的分量形式為的分量形式為 )4 . 5 . 7(. ), 2, 1(),(, 2, 1,1,miyyyxhfyynmnnninini),(),(),(, 2, 1, 2, 12, 2, 11, 2, 11,1, 21, 1nmnnnmnmnnnnmnnnnmnnnmnnyyyxfyyyxfyyyxfhyyyyyy或或 7: ODE45四階標(biāo)準(zhǔn)的四階標(biāo)準(zhǔn)的Runge-Kutta公式公式 )5 . 5 . 7()22(643211KKKKhYYnn),(1nnYxFK )2,2(12KhYhxFKnn)2,2(23

37、KhYhxFKnn),(34hKYhxFKnn設(shè)設(shè) 則則 (7.5.5)的分量形式為的分量形式為 ,)4, 3, 2, 1(),(T,2, 1jkkkKjmjjj7: ODE46四階標(biāo)準(zhǔn)的四階標(biāo)準(zhǔn)的Runge-Kutta公式的分量形式公式的分量形式 ,1,1,2,3,4,11,2,21,12,1,1,31,22,2,2,41,32,(22),6(,),(,),2222(,),2222(,i ni niiiiiinnnm niinninim niiinninim niiinnihyykkkkkf xyyyhhhhkfxykykykhhhhkfxykykykkfxh yhky,3,3,),nim

38、nihkyhk., 2, 1mi其中其中7: ODE47例：試寫出用中點(diǎn)公式解下列初值問題的例：試寫出用中點(diǎn)公式解下列初值問題的計算公式計算公式. 3)0(, 1)0(,35,643zyzyzzyxy解解:令令 ,zyY則則 .nnnzyY.35643),(zyzyxYxF,)2 , 1(, 2, 1jkkKjjj再取再取由向量形式的中點(diǎn)公式由向量形式的中點(diǎn)公式7: ODE48中點(diǎn)公式的向量形式中點(diǎn)公式的向量形式 )2,2(),(12121KhYhxFKYxFKhKYYnnnnnn上述中點(diǎn)公式的分量計算形式為上述中點(diǎn)公式的分量計算形式為2, 22, 111kkhzyzynnnnnnnnnzyzyxkk356431 , 21 , 1)2( 3)2( 5)2(6)2(4)2( 31 , 21 , 11 , 21 , 12, 22, 1khzkhykhzkhyhxkknnnnn分量計算形式為分量計算形式為7: ODE49整理得分量計算格式整理得分量計算格式. 3, 1, )2( 3)2(5, )2(6)2(4)2( 3,35,643,001 , 21 , 12, 21 , 21 , 12, 11 , 21 , 12, 212, 11zyhkzhkykhkzhkyhxkzykzyxkkhzzkhyynnnn