chap常微分方程的數(shù)值解法PPT課件

1、1第第7章章 常微分方程（組）的常微分方程（組）的數(shù)值解法數(shù)值解法劉東毅劉東毅天津大學(xué)理學(xué)院數(shù)學(xué)系天津大學(xué)理學(xué)院數(shù)學(xué)系第1頁(yè)/共56頁(yè)2第7章 常微分方程（組）的數(shù)值解法主要目標(biāo)： 掌握常微分方程初值問(wèn)題數(shù)值解法的基本理論掌握計(jì)算機(jī)上的常用算法主要內(nèi)容：初值問(wèn)題計(jì)算格式的建立Runge-Kutta 方法一階常微分方程組與高階方程的數(shù)值解法第2頁(yè)/共56頁(yè)3第7章 常微分方程（組）的數(shù)值解法 在科學(xué)研究和工程實(shí)踐中會(huì)遇到很多微分方程，雖然從理論上可以證明其解的存在性，但其解的解析表達(dá)式往往是很難求解的，或者即使可以寫(xiě)出來(lái)，但也難于計(jì)算，此時(shí)，只能借助數(shù)值解來(lái)解決問(wèn)題. 常微分方程（組）定解問(wèn)題是

2、自然科學(xué)和工程技術(shù)領(lǐng)域中常見(jiàn)的數(shù)學(xué)模型. 本章介紹求解此類問(wèn)題的基本理論和數(shù)值解法。 第3頁(yè)/共56頁(yè)4yR定義定義7.0.1 若存在常數(shù)若存在常數(shù) L 0, 使得對(duì)一切的使得對(duì)一切的 xa , b 及及 y , , 均有均有 |( , )( , )|,f x yf x yL yy則稱則稱 f (x, y) 在在 D 上關(guān)于上關(guān)于 y 滿足滿足 Lipschitz 條件條件, 其中其中 L 稱為稱為 Lipschitz 常數(shù)常數(shù). 我們首先考慮一階常微分方程我們首先考慮一階常微分方程初值問(wèn)題初值問(wèn)題0( , ),( ).yf x yaxby ay 其中其中 f (x , y) 是區(qū)域是區(qū)域(

3、, )|,Dx yaxbyR上的實(shí)值函數(shù)上的實(shí)值函數(shù).(7.0.1)第4頁(yè)/共56頁(yè)5 我們首先給出常微分方程初值問(wèn)題解的存在惟一性定理。定理 假設(shè) f (x, y)C(D), ,且關(guān)于 y 滿足 Lipschitz 條件, 則一階常微分方程初值問(wèn)題(7.0.1) 存在唯一解. 下面在此前提下, 我們討論上述初值問(wèn)題 (7.0.1) 的數(shù)值解法。0( , ),( ).yf x yaxby ay 第5頁(yè)/共56頁(yè)6012,Naxxxxb然后在節(jié)點(diǎn)上建立逼近于原初值問(wèn)題的計(jì)算格式然后在節(jié)點(diǎn)上建立逼近于原初值問(wèn)題的計(jì)算格式 (或差分格式或差分格式), 由此計(jì)算出由此計(jì)算出原問(wèn)題原問(wèn)題的解的解 y(

4、x ) 在節(jié)點(diǎn)在節(jié)點(diǎn) x1 , x2 , . . . , xN 處的近似值：處的近似值： y1 , y2 , . . ., yN， 稱它們?yōu)榉Q它們?yōu)槌Ｎ⒎址匠坛醭Ｎ⒎址匠坛踔祮?wèn)題的數(shù)值解值問(wèn)題的數(shù)值解. 相鄰兩個(gè)節(jié)點(diǎn)的距離相鄰兩個(gè)節(jié)點(diǎn)的距離 hn= xn+1 - xn 稱為稱為步長(zhǎng)步長(zhǎng), 通常取定步長(zhǎng)通常取定步長(zhǎng)h 0, 即即節(jié)點(diǎn)節(jié)點(diǎn) xn = x0 + nh , n = 0,1, , N. 其基本思想是在區(qū)間 a , b 上引入一系列節(jié)點(diǎn)第6頁(yè)/共56頁(yè)7 初值問(wèn)題計(jì)算格式的建立1. 數(shù)值微分方法數(shù)值微分方法在等距節(jié)點(diǎn)下討論問(wèn)題在等距節(jié)點(diǎn)下討論問(wèn)題. 利用兩點(diǎn)數(shù)值微分公式利用兩點(diǎn)數(shù)值微分公式

5、.1 計(jì)算格式的建立計(jì)算格式的建立將上式代入初值問(wèn)題將上式代入初值問(wèn)題 (7.0.1),)()(0yayx,yfy1()()(, ()(),2nnnnny xy xhf xy xyh有有11()()()(),()2nnnnnnny xy xhy xyxxh(7.1.1) 第7頁(yè)/共56頁(yè)8略去余項(xiàng), 并以數(shù)值解 yn , yn+1 代替 y (xn) 及 y (xn+1), 則得差分方程差分方程上式稱為上式稱為 Euler 公式公式. 利用此式可由利用此式可由初值初值 y0 出發(fā)按出發(fā)按“步進(jìn)式步進(jìn)式” 方法方法, 逐步逐步求得數(shù)值解求得數(shù)值解y1 , y2 , . . . , yN .)2

6、. 1 . 7(),(1nnnnyxhfyy由于計(jì)算由于計(jì)算yn+1時(shí)時(shí), 只只用到它前一步的結(jié)用到它前一步的結(jié)果果yn , 這類公式稱為這類公式稱為單步法單步法. 又因?yàn)槠潢P(guān)又因?yàn)槠潢P(guān)于于yn+1是顯式形式是顯式形式, 故稱該故稱該Euler公式為公式為顯格式顯格式.第8頁(yè)/共56頁(yè)9如果利用下列數(shù)值微分公式111()()()(),()2nnnnnnny xy xhy xyxxh由由 類似的可導(dǎo)出類似的可導(dǎo)出)()(xx,yfxy),(111nnnnyxhfyy上述公式稱為上述公式稱為后退的后退的 Euler 公式公式, 此公式為單步法公式此公式為單步法公式. 又因?yàn)樗P(guān)于又因?yàn)樗P(guān)于 yn

7、+1 成成隱式形式隱式形式, 所以該公式為所以該公式為隱式公式隱式公式，簡(jiǎn)稱，簡(jiǎn)稱隱格式隱格式.第9頁(yè)/共56頁(yè)102111()()()(),()26nnnnnnny xy xhy xyxxh類似地，可導(dǎo)出類似地，可導(dǎo)出上述公式稱為上述公式稱為Euler兩步法兩步法公式公式. 這因?yàn)?，?dāng)計(jì)算這因?yàn)椋?dāng)計(jì)算 yn+1 時(shí)時(shí), 要用到要用到 yn -1 與與 yn . 顯然它顯然它也是顯格式也是顯格式.如果利用下列三點(diǎn)數(shù)值微分公式112).(,nnnnyyhf xy(7.1.3) 第10頁(yè)/共56頁(yè)11設(shè)設(shè) y (x)C2a , b, 由由 Taylor 公式有公式有由于由于 故上式即為故上式即為

8、 ( )( , ( ) ,y xf x y x略去余項(xiàng)略去余項(xiàng), 并以并以 yn , yn+1 代替代替 y (xn) 及及 y (xn+1), 得到的差分方程正是得到的差分方程正是Euler 公式公式.211()()()()().2nnnnnnnhy xy xhy xyxx(7.1.4) )5 . 1 . 7()(! 2)(,()()(21nnnnnyhxyxhfxyxy 2. Taylor 展開(kāi)展開(kāi)法法第11頁(yè)/共56頁(yè)123. 數(shù)值積分方法數(shù)值積分方法對(duì)對(duì) ，在區(qū)間，在區(qū)間 xn , xn+1 上積分，得上積分，得( )( , )y xf x y11( )( , ( ),nnnnxxxx

9、y x dxf x y x dx則有則有11()()( , ( ).nnxnnxy xy xf x y x dx對(duì)上式中的積分采用不同的數(shù)值積分公式可得到不同的差分方程對(duì)上式中的積分采用不同的數(shù)值積分公式可得到不同的差分方程. 例如例如, 對(duì)上式對(duì)上式的積分采用左矩形公式的積分采用左矩形公式, 可得到可得到 Euler 公式公式. 第12頁(yè)/共56頁(yè)1311131()() (, ()(, ()2(, ()().12nnnnnnnnnnnhy xy xf xy xf xy xhfyxx若對(duì)此式的積分采用梯形公式, 11()()( , ( ).nnxnnxy xy xf x y x dx則有則有若

10、略去余項(xiàng)若略去余項(xiàng), 以以 yn , yn+1 代替代替 y (xn) 及及 y (xn+1), 得到的差分方程得到的差分方程111 (,)(,2.)nnnnnnhyyf xyf xy第13頁(yè)/共56頁(yè)14上式稱為上式稱為梯形公式梯形公式. 由于它關(guān)于由于它關(guān)于 yn+1 成隱式形式成隱式形式, 故其為故其為隱格式隱格式. 隱格式求解比較困難隱格式求解比較困難, 當(dāng)當(dāng) yn 已知時(shí)已知時(shí), 要求要求 yn+1 ，需解關(guān)于需解關(guān)于 yn+1 的非線性方程的非線性方程. 在實(shí)際應(yīng)用時(shí)在實(shí)際應(yīng)用時(shí), 上式常與上式常與 Euler 公式聯(lián)合使用公式聯(lián)合使用, 構(gòu)成如下計(jì)算格式構(gòu)成如下計(jì)算格式:111

11、(,)(,2.)nnnnnnhyyf xyf xy(7.1.6)(0)1(1)( )111(,), (,)(,)(0,1,)2nnnnkknnnnnnyyh f xyhyyf xyf xyk第14頁(yè)/共56頁(yè)15隱式梯形公式的迭代格式(0)1(1)( )111(,), (,)(,)(0,1,)2nnnnkknnnnnnyyh f xyhyyf xyf xyk(7.1.7)由由上式可以得到一個(gè)序列上式可以得到一個(gè)序列: , k = 0,1, , 關(guān)于此序列的收斂性關(guān)于此序列的收斂性, 有如下的定理有如下的定理.( )1kny第15頁(yè)/共56頁(yè)16 設(shè) f (x , y) 在區(qū)域 D 上關(guān)于 y

12、滿足 Lipschitz 條件, 即|( , )( , )|()|.f x yf x yLyy其中其中 L 為為 Lipschitz 常數(shù)常數(shù), 當(dāng)步長(zhǎng)當(dāng)步長(zhǎng) 時(shí)時(shí), 對(duì)任意的初值對(duì)任意的初值 按格式（按格式（7.1.7）2hL(0)1ny(0)1(1)( )111(,), (,)(,)(0,1,)2nnnnkknnnnnnyyh f xyhyyf xyf xyk生成的序列生成的序列 收斂于梯形公式（收斂于梯形公式（7.1.6）( )1kny111 (,)(,)2nnnnnnhyyf xyf xy1ny的解的解 .第16頁(yè)/共56頁(yè)17為了減少計(jì)算量為了減少計(jì)算量, 可采用可采用預(yù)測(cè)預(yù)測(cè)- -

13、校正格式校正格式. 方法是先用方法是先用Euler公式求得一個(gè)初公式求得一個(gè)初始近似值始近似值 稱為稱為預(yù)測(cè)值預(yù)測(cè)值, 再把再把 帶入梯形公式右端計(jì)算一次求得帶入梯形公式右端計(jì)算一次求得 yn+1 稱之為稱之為校正校正值值, 即即1ny1ny預(yù)測(cè)預(yù)測(cè):校正校正:1(,),nnnnyyhf xy111 (,)(,).2nnnnnnhyyf xyf xy上式稱為上式稱為預(yù)測(cè)預(yù)測(cè) - - 校正公式校正公式或或改進(jìn)的改進(jìn)的 EulerEuler 公式公式. 上式也可寫(xiě)成如下形式上式也可寫(xiě)成如下形式:11 (,)(,(,)2.nnnnnnnnhyyf xyf xyhf xy第17頁(yè)/共56頁(yè)18例：利用

14、Euler公式與改進(jìn)的Euler公式求解初值問(wèn)題（步長(zhǎng)）. 1)0(, 10,2yxyxyy解：由步長(zhǎng)解：由步長(zhǎng)h=0.1，知節(jié)點(diǎn)知節(jié)點(diǎn) 設(shè)數(shù)值解為設(shè)數(shù)值解為 利用利用Euler公式得公式得,1 . 00nnhxxn.10, 2 , 1 , 0n,nynnnnnnnnnnnyxyyxyyyxhfyy2 . 01 . 121 . 0),(1第18頁(yè)/共56頁(yè)19計(jì)算結(jié)果見(jiàn)下表(見(jiàn)書(shū)P227表7.1) 此初值問(wèn)題的解析解為此初值問(wèn)題的解析解為 , 從上表可以看出從上表可以看出, 數(shù)值解數(shù)值解 yn與解析與解析解解 y(xn) 比較比較, yn精度較差精度較差. xy21第19頁(yè)/共56頁(yè)20解此問(wèn)

15、題的改進(jìn)的Euler公式為 ),(1nnnnyxhfyy),(),(2111nnnnnnyxfyxfhyy.1 . 011. 0105. 12221 . 0,2 . 01 . 121 . 01111111nnnnnnnnnnnnnnnnnnnnnyxyxyyxyyxyyyyxyyxyyy同同Euler公式比較公式比較, 改進(jìn)的改進(jìn)的Euler法顯然精度提高了法顯然精度提高了.由于誤差大小是評(píng)價(jià)計(jì)算格式優(yōu)劣的重要依據(jù)由于誤差大小是評(píng)價(jià)計(jì)算格式優(yōu)劣的重要依據(jù), 故需要給出有關(guān)誤差的概念故需要給出有關(guān)誤差的概念. 計(jì)算結(jié)計(jì)算結(jié)果見(jiàn)下果見(jiàn)下表表(見(jiàn)書(shū)見(jiàn)書(shū)P228表表7.2)第20頁(yè)/共56頁(yè)21 截?cái)?/p>

16、誤差與方法的精度定義定義 7 稱誤差稱誤差 en+1 = y ( xn+1 ) - - yn+1為為數(shù)值方法在點(diǎn)數(shù)值方法在點(diǎn) xn+1 的截?cái)嗾`差的截?cái)嗾`差, 又稱又稱整體截?cái)嗾`差整體截?cái)嗾`差. 設(shè)設(shè) yk= y ( xk ) (k = 0,1,. . . , n)，則則 為數(shù)值方法在點(diǎn)為數(shù)值方法在點(diǎn) xn+1 的的局部截?cái)嗾`差局部截?cái)嗾`差.111()nnny xy第21頁(yè)/共56頁(yè)22整體截?cái)嗾`差 en+1 是在沒(méi)有引進(jìn)舍入誤差的情況下, 純粹因?yàn)椴粶?zhǔn)確的計(jì)算格式造成的, 故又稱為方法誤差.它不僅與 x = xn+1 這一步的計(jì)算有關(guān), 而且和 xn , xn-1 ,. . . , x1 這

17、幾步的計(jì)算都有關(guān)系. 局部截?cái)嗾`差是假設(shè)局部截?cái)嗾`差是假設(shè) xn 之前各數(shù)值解沒(méi)有誤差之前各數(shù)值解沒(méi)有誤差, 僅由僅由 xn 到到 xn+1 這一步這一步計(jì)算由計(jì)算計(jì)算由計(jì)算格式引起的誤差格式引起的誤差.第22頁(yè)/共56頁(yè)23如Euler公式1(,)nnnnyyhf xy在點(diǎn)在點(diǎn) xn+1 的整體截?cái)嗾`差的整體截?cái)嗾`差 en+1 = y (xn+1)- - yn+122()(, ()()(,)2(, ()(,)()2nnnnnnnnnnnnnhy xhf xy xyyhf xyheh f xy xf xyy局部截?cái)嗾`差局部截?cái)嗾`差2111()()2nnnnhy xyy第23頁(yè)/共56頁(yè)24定義

18、定義7.1.2 若某數(shù)值方法的局部截若某數(shù)值方法的局部截?cái)鄶嗾`差為誤差為 則稱該方法具有則稱該方法具有 P 階精度階精度, 或稱其為或稱其為 P 階方法階方法.11(),pnO h可以證明可以證明:Euler 方法的局部截?cái)嗾`差方法的局部截?cái)嗾`差 其具有其具有一階一階精度精度. 梯形方法的局部截?cái)嗾`差梯形方法的局部截?cái)嗾`差 其具有其具有二階二階精度精度.改進(jìn)的改進(jìn)的 Euler 方法方法的局部截?cái)嗾`差的局部截?cái)嗾`差 具有具有二階二階精度精度.21(),nO h31(),nO h31(),nO h第24頁(yè)/共56頁(yè)257.2 Runge-Kutta 方法 繼續(xù)討論前面的繼續(xù)討論前面的 Taylo

19、r 展開(kāi)展開(kāi)法。法。設(shè)設(shè) y (x)C2a , b, 由由 Taylor 公式有公式有211()()()()().2nnnnnnnhy xy xhy xyxx 由由 故上式即為故上式即為 ( )( , ( ) ,y xf x y x)(! 2)(,()()(21nnnnnyhxyxhfxyxy 略去余項(xiàng)略去余項(xiàng), 并以并以 yn , yn+1 代替代替 y (xn) 及及 y (xn+1), 得到得到Euler 公式公式.第25頁(yè)/共56頁(yè)26進(jìn)一步假設(shè)設(shè) y (x)Cp+1a , b,由 Taylor 公式有) 1 . 2 . 7(,)(!)(! 2)()()()(21nnppnnnnRxy

20、phxyhxyhxyxy )2 . 2 . 7.(),()()!1(11)1(1nnnpnppnxxhOyphR其中其中 由由 故式故式(7.2.1)即為即為 ( )( , ( ) ,y xf x y x第26頁(yè)/共56頁(yè)277.2 Runge-Kutta 方法略去余項(xiàng)略去余項(xiàng), 并以數(shù)值解并以數(shù)值解 yn , yn+1 代替代替 (7.2.3) 中的中的 解析解解析解y (xn) 及及 y (xn+1), 可得到一個(gè)可得到一個(gè)差分方程，即差分方程，即nnppnnnnRxyphxyhxyhxyxy )(!)(! 2)()()()(21) 3 . 2 . 7(,)(,(!)(,(! 2)(,()

21、()()1(21nnnppnnnnnnRxyxfphxyxfhxyxhfxyxy., )()(,()!1(11)(1nnnpnnppnxxhOyfphR其中余項(xiàng)可寫(xiě)成其中余項(xiàng)可寫(xiě)成)(,(),()()(xyxfdxdyxfkkk注：這里注：這里第27頁(yè)/共56頁(yè)28在(7.2.3)中略去余項(xiàng),用yn , yn+1 代替 y (xn) 及 y (xn+1)4 . 2 . 7(.! 2! 2)1(1)1(21pnpnnnpnpnnnnfphfhfhyfphfhfhyy, ),(nnnyxff . ) 1, 2 , 1(),()()(piyxffnniin其中其中稱稱 (7.2.4) 式為求解常微分

22、方程初值問(wèn)題數(shù)值解式為求解常微分方程初值問(wèn)題數(shù)值解Taylor的格式的格式 .21(1)()()(,()(,()2!(,(),(7.2.3)!nnnnnnppnnnhy xy xhf xy xfxy xhfxy xRp第28頁(yè)/共56頁(yè)29由于局部截?cái)嗾`差 ),(O)(1111pnnnhyxy可知它是一個(gè)可知它是一個(gè) p 階方法。當(dāng)階方法。當(dāng)p=1時(shí)，上式正是時(shí)，上式正是Euler 公式。但當(dāng)公式。但當(dāng) p 2 時(shí)，需要計(jì)算時(shí)，需要計(jì)算f (x, y(x) ) 的高階導(dǎo)數(shù)，特別是對(duì)于復(fù)雜函數(shù)的高階導(dǎo)數(shù)，特別是對(duì)于復(fù)雜函數(shù) f (x, y(x) 的求導(dǎo)，這無(wú)疑是大大增加計(jì)的求導(dǎo)，這無(wú)疑是大大增加

23、計(jì)算量，這是它最大的缺點(diǎn)。因此高階的算量，這是它最大的缺點(diǎn)。因此高階的Taylor方法是不實(shí)用的。方法是不實(shí)用的。 德國(guó)數(shù)學(xué)家德國(guó)數(shù)學(xué)家C.Runge 及提出了一種改進(jìn)策略，得到了至今還被作為高精度的單及提出了一種改進(jìn)策略，得到了至今還被作為高精度的單步法廣泛使用龍格步法廣泛使用龍格-庫(kù)塔法庫(kù)塔法 ( Runge- Kutta method )。 第29頁(yè)/共56頁(yè)307.2.1 Runge-Kutta方法的基本思想.! 2)1(11pnpnnnnfphfhfhyyRunge-Kutta方法是方法是利用利用 f 在在某些點(diǎn)處函數(shù)值的線性組合替代（某些點(diǎn)處函數(shù)值的線性組合替代（7.2.4）步長(zhǎng)）

24、步長(zhǎng) h 后面括號(hào)中的因子來(lái)構(gòu)造差分方程后面括號(hào)中的因子來(lái)構(gòu)造差分方程, 從而避免了高階導(dǎo)數(shù)的計(jì)算從而避免了高階導(dǎo)數(shù)的計(jì)算, 這就是這就是 Runge-Kutta 方法的基本思想方法的基本思想. )1(1! 2pnpnnfphfhf用用f 在在某些點(diǎn)處函數(shù)值的線性組合替代這一部分某些點(diǎn)處函數(shù)值的線性組合替代這一部分第30頁(yè)/共56頁(yè)31其一般形式為:11111,(,),(,)(2,3, )rnniiinniininijjjyyhkkf xykf xh yhkir其中其中 r 是上式中調(diào)用是上式中調(diào)用 f 的個(gè)數(shù)的個(gè)數(shù), r 稱為級(jí)數(shù)，稱為級(jí)數(shù)， 為待定參數(shù)為待定參數(shù), 適當(dāng)確定這些參數(shù)適當(dāng)確定

25、這些參數(shù), 可使之具有盡可能高的精度可使之具有盡可能高的精度. 如局部截?cái)嗳缇植拷財(cái)嗾`差滿足誤差滿足,iii j . )(O)(1111rnnnhyxy第31頁(yè)/共56頁(yè)32 二階 Runge-Kutta 方法考慮考慮 r = 2 的情況的情況, 此時(shí)有此時(shí)有11 122121(),(,),(,).nnnnnnyyhkkkf xykf xh yh k利用二元函數(shù)的利用二元函數(shù)的 一階一階Taylor 公式，即全微分公式公式，即全微分公式 希望適當(dāng)選擇參數(shù)希望適當(dāng)選擇參數(shù) 使上式的局部截?cái)嗾`差為使上式的局部截?cái)嗾`差為12, , 3111()()nnny xyO h即為二階方法即為二階方法.,)(

26、)()(,()(,(),(),(22nnnnnynnnxnnyyxxOyyyxfxxyxfyxfyxf下面將下面將yn+1與與y(xn+1)作比較作比較第32頁(yè)/共56頁(yè)332222( ,)(,)(,)()(,)()()()() ()() ()()()nnxnnnynnnnndefnxnnynnnnf x yf xyfxyxxfxyyyOxxyyffxxfyyOxxyy. )O()()()(212hkfhfhfknynxn從而有, )(22111kkhyynn將上式代入),(12khyhxfknn再由第33頁(yè)/共56頁(yè)34得到. )()()()()O()()()()(3222121211221

27、11hOfffhfhyhkfhfhfhkhykkhyynynnxnnnynxnnnn在下面要將yn+1與y(xn+1)作比較，使它們的局部截?cái)嗾`差滿足. )O()()()(212hkfhfhfknynxn, )(O)(3111hyxynnn為此考慮y(xn+1)。第34頁(yè)/共56頁(yè)35再根據(jù) y(xn+1) 在點(diǎn) xn的一元 3 階 Taylor 展開(kāi)式. )()()(! 2)(),(! 2)()(! 2)()()(323),(2321hOfffhhfyhOyxfhhfyhOxyhxyhxyxynynnxnnyxnnnnnnnn )()()()(322211hOfffhfhyynynnxnnn

28、由剛才已得到的, )(O)(3111hyxynnn讓它們滿足第35頁(yè)/共56頁(yè)36即由122211212左式含有四個(gè)未知元三個(gè)左式含有四個(gè)未知元三個(gè)方程方程, 因此因此解不唯一解不唯一. 參數(shù)參數(shù)滿足左式的滿足左式的一族公式一族公式統(tǒng)稱統(tǒng)稱二階二階 Runge-Kutta 公式公式. 可得參數(shù)應(yīng)滿足下列方程組:. )()()(! 2)(321hOfffhhfyxynynnxnnn)()()()(322211hOfffhfhyynynnxnnn第36頁(yè)/共56頁(yè)37 取上式稱為上式稱為中點(diǎn)公式中點(diǎn)公式 .1210,1,212,nnyyhk1(,),nnkf xy21(,).22nnhhkf xy

29、k 取上式稱為上式稱為 Heun 公式公式 .112(3),4nnhyykk1(,),nnkf xy2122(,).33nnkf xh yhk12132,443 可見(jiàn)可見(jiàn), 二階二階 Runge-Kutta 公式公式, 每計(jì)算一步需要每計(jì)算一步需要 兩次調(diào)用兩次調(diào)用 f 的函數(shù)值的函數(shù)值. 取121,1,2112121(),2(,),(,).nnnnnnhyykkkf xykf xh yhk得這正是這正是改進(jìn)的改進(jìn)的 Euler 公式公式. 第37頁(yè)/共56頁(yè)38 四階Runge-Kutta方法當(dāng)當(dāng) r = 4 時(shí)時(shí), 類似地可導(dǎo)出四階類似地可導(dǎo)出四階 Runge-Kutta 公式公式, 這種

30、公式也有一族這種公式也有一族, 其中常用地有其中常用地有:標(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 yk第38頁(yè)/共56頁(yè)3911234(22)(22),6nnhyykkkk1(,),nnkf xy21(,),22nnhhkf xyk3122122(,),222nnhkf xyhkhk423222(,).22nnkf xh yhkhk Gill 公式公式 Gill 公式是標(biāo)準(zhǔn)的公式是標(biāo)準(zhǔn)的 Runge-Kutta

31、 公式的改進(jìn)形式公式的改進(jìn)形式, 這種算法可節(jié)省存儲(chǔ)單元這種算法可節(jié)省存儲(chǔ)單元, 并能控并能控制舍入誤差的增長(zhǎng)制舍入誤差的增長(zhǎng).四階四階 Runge-Kutta 公式公式, 每一步計(jì)算需四次調(diào)用每一步計(jì)算需四次調(diào)用 f 的函數(shù)值的函數(shù)值, 計(jì)算量較大計(jì)算量較大, 但其局部截?cái)嗟渚植拷財(cái)嗾`差可達(dá)誤差可達(dá) O(h5), 精度較高精度較高.第39頁(yè)/共56頁(yè)40例例7.2.1 用標(biāo)準(zhǔn)的用標(biāo)準(zhǔn)的 四階四階Rung-Kutta 法解初值問(wèn)題法解初值問(wèn)題,取取步長(zhǎng)步長(zhǎng)h=0.2.2, (0)1,(01).xyyyxy 解解: 解此問(wèn)題的計(jì)算公式為解此問(wèn)題的計(jì)算公式為112340.2(22),6nnyyk

32、kkk12,nnnxkyy21122,22()nnnhxhkykhyk32222,22()nnnhxhkykhyk4332,2()nnnxhhkykyhkxn0.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計(jì)算結(jié)果如下計(jì)算結(jié)果如下:顯然在計(jì)算量大致相同的情顯然在計(jì)算量大致相同的情況下況下, 標(biāo)準(zhǔn)的標(biāo)準(zhǔn)的 Runge-Kutta方方法比改進(jìn)的法比改進(jìn)的 Euler 方法精確方法精確度更高度更高. （參見(jiàn)（參見(jiàn)p227和和p228的表）的表）第40

33、頁(yè)/共56頁(yè)41一階常微分方程組與高階方程初值問(wèn)題的數(shù)值解法 7.5.1 一階常微分方程組初值問(wèn)題 )y,y,y(x,fdxdy)y,y,y(x,fdxdy)y,y,y(x,fdxdym21mmm2122m2111( )( )( )1122mmy a = s ,ya = s , ya = s ,xa,b第41頁(yè)/共56頁(yè)42寫(xiě)成向量形式： 0( ,)( ) (7.5.2)dYF x YdxY a =Y,)()()(Yn21n21xyxyxyyyy.s ,s ,sYTn210,)()()()(21212211mmmm,y,yx,yf,y,yx,yf,y,yx,yfx,YF其中其中注意注意: 在形

34、式上在形式上 (7.5.2) 與與 (7.0.1) 一樣一樣, 所以可以把求解常微分方程初值問(wèn)題的各所以可以把求解常微分方程初值問(wèn)題的各種數(shù)值方法推廣到方程組上來(lái)種數(shù)值方法推廣到方程組上來(lái). 第42頁(yè)/共56頁(yè)43利用向量值函數(shù)的微積分理論, 很容易推導(dǎo)出一階常微分方程組初值問(wèn)題的數(shù)值解法. )3 . 5 . 7(),(1nnnnYxhFYY如如Euler公式公式,),(T,2, 1nmnnnyyyY其中其中 .),(),(),(),(, 2, 1, 2, 12, 2, 11nmnnnmnmnnnnmnnnnnyyyxfyyyxfyyyxfYxF第43頁(yè)/共56頁(yè)44(7.5.3)的分量形式為

35、 )4 . 5 . 7(. ), 2, 1(),(, 2, 1,1,miyyyxhfyynmnnninini),(),(),(, 2, 1, 2, 12, 2, 11, 2, 11,1, 21, 1nmnnnmnmnnnnmnnnnmnnnmnnyyyxfyyyxfyyyxfhyyyyyy或或 第44頁(yè)/共56頁(yè)45四階標(biāo)準(zhǔn)的Runge-Kutta公式 )5 . 5 . 7()22(643211KKKKhYYnn),(1nnYxFK )2,2(12KhYhxFKnn)2,2(23KhYhxFKnn),(34hKYhxFKnn設(shè)設(shè) 則則 (7.5.5)的分量形式為的分量形式為 ,)4, 3, 2

36、, 1(),(T,2, 1jkkkKjmjjj第45頁(yè)/共56頁(yè)46四階標(biāo)準(zhǔn)的Runge-Kutta公式的分量形式 , )2,2,2,2(),2,2,2,2(),2,2,2,2(),(),22(63 ,3 , 23 , 14,2,2, 22, 13 ,1 ,1 , 21 , 12, 2, 11 ,4,3 ,2,1 ,1,inmininniiinmininniiinmininniinmnnniiiiiininikhykhykhyhxfkkhykhykhyhxfkkhykhykhyhxfkyyyxfkkkkkhyy., 2, 1mi其中其中第46頁(yè)/共56頁(yè)47例：試寫(xiě)出用中點(diǎn)公式解下列初值問(wèn)題的

37、計(jì)算公式. 3)0(, 1)0(,35,643zyzyzzyxy解解:令令 ,zyY則則 .nnnzyY.35643),(zyzyxYxF,)2 , 1(, 2, 1jkkKjjj再取再取由向量形式的中點(diǎn)公式由向量形式的中點(diǎn)公式第47頁(yè)/共56頁(yè)48中點(diǎn)公式的向量形式 )2,2(),(12121KhYhxFKYxFKhKYYnnnnnn上述中點(diǎn)公式的分量計(jì)算形式為上述中點(diǎn)公式的分量計(jì)算形式為2, 22, 111kkhzyzynnnnnnnnnzyzyxkk356431 , 21 , 1)2( 3)2( 5)2(6)2(4)2( 31 , 21 , 11 , 21 , 12, 22, 1khzkhykhzkhyhxkknnnnn分量計(jì)算形式為分量