数学 - 微分方程数值解 - 第 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 方法