龍格庫塔方法基本原理

《龍格庫塔方法基本原理》由會員分享，可在線閱讀，更多相關(guān)《龍格庫塔方法基本原理（54頁珍藏版）》請在裝配圖網(wǎng)上搜索。

1、會計學1龍格庫塔方法基本原理龍格庫塔方法基本原理第1頁/共54頁2022-6-163得到高精度方法的一個直接想法是利用Taylor展開假設式 y =f(x,y) (axb) 中的 f(x,y) 充分光滑,將y(xi+1)在x i點作Taylor展開，若取右端不同的有限項作為y(xi+1)的近似值，就可得到計算y(xi+1)的各種不同截斷誤差的數(shù)值公式。例如：取前兩項可得到例如：取前兩項可得到)()()()(21hOxyhxyxyiii)(),()()(,()(22hOyxhfyhOxyxhfxyiiiiii3.2 龍格庫塔方法龍格庫塔方法第2頁/共54頁2022-6-164)(21! 2PiP

2、iiiiyPhyhyhyy 其中ffffffffff ffyf ffyxfyfyyyxyyxyxxxyxiyxxiiii22)(2)(),(, P階泰勒方法若取前三項，可得到截斷誤差為若取前三項，可得到截斷誤差為O(h3)的公式的公式 )()(2)()()(321hOxyhxyhxyxyiiii )(),(),(),(2),(32hOyxfyxfyxfhyxhfyiiyiiiixiii類似地，若取前類似地，若取前P+1項作為項作為y(xi+1)的近似值，便得到的近似值，便得到第3頁/共54頁2022-6-165顯然p=1時, y i+1=y i+hf(xi,y i)它即為我們熟悉的Euler方

3、法。當p2時,要利用泰勒方法就需要計算f(x,y)的高階微商。這個計算量是很大的，尤其當f(x,y)較復雜時，其高階導數(shù)會很復雜。因此，利用泰勒公式構(gòu)造高階公式是不實用的。但是泰勒級數(shù)展開法的基本思想是許多數(shù)值方法的基礎。R-K方法不是直接使用Taylor級數(shù),而是利用它的思想第4頁/共54頁2022-6-166龍格龍格- -庫塔庫塔( (R-K) )法的基本思想法的基本思想Euler公式可改寫成 ),(1iiiiyxhfKKyy則yi+1的表達式與y(xi+1)的Taylor展開式的前兩項完全相同,即局部截斷誤差局部截斷誤差為O(h2)。Runge-Kutta 方法是一種高精度的單步法方法是

4、一種高精度的單步法, ,簡稱簡稱R-K法法第5頁/共54頁2022-6-167同理，改進Euler公式可改寫成 ),(),(2121121211KyhxhfKyxhfKKKyyiiiiii 上述兩組公式在形式上共同點:都是用f(x,y)在某些點上值的線性組合得出y(xi+1)的近似值yi+1, 且增加計算的次數(shù)f(x,y)的次數(shù),可提高截斷誤差的階。如歐拉法:每步計算一次f(x,y)的值,為一階方法。改進歐拉法需計算兩次f(x,y)的值，為二階方法。局部截斷誤差局部截斷誤差為O(h3)第6頁/共54頁2022-6-168 于是可考慮用函數(shù)f(x,y)在若干點上的函數(shù)值的線性組合來構(gòu)造近似公式，

5、構(gòu)造時要求近似公式在(xi,yi)處的Taylor展開式與解y(x)在xi處的Taylor展開式的前面幾項重合，從而使近似公式達到所需要的階數(shù)。既避免求高階導數(shù),又提高了計算方法精度的階數(shù)?；蛘哒f,在xi,xi+1這一步內(nèi)多計算幾個點的斜率值，然后將其進行加權(quán)平均作為平均斜率，則可構(gòu)造出更高精度的計算格式，這就是龍格龍格庫塔庫塔（Runge-Kutta）法的基本思想法的基本思想。 第7頁/共54頁),(),(),(11,1112122122111ppppipipiiiippiiKbKbyhaxhfK KbyhaxhfKyxhfKKcKcKcyy一般龍格庫塔方法的形式為2022-6-169其中a

6、i,bij,ci為待定參數(shù)，要求上式y(tǒng)i+1在點(xi,yi)處作Tailor展開，通過相同項的系數(shù)確定參數(shù)。稱為P階龍格庫塔方法。第8頁/共54頁10Runge-Kutta方法的推導思想0)(),(yaybxayxfy對于常微分方程的初值問題的解y=y(x)，在區(qū)間xi, xi+1上使用微分中值定理，有)()()(1iiiyhxyxy即)()()(11iiiiixxyxyxy),(1iiixx其中2022-6-16第9頁/共54頁11上的平均斜率在區(qū)間可以認為是,)(1iixxxyyKKxyxyii)()(1引入記號)(,)(iiiyhfyhKKxxxyii上平均斜率的近似值間在區(qū)出只要使用

7、適當?shù)姆椒ㄇ?)(1就可得到相應的Runge-Kutta方法ix1ixxy)(xyy Kyyii1K2022-6-16第10頁/共54頁12ix1ixxy)(xyy 如下圖Kxxxyxxyiii上的平均斜率在處的斜率作為在如果以,)()(1即則上式化為),(1iiiiyxhfyy)(ixyhK)(,iixyxhf),(iiyxhf即Euler方法Euler方法也稱為一階一階Runge-Kutta方法方法KK2022-6-16第11頁/共54頁二階龍格二階龍格庫塔法庫塔法 在xi, xi+1上取兩點xi和xi+a2= xi +a2h,以該兩點處的斜率值K1和K2的加權(quán)平均(或稱為線性組合)來求取

8、平均斜率k*的近似值K，即 2211KcKcK式中:K1 1為xi點處的切線斜率值 K1 =hf(xi, yi)=hy(xi) K2 2為xi +a2h點處的切線斜率值,比照改進的歐拉法,將xi+a2視為xi+1，即可得 ),(12122KbyhaxhfKii2022-6-1613確定系數(shù) c1、c2、a2、b21 ，可得到有2階精度的算法格式第12頁/共54頁2022-6-1614因此 Kxyxyii)()(1)()(2211KcKcxyi將y(xi+1)在x=xi處進行Taylor展開： )()(! 2)()()(321hOxyhxyhxyxyiiii )(! 2),()(32hOf ff

9、hyxhfxyyxiii將 在x=xi處進行Taylor展開： ),()(121222KbyhaxhfhaxyhKiii第13頁/共54頁2022-6-1615)(),(212122hOfKbf hayxfhKyxiiK1 =hf(xi, yi)(),(3212hOfhfbf hayxfhyxii)(22111KcKcyyii)(),(32122hOfhfbf hayxfhcyxii),()(1iiiyxhfcxy),()()(21iiiyxhfccxy)(32221222hOf fhcbfhcayx第14頁/共54頁2022-6-161621,21,12212221cbcacc這里有這里有

10、4 個未知個未知數(shù)，數(shù)，3 個方程。個方程。存在無窮多個解無窮多個解。所有滿足上式的格式統(tǒng)稱為2階階龍格龍格 - - 庫塔格式庫塔格式。令 11)(iiyxy對應項的系數(shù)相等，得到 第15頁/共54頁2022-6-1617注意到，注意到， 就是二階就是二階龍格龍格 - - 庫塔庫塔公式，也就是公式，也就是改進的歐拉法。改進的歐拉法。 21, 121212ccba),(),(21121211KyhxhfKyxhfKKKyyiiiiiiix1ixxy)(xyyK1K2K 因此，凡滿足條件式有一簇形如上式的計算格式，這些格式統(tǒng)稱為二階龍格庫塔格式。因此改進的歐拉格式是眾多的二階龍格庫塔法中的一種特殊

11、格式。 第16頁/共54頁若取若取 ，就是另一種形式的二，就是另一種形式的二階階龍格龍格 - - 庫塔公式庫塔公式。 1, 0,2121212ccba2022-6-1618)21,21(),(12121KyhxhfKyxhfKKyyiiiiii此計算公式稱為變形的二階龍格庫塔法。式中 為區(qū)間 的中點。也稱中點公式也稱中點公式。 hxi211,iixxQ:為獲得更高的精度，應該如何進一步推廣？為獲得更高的精度，應該如何進一步推廣？第17頁/共54頁2022-6-1619 二級R-K方法是顯式單步式,每前進一步需要計算兩個函數(shù)值。由上面的討論可知,適當選擇四個參數(shù)c1,c2,a2, b21，可使每

12、步計算兩次函數(shù)值的二階R-K方法達到二階精度。能否在計算函數(shù)值次數(shù)不變的情況下,通過選擇不同的參數(shù)值,使得二階R-K方法的精度再提高呢? 答案是否定的！無論四個參數(shù)怎樣選擇,都不能使公式的局部截斷誤差提高到三階。 這說明每一步計算兩個函數(shù)值的二階R-K方法最高階為二階。若要獲得更高階得數(shù)值方法若要獲得更高階得數(shù)值方法, ,就必須增加計算就必須增加計算函數(shù)值的次數(shù)。函數(shù)值的次數(shù)。第18頁/共54頁),(),(),(232131331212213322111KbKbyhaxhfKKbyhaxhfKyxhfKKcKcKcyyiiiiiiii三階龍格三階龍格庫塔法庫塔法2022-6-1620 為進一步

13、提高精度，在區(qū)間xi, xi+1上除兩點xi和xi+a2= xi +a2h,以外，再增加一點xi+a3= xi +a3h ，用這三點處的斜率值K1、K2和K3的加權(quán)平均得出平均斜率K*的近似值K，這時計算格式具有形式: ix3aix2aixxy)(xyy K1K2K3K第19頁/共54頁2022-6-1621 同理推導二階公式，將y(xi+1)和yi+1在x=xi處進行Taylor展開，使局部截斷誤差達到O(h4)，使對應項的系數(shù)相等，得到系數(shù)方程組：61,6131)(31)(,3121)(,211332212322323231222132313322123232223323122133223

14、21cbbcbacbbcbbbcacbacacacbbcbcacaccc第20頁/共54頁參數(shù)的選擇不唯一，從而構(gòu)成一類不同的三階R-K公式，下面給出一種常用的三階R-K公式，形似simpson公式：)2,()21,2(),()4(612131213211KKyhxhfKKyhxhfKyxhfKKKKyyiiiiiiii2022-6-1622第21頁/共54頁2022-6-1623四階四階( (經(jīng)典經(jīng)典) )龍格龍格庫塔法庫塔法 如果需要再提高精度，用類似上述的處理方法，只需在區(qū)間xi,xi+1上用四個點處的斜率加權(quán)平均作為平均斜率K*的近似值，構(gòu)成一系列四階龍格庫塔公式。具有四階精度，即局部

15、截斷誤差是O(h5)。 推導過程與前面類似，由于過程復雜，這里從略，只介紹最常用的一種四階經(jīng)典龍格四階經(jīng)典龍格庫塔公式庫塔公式。 第22頁/共54頁2022-6-1624 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、c3、c4、a2、a3、a4、b21、b31、b32、b41、b42、b43均為待定系數(shù)。這里K1、K2、K3、K4為四個不同點上的函數(shù)值，分別設其為 設yi+1=yi+c1K1+c2K2+c3K

16、3+c4K4第23頁/共54頁2022-6-1625 類似于前面的討論，把K2、K3、K4分別在xi點展成h的冪級數(shù)，代入線性組合式中，將得到的公式與y(xi+1)在xi點上的泰勒展開式比較，使其兩式右端直到h4的系數(shù)相等，經(jīng)過較復雜的解方程過程便可得到關(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頁/共54頁26 四階(經(jīng)典)Runge-Kutta方法),()21,2()21,2(),()22(61342312143211KyhxhfKKyhxhfKKyhxhfKyxhfKKK

17、KKyyiiiiiiiiii2022-6-16第25頁/共54頁27例1. 使用高階R-K方法計算初值問題1)0(5 . 002yxyy. 1 . 0h取解:(1) 使用三階R-K方法時0i1 . 0201 hyK1103. 0)21(2102KyhK1256. 0)2(22103KKyhK1111. 1)4(6132101KKKyy2022-6-16第26頁/共54頁28其余結(jié)果如下:(2) 如果使用四階R-K方法時0i1 . 0201 hyK1103. 0)21(2102KyhK i xi k1 k2 k3 yi 1.0000 0.1000 0.1000 0.1103 0.1256 1.1

18、111 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-16第27頁/共54頁291113. 0)21(2203KyhK1235. 0)(2304KyhK)22(61 . 0432101KKKKyy1111. 1其余結(jié)果如下: i xi k1 k2 k3 k4 yi 1.0000 0.1000 0.

19、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-16第28頁/共54頁【例】 已知一階系統(tǒng)的微分方程為： ，初始條件 ，取仿真步長h=0.1，分別用歐拉法、梯形法和龍格庫塔法計算

20、該系統(tǒng)仿真第一步的值。102ydtdy1)(00 yty解：原方程可變?yōu)? ),(210kkytfydtdy即 1210),(0yyytfkkk數(shù)值積分公式應用數(shù)值積分公式應用 第29頁/共54頁 （1）用歐拉法計算 根據(jù)歐拉公式，將函數(shù)表達式及其初始值代入后，可得該系統(tǒng)仿真第一步的值：8 .1)1210(1 .01),(0001ythfyy數(shù)值積分公式應用數(shù)值積分公式應用 ),(1kkkkythfyy第30頁/共54頁（2）用梯形法計算： 根據(jù)預報校正公式，將函數(shù)表達式及其初始值代入后，可得仿真第一步的值。 用預報公式求起始值：8 . 1) 1210(1 . 01),(0001)0(ythf

21、yy數(shù)值積分公式應用數(shù)值積分公式應用 ),(),(21),(1)0(111)0(kkkkkkkkkkytfytfhyyythfyy第31頁/共54頁再用校正公式得到系統(tǒng)仿真第一步的值：72. 1)8 . 1210() 1210(1 . 0211),(),(211)0(10001ytfytfhyy數(shù)值積分公式應用數(shù)值積分公式應用 第32頁/共54頁)(2)(2111kkhyytynnn),(),(121hkyhtfkytfknnnn第33頁/共54頁（3）用二階龍格庫塔法計算 根據(jù)公式先計算出兩個系數(shù)，再計算仿真第一步的值： 812100),(001yytfk4 . 6)81 . 01 (210

22、)(210),(101002hkyhkyhtfk數(shù)值積分公式應用數(shù)值積分公式應用 第34頁/共54頁則系統(tǒng)仿真第一步的值為：72. 1)4 . 68(1 . 0211)(22101kkhyy數(shù)值積分公式應用數(shù)值積分公式應用 第35頁/共54頁),()2,2()2,2(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyykkkkkkkkkk四階龍格四階龍格庫塔（庫塔（RungeKutta）法）法5h第36頁/共54頁（4）用四階龍格庫塔公式計算根據(jù)公式先計算出4個系數(shù)，再計算仿真第一步的值： 81210210),(0001yytfk2 . 7)

23、81 . 0211 (210)2(210)2,2(101002khykhyhtfk數(shù)值積分公式應用數(shù)值積分公式應用 第37頁/共54頁28. 7)2 . 71 . 0211 (210)2(210)2,2(202003khykhyhtfk544. 6)28. 71 . 01 (210)(210),(303004hkyhkyhtfk數(shù)值積分公式應用數(shù)值積分公式應用 第38頁/共54頁則系統(tǒng)仿真第一步的值為：725067. 1)544. 628. 722 . 728(1 . 0611)22(6432101kkkkhyy數(shù)值積分公式應用數(shù)值積分公式應用 第39頁/共54頁從上述結(jié)果可以看出從上述結(jié)果可

24、以看出: 對于同一個系統(tǒng)進行仿真計算時，其值的精度是隨對于同一個系統(tǒng)進行仿真計算時，其值的精度是隨著數(shù)值積分公式的變化而改變的，其中歐拉法計算精著數(shù)值積分公式的變化而改變的，其中歐拉法計算精度最低，其次為梯形法和二階龍格度最低，其次為梯形法和二階龍格庫塔法，四階龍庫塔法，四階龍格格庫塔法計算精度最高。庫塔法計算精度最高。數(shù)值積分公式應用數(shù)值積分公式應用 例例2：用：用matlab演示演示第40頁/共54頁第41頁/共54頁)k(k512hyyE)255k(k256hyy)255k(k256hyh,ftk)k2hy,2hf(tk)y,f(tk)k510k(k512hyy311n1nk21k1n2

25、1nn31nn2nn1321n1n第42頁/共54頁)k8k9k(2k6hE)3k4k(k2hyh,ftk)3k(k8hy,2hftk)k(k6hy,3hftk)k3hy,3hf(tk)y,f(tk)k4k(k6hyy5431n341nn531nn421nn31nn2nn1541k1n第43頁/共54頁第44頁/共54頁)(,(tuyfdtdy),(),()(211121211hkyutfkyutfkkkhyynnnnnnnn第45頁/共54頁1nu采入2k計算1k計算nt2htn1nntht21htn21nntht?第46頁/共54頁圖圖 RK-2 2的計算流程的計算流程nt2htn1nu采

26、入并輸出計算1ny1k計算下一個2k計算1k計算1nntht21htn21nntht第47頁/共54頁)2,(),(12/12/12121khyutfkyutfkhkyynnnnnnnn 其計算流程：假定在h/2的時間內(nèi)計算機也剛好計算一次右端函數(shù)f,則計算也分為兩步： 1 在tn時刻利用當前的un、yn計算K1； 2 在tn+h/2時刻計算K2,此時un+1/2可以得到，不會引入新的誤差，可實時得到y(tǒng)n+1圖圖 實時實時RK-2RK-2公式計算流程公式計算流程1k計算2k計算1k計算下一個)2/(htun采入并輸出計算1nynt2htn1nntht21htn21nntht第48頁/共54頁第

27、49頁/共54頁2022-6-1651由上節(jié)分析常微分方程數(shù)值解法穩(wěn)定性問題的方法，可得到各階Runge-Kutta公式的穩(wěn)定性條件：二階二階121122hh與歐拉預估校正公式一致三階三階1)(! 31)(! 21132hhh四階四階1)(! 41)(! 31)(! 211432hhhh龍格庫塔方法的穩(wěn)定性條件龍格庫塔方法的穩(wěn)定性條件第50頁/共54頁2022-6-1652 龍格庫塔方法的推導基于Taylor展開方法，因而它要求所求的解具有較好的光滑性。如果解的光滑性差，那么，使用四階龍格庫塔方法求得的數(shù)值解，其精度可能反而不如改進的歐拉方法。在實際計算時，應當針對問題的具體特點選擇合適的算法

28、。對于光滑性不太好的解，最好采用對于光滑性不太好的解，最好采用低階算低階算法法而將步長而將步長h 取小取小。 第51頁/共54頁2022-6-1653 前面已經(jīng)看到,二階、四階R-K方法可分別達到最高階數(shù)2 2階、4 4階，但是N階R-K方法的最高階卻不一定是N階。R-K方法的級數(shù)表示公式中計算函數(shù)值f 的次數(shù)。Butcher于1965年給出了R-K方法計算函數(shù)值f 的次數(shù)與可達到的最高精度階數(shù)之間的關(guān)系表，如下： 由表可見，四級以下R-K的方法其最高階數(shù)與計算f 的次數(shù)一致，對m階R-K公式，當m4，雖然計算f 的次數(shù)增加，但是方法階數(shù)不一定增加。因此四階R-K公式是應用最為廣泛的公式。753可達到的最高精度可達到的最高精度642每步須算每步須算Ki 的個數(shù)的個數(shù))(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n第52頁/共54頁2022-6-1654 順便指出，當常微分方程中的f (t,u)與u無關(guān)時，常微分方程初值問題便簡化為計算定積分的問題。這時，第一節(jié)介紹的Euler方法就是求定積分的矩形公式，改進的Euler方法就是求定積分的梯形公式，而三階Runge-Kutta公式就是計算定積分的Simpson公式，它們的精度也是一致的。第53頁/共54頁