数学 - 微分方程数值解 - 第 1 章 一阶常微分方程初值问题 - 1.3 龙格法

1.3 Runge - Kutta 法

1.3.1 RK 方法的构造

一般地，RK 方法设近似公式为：

\[\left\{ \begin{align*} & y_{n+1} = y_n + h \sum_{i=1}^{p} c_i K_i \\ & K_1 = f(x_n, y_n) \\ & K_i = f(x_n + h a_i, y_n + h \sum_{j=1}^{i-1} b_{ij} K_j) \quad (i=1,2,\cdots,p) \end{align*} \right. \tag{1.3.1} \]

上式中 \(a_i\)、\(b_{ij}\)、\(c_i\) 都是参数。确定它们的原则是使近似公式在 \(x_n\) 处的 Taylor 展开式与 \(y(x)\) 在 \(x_n\) 处的 Taylor 展开式的前面的项尽可能多地重合，使之具有尽可能高阶的局部截断误差。

1.3.1 二阶 RK 公式

当 \(p=2\) 时，近似公式为：

\[\left\{ \begin{align*} & y_{n+1} = y_n + h（c_1 K_1 + c_2 K_2） \\ & K_1 = f(x_n, y_n) \\ & K_2 = f(x_n + h a_2, y_n + h b_{21} K_1) \end{align*} \right. \tag{1.3.2} \]

计算局部截断误差时，为方便书写，此处令 \(y_{n+1}=y(x_{n+1})\) 和 \(y_{n}=y(x_{n})\)。上式在 \((x_n, y_n)\) 处的 Taylor 展开式为

\[\begin{align*} y_{n+1} & = y_n + h \left[ c_1 f(x_n, y_n) + c_2 f(x_n + h a_2, y_n + h b_{21} K_1) \right] \\ & = y_n + h \left\{ c_1 f(x_n, y_n) + c_2 \left[ f(x_n, y_n) + a_2 h f^{'}_{x}(x_n, y_n) + b_{21} h f^{'}_{y}(x_n, y_n) f(x_n, y_n)\right] \right\} + O(h^3) \\ & = y_n + (c_1 + c_2) f(x_n, y_n) h + c_2 \left[ a_2 h f^{'}_{x}(x_n, y_n) + b_{21} h f^{'}_{y}(x_n, y_n) f(x_n, y_n)\right] h^2 + O(h^3) \\ \end{align*} \tag{1.3.3} \]

\(y(x_{n+1})\) 在 \(x_n\) 处的 Taylor 展开式为：

\[\begin{align*} y(x_{n+1}) & = y(x_n) + h y'(x_n) + \frac{h^2}{2} y''(x_n) + O(h^3) \\ & = y_n + f(x_n, y_n) h + \frac{h^2}{2} \left[ f^{'}_{x}(x_n, y_n) + f^{'}_{y}(x_n, y_n) f(x_n, y_n)\right] + O(h^3) \end{align*} \tag{1.3.4} \]

注意，式 \((1.3.3)\) 和式 \((1.3.4)\) 中的 \(f^{'}_{x}\) 表示函数 \(f(x,y)\) 对 \(x\) 的一阶导数。

要使近似公式 \((1.3.2)\) 的局部截断误差为 \(O(h^3)\)，则应要求式 \((1.3.3)\) 和式 \((1.3.4)\) 的前三项相同，于是有

\[\left\{ \begin{align*} & c_1 + c_2 = 1 \\ & c_2 a_2 = \frac{1}{2} \\ & c_2 b_{21} = \frac{1}{2} \end{align*} \right. \tag{1.3.5} \]

式 \((1.3.5)\) 有无穷多组解，它的每一组解代入式 \((1.3.2)\) 中得到的近似公式，局部截断误差均为 \(O(h^3)\)，故这些方法都为二阶方法。

如果我们取 \(c_1 = c_2 = 1/2\)，\(a_2 = b_{21} = 1\)，可得改进 Euler 公式：

\[\left\{ \begin{align*} & y_{n+1} = y_n + \frac{h}{2} (K_1 + K_2) \\ & K_1 = f(x_n, y_n) \\ & K_2 = f(x_n + h, y_n + h K_1) \end{align*} \right. \tag{1.3.6} \]

如果我们取 \(c_1 = 0\)，\(c_2 = 1\)，\(a_2 = b_{21} = 1/2\)，可得中点公式：

\[\left\{ \begin{align*} & y_{n+1} = y_n + h K_2 \\ & K_1 = f(x_n, y_n) \\ & K_2 = f(x_n + h / 2, y_n + K_1 h / 2) \end{align*} \right. \tag{1.3.7} \]

理论上可以证明，无论怎样选取参数，式 \((1.3.2)\) 都不可能具有更高的精度，最多只能得到二阶公式。

1.3.2 三阶 RK 公式和四阶 RK 公式

类似地，对 \(p=3\) 和 \(p=4\) 可以类似地推导，得到三阶和四阶 RK 公式，其中常用的三阶 RK 公式为

\[\left\{ \begin{align*} & y_{n+1} = y_n + \frac{h}{6} (K_1 + 4 K_2 + K_3) \\ & K_1 = f(x_n, y_n) \\ & K_2 = f(x_n + \frac{h}{2}, y_n + \frac{h}{2} K_1) \\ & K_3 = f(x_n + h, y_n - h K_1 + 2 h K_2) \end{align*} \right. \tag{1.3.8} \]

常用的四阶 RK 公式为

\[\left\{ \begin{align*} & y_{n+1} = y_n + \frac{h}{6} (K_1 + 2 K_2 + 2 K_3 + K_4) \\ & K_1 = f(x_n, y_n) \\ & K_2 = f(x_n + \frac{h}{2}, y_n + \frac{h}{2} K_1) \\ & K_3 = f(x_n + \frac{h}{2}, y_n + \frac{h}{2} K_2) \\ & K_4 = f(x_n + h, y_n + h K_3) \end{align*} \right. \tag{1.3.9} \]

式 \((1.3.9)\) 也被称为经典形式四阶 RK 公式。

1.3.3 RK 公式说明

当 \(p=1,2,3,4\) 时，RK 公式的局部截断误差的最高阶数恰好是 \(p\)，当 \(p>4\) 时，RK 公式的局部截断误差的最高阶数不是 \(p\)，如 \(p=5\) 时，RK 公式的最高阶数仍为 \(4\)，如 \(p=6\) 时，RK 公式的最高阶数为 \(5\)。

值得注意的是，RK 方法的导出基于 Taylor 展开，故它要求所求问题的解具有较高的光滑度。当解充分光滑时，四阶 RK 方法优于改进 Euler 方法；如果解的光滑性较差，则用四阶 RK 方法求数值解的效果可能不如改进 Euler 方法。

1.3.4 变步长 RK 方法