龙格－库塔（Runge-Kutta methods）

龙格－库塔（Runge-Kutta methods）是用于非线性常微分方程的解的重要的一类迭代法。 由数学家卡尔·龙格和马丁·威尔海姆·库塔于 1900 年左右发明。

• 原理
• 显式龙格-库塔
  • 四阶
    • 经典龙格-库塔方法
    • 基尔方法
    • 其他
  • 高阶
• 隐式
  • 四阶龙格-库塔
  • 隐式中点法
  • Hammer-Hollingsworth 方法
  • 一种三级六阶的隐式龙格库塔方法
• 变步长龙格-库塔
• 单步法的绝对稳定性
• 代码实现

原理

龙格-库塔方法实质上是简介的使用泰勒级数法的一种技术。

考虑插商 $\displaystyle \frac{y(x_{n+1}-y(x_n))}{h}$，根据微分中值定理，存在 $0<\theta<1$，使得

$$\frac{y(x_{n+1}-y(x_n))}{h} = y'(x_n + \theta h)$$

于是，利用所给方程 $y'=f(x, y)$ 得到

$$y(x_{n+1}) = y(x_n) + hf(x_n + \theta h, y(x_n + \theta h))$$

这里 $K^* = f(x_n + \theta h, y(x_n + \theta h))$ 称作区间 $[x_n,x_{n+1}]$ 上的平均斜率。

如果设法在区间 $[x_n,x_{n+1}]$ 内多预测几个点的斜率值，然后将他们加权平均作为平均斜率 $K^*$，则有可能构造出更高精度的计算公式，这就是龙格-库塔的基本思想。

显式龙格-库塔

四阶

经典龙格-库塔方法

在各种 Runge-Kutta 法当中有一个方法十分常用，以至于经常被称为 “RK4” （四阶四级）。

令初值问题表述如下。

$$y'=f(t,y),\quad y(t_{0})=y_{0}$$

则，对于该问题的 RK4 由如下方程给出：

$$y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$

其中

$$\begin{aligned} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_1) \\ k_3 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_2) \\ k_4 &= f(t_n + h, y_n + hk_3) \end{aligned}$$

这样，下一个值（$y_n+1$）由现在的值（$y_n$）加上时间间隔（$h$）和一个估算的斜率的乘积所决定。 该斜率是以下斜率的加权平均：

• $k_1$ 是时间段开始时的斜率；
• $k_2$ 是时间段中点的斜率，通过欧拉法采用斜率 $k_1$ 来决定 $y$ 在点 $t_n + \frac{h}{2}$ 的值；
• $k_3$ 也是中点的斜率，但是这次采用斜率 $k_2$ 决定 $y$ 值；
• $k_4$ 是时间段终点的斜率，其 $y$ 值用 $k_3$ 决定。

当四个斜率取平均时，中点的斜率有更大的权值：

$$slope = \frac{k_1 + 2k_2 + 2k_3 + k_4}{6}.$$

RK4 法是四阶方法，也就是说每步的误差是 $h^5$ 阶，而总积累误差为 $h^4$ 阶。

注意上述公式对于标量或者向量函数（$y$ 可以是向量）都适用。

::: tips 值得注意的是龙格-库塔方法基于泰勒展开，因而他要求所具有的解具有较好的光滑性质。 反之，如果解的光滑性较差，那么使用龙格-库塔得到的解，其精度可能反而不如低阶方法。 因此应结合具体问题的特点选择合适的算法。 :::

基尔方法

四阶四级。

$$y_{n+1} = y_n + \frac{h}{6}(k_1 + (2-\sqrt{2})k_2 + (2+\sqrt{2})k_3 + k_4)$$

其中

$$\begin{aligned} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_1) \\ k_3 &= f(t_n + \frac{h}{2}, y_n + \frac{\sqrt{2}-1}{2}hk_1+(1-\frac{\sqrt{2}}{2})hk_2) \\ k_4 &= f(t_n + h, y_n - \frac{\sqrt{2}}{2}hk_2 + (1+\frac{\sqrt{2}}{2})hk_3) \end{aligned}$$

其他

四阶四级。

$$y_{n+1} = y_n + \frac{h}{6}(k_1 + 3k_2 + 3k_3 + k_4)$$

其中

$$\begin{aligned} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + \frac{h}{3}, y_n + \frac{h}{3}k_1) \\ k_3 &= f(t_n + \frac{2}{3}h, y_n - \frac{1}{3}h(k_1-3k_2)) \\ k_4 &= f(t_n + h, y_n + h(k_1-k_2+k_3)) \end{aligned}$$

高阶

高阶相对于四阶而言，多了更多的计算量，却没有相对来说性价比较高的精度进步，故而一般还是采用四阶龙格-库塔，例如，6 级龙格-库塔拥有 5 阶精度，8 级龙格库塔拥有 6 级精度。

隐式

四阶龙格-库塔

显式 Runge-Kutta 法一般来讲不适用于求解刚性方程。 这是因为显式 Runge-Kutta 方法的稳定区域被局限在一个特定的区域里。 显式 Runge-Kutta 方法的这种缺陷使得人们开始研究隐式 Runge-Kutta 方法，一般而言，隐式 Runge-Kutta 方法具有以下形式：

$$y_{n+1} = y_n + h \sum_{i=1}^{s} b_i k_i$$

其中，

$$k_i = f(t_n + c_i h, y_n + h \sum_{j=1}^{s} a_{ij} k_j), \qquad i = 1, \dotsb, s.$$

在显式 Runge-Kutta 方法的框架里，定义参数 ${a_{ij}}$ 的矩阵是一个下三角矩阵，而隐式 Runge-Kutta 方法并没有这个性质，这是两个方法最直观的区别。

需要注意的是，与显式 Runge-Kutta 方法不同，隐式 Runge-Kutta 方法在每一步的计算里需要求解一个线性方程组，这相应的增加了计算的成本。

隐式中点法

该方法为一级二阶方法。

$$\begin{cases} y_{n+1} = y_n + hK_1 \\ K_1 = f(x_n+\frac{h}{2}, y_n + \frac{h}{2} K_1) \end{cases}$$

Hammer-Hollingsworth 方法

这是一种二级四阶的方法，其系数表有

$\frac{1}{2}-\frac{\sqrt{3}}{6}$	$\frac{1}{4}$	$\frac{1}{4}-\frac{\sqrt{3}}{6}$
$\frac{1}{2}+\frac{\sqrt{3}}{6}$	$\frac{1}{4}+\frac{\sqrt{3}}{6}$	$\frac{1}{4}$
	$\frac{1}{2}$	$\frac{1}{2}$

一种三级六阶的隐式龙格库塔方法

其系数表为

$\frac{1}{10}(5-\sqrt{15})$	$\frac{5}{36}$	$\frac{1}{45}(10-3\sqrt{15})$	$\frac{1}{180}(25-6\sqrt{15})$
$\frac{1}{2}$	$\frac{1}{72}(10+3\sqrt{15})$	$\frac{2}{9}$	$\frac{1}{72}(10-3\sqrt{15})$
$\frac{1}{10}(5+\sqrt{15})$	$\frac{1}{180}(25+6\sqrt{15})$	$\frac{1}{45}(10+3\sqrt{15})$	$\frac{5}{36}$
	$\frac{5}{18}$	$\frac{4}{9}$	$\frac{5}{18}$

变步长龙格-库塔

步长越小，截断误差也越小，但随着步长的缩小，在一定范围内所要完成的步数就增加了，步数的增加不但引起了计算量的增大，而且可能导致舍入误差的严重积累，因此同积分的数值计算一样，微分方程的数值解法也有一个选者步长的问题。

这就引入了两个问题：

1. 怎样衡量和检验计算结果的精度？
2. 如何依据的精度处理步长？

从经典的四阶龙格-库塔方法出发，由于公式的局部误差为 $O(h^5)$，故有

$$y(x_{n+1}) - y_{n+1}^{(h)} \approx ch^5 \tag{1}$$

然后将步长折半， 即步长为 $\displaystyle \frac{h}{2}$，从 $x_n$ 跨两步到 $x_{n+1}$，再求得一个近似值 $y_{n+1}^{\frac{h}{2}}$，每跨一步的截断误差为 $\displaystyle c(\frac{h}{2})^5$，因此有

$$y(x_{n+1}) - y_{n+1}^{(\frac{h}{2})} = 2c(\frac{h}{2})^5 \tag{2}$$

比较 $(1)$ 和 $(2)$ 式，我们可以看到，步长折半以后误差大约只有原误差的 $\displaystyle \frac{1}{16}$，既有

$$\frac{y(x_{n+1}) - y_{n+1}^{(\frac{h}{2})}}{y(x_{n+1}) - y_{n+1}^{(h)}} \approx \frac{1}{16}$$

于是得到事后估计式：

$$y(x_{n+1}) - y_{n+1}^{(\frac{h}{2})} \approx \frac{1}{15}[y_{n+1}^{(\frac{h}{2})} - y_{n+1}^{(h)}]$$

随后我们检查步长，折半前后的两次计算结果的偏差：

$$\Delta = |y_{n+1}^{(\frac{h}{2})} - y_{n+1}^{(h)}|$$

来判断所选的步长是否合适，即：

对于给定的精度 $\varepsilon$，

1. 如果 $\Delta > \varepsilon$，我们将反复将步长折半进行计算，直到 $\Delta < \varepsilon$ 为止，这时取最终的 $y_{n+1}^{(h)}$ 为结果。
2. 如果 $\Delta < \varepsilon$，我们将反复将步长加倍进行计算，直到 $\Delta > \varepsilon$ 为止，这时在将步长折半一次， 这时取最终的 $y_{n+1}^{(h)}$ 为结果。

当然，必不可少的引入了更多的计算量，但总体考虑有可能是合算的。

单步法的绝对稳定性

可将方程局部线性化为

$$\frac{dy}{dx} = \lambda y$$

从 $y_n$ 计算一步得

$$y_{n+1} = E(\lambda h) y_n$$

$E(\lambda h)$ 依赖于所选的方法，而且是 $e^{\lambda h}$ 的一个近似值。 如果在 $y_n$ 的计算有误差 $\varepsilon$ 则会引起 $y_{n+1}$ 的误差 $E(\lambda h) \varepsilon$。 如果在 $y_0$ 的计算有误差 $\varepsilon_0$ 则会引起 $y_{n+1}$ 的误差 $(E(\lambda h))^{n+1} \varepsilon_0$。

若某方法满足 $|E(\lambda h) < 1|$，称此方法为绝对稳定；在 $\mu$ 复平面上复变量 $\mu$ 满足 $|E(\mu)| < 1$ 的区域称之为绝对稳定区域，它与实轴的交称之为绝对稳定区间。

例子：

1. 将欧拉法用于 $\displaystyle \frac{dy}{dx} = \lambda y$，有

$$y_{n+1} = y_n + \lambda h y_n = (1+\mu) y_n$$

其中 $\mu = \lambda h, E(\mu) = 1 + \mu$。 $E(\mu)$ 是 $e^\mu$ 的一阶泰勒近似，当 $|1+\mu| < 1$ 时，欧拉法是绝对稳定的，即欧拉法是一个以 $1-$ 为中心的单位圆，绝对稳定区间为 $(-2, 0)$。

2. 将四阶经典龙格-库塔用于 $\displaystyle \frac{dy}{dx} = \lambda y$，有

$$y_{n+1} = 1 + \lambda h + \frac{(\lambda h)^2}{2} + \frac{(\lambda h)^3}{6} + \frac{(\lambda h)^4}{24} = e^{i\theta}$$

有 $\displaystyle E(\mu) = 1 + \mu + \frac{\mu^2}{2} + \frac{\mu^3}{6} + \frac{\mu^4}{24}$ 为 $e^{\mu}$ 的四阶泰勒近似。

代码实现

TODO