Runge-Kutta方法

1、1第三节第三节 Runge-KuttaRunge-Kutta方法方法2 2009, Henan Polytechnic University23 龙格龙格 - 库塔法库塔法 /* Runge-Kutta Method */ 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii 斜率斜率一定取一定取K1 K2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，以点出发，以某一斜某一斜率率沿直线达到沿直线达到 (

2、 xi+1 , yi+1 ) 点。欧拉法及其各种变形所点。欧拉法及其各种变形所能达到的最高精度为能达到的最高精度为2阶阶。建立高精度的单步递推格式。建立高精度的单步递推格式。3 2009, Henan Polytechnic University3首先希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 )(iixyy )()(311hOyxyRiii Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyx

3、phfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 将改进欧拉法推广为：将改进欧拉法推广为：),(),(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 这里有这里有 个未知个未知数，数， 个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格 - 库库塔格式塔格式。21, 121 p注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广

5、？5 2009, Henan Polytechnic University5其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为待定均为待定系数，确定这些系数的系数，确定这些系数的步骤与前面相似。步骤与前面相似。 ).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法 /* Classical

6、Runge-Kutta Method */ ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 6 2009, Henan Polytechnic University6注：注： 龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度642每步须算Ki 的个数)(2hO)(3hO)(4hO)(5hO)(6

7、hO)(4hO)(2nhO8n 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。深入研究龙格深入研究龙格-库塔法请看库塔法请看此处此处！7 2009, Henan Polytechnic University77.2 RungeKutta方法 7.2.1 构造高阶单步法的直接方法 由Taylor公式： 当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 当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程： (7.2-1) 其局部截断误差为： 即(7.2-1)为p阶方式，上述方式称为Taylor方式。 注：利用Tayl

9、er公式构造，不实用，高阶导数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公式 一阶精度 取 梯形公式 二阶精度 猜想：若能多预测几个点的斜率，再取其加权平均作为K，可望得到较高精度的数值解，从而避免求f 的高阶导数。1)(,()()(1iixxiidxxyxfxyxy)(2121KK 10 2009, Henan Polytechnic University10 2. RK公式 (7.2-4) 其中Kj为y = y(x)在xi + ajh (0 aj 1)处的斜率预测值。aj，bjs，cj为特定常数。111112),(),(jssjsijijiipjjjiipjKbh

11、yhaxfKyxfKKchyy11 2009, Henan Polytechnic University11 3. 常数的确定 确定的原则是使精度尽可能高。 以二阶为例： (7.2-5) (希望y(xi+1) yi+1 = O(hp)的阶数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). 则应： 特例：a2 = 1 c1 = c2 = 1/2，b21 = 1，得2阶R-K公式： 改进欧拉公式。12112212221abaccc),(),()(2121211hKyhxfKyxfKKKhyyiiiiii14 2009, Henan Polytechnic University14 希望：ei+1 = y(xi+1) yi+1 = O(h3). 则应： 特例：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公式 标准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