Runge-Kutta方法

1、1第三節(jié)第三節(jié) Runge-KuttaRunge-Kutta方法方法2 2009, Henan Polytechnic University23 龍格龍格 - 庫塔法庫塔法 /* Runge-Kutta Method */ 考察改進(jìn)的歐拉法，可以將其改寫為：考察改進(jìn)的歐拉法，可以將其改寫為：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii 斜率斜率一定取一定取K1 K2 的的平均值平均值嗎？嗎？步長一定是一個步長一定是一個h 嗎？嗎？單步遞推法的單步遞推法的基本思想基本思想是從是從 ( xi , yi ) 點出發(fā)，以點出發(fā)，以某一斜某一斜率率沿直線達(dá)到沿直線達(dá)到 (

2、 xi+1 , yi+1 ) 點。歐拉法及其各種變形所點。歐拉法及其各種變形所能達(dá)到的最高精度為能達(dá)到的最高精度為2階階。建立高精度的單步遞推格式。建立高精度的單步遞推格式。3 2009, Henan Polytechnic University3首先希望能確定系數(shù)首先希望能確定系數(shù) 1、 2、p，使得到的算法格式有，使得到的算法格式有2階階精度，即在精度，即在 的前提假設(shè)下，使得的前提假設(shè)下，使得 )(iixyy )()(311hOyxyRiii Step 1: 將將 K2 在在 ( xi , yi ) 點作點作 Taylor 展開展開)(),(),(),(),(2112hOyxfphKyx

3、phfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 將改進(jìn)歐拉法推廣為：將改進(jìn)歐拉法推廣為：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2: 將將 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii 4 2009, Henan Polytechnic University4Step 3: 將將 yi+1 與與 y(

4、 xi+1 ) 在在 xi 點的點的泰勒泰勒展開作比較展開作比較)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，則必須有：，則必須有：)()(311hOyxyRiii21,1221 p 這里有這里有 個未知個未知數(shù)，數(shù)， 個方程。個方程。32存在存在無窮多個解無窮多個解。所有滿足上式的格式統(tǒng)稱為。所有滿足上式的格式統(tǒng)稱為2階龍格階龍格 - 庫庫塔格式塔格式。21, 121 p注意到，注意到， 就是改進(jìn)的歐拉法。就是改進(jìn)的歐拉法。 Q: 為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣

5、？5 2009, Henan Polytechnic University5其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均為待定均為待定系數(shù)，確定這些系數(shù)的系數(shù)，確定這些系數(shù)的步驟與前面相似。步驟與前面相似。 ).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用為四級最常用為四級4階階經(jīng)典龍格經(jīng)典龍格-庫塔法庫塔法 /* Classical

6、Runge-Kutta Method */ ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 6 2009, Henan Polytechnic University6注：注： 龍格龍格-庫塔法庫塔法的主要運算在于計算的主要運算在于計算 Ki 的值，即計算的值，即計算 f 的的值。值。Butcher 于于1965年給出了計算量與可達(dá)到的最高精年給出了計算量與可達(dá)到的最高精度階數(shù)的關(guān)系：度階數(shù)的關(guān)系：753可達(dá)到的最高精度642每步須算Ki 的個數(shù))(2hO)(3hO)(4hO)(5hO)(6

7、hO)(4hO)(2nhO8n 由于龍格由于龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對于光滑性不太好的解，最好解函數(shù)的光滑性影響。對于光滑性不太好的解，最好采用采用低階算法低階算法而將步長而將步長h 取小取小。深入研究龍格深入研究龍格-庫塔法請看庫塔法請看此處此處！7 2009, Henan Polytechnic University77.2 RungeKutta方法 7.2.1 構(gòu)造高階單步法的直接方法 由Taylor公式： 當(dāng)h充分小時，略去Taylor公式余項，并以yi、yi+1分別代替y(xi)、y(xi+1)，得到差分方

8、程： (7.2-1) 其局部截斷誤差為：)()!1()(!.)(! 2)( )()()()1(1)(21ppippiiiiiyphxyphxyhxhyxyhxyxy),(!.),( ! 2),()1(21iippiiiiiiyxfphyxfhyxhfyy)()!1()()1(111ppiiyphyxy8 2009, Henan Polytechnic University8 當(dāng)h充分小時，略去Taylor公式余項，并以yi、yi+1分別代替y(xi)、y(xi+1)，得到差分方程： (7.2-1) 其局部截斷誤差為： 即(7.2-1)為p階方式，上述方式稱為Taylor方式。 注：利用Tayl

9、er公式構(gòu)造，不實用，高階導(dǎo)數(shù)f (i)不易計算。)()!1()(!.)(! 2)( )()()()1(1)(21ppippiiiiiyphxyphxyhxhyxyhxyxy),(!.),( ! 2),()1(21iippiiiiiiyxfphyxfhyxhfyy)()!1()()1(111ppiiyphyxy9 2009, Henan Polytechnic University9 7.2.2 RungeKutta方法 1. 基本思想 因為 = y(xi) + hf (，y() = y(xi) + hK 其中K = f (，y()稱為y(x)在xi，xi+1上的平均斜率。 若取 K1 = f

10、 (xi，y(xi) Euler公式 取 K2 = f (xi+1，y(xi+1) 向后Euler公式 一階精度 取 梯形公式 二階精度 猜想：若能多預(yù)測幾個點的斜率，再取其加權(quán)平均作為K，可望得到較高精度的數(shù)值解，從而避免求f 的高階導(dǎo)數(shù)。1)(,()()(1iixxiidxxyxfxyxy)(2121KK 10 2009, Henan Polytechnic University10 2. RK公式 (7.2-4) 其中Kj為y = y(x)在xi + ajh (0 aj 1)處的斜率預(yù)測值。aj，bjs，cj為特定常數(shù)。111112),(),(jssjsijijiipjjjiipjKbh

11、yhaxfKyxfKKchyy11 2009, Henan Polytechnic University11 3. 常數(shù)的確定 確定的原則是使精度盡可能高。 以二階為例： (7.2-5) (希望y(xi+1) yi+1 = O(hp)的階數(shù)p盡可能高) 首先：),(),()(12122122111hKbyhaxfKyxfKKcKchyyiiiiii)()(! 21)( )()(321hOxyhxhyxyxyiiii12 2009, Henan Polytechnic University12 另一方面： 將K2在(xi，yi)處展開。K2 = f (xi，yi) + a2hf x(xi，yi)

12、 + b21hK1 f y(xi，yi) + O(h2). 代入(7.2-5)得：yi+1 = yi + hc1 f (xi，yi) + hc2 f (xi，yi) + h2c2a2 f x(xi，yi) + b21K1 f y(xi，yi) + O(h3) = yi + h(c1 + c2) f (xi，yi) + c2a2h2f x(xi，yi) + (b21/a12) f (xi，yi) f y(xi，yi)+O(h3) （希望）.),()(),(),(000000yxfyyxxyxfyyxxf)()()( )(322221hOxyhacxycchyiii13 2009, Henan P

13、olytechnic University13 希望：ei+1 = y(xi+1) yi+1 = O(h3). 則應(yīng)： 特例：a2 = 1 c1 = c2 = 1/2，b21 = 1，得2階R-K公式： 改進(jìn)歐拉公式。12112212221abaccc),(),()(2121211hKyhxfKyxfKKKhyyiiiiii14 2009, Henan Polytechnic University14 希望：ei+1 = y(xi+1) yi+1 = O(h3). 則應(yīng)： 特例：c1 = 0 c2 = 1，a2 = 1/2，b21 = 1/2，得： (7.2-7) 稱為中點公式。1211221

14、2221abaccc)2,2(),(12121KhyhxfKyxfKhKyyiiiiii15 2009, Henan Polytechnic University15 4. 最常用的R-K公式 標(biāo)準(zhǔn)4階R-K公式 (7.2-8) 算法：),()2,2()2,2(),()22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyyiiiiiiiiii輸入輸入a，b，n，y0h=(b-a)n，x0 = afor i = 1, i=n, i+K1 = f(x0, y0)K2 = f(x0+h/2, y0+h*K1/2)K3 = f(x0+h/2, y0+h*K2/

15、2)K4 = f(x0+h, y0+h*K3)x0 = x0+hy0 = y0 + h*(K1+2*K2+2*K3+K4)/6輸出輸出x0，y016 2009, Henan Polytechnic University16 Matlab代碼： function Runge_Kutta4(a,b,h,y0) n=(b-a)/h;x0=a; for i=1:n K1=f(x0,y0) K2=f(x0+h/2,y0+h*K1/2) K3=f(x0+h/2,y0+h*K2/2) K4=f(x0+h,y0+h*K3) x0=x0+h y1=y0+h*(K1+2*K2+2*K3+K4)/6; y0=y1 end; end; function f=f(x,y) f=x+y; end;輸入輸入a，b，n，y0h=(b-a)n，x0 = afor i = 1, i 時，說明步長h