計(jì)算方法龍格庫(kù)塔方法

1、2022-6-31得到高精度方法的一個(gè)直接想法是利用Taylor展開(kāi)假設(shè)式 y =f(x,y) (axb) 中的 f(x,y) 充分光滑,將y(xi+1)在x i點(diǎn)作Taylor展開(kāi)，若取右端不同的有限項(xiàng)作為y(xi+1)的近似值，就可得到計(jì)算y(xi+1)的各種不同截?cái)嗾`差的數(shù)值公式。例如：取前兩項(xiàng)可得到例如：取前兩項(xiàng)可得到)()()()(21hOxyhxyxyiii)(),()()(,()(22hOyxhfyhOxyxhfxyiiiiii9.4 龍格庫(kù)塔方法龍格庫(kù)塔方法2022-6-32)(21! 2PiPiiiiyPhyhyhyy 其中ffffffffff ffyf ffyxfyfyyy

2、xyyxyxxxyxiyxxiiii22)(2)(),(, P階泰勒方法若取前三項(xiàng)，可得到截?cái)嗾`差為若取前三項(xiàng)，可得到截?cái)嗾`差為O(h3)的公式的公式 )()(2)()()(321hOxyhxyhxyxyiiii )(),(),(),(2),(32hOyxfyxfyxfhyxhfyiiyiiiixiii類似地，若取前類似地，若取前P+1項(xiàng)作為項(xiàng)作為y(xi+1)的近似值，便得到的近似值，便得到2022-6-33顯然p=1時(shí), y i+1=y i+hf(xi,y i)它即為我們熟悉的Euler方法。當(dāng)p2時(shí),要利用泰勒方法就需要計(jì)算f(x,y)的高階微商。這個(gè)計(jì)算量是很大的，尤其當(dāng)f(x,y)較

3、復(fù)雜時(shí)，其高階導(dǎo)數(shù)會(huì)很復(fù)雜。因此，利用泰勒公式構(gòu)造高階公式是不實(shí)用的。但是泰勒級(jí)數(shù)展開(kāi)法的基本思想是許多數(shù)值方法的基礎(chǔ)。R-K方法不是直接使用Taylor級(jí)數(shù),而是利用它的思想2022-6-349.4.1 9.4.1 龍格龍格- -庫(kù)塔庫(kù)塔( (R-K) )法的基本思想法的基本思想Euler公式可改寫(xiě)成 ),(1iiiiyxhfKKyy則yi+1的表達(dá)式與y(xi+1)的Taylor展開(kāi)式的前兩項(xiàng)完全相同,即局部截?cái)嗾`差局部截?cái)嗾`差為O(h2)。Runge-Kutta 方法是一種高精度的單步法方法是一種高精度的單步法, ,簡(jiǎn)稱簡(jiǎn)稱R-K法法2022-6-35同理，改進(jìn)Euler公式可改寫(xiě)成 )

4、,(),(2121121211KyhxhfKyxhfKKKyyiiiiii 上述兩組公式在形式上共同點(diǎn):都是用f(x,y)在某些點(diǎn)上值的線性組合得出y(xi+1)的近似值yi+1, 且增加計(jì)算的次數(shù)f(x,y)的次數(shù),可提高截?cái)嗾`差的階。如歐拉法:每步計(jì)算一次f(x,y)的值,為一階方法。改進(jìn)歐拉法需計(jì)算兩次f(x,y)的值，為二階方法。局部截?cái)嗾`差局部截?cái)嗾`差為O(h3)2022-6-36 于是可考慮用函數(shù)f(x,y)在若干點(diǎn)上的函數(shù)值的線性組合來(lái)構(gòu)造近似公式，構(gòu)造時(shí)要求近似公式在(xi,yi)處的Taylor展開(kāi)式與解y(x)在xi處的Taylor展開(kāi)式的前面幾項(xiàng)重合，從而使近似公式達(dá)到所

5、需要的階數(shù)。既避免求高階導(dǎo)數(shù),又提高了計(jì)算方法精度的階數(shù)。或者說(shuō),在xi,xi+1這一步內(nèi)多計(jì)算幾個(gè)點(diǎn)的斜率值，然后將其進(jìn)行加權(quán)平均作為平均斜率，則可構(gòu)造出更高精度的計(jì)算格式，這就是龍格龍格庫(kù)塔庫(kù)塔（Runge-Kutta）法的基本思想法的基本思想。 ),(),(),(11,1112122122111ppppipipiiiippiiKbKbyhaxhfK KbyhaxhfKyxhfKKcKcKcyy一般龍格庫(kù)塔方法的形式為2022-6-37其中ai,bij,ci為待定參數(shù)，要求上式y(tǒng)i+1在點(diǎn)(xi,yi)處作Tailor展開(kāi)，通過(guò)相同項(xiàng)的系數(shù)確定參數(shù)。稱為P階龍格庫(kù)塔方法。8Runge-Ku

6、tta方法的推導(dǎo)思想0)(),(yaybxayxfy對(duì)于常微分方程的初值問(wèn)題的解y=y(x)，在區(qū)間xi, xi+1上使用微分中值定理，有)()()(1iiiyhxyxy即)()()(11iiiiixxyxyxy),(1iiixx其中2022-6-39上的平均斜率在區(qū)間可以認(rèn)為是,)(1iixxxyyKKxyxyii)()(1引入記號(hào))(,)(iiiyhfyhKKxxxyii上平均斜率的近似值間在區(qū)出只要使用適當(dāng)?shù)姆椒ㄇ?)(1就可得到相應(yīng)的Runge-Kutta方法ix1ixxy)(xyy Kyyii1K2022-6-310ix1ixxy)(xyy 如下圖Kxxxyxxyiii上的平均斜率在

7、處的斜率作為在如果以,)()(1即則上式化為),(1iiiiyxhfyy)(ixyhK)(,iixyxhf),(iiyxhf即Euler方法Euler方法也稱為一階一階Runge-Kutta方法方法KK2022-6-39.4.2 9.4.2 二階龍格二階龍格庫(kù)塔法庫(kù)塔法 在xi, xi+1上取兩點(diǎn)xi和xi+a2= xi +a2h,以該兩點(diǎn)處的斜率值K1和K2的加權(quán)平均(或稱為線性組合)來(lái)求取平均斜率k*的近似值K，即 2211KcKcK式中:K1 1為xi點(diǎn)處的切線斜率值 K1 =hf(xi, yi)=hy(xi) K2 2為xi +a2h點(diǎn)處的切線斜率值,比照改進(jìn)的歐拉法,將xi+a2視為

8、xi+1，即可得 ),(12122KbyhaxhfKii2022-6-311確定系數(shù) c1、c2、a2、b21 ，可得到有2階精度的算法格式2022-6-312因此 Kxyxyii)()(1)()(2211KcKcxyi將y(xi+1)在x=xi處進(jìn)行Taylor展開(kāi)： )()(! 2)()()(321hOxyhxyhxyxyiiii )(! 2),()(32hOf ffhyxhfxyyxiii將 在x=xi處進(jìn)行Taylor展開(kāi)： ),()(121222KbyhaxhfhaxyhKiii2022-6-313)(),(212122hOfKbf hayxfhKyxiiK1 =hf(xi, yi)

9、(),(3212hOfhfbf hayxfhyxii)(22111KcKcyyii)(),(32122hOfhfbf hayxfhcyxii),()(1iiiyxhfcxy),()()(21iiiyxhfccxy)(32221222hOf fhcbfhcayx2022-6-31421,21,12212221cbcacc這里有這里有 4 個(gè)未知個(gè)未知數(shù)，數(shù)，3 個(gè)方程。個(gè)方程。存在無(wú)窮多個(gè)解無(wú)窮多個(gè)解。所有滿足上式的格式統(tǒng)稱為2階階龍格龍格 - - 庫(kù)塔格式庫(kù)塔格式。令 11)(iiyxy對(duì)應(yīng)項(xiàng)的系數(shù)相等，得到 2022-6-315注意到，注意到， 就是二階就是二階龍格龍格 - - 庫(kù)塔庫(kù)塔公式

10、，也就是公式，也就是改進(jìn)的歐拉法。改進(jìn)的歐拉法。 21, 121212ccba),(),(21121211KyhxhfKyxhfKKKyyiiiiiiix1ixxy)(xyyK1K2K 因此，凡滿足條件式有一簇形如上式的計(jì)算格式，這些格式統(tǒng)稱為二階龍格庫(kù)塔格式。因此改進(jìn)的歐拉格式是眾多的二階龍格庫(kù)塔法中的一種特殊格式。 若取若取 ，就是另一種形式的二，就是另一種形式的二階階龍格龍格 - - 庫(kù)塔公式庫(kù)塔公式。 1, 0,2121212ccba2022-6-316)21,21(),(12121KyhxhfKyxhfKKyyiiiiii此計(jì)算公式稱為變形的二階龍格庫(kù)塔法。式中 為區(qū)間 的中點(diǎn)。也稱

11、中點(diǎn)公式也稱中點(diǎn)公式。 hxi211,iixxQ:為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？2022-6-317 二級(jí)R-K方法是顯式單步式,每前進(jìn)一步需要計(jì)算兩個(gè)函數(shù)值。由上面的討論可知,適當(dāng)選擇四個(gè)參數(shù)c1,c2,a2, b21，可使每步計(jì)算兩次函數(shù)值的二階R-K方法達(dá)到二階精度。能否在計(jì)算函數(shù)值次數(shù)不變的情況下,通過(guò)選擇不同的參數(shù)值,使得二階R-K方法的精度再提高呢? 答案是否定的！無(wú)論四個(gè)參數(shù)怎樣選擇,都不能使公式的局部截?cái)嗾`差提高到三階。 這說(shuō)明每一步計(jì)算兩個(gè)函數(shù)值的二階R-K方法最高階為二階。若要獲得更高階得數(shù)值方法若要獲得更高階得數(shù)值方法, ,就

12、必須增加計(jì)算函就必須增加計(jì)算函數(shù)值的次數(shù)。數(shù)值的次數(shù)。),(),(),(232131331212213322111KbKbyhaxhfKKbyhaxhfKyxhfKKcKcKcyyiiiiiiii9.4.3 三階龍格三階龍格庫(kù)塔法庫(kù)塔法2022-6-318 為進(jìn)一步提高精度，在區(qū)間xi, xi+1上除兩點(diǎn)xi和xi+a2= xi +a2h,以外，再增加一點(diǎn)xi+a3= xi +a3h ，用這三點(diǎn)處的斜率值K1、K2和K3的加權(quán)平均得出平均斜率K*的近似值K，這時(shí)計(jì)算格式具有形式: ix3aix2aixxy)(xyy K1K2K3K2022-6-319 同理推導(dǎo)二階公式，將y(xi+1)和yi+

13、1在x=xi處進(jìn)行Taylor展開(kāi)，使局部截?cái)嗾`差達(dá)到O(h4)，使對(duì)應(yīng)項(xiàng)的系數(shù)相等，得到系數(shù)方程組：61,6131)(31)(,3121)(,21133221232232323122213231332212323222332312213322321cbbcbacbbcbbbcacbacacacbbcbcacaccc參數(shù)的選擇不唯一，從而構(gòu)成一類不同的三階R-K公式，下面給出一種常用的三階R-K公式，形似simpson公式：)2,()21,2(),()4(612131213211KKyhxhfKKyhxhfKyxhfKKKKyyiiiiiiii2022-6-3202022-6-3219.4.4

14、 四階四階( (經(jīng)典經(jīng)典) )龍格龍格庫(kù)塔法庫(kù)塔法 如果需要再提高精度，用類似上述的處理方法，只需在區(qū)間xi,xi+1上用四個(gè)點(diǎn)處的斜率加權(quán)平均作為平均斜率K*的近似值，構(gòu)成一系列四階龍格庫(kù)塔公式。具有四階精度，即局部截?cái)嗾`差是O(h5)。 推導(dǎo)過(guò)程與前面類似，由于過(guò)程復(fù)雜，這里從略，只介紹最常用的一種四階經(jīng)典龍格四階經(jīng)典龍格庫(kù)塔公式庫(kù)塔公式。 2022-6-322 K1=hf (xi, yi) K2=hf (xi+a2h, yi+b21K1) K3=hf (xi+a3h, yi+b31K1+b32K2) K4=hf (xi+a4h, yi+b41K1+b42K2+b43K3) 其中c1、c2

15、、c3、c4、a2、a3、a4、b21、b31、b32、b41、b42、b43均為待定系數(shù)。這里K1、K2、K3、K4為四個(gè)不同點(diǎn)上的函數(shù)值，分別設(shè)其為 設(shè)yi+1=yi+c1K1+c2K2+c3K3+c4K42022-6-323 類似于前面的討論，把K2、K3、K4分別在xi點(diǎn)展成h的冪級(jí)數(shù)，代入線性組合式中，將得到的公式與y(xi+1)在xi點(diǎn)上的泰勒展開(kāi)式比較，使其兩式右端直到h4的系數(shù)相等，經(jīng)過(guò)較復(fù)雜的解方程過(guò)程便可得到關(guān)于ci，ai，bij的一組特解 a2=a3=b21=b32=1/2 b31=b41=b42=0 a4=b43=1 c1=c4=1/6 c2=c3=1/3 24 四階(

16、經(jīng)典)Runge-Kutta方法),()21,2()21,2(),()22(61342312143211KyhxhfKKyhxhfKKyhxhfKyxhfKKKKKyyiiiiiiiiii2022-6-325例1. 使用高階R-K方法計(jì)算初值問(wèn)題1)0(5 . 002yxyy. 1 . 0h取解:(1) 使用三階R-K方法時(shí)0i1 . 0201 hyK1103. 0)21(2102KyhK1256. 0)2(22103KKyhK1111. 1)4(6132101KKKyy2022-6-326其余結(jié)果如下:(2) 如果使用四階R-K方法時(shí)0i1 . 0201 hyK1103. 0)21(2102

17、KyhK i xi k1 k2 k3 yi 1.0000 0.1000 0.1000 0.1103 0.1256 1.1111 2.0000 0.2000 0.1235 0.1376 0.1595 1.2499 3.0000 0.3000 0.1562 0.1764 0.2092 1.4284 4.0000 0.4000 0.2040 0.2342 0.2866 1.6664 5.0000 0.5000 0.2777 0.3259 0.4163 1.99932022-6-3271113. 0)21(2203KyhK1235. 0)(2304KyhK)22(61 . 0432101KKKKyy1

18、111. 1其余結(jié)果如下: i xi k1 k2 k3 k4 yi 1.0000 0.1000 0.1000 0.1103 0.1113 0.1235 1.1111 2.0000 0.2000 0.1235 0.1376 0.1392 0.1563 1.2500 3.0000 0.3000 0.1562 0.1764 0.1791 0.2042 1.4286 4.0000 0.4000 0.2040 0.2342 0.2389 0.2781 1.6667 5.0000 0.5000 0.2777 0.3259 0.3348 0.4006 2.00002022-6-32022-6-328由上節(jié)分析常微分方程數(shù)值解法穩(wěn)定性問(wèn)題的方法，可得到各階Runge-Kutta公式的穩(wěn)定性條件：二階二階121122hh與歐拉預(yù)估校正公式一致三階三階1)(! 31)(! 21132hhh四階四階1)(! 41)(! 31)(! 211432hhhh9.4.5 龍格庫(kù)塔方法的穩(wěn)定性條件龍格庫(kù)塔方法的穩(wěn)定性條件2022-6-329 龍格庫(kù)塔方法的推導(dǎo)基于Taylor展開(kāi)方法，因而它要求所求的解