《数值分析Runge-Kutta方法解析》由会员分享，可在线阅读，更多相关《数值分析Runge-Kutta方法解析（24页珍藏版）》请在金锄头文库上搜索。
1、1 第三节 Runge-Kutta方法 2 2009, Henan Polytechnic University 23 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 3 龙格 - 库塔法 /* Runge-Kutta Method */ 考察改进的欧拉法，可以将其改写为： 斜率 一定取K1 K2 的平均值吗？ 步长一定是一个h 吗 ？ 单步递推法的基本思想是从 ( xi , yi ) 点出发，以某一斜 率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所 能达到的最高精度为2阶。 建立高精度的单步递推格式。 3
2、2009, Henan Polytechnic University 33 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 首先希望能确定系数 1、2、p，使得到的算法格式有2阶 精度，即在 的前提假设下，使得 Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开 将改进欧拉法推广为： ),( ),( 12 1 22111 phKyphxfK yxfK KKhyy ii ii ii += = += + Step 2: 将 K2 代入第1式，得到 4 2009, Henan Polytechnic Unive
3、rsity 43 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较 要求 ，则必须有： 这里有 个未知 数， 个方程。 3 2 存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库 塔格式。注意到， 就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？ 5 2009, Henan Polytechnic University 53 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解
4、法 其中i ( i = 1, , m )，i ( i = 2, , m ) 和 ij ( i = 2, , m; j = 1, , i1 ) 均为待定 系数，确定这些系数的 步骤与前面相似。 ).,( . ),( ),( ),( . 112211 23213133 12122 1 22111 + += += += = += mm mmmmim ii ii ii mmii hKhKhKyhxfK hKhKyhxfK hKyhxfK yxfK KKKhyy 最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ ： 6 2009, Henan Pol
5、ytechnic University 63 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 注： 龙格-库塔法的主要运算在于计算 Ki 的值，即计算 f 的 值。Butcher 于1965年给出了计算量与可达到的最高精 度阶数的关系： 753 可达到的最高精度 642 每步须算Ki 的个数 由于龙格-库塔法的导出基于泰勒展开，故精度主要受 解函数的光滑性影响。对于光滑性不太好的解，最好 采用低阶算法而将步长h 取小。 深入研究龙格- 库塔法请看此处！ 7 2009, Henan Polytechnic University 73
6、3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 7.2 RungeKutta方法 7.2.1 构造高阶单步法的直接方法 由Taylor公式： 当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)， 得到差分方程： (7.2-1) 其局部截断误差为： 8 2009, Henan Polytechnic University 83 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 当h充分小时，略去Taylor公式余项，并以yi、yi+1
7、分别代替y(xi)、y(xi+1)， 得到差分方程： (7.2-1) 其局部截断误差为： 即(7.2-1)为p阶方式，上述方式称为Taylor方式。 注：利用Tayler公式构造，不实用，高阶导数f (i)不易计算。 9 2009, Henan Polytechnic University 93 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 7.2.2 RungeKutta方法 1. 基本思想 因为 = y(xi) + hf (，y() = y(xi) + hK 其中K = f (，y()称为y(x)在xi，xi+1上的平均斜率。
8、 若取 K1 = f (xi，y(xi) Euler公式 取 K2 = f (xi+1，y(xi+1) 向后Euler公式 一阶精度 取 梯形公式 二阶精度 猜想：若能多预测几个点的斜率，再取其加权平均作为K，可望得到较高精 度的数值解，从而避免求f 的高阶导数。 10 2009, Henan Polytechnic University 103 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 2. RK公式 (7.2-4) 其中Kj为y = y(x)在xi + ajh (0 aj 1)处的斜率预测值。 aj，bjs，cj为特定常数
9、。 11 2009, Henan Polytechnic University 113 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 3. 常数的确定 确定的原则是使精度尽可能高。 以二阶为例： (7.2-5) (希望y(xi+1) yi+1 = O(hp)的阶数p尽可能高) 首先： 12 2009, Henan Polytechnic University 123 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 另一方面： 将K2在(xi，yi)处展开。 K2 = f
10、 (xi，yi) + a2hf x(xi，yi) + 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) （希望） 13 2009, Henan Polytechnic University 133 3 Runge-Kutt
11、aRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 希望：ei+1 = y(xi+1) yi+1 = O(h3). 则应： 特例：a2 = 1 c1 = c2 = 1/2，b21 = 1，得2阶R-K公式： 改进欧拉公式。 14 2009, Henan Polytechnic University 143 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 希望：ei+1 = y(xi+1) yi+1 = O(h3). 则应： 特例：c1 = 0 c2 = 1，a2 = 1/2，b21 = 1/2，得：
12、 (7.2-7) 称为中点公式。 15 2009, Henan Polytechnic University 153 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 4. 最常用的R-K公式 标准4阶R-K公式 (7.2-8) 算法： 输入a，b，n，y0 h=(b-a)n，x0 = a for 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/2) K4 = f(x0+h, y0+h*K3) x0 = x0+h y0 = y
13、0 + h*(K1+2*K2+2*K3+K4)/6 输出x0，y0 16 2009, Henan Polytechnic University 163 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 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，y0 h=(b-a)n，x0 = a for i = 1, i 时，说明步长h / 2仍然偏大，须将步长减半，继续计算； 当 / 2p时，说明已有 ，步长h / 2偏小，应取 步长h。 24 2009, Henan Polytechnic University 243 3 Runge-KuttaRunge-Kutta方法方法 第五章第五章 常微分方程数值解法常微分方程数值解法 例2(P481) 用变步长的标准4阶R-K方法求初值问题： 的数值解，要求精度为
