查看“︁龙格-库塔法”︁的源代码

[[数值分析]]中，'''龙格-库塔法'''（英文：Runge-Kutta methods）是用于[[非線性系統|非线性常微分方程]]的解的重要的一类隐式或显式迭代法。这些技术由数学家[[卡尔·龙格]]和[[马丁·威尔海姆·库塔]]于1900年左右发明。<!--背景知识和其它方法请参看[[数值常微分方程]]条目。-->

== 四阶龙格-库塔法 ==
在各種龙格－库塔法當中有一個方法十分常用，以至于经常被称为“RK4”或者就是“龙格－库塔法”。该方法主要是在已知方程导数和初始值时，利用计算机的仿真应用，省去求解微分方程的复杂过程。

令[[初值问题]]表述如下。

:<math> y' = f(t, y), \quad y(t_0) = y_0 </math>

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

:<math> y_{n+1} = y_n + {h \over 6} (k_1 + 2k_2 + 2k_3 + k_4) </math>

其中

:<math> k_1 = f \left( t_n, y_n \right) </math>

:<math> k_2 = f \left( t_n + {h \over 2}, y_n + {h \over 2} k_1 \right) </math>

:<math> k_3 = f \left( t_n + {h \over 2}, y_n + {h \over 2} k_2 \right) </math>

:<math> k_4 = f \left( t_n + h, y_n + hk_3 \right) </math>

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

*''k''<sub>1</sub>是时间段开始时的斜率；
*''k''<sub>2</sub>是时间段中点的斜率，通过[[欧拉法]]采用斜率''k''<sub>1</sub>来决定''y''在点''t''<sub>''n''</sub> + ''h''/2的值；
*''k''<sub>3</sub>也是中点的斜率，但是这次采用斜率''k''<sub>2</sub>决定''y''值；
*''k''<sub>4</sub>是时间段终点的斜率，其''y''值用''k''<sub>3</sub>决定。

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

:<math>\mbox{slope} = \frac{k_1 + 2k_2 + 2k_3 + k_4}{6}.</math>

RK4法是四阶方法，也就是说每步的误差是''h''<sup>5</sup>[[大O记法|阶]]，而总积累误差为''h''<sup>4</sup>阶。

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

== 显式龙格-库塔法 ==
显式龙格-库塔法是上述RK4法的一个推广。它由下式给出
:<math> y_{n+1} = y_n + h\sum_{i=1}^s b_i k_i, </math> 
其中
:<math> k_1 = f(t_n, y_n), \, </math>
:<math> k_2 = f(t_n+c_2h, y_n+a_{21}hk_1), \, </math>
:<math> k_3 = f(t_n+c_3h, y_n+a_{31}hk_1+a_{32}hk_2), \, </math>
:::<math> \vdots </math>
:<math> k_s = f(t_n+c_sh, y_n+a_{s1}hk_1+a_{s2}hk_2+\cdots+a_{s,s-1}hk_{s-1}). </math>
（注意：上述方程在不同著述中有不同但等价的定义）。

要给定一个特定的方法，必须提供整数''s''（级数），以及系数 ''a''<sub>''ij''</sub>（对于1 &le; ''j'' < ''i'' &le; ''s''）, ''b''<sub>''i''</sub>（对于''i'' = 1, 2, ..., ''s''）和''c''<sub>''i''</sub>（对于''i'' = 2, 3, ..., ''s''）。这些数据通常排列在一个助记工具中，称为''Butcher tableau''（得名自{{le|約翰·C·布徹|John C. Butcher}}）：

{| cellpadding=3px cellspacing=0px
|width="20px"| || style="border-right:1px solid;" | 0
|- 
||| style="border-right:1px solid;" | <math> c_2 </math> || <math> a_{21} </math>
|- 
||| style="border-right:1px solid;" | <math> c_3 </math> || <math> a_{31} </math> || <math> a_{32} </math>
|- 
||| style="border-right:1px solid;" | <math> \vdots </math> || <math> \vdots </math> || || <math> \ddots </math>
|- 
||| style="border-right:1px solid; border-bottom:1px solid;" | <math> c_s </math> 
| style="border-bottom:1px solid;" | <math> a_{s1} </math> 
| style="border-bottom:1px solid;" | <math> a_{s2} </math> 
| style="border-bottom:1px solid;" | <math> \cdots </math>
| style="border-bottom:1px solid;" | <math> a_{s,s-1} </math> || style="border-bottom:1px solid;" | 
|-
||| style="border-right:1px solid;" | || <math> b_1 </math> || <math> b_2 </math> || <math> \cdots </math> || <math> b_{s-1} </math> || <math> b_s </math>
|}

龙格－库塔法是自洽的，如果
:<math>\sum_{j=1}^{i-1} a_{ij} = c_i\ \mathrm{for}\ i=2, \ldots, s.</math>
如果要求方法的精度為''p''階，即截斷誤差為O(''h''<sup>''p''+1</sup>)的，则还有相应的条件。这些可以从截斷誤差本身的定义中导出。例如，一个2级2阶方法要求''b''<sub>1</sub> + ''b''<sub>2</sub> = 1, ''b''<sub>2</sub>''c''<sub>2</sub> = 1/2, 以及''b''<sub>2</sub>''a''<sub>21</sub> = 1/2。

=== 例子 ===
RK4法处于这个框架之内。其表为：

{| cellpadding=3px cellspacing=0px
|width="20px"| || style="border-right:1px solid;" | 0
|- 
||| style="border-right:1px solid;" | 1/2 || 1/2
|- 
||| style="border-right:1px solid;" | 1/2 || 0 || 1/2
|- 
||| style="border-right:1px solid; border-bottom:1px solid;" | 1 || style="border-bottom:1px solid;" | 0
| style="border-bottom:1px solid;" | 0 || style="border-bottom:1px solid;" | 1
| style="border-bottom:1px solid;" |
|-
||| style="border-right:1px solid;" | || 1/6 || 1/3 || 1/3 || 1/6
|}

然而，最简单的龙格－库塔法是（更早发现的）[[欧拉方法]]，其公式為<math> y_{n+1} = y_n + hf(t_n,y_n) </math>。这是唯一自洽的一级显式龙格－库塔法。相应的表为：

{| cellpadding=3px cellspacing=0px
| width="20px" | || style="border-right:1px solid; border-bottom:1px solid;" width="10px" | 0 
| style="border-bottom:1px solid;" width="10px" | 
|-
||| style="border-right:1px solid;" | || 1
|}

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

:<math> y_{n+1} = y_n + h \sum_{i=1}^s b_i k_i, </math>

其中

:<math> 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.</math>

在显式龙格－库塔法的框架里，定义参数<math>a_{ij} </math>的矩阵是一个[[下三角矩阵]]，而隐式龙格－库塔法并没有这个性质，这是两个方法最直观的区别：

:<math>
\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}
</math>

需要注意的是，与显式龙格－库塔法不同，隐式龙格－库塔法在每一步的计算里需要求解一个线性方程組，这相应的增加了计算的成本。

== 参考 ==

* 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.)''
{{常微分方程数值方法}}
[[Category:数值微分方程]]