常微分方程的龍格-庫塔方法

1、第八章常微分方程數(shù)值解法8.2 Runge-Kutta 方法方法8.2.2 幾類顯式幾類顯式Runge-Kutta方法方法8.2.1 Runge-Kutta方法的基本思想方法的基本思想第八章常微分方程數(shù)值解法8.2.1 Runge-Kutta方法的基本思想方法的基本思想 顯式顯式Euler方法是最簡單的單步法，它是一階的，它可以看作方法是最簡單的單步法，它是一階的，它可以看作Talylor展開后展開后取前兩項(xiàng)。因此，得到高階方法的一個直接想法是用取前兩項(xiàng)。因此，得到高階方法的一個直接想法是用Talylor展開，如果能計(jì)算展開，如果能計(jì)算 的高階導(dǎo)數(shù)，則可寫出的高階導(dǎo)數(shù)，則可寫出p階方法的計(jì)算方

2、法階方法的計(jì)算方法 xy ，pnpnnnnyphyhyhyy!221 其中其中 是是 的近似值，的近似值， 若將若將 分別記成分別記成 則對于二階和三階導(dǎo)數(shù)可表示為則對于二階和三階導(dǎo)數(shù)可表示為 jny njxy。，pj210，yfxfyxf，yxfff。，fffffffffyfffyyyyyxxyxxyx222 第八章常微分方程數(shù)值解法這個方法并不實(shí)用，因?yàn)橐话闱闆r下，求這個方法并不實(shí)用，因?yàn)橐话闱闆r下，求 的導(dǎo)數(shù)相當(dāng)麻煩。從計(jì)算高的導(dǎo)數(shù)相當(dāng)麻煩。從計(jì)算高階導(dǎo)數(shù)的公式知道，方法的截?cái)嗾`差提高一階，需要增加的計(jì)算量很大。但階導(dǎo)數(shù)的公式知道，方法的截?cái)嗾`差提高一階，需要增加的計(jì)算量很大。但是是由此

3、啟發(fā)我們用區(qū)間上若干個點(diǎn)的導(dǎo)數(shù)由此啟發(fā)我們用區(qū)間上若干個點(diǎn)的導(dǎo)數(shù) ，而不是高階導(dǎo)數(shù)，將它們作線性，而不是高階導(dǎo)數(shù)，將它們作線性組合得到平均斜率，將其與解的組合得到平均斜率，將其與解的Taylor展開相比較，使前面若干項(xiàng)吻合，從而展開相比較，使前面若干項(xiàng)吻合，從而得到具有一定階的方法。這就是得到具有一定階的方法。這就是Runge-Kutta方法的基本思想，其一般形式為方法的基本思想，其一般形式為yxf，f，LiKahcyhcxfKyxfKKhyyijjijinininniLiinn32,11111(8.2.1)第八章常微分方程數(shù)值解法其中，其中， 與與 的區(qū)別在于：用微分方程準(zhǔn)確解的區(qū)別在于：用

4、微分方程準(zhǔn)確解 代替代替 中的中的 就得就得到到 。參數(shù)。參數(shù) 和和 待定，確定它們的原則和方法是：將（待定，確定它們的原則和方法是：將（8.2.2）式中）式中的的 在在 處作處作Taylor展開，將展開，將 在在 處作二元處作二元Taylor展開，展開，將展開式按將展開式按H的冪次整理后，令的冪次整理后，令 中中h的低次冪的系數(shù)為零，使的低次冪的系數(shù)為零，使 首項(xiàng)中首項(xiàng)中h的冪次盡量高，比如使的冪次盡量高，比如使 ，則稱（，則稱（8.2.1）式為）式為L級級p階階Runge-Kutta方法方法（簡稱（簡稱R-K法法）。）。 *iKiK nxyiKny*iKiic， ija 1 nxynx n

5、nxyx ，1 nT1 nT 11 pnhOT*iK其中，其中， 。它的局部截?cái)嗾`差是。它的局部截?cái)嗾`差是111111 ijijLiiiac， *111iLiinnnKhxyxyT ， （8.2.2）第八章常微分方程數(shù)值解法它與顯式它與顯式R-K公式的區(qū)別在于：顯式公式中，對系數(shù)公式的區(qū)別在于：顯式公式中，對系數(shù) 求和的上限是求和的上限是 ，從，從而而 構(gòu)成的矩陣是一個嚴(yán)格下三角陣。而在隱式公式中，對系數(shù)構(gòu)成的矩陣是一個嚴(yán)格下三角陣。而在隱式公式中，對系數(shù) 求和的上求和的上限是限是L，從而，從而 構(gòu)成的矩陣是方陣，需要用迭代法求出近似斜率構(gòu)成的矩陣是方陣，需要用迭代法求出近似斜率 推導(dǎo)隱式公式

6、的思路和方法與顯式公式法類似。推導(dǎo)隱式公式的思路和方法與顯式公式法類似。ija1 iijaijaija LiKi，21 類似于顯式類似于顯式R-K公式（公式（8.2.1），稍加改變，就得到隱式），稍加改變，就得到隱式R-K公式公式，iLiinnKhyy11。，LiKahcyhcxfKLjjijinini211第八章常微分方程數(shù)值解法8.2.2 幾類顯式幾類顯式Runge-Kutta方法方法 對于對于L=2，則，則。，1222122111hKcyhcxfKyxfKKKhyynnnnnn其局部截?cái)嗾`差是其局部截?cái)嗾`差是 *22*1111KKhxyxyTnnn（8.2.3）第八章常微分方程數(shù)值解法

7、將將 中的各項(xiàng)作中的各項(xiàng)作Taylor展開，并利用展開，并利用 則有則有1nT ，fffyxyxfxyyxnnn ，nnnnnnnnnnnxyhcxyhcxfKxyxyxfKhOxyhxyhxyhxyxy 22*2*1432162 ，32222 222hOfffffhcxhycxyyyxyxxnn將它們代入（將它們代入（8.2.3）式，整理后得）式，整理后得 nnnxyhcxyhT 222211211 4222232261hOfffffcxyhyyxyxxn 第八章常微分方程數(shù)值解法選取選取 和和 ，使方法的階盡可能高，就是使，使方法的階盡可能高，就是使 h 和和 的系數(shù)為零，因?yàn)榈南禂?shù)為零，

8、因?yàn)?的系數(shù)一般不為零。于是得到方程組的系數(shù)一般不為零。于是得到方程組21，2c2h3h。，2112221c顯然，該方程組有無窮多組解，從而得到一族顯然，該方程組有無窮多組解，從而得到一族二級二階二級二階R-K方法方法。2c若以若以 為自由參數(shù)，取為自由參數(shù)，取 得得中點(diǎn)公式中點(diǎn)公式212c，nnnnnnyxfhyhxhfyy221 （8.2.4）取取c=2/3得得Heun公式公式nnnnnnnnyxhfyhxfyxfhyy，3232341， （8.2.5）第八章常微分方程數(shù)值解法取取c=1得改進(jìn)的得改進(jìn)的Euler公式（公式（8.1.6）。）。對于對于L=3的情形，要計(jì)算三個斜率的近似值：的

9、情形，要計(jì)算三個斜率的近似值： 。，23213133312221KaKahcyhcxfKhKcyhcxfKyxfKnnnnnn 類似于二階方法的推導(dǎo)，可以得三階的方法，所得系數(shù)應(yīng)滿足的方程組是類似于二階方法的推導(dǎo)，可以得三階的方法，所得系數(shù)應(yīng)滿足的方程組是。，161312111323132323233222332221321aaacccccca該方程組的解也是該方程組的解也是 不唯一的。常見的一種不唯一的。常見的一種三級三階方法三級三階方法是是第八章常微分方程數(shù)值解法 。，2131221321122246hKhKyhxfKKhyhxfKyxfKKKKhyynnnnnnn 對于對于L=4的情形，

10、可進(jìn)行類似推導(dǎo)。最常用的四級四階方法是如下的情形，可進(jìn)行類似推導(dǎo)。最常用的四級四階方法是如下經(jīng)典經(jīng)典R-K方法方法 。，3423121432112222226hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn （8.2.6）第八章常微分方程數(shù)值解法為了分析經(jīng)典為了分析經(jīng)典R-K公式的計(jì)算量和計(jì)算精度，將四階經(jīng)典公式的計(jì)算量和計(jì)算精度，將四階經(jīng)典R-K公式（公式（8.2.6）與一階顯式與一階顯式Euler公式（公式（8.1.2）及二階改進(jìn)的）及二階改進(jìn)的Euler公式相比較。一般說來，公公式相比較。一般說來，公式的級數(shù)越大，計(jì)算右端項(xiàng)式的級數(shù)越大，計(jì)算右端項(xiàng) f

11、 的次數(shù)越多，計(jì)算量越大。在同樣步長的情況的次數(shù)越多，計(jì)算量越大。在同樣步長的情況下，下，Euler方法每步只計(jì)算一個函數(shù)值，而經(jīng)典方法要計(jì)算方法每步只計(jì)算一個函數(shù)值，而經(jīng)典方法要計(jì)算4個函數(shù)值。四階個函數(shù)值。四階R-K法的計(jì)算量差不多是改進(jìn)的法的計(jì)算量差不多是改進(jìn)的Euler公式的公式的2倍，是顯式倍，是顯式Euler公式的公式的4倍。下面倍。下面的例子中的例子中Euler方法用步長方法用步長 ，二階改進(jìn)的，二階改進(jìn)的Euler法用步長法用步長 ，而四階經(jīng)典公，而四階經(jīng)典公式用步長式用步長 。這樣，從。這樣，從 到到 三種方法都計(jì)算了三種方法都計(jì)算了4個函數(shù)制，計(jì)算個函數(shù)制，計(jì)算量大體相當(dāng)。

12、量大體相當(dāng)。12h14h1hnx14hxn例例8.3 考慮初值問題考慮初值問題 。，001 yyy第八章常微分方程數(shù)值解法其解析解為其解析解為 。分別用。分別用h=0.025的顯式的顯式Euler方法，方法，h=0.05改進(jìn)改進(jìn)Euler法和法和h=0.1的經(jīng)典的經(jīng)典R-K方法計(jì)算到方法計(jì)算到x=0.5。三種方法在。三種方法在x方方向每前進(jìn)向每前進(jìn)0.1都要計(jì)算都要計(jì)算4個右端函數(shù)值，計(jì)算量相當(dāng)。計(jì)算結(jié)果列于個右端函數(shù)值，計(jì)算量相當(dāng)。計(jì)算結(jié)果列于表表8-3。從計(jì)算結(jié)果看，在工作量大致相同的情況下，還是經(jīng)典方法。從計(jì)算結(jié)果看，在工作量大致相同的情況下，還是經(jīng)典方法比其他兩種方法的結(jié)果好得多。在比

13、其他兩種方法的結(jié)果好得多。在x=0.5處，三種方法的誤差分別處，三種方法的誤差分別是是 。經(jīng)典。經(jīng)典R-K法對多數(shù)好條件問法對多數(shù)好條件問題（題（ ，參見下節(jié)單步法的穩(wěn)定性），能獲得好的效果。，參見下節(jié)單步法的穩(wěn)定性），能獲得好的效果。 xexy 1743108 .2103 .1108 .3 ，0yf第八章常微分方程數(shù)值解法表表 8-3Euler法法 改進(jìn)改進(jìn)Euler法法 經(jīng)典經(jīng)典R-K法法 準(zhǔn)確解準(zhǔn)確解h=0.025 h=0.05 h=0.1 nx 0.1 0.096312 0.095123 0.09516250 0.09516258 0.2 0.183348 0.181193 0.181

14、26910 0.18126925 0.3 0.262001 0.259085 0.25918158 0.25918178 0.4 0.333079 0.329563 0.32967971 0.32967995 0.5 0.397312 0.393337 0.39346906 0.39346934 nxy 在微分方程數(shù)值解法的實(shí)際計(jì)算中，有個如何選擇步長的問題。在微分方程數(shù)值解法的實(shí)際計(jì)算中，有個如何選擇步長的問題。因?yàn)閱螐拿恳徊娇矗介L越小，截?cái)嗾`差越小。但隨著步長的縮小，因?yàn)閱螐拿恳徊娇?，步長越小，截?cái)嗾`差越小。但隨著步長的縮小，在一定求解范圍內(nèi)所要完成的步數(shù)就增加了。步數(shù)的增加不但引起計(jì)在

15、一定求解范圍內(nèi)所要完成的步數(shù)就增加了。步數(shù)的增加不但引起計(jì)算量的增大，而且可能導(dǎo)致舍入誤差的嚴(yán)重積累。算量的增大，而且可能導(dǎo)致舍入誤差的嚴(yán)重積累。第八章常微分方程數(shù)值解法 在選擇步長時，我們需要衡量和檢驗(yàn)計(jì)算結(jié)果的精度，并依據(jù)所獲得的精度在選擇步長時，我們需要衡量和檢驗(yàn)計(jì)算結(jié)果的精度，并依據(jù)所獲得的精度處理步長。下面以經(jīng)典處理步長。下面以經(jīng)典R-K方法為例進(jìn)行說明。方法為例進(jìn)行說明。 從節(jié)點(diǎn)從節(jié)點(diǎn) 出發(fā)，先以出發(fā)，先以h為步長求出一個近似值為步長求出一個近似值 ，由于公式的局部截?cái)?，由于公式的局部截?cái)嗾`差為誤差為 ，故有，故有nx hny1 5hOnx然后將步長折半，既然后將步長折半，既 h/

16、2 為步長，從為步長，從 跨兩步到跨兩步到 ，再求得一個近似，再求得一個近似值值 ，每跨一步的截?cái)嗾`差約為，每跨一步的截?cái)嗾`差約為 ，因此有，因此有1nx 21hnh 52hc 511chyxyhnn第八章常微分方程數(shù)值解法521122hcyxyhnn比較上述二式，有比較上述二式，有。16111211hnnhnnyxyyxy由此易得下列事后估計(jì)式由此易得下列事后估計(jì)式 。hnhnhnnyyyxy121211151這樣，我們可以通過檢查步長折半前后兩次計(jì)算結(jié)果的偏差這樣，我們可以通過檢查步長折半前后兩次計(jì)算結(jié)果的偏差 hnhnyy121來判定所選的步長是否合適。來判定所選的步長是否合適。第八章常微分方程數(shù)值解法 具體地說，對于給定的精度具體地說，對于給定的精度 ，將按兩種情況處理。如果，將按兩種情況處理。如果 ，我，我們反復(fù)將步