030 The Runge-Kutta Method 龙格-库塔法

在第 8.1 和 8.2 节中，我们讨论了欧拉公式、后向欧拉公式和改进的欧拉公式，作为数值近似求解初值问题 $$y'=f(t,y),\quad y(t_0)=y_0\tag{1}$$ 的方法。这些方法的局部截断误差分别正比于 $h^2,h^2,h^3$ 成正比。欧拉法和改进欧拉法都属于现在所称的龙格-库塔类方法。

在本节中，要讨论最初由龙格（Runge）和库塔（Kutta）发明的方法。该方法现在被称为经典的四阶四段龙格-库塔法（fourth-order four-stage Runge-Kutta method），简便起见，通常将其简称为龙格-库塔法。该方法的局部截断误差正比于 $h^5$。因此它的精度比改进欧拉法高出两个数量级，比普通欧拉法高出三个数量级。它使用起来相对简单并且精度足以高效地处理许多问题。对于自适应龙格-库塔法（adaptive Runge-Kutta methods）来说更是如此，可以根据需要调整步长。本节末尾会讨论这一话题。

龙格-库塔法涉及 $t_n\leq t\leq t_{n+1}$ 区间上四个不同点处 $f(t,y)$ 的加权平均值。具体公式是 $$y_{n+1}=y_n+h(\frac{k_{n1}+2k_{n2}+2k_{n3}+k_{n4}}{6})\tag{2}$$ 其中 $$\begin{aligned} k_{n1}&=f(t_n,y_n)\\ k_{n2}&=f(t_n+\frac{1}{h},y_n+\frac{1}{2}hk_{n1})\\ k_{n3}&=f(t_n+\frac{1}{h},y_n+\frac{1}{2}hk_{n2})\\ k_{n4}&=f(t_n+h,y_n+hk_{n3}) \end{aligned}\tag{3}$$ $\frac{1}{6}(k_{n1}+2k_{n2}+2k_{n3}+k_{n4})$ 可以解释为平均斜率。$k_{n1}$ 是区间左端点处的斜率；$k_{n2}$ 是中点处的斜率，它是利用欧拉公式从 $t_n$ 到 $t_n+h/2$ 得到的；$k_{n3}$ 是对中点处斜率的第二次近似；而 $k_{n4}$ 则是 $t_n+h$ 处的斜率，是用欧拉公式以及斜率 $k_{n3}$ 从 $t_n$ 到 $t_n+h$ 得到的。

虽然从原理上证明公式 $(2)$ 与解 $\phi$ 的泰勒展开式之间仅存在与 $h^5$ 成正比并不困难，但其代数推导过程相当冗长。因此这里直接给出结论而忽略证明。使用公式 $(2)$ 的局部截断误差正比于 $h^5$，在有限区间内，其全局截断误差最多是一个常数乘以 $h^4$。此前将该方法称为四阶四段法，正是以为其全局截断误差相对于步长 $h$ 是四阶的，并且计算过程中包含四个中间阶段（即 $k_{n1},\dots,k_{n4}$ 的计算）。

显然，公式 $(2)$ 和 $(3)$ 的龙格-库塔公式更为复杂。但在实际中并不重要，因为写计算机程序来实现该方法并不困难。

如果 $f$ 不依赖于 $y$，那么 $$k_{n1}=f(t_n),k_{n2}=k_{n3}=f\left(t_n+\frac{h}{2}\right),k_{n4}=f(t_n+h)\tag{4}$$ 公式 $(2)$ 化简为 $$y_{n+1}-y_n=\frac{h}{6}\left[f(t_n)+4f\left(t_n+\frac{h}{2}\right)+f(t_n+h)\right]\tag{5}$$ 公式 $(5)$ 可以用于近似计算 $y'=f(t)$ 的积分，称为辛普森法则（Simpson's rule）。辛普森法则的误差正比于 $h^5$，与龙格-库塔公式的局部截断误差是一致的。

例 1 使用龙格-库塔法求解初值问题 $$y'=1-t+4y,\quad y(0)=1\tag{6}$$ 的解 $y=\phi(t)$ 的近似值。

解：取步长 $h=0.2$，有 $$\begin{aligned} k_{01}&=f(0,1)=5&&\frac{1}{2}hk_{01}=0.5\\ k_{02}&=f(0+0.1,1+0.5)=6.9&&\frac{1}{2}hk_{02}=0.69\\ k_{03}&=f(0+0.1,1+0.69)=7.66&&hk_{03}=1.532\\ k_{04}&=f(0+0.2,1+1.532)=10.928 \end{aligned}$$ 那么 $$\begin{aligned} y_1&=1+\frac{0.2}{6}\left(5+2(6.9)+2(7.66)+10.928\right)\\ &=1+1.5016=2.5016 \end{aligned}$$ 下表的三、四、五列是使用步长 $h=0.2,h=0.1,h=0.05$ 的龙格-库塔法近似结果。第二列是改进的欧拉法的近似结果。当步长 $h=0.1$ 时，龙格-库塔法在 $t=2$ 处的近似值与精确值的误差仅为 0.122%，而当 $h=0.05$ 时，误差仅为 0.00903%。对于后者，误差小于万分之一，计算出的 $t=2$ 处的值前四位都是准确的。

$t$	$h=0.025$	$h=0.2$	$h=0.1$	$h=0.05$	$Exact$
0.0	1.0000000	1.0000000	1.0000000	1.0000000	1.0000000
0.1		1.6079462	1.6089333	1.6090338	1.6090418
0.2	2.5020618	2.5016000	2.5050062	2.5053060	2.5053299
0.3		3.8228282	3.8294145	3.8300854	3.8301388
0.4	5.7796888	5.7776358	5.7927853	5.7941197	5.7942260
0.5		8.6849039	8.7093175	8.7118060	8.7120041
1.0	64.497931	64.441579	64.858107	64.894875	64.897803
1.5		474.83402	478.81928	479.22674	479.25919
2.0	3496.6702	3490.5574	3535.8667	3539.8804	3540.2001

下面比较这些方法的计算量，使用 $h=0.05$ 的龙格-库塔法和使用 $h=0.025$ 的改进欧拉法，计算到 $t=2$ 都需要计算 160 次 $f$。然而，改进欧拉法在 $t=2$ 处的结果误差为 1.23%。虽然在某些情况下这个误差是可以接受的，但在相同的计算工作量下，它的误差是龙格-库塔法的 135 倍以上。另外可以看到，使用 $h=0.2$ 的龙格-库塔法（仅需计算 40 次 $f$），在 $t=2$ 处产生的误差为 1.40%，这仅比计算了 160 次的改进欧拉法的误差略大一点。由此可以得出更精确的算法效率更高，它能在相似的工作量下产生更好的结果，或者以更少的工作量产生相似的结果。

龙格-库塔法与其它固定步长方法有着同样的问题，就是当局部截断误差在研究区间内波动很大时它们的表现往往不是很好。为了在区间某些部分达到可以接受的精度而设定的很小步长，在区间的其他部分可能远小于实际所需的步长。因此发展出自适应龙格-库塔法（adaptive Runge-Kutta method）。可以在计算过程中自动调整步长，从而使局部截断误差保持在指定的容差附近。如 8.2 节所述，在每一步都要对局部截断误差进行估算。实现的一种方法是使用五阶方法（其局部截断误差正比于 $h^6$）重复计算，然后将两种方法的结果之差作为误差的估算值。如果直接这样做，除了原有的四阶方法所需的求值次数外，五阶方法在每一步至少还需要额外计算五次 $f$ 的值。如果我们对四阶龙格-库塔公式中的中间点和权重系数进行恰当的选择，那么计算 $k_{n1},\dots,k_{n4}$ 的这些值就可以在相应的五阶方法中被重复使用，且仅需额外增加一个阶段（即多算一个 $k$ 值）。这极大地提高了效率。事实证明，实现这一点的方法不止一种。

第一个四阶和五阶龙格-库塔方法是由埃尔温·费尔伯格（Erwin Fehlberg）在 20 世纪 60 年代后期研究得到的，现在被称为龙格-库塔-费尔伯格法（Runge-Kutta-Fehlberg method），简称 RKF 方法。RKF 方法及其他自适应龙格-库塔法是为极其广泛的初值问题寻找数值近似解的强大且高效的方法。这些方法在很多在商业软件库中都有实现。