龍格-庫塔法

數值分析中，龍格-庫塔法（英文：Runge-Kutta methods）是用於非線性常微分方程的解的重要的一類隱式或顯式迭代法。這些技術由數學家卡爾·龍格和馬丁·威爾海姆·庫塔於1900年左右發明。

四階龍格-庫塔法

在各種龍格－庫塔法當中有一個方法十分常用，以至於經常被稱為「RK4」或者就是「龍格－庫塔法」。該方法主要是在已知方程導數和初始值時，利用計算機的仿真應用，省去求解微分方程的複雜過程。

令初值問題表述如下。

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

則，對於該問題的RK4由如下方程給出：

${\displaystyle y_{n+1}=y_{n}+{h \over 6}(k_{1}+2k_{2}+2k_{3}+k_{4})}$

其中

${\displaystyle k_{1}=f\left(t_{n},y_{n}\right)}$

${\displaystyle k_{2}=f\left(t_{n}+{h \over 2},y_{n}+{h \over 2}k_{1}\right)}$

${\displaystyle k_{3}=f\left(t_{n}+{h \over 2},y_{n}+{h \over 2}k_{2}\right)}$

${\displaystyle k_{4}=f\left(t_{n}+h,y_{n}+hk_{3}\right)}$

這樣，下一個值（y_n+1）由現在的值（y_n）加上時間間隔（h）和一個估算的斜率的乘積所決定。該斜率是以下斜率的加權平均：

• k₁是時間段開始時的斜率；
• k₂是時間段中點的斜率，通過歐拉法採用斜率k₁來決定y在點t_n + h/2的值；
• k₃也是中點的斜率，但是這次採用斜率k₂決定y值；
• k₄是時間段終點的斜率，其y值用k₃決定。

當四個斜率取平均時，中點的斜率有更大的權值：

${\displaystyle {\mbox{slope}}={\frac {k_{1}+2k_{2}+2k_{3}+k_{4}}{6}}.}$

RK4法是四階方法，也就是說每步的誤差是h⁵階，而總積累誤差為h⁴階。

注意上述公式對於純量或者向量函數（y可以是向量）都適用。

顯式龍格-庫塔法

顯式龍格-庫塔法是上述RK4法的一個推廣。它由下式給出

${\displaystyle y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},}$

其中

${\displaystyle k_{1}=f(t_{n},y_{n}),\,}$

${\displaystyle k_{2}=f(t_{n}+c_{2}h,y_{n}+a_{21}hk_{1}),\,}$

${\displaystyle k_{3}=f(t_{n}+c_{3}h,y_{n}+a_{31}hk_{1}+a_{32}hk_{2}),\,}$

${\displaystyle \vdots }$

${\displaystyle k_{s}=f(t_{n}+c_{s}h,y_{n}+a_{s1}hk_{1}+a_{s2}hk_{2}+\cdots +a_{s,s-1}hk_{s-1}).}$

（注意：上述方程在不同著述中有不同但等價的定義）。

要給定一個特定的方法，必須提供整數s（級數），以及係數 a_ij（對於1 ≤ j < i ≤ s）, b_i（對於i = 1, 2, ..., s）和c_i（對於i = 2, 3, ..., s）。這些數據通常排列在一個助記工具中，稱為Butcher tableau（得名自約翰·C·布徹（英語：John C. Butcher））：

	0
	${\displaystyle c_{2}}$	${\displaystyle a_{21}}$
	${\displaystyle c_{3}}$	${\displaystyle a_{31}}$	${\displaystyle a_{32}}$
	${\displaystyle \vdots }$	${\displaystyle \vdots }$		${\displaystyle \ddots }$
	${\displaystyle c_{s}}$	${\displaystyle a_{s1}}$	${\displaystyle a_{s2}}$	${\displaystyle \cdots }$	${\displaystyle a_{s,s-1}}$
		${\displaystyle b_{1}}$	${\displaystyle b_{2}}$	${\displaystyle \cdots }$	${\displaystyle b_{s-1}}$	${\displaystyle b_{s}}$

龍格－庫塔法是自洽的，如果

${\displaystyle \sum _{j=1}^{i-1}a_{ij}=c_{i}\ \mathrm {for} \ i=2,\ldots ,s.}$

如果要求方法的精度為p階，即截斷誤差為O(h^p+1)的，則還有相應的條件。這些可以從截斷誤差本身的定義中導出。例如，一個2級2階方法要求b₁ + b₂ = 1, b₂c₂ = 1/2, 以及b₂a₂₁ = 1/2。

例子

RK4法處於這個框架之內。其表為：

	0
	1/2	1/2
	1/2	0	1/2
	1	0	0	1
		1/6	1/3	1/3	1/6

然而，最簡單的龍格－庫塔法是（更早發現的）歐拉方法，其公式為${\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n})}$。這是唯一自洽的一級顯式龍格－庫塔法。相應的表為：

	0
		1

隱式龍格-庫塔法

以上提及的顯式龍格－庫塔法一般來講不適用於求解剛性方程。這是因為顯式龍格－庫塔法的穩定區域被局限在一個特定的區域裏。顯式龍格－庫塔法的這種缺陷使得人們開始研究隱式龍格－庫塔法，一般而言，隱式龍格－庫塔法具有以下形式：

${\displaystyle y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},}$

其中

${\displaystyle k_{i}=f\left(t_{n}+c_{i}h,y_{n}+h\sum _{j=1}^{s}a_{ij}k_{j}\right),\quad i=1,\ldots ,s.}$

在顯式龍格－庫塔法的框架里，定義參數${\displaystyle a_{ij}}$的矩陣是一個下三角矩陣，而隱式龍格－庫塔法並沒有這個性質，這是兩個方法最直觀的區別：

${\displaystyle {\begin{array}{c|cccc}c_{1}&a_{11}&a_{12}&\dots &a_{1s}\\c_{2}&a_{21}&a_{22}&\dots &a_{2s}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{s1}&a_{s2}&\dots &a_{ss}\\\hline &b_{1}&b_{2}&\dots &b_{s}\\\end{array}}={\begin{array}{c|c}\mathbf {c} &A\\\hline &\mathbf {b^{T}} \\\end{array}}}$

需要注意的是，與顯式龍格－庫塔法不同，隱式龍格－庫塔法在每一步的計算里需要求解一個線性方程組，這相應的增加了計算的成本。

參考

• George E. Forsythe, Michael A. Malcolm, and Cleve B. Moler. Computer Methods for Mathematical Computations. Englewood Cliffs, NJ: Prentice-Hall, 1977. (See Chapter 6.)
• Ernst Hairer, Syvert Paul Nørsett, and Gerhard Wanner. Solving ordinary differential equations I: Nonstiff problems, second edition. Berlin: Springer Verlag, 1993. ISBN 3-540-56670-8.
• William H. Press, Brian P. Flannery, Saul A. Teukolsky, William T. Vetterling. Numerical Recipes in C. Cambridge, UK: Cambridge University Press, 1988. (See Sections 16.1 and 16.2.)