§53Runge-Kutta方法技巧.ppt

第三节 Runge-Kutta方法 7.2 Runge?Kutta方法 7.2.1 构造高阶单步法的直接方法 由Taylor公式： 当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程： (7.2-1) 其局部截断误差为： 当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程： (7.2-1) 其局部截断误差为： 即(7.2-1)为p阶方式，上述方式称为Taylor方式。 注：利用Tayler公式构造，不实用，高阶导数f (i)不易计算。 7.2.2 Runge?Kutta方法 1. 基本思想 因为 = y(xi) + hf (?，y(?)) = y(xi) + hK? 其中K? = f (?，y(?))称为y(x)在[xi，xi+1]上的平均斜率。 若取 K1 = f (xi，y(xi)) ——Euler公式 取 K2 = f (xi+1，y(xi+1)) —— 向后Euler公式 一阶精度 取 —— 梯形公式 二阶精度 猜想：若能多预测几个点的斜率，再取其加权平均作为K?，可望得到较高精度的数值解，从而避免求f 的高阶导数。 2. R?K公式 (7.2-4) 其中Kj为y = y(x)在xi + ajh (0 ? aj ? 1)处的斜率预测值。 aj，bjs，cj为特定常数。 3. 常数的确定 确定的原则是使精度尽可能高。 以二阶为例： (7.2-5) (希望y(xi+1) – yi+1 = O(hp)的阶数p尽可能高) 首先： 另一方面： 将K2在(xi，yi)处展开。 K2 = f (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) + h2c2[a2 f x(xi，yi) + b21K1 f y(xi，yi)] + O(h3) = yi + h(c1 + c2) f (xi，yi) + c2a2h2[f x(xi，yi) + (b21/a12) f (xi，yi) f y(xi，yi)]+O(h3) （希望） 希望：ei+1 = y(xi+1) – yi+1 = O(h3). 则应： 特例：a2 = 1 ? c1 = c2 = 1/2，b21 = 1，得2阶R-K公式： ??改进欧拉公式。 希望：ei+1 = y(xi+1) – yi+1 = O(h3). 则应： 特例：c1 = 0 ? c2 = 1，a2 = 1/2，b21 = 1/2，得： (7.2-7) 称为中点公式。 4. 最常用的R-K公式 —— 标准4阶R-K公式 (7.2-8) 算法： 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; 例1：（P478）用标准4阶R-K公式求： 的数值解。取h = 0.2，并与标准解y = 2ex – x – 1比较。 解：因为f(x，y) = x + y，从而由(7.2-8)得： 例1：（P478）用标准4阶R-K公式求： 的数值解。取h = 0.2，并与标准解y = 2ex – x – 1比较。 解： 注：步长h的选择 (P479) 使用数值解法求解初值问题，选择步长h是一个重要问题。 从每一步看，h小局部截断误差小，整体截断误差也就小； 但从整个区间看，h小则节点多