# Runge-Kutta法

## 經典RK法

Runge-Kutta系列裏便最廣為人知嘅一種通常着喊做「 RK4」、「經典Runge-Kutta法」抑或直程「 Runge-Kutta法」，「4」係因爲有四階。

${\displaystyle {\frac {dy}{dt}}=f(t,y),\quad y(t_{0})=y_{0}.}$

{\displaystyle {\begin{aligned}y_{n+1}&=y_{n}+{\frac {1}{6}}h\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right),\\t_{n+1}&=t_{n}+h\\\end{aligned}}}

{\displaystyle {\begin{aligned}k_{1}&=\ f(t_{n},y_{n}),\\k_{2}&=\ f\left(t_{n}+{\frac {h}{2}},y_{n}+h{\frac {k_{1}}{2}}\right),\\k_{3}&=\ f\left(t_{n}+{\frac {h}{2}},y_{n}+h{\frac {k_{2}}{2}}\right),\\k_{4}&=\ f\left(t_{n}+h,y_{n}+hk_{3}\right).\end{aligned}}}
（注意：以上等式喺唔同嘅文本裏便具有唔同但等效嘅定義）。[3]

• ${\displaystyle k_{1}}$係區間開頭嘅斜率，用到${\displaystyle y}$歐拉法）;
• ${\displaystyle k_{2}}$係區間中點嘅斜率，用到${\displaystyle y}$${\displaystyle k_{1}}$ ;
• ${\displaystyle k_{3}}$又係中點度嘅斜率，之用到${\displaystyle y}$${\displaystyle k_{2}}$ ;
• ${\displaystyle k_{4}}$係區間孻尾嘅斜率，用到${\displaystyle y}$${\displaystyle k_{3}}$

## 顯式RK法

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

{\displaystyle {\begin{aligned}k_{1}&=f(t_{n},y_{n}),\\k_{2}&=f(t_{n}+c_{2}h,y_{n}+h(a_{21}k_{1})),\\k_{3}&=f(t_{n}+c_{3}h,y_{n}+h(a_{31}k_{1}+a_{32}k_{2})),\\&\ \ \vdots \\k_{s}&=f(t_{n}+c_{s}h,y_{n}+h(a_{s1}k_{1}+a_{s2}k_{2}+\cdots +a_{s,s-1}k_{s-1})).\end{aligned}}}
（注意：以上公式喺某些文本裏便可能具有唔同但等效嘅定義）。

 ${\displaystyle 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 _{i=1}^{s}b_{i}=1.}$

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

${\displaystyle {\begin{array}{c|cccccccc}p&1&2&3&4&5&6&7&8\\\hline \min s&1&2&3&4&6&7&9&11\end{array}}}$

### 例子

RK4法屬於呢個框架。佢個Butcher表係：[12]

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

 0 1/3 1/3 2/3 -1/3 1 1 1 −1 1 1/8 3/8 3/8 1/8

 0 1

### 兩級二階方法

${\displaystyle y_{n+1}=y_{n}+hf\left(t_{n}+{\frac {1}{2}}h,y_{n}+{\frac {1}{2}}hf(t_{n},\ y_{n})\right).}$

 0 1/2 1/2 0 1

${\displaystyle y_{n+1}=y_{n}+h{\bigl (}(1-{\tfrac {1}{2\alpha }})f(t_{n},y_{n})+{\tfrac {1}{2\alpha }}f(t_{n}+\alpha h,y_{n}+\alpha hf(t_{n},y_{n})){\bigr )}.}$

 0 ${\displaystyle \alpha }$ ${\displaystyle \alpha }$ ${\displaystyle (1-{\tfrac {1}{2\alpha }})}$ ${\displaystyle {\tfrac {1}{2\alpha }}}$

## 使法

 0 2/3 2/3 1/4 3/4

{\displaystyle {\begin{aligned}k_{1}&=f(t_{n},\ y_{n}),\\k_{2}&=f(t_{n}+{\tfrac {2}{3}}h,\ y_{n}+{\tfrac {2}{3}}hk_{1}),\\y_{n+1}&=y_{n}+h\left({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2}\right).\end{aligned}}}

${\displaystyle {\frac {dy}{dt}}=\tan(y)+1,\quad y_{0}=1,\ t\in [1,1.1]}$

 ${\displaystyle t_{0}=1\colon }$ ${\displaystyle y_{0}=1}$ ${\displaystyle t_{1}=1.025\colon }$ ${\displaystyle y_{0}=1}$ ${\displaystyle k_{1}=2.557407725}$ ${\displaystyle k_{2}=f(t_{0}+{\tfrac {2}{3}}h,\ y_{0}+{\tfrac {2}{3}}hk_{1})=2.7138981400}$ ${\displaystyle y_{1}=y_{0}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.066869388}}}$ ${\displaystyle t_{2}=1.05\colon }$ ${\displaystyle y_{1}=1.066869388}$ ${\displaystyle k_{1}=2.813524695}$ ${\displaystyle k_{2}=f(t_{1}+{\tfrac {2}{3}}h,\ y_{1}+{\tfrac {2}{3}}hk_{1})}$ ${\displaystyle y_{2}=y_{1}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.141332181}}}$ ${\displaystyle t_{3}=1.075\colon }$ ${\displaystyle y_{2}=1.141332181}$ ${\displaystyle k_{1}=3.183536647}$ ${\displaystyle k_{2}=f(t_{2}+{\tfrac {2}{3}}h,\ y_{2}+{\tfrac {2}{3}}hk_{1})}$ ${\displaystyle y_{3}=y_{2}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.227417567}}}$ ${\displaystyle t_{4}=1.1\colon }$ ${\displaystyle y_{3}=1.227417567}$ ${\displaystyle k_{1}=3.796866512}$ ${\displaystyle k_{2}=f(t_{3}+{\tfrac {2}{3}}h,\ y_{3}+{\tfrac {2}{3}}hk_{1})}$ ${\displaystyle y_{4}=y_{3}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.335079087}}.}$

## 自適應RK法

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

${\displaystyle e_{n+1}=y_{n+1}-y_{n+1}^{*}=h\sum _{i=1}^{s}(b_{i}-b_{i}^{*})k_{i},}$

 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 b_{1}^{*}}$ ${\displaystyle b_{2}^{*}}$ ${\displaystyle \cdots }$ ${\displaystyle b_{s-1}^{*}}$ ${\displaystyle b_{s}^{*}}$

Runge-Kutta-Fehlberg法具有5跟4階兩種方法。佢嘅擴展後嘅Butcher表係：

 0 1/4 1/4 3/8 3/32 9/32 12/13 1932/2197 −7200/2197 7296/2197 1 439/216 -8 3680/513 -845/4104 1/2 −8/27 2 −3544/2565 1859/4104 −11/40 16/135 0 6656/12825 28561/56430 −9/50 2/55 25/216 0 1408/2565 2197/4104 −1/5 0

 0 1 1 1/2 1/2 1 0

## Runge-Kutta-Nyström法

Runge-Kutta-Nyström法係特定嘅Runge-Kutta法，之針對二階微分方程進行唨改良[16] [17]

${\displaystyle {\frac {d^{2}y}{dt^{2}}}=f(y,{\dot {y}},t).}$

## 隱式RK法

${\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.}$ [19]

${\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}\\&b_{1}^{*}&b_{2}^{*}&\dots &b_{s}^{*}\\\end{array}}={\begin{array}{c|c}\mathbf {c} &A\\\hline &\mathbf {b^{T}} \\\end{array}}}$

### 例子

${\displaystyle y_{n+1}=y_{n}+hf(t_{n}+h,\ y_{n+1}).\,}$

${\displaystyle {\begin{array}{c|c}1&1\\\hline &1\\\end{array}}}$

${\displaystyle k_{1}=f(t_{n}+h,\ y_{n}+hk_{1})\quad {\text{and}}\quad y_{n+1}=y_{n}+hk_{1},}$

${\displaystyle {\begin{array}{c|cc}0&0&0\\1&{\frac {1}{2}}&{\frac {1}{2}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\&1&0\\\end{array}}}$

Gauss-Legendre法構成唨一系列基於Gauss正交嘅搭配方法。s級嘅Gauss-Legendre法嘅階數為2 s （因此，可以構造任意高階嘅方法）。 [22]兩級（即有四階）嘅方法會有噉樣嘅Butcher表：

${\displaystyle {\begin{array}{c|cc}{\frac {1}{2}}-{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}&{\frac {1}{4}}-{\frac {1}{6}}{\sqrt {3}}\\{\frac {1}{2}}+{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}+{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\&{\frac {1}{2}}+{\frac {1}{2}}{\sqrt {3}}&{\frac {1}{2}}-{\frac {1}{2}}{\sqrt {3}}\end{array}}}$ [20]

### 穩定

${\displaystyle r(z)=1+zb^{T}(I-zA)^{-1}e={\frac {\det(I-zA+zeb^{T})}{\det(I-zA)}},}$ [23]

s級嘅Gauss-Legendre法嘅階數為2 s ，因此佢個穩定性函數為m = n = s嘅Padé近似值。因此，呢個方法係A穩定嘅。 [27]噉樣表明A穩定嘅Runge-Kutta可以有任意高嘅階。相反，A穩定線性多步法嘅階數唔超得過2。

## B穩定性

${\displaystyle B}$${\displaystyle M}$${\displaystyle Q}$係三個${\displaystyle s\times s}$矩陣，定義爲：

${\displaystyle B=\operatorname {diag} (b_{1},b_{2},\ldots ,b_{s}),\,M=BA+A^{T}B-bb^{T},\,Q=BA^{-1}+A^{-T}B-A^{-T}bb^{T}A^{-1}.}$

## 四階RK法嘅推導

${\displaystyle y_{t+h}=y_{t}+h\cdot \sum _{i=1}^{s}a_{i}k_{i}+{\mathcal {O}}(h^{s+1}),}$

${\displaystyle k_{i}=y_{t}+h\cdot \sum _{j=1}^{s}\beta _{ij}f\left(k_{j},\ t_{n}+\alpha _{i}h\right)}$

{\displaystyle {\begin{aligned}&\alpha _{i}&&\beta _{ij}\\\alpha _{1}&=0&\beta _{21}&={\frac {1}{2}}\\\alpha _{2}&={\frac {1}{2}}&\beta _{32}&={\frac {1}{2}}\\\alpha _{3}&={\frac {1}{2}}&\beta _{43}&=1\\\alpha _{4}&=1&&\\\end{aligned}}}

{\displaystyle {\begin{aligned}y_{t+h}^{1}&=y_{t}+hf\left(y_{t},\ t\right)\\y_{t+h}^{2}&=y_{t}+hf\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)\\y_{t+h}^{3}&=y_{t}+hf\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)\end{aligned}}}

{\displaystyle {\begin{aligned}k_{1}&=f(y_{t},\ t)\\k_{2}&=f\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right)\\k_{3}&=f\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{2},\ t+{\frac {h}{2}}\right)\\k_{4}&=f\left(y_{t+h}^{3},\ t+h\right)=f\left(y_{t}+hk_{3},\ t+h\right)\end{aligned}}}

{\displaystyle {\begin{aligned}k_{2}&=f\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right)\\&=f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\\k_{3}&=f\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right),\ t+{\frac {h}{2}}\right)\\&=f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\right]\\k_{4}&=f\left(y_{t+h}^{3},\ t+h\right)=f\left(y_{t}+hf\left(y_{t}+{\frac {h}{2}}k_{2},\ t+{\frac {h}{2}}\right),\ t+h\right)\\&=f\left(y_{t}+hf\left(y_{t}+{\frac {h}{2}}f\left(y_{t}+{\frac {h}{2}}f\left(y_{t},\ t\right),\ t+{\frac {h}{2}}\right),\ t+{\frac {h}{2}}\right),\ t+h\right)\\&=f\left(y_{t},\ t\right)+h{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\right]\right]\end{aligned}}}

${\displaystyle {\frac {d}{dt}}f(y_{t},\ t)={\frac {\partial }{\partial y}}f(y_{t},\ t){\dot {y}}_{t}+{\frac {\partial }{\partial t}}f(y_{t},\ t)=f_{y}(y_{t},\ t){\dot {y}}+f_{t}(y_{t},\ t):={\ddot {y}}_{t}}$

{\displaystyle {\begin{aligned}y_{t+h}={}&y_{t}+h\left\lbrace a\cdot f(y_{t},\ t)+b\cdot \left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right.+\\&{}+c\cdot \left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right]+\\&{}+d\cdot \left[f(y_{t},\ t)+h{\frac {d}{dt}}\left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}\left[f(y_{t},\ t)+\left.{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right]\right]\right\rbrace +{\mathcal {O}}(h^{5})\\={}&y_{t}+a\cdot hf_{t}+b\cdot hf_{t}+b\cdot {\frac {h^{2}}{2}}{\frac {df_{t}}{dt}}+c\cdot hf_{t}+c\cdot {\frac {h^{2}}{2}}{\frac {df_{t}}{dt}}+\\&{}+c\cdot {\frac {h^{3}}{4}}{\frac {d^{2}f_{t}}{dt^{2}}}+d\cdot hf_{t}+d\cdot h^{2}{\frac {df_{t}}{dt}}+d\cdot {\frac {h^{3}}{2}}{\frac {d^{2}f_{t}}{dt^{2}}}+d\cdot {\frac {h^{4}}{4}}{\frac {d^{3}f_{t}}{dt^{3}}}+{\mathcal {O}}(h^{5})\end{aligned}}}

{\displaystyle {\begin{aligned}y_{t+h}&=y_{t}+h{\dot {y}}_{t}+{\frac {h^{2}}{2}}{\ddot {y}}_{t}+{\frac {h^{3}}{6}}y_{t}^{(3)}+{\frac {h^{4}}{24}}y_{t}^{(4)}+{\mathcal {O}}(h^{5})=\\&=y_{t}+hf(y_{t},\ t)+{\frac {h^{2}}{2}}{\frac {d}{dt}}f(y_{t},\ t)+{\frac {h^{3}}{6}}{\frac {d^{2}}{dt^{2}}}f(y_{t},\ t)+{\frac {h^{4}}{24}}{\frac {d^{3}}{dt^{3}}}f(y_{t},\ t)\end{aligned}}}

${\displaystyle {\begin{cases}&a+b+c+d=1\\[6pt]&{\frac {1}{2}}b+{\frac {1}{2}}c+d={\frac {1}{2}}\\[6pt]&{\frac {1}{4}}c+{\frac {1}{2}}d={\frac {1}{6}}\\[6pt]&{\frac {1}{4}}d={\frac {1}{24}}\end{cases}}}$

## 註

1. DEVRIES, Paul L. ; HASBUN, Javier E. A first course in computational physics. Second edition. Jones and Bartlett Publishers: 2011. p. 215.
2. Press et al. 2007, p. 908; Süli & Mayers 2003, p. 328
3. Atkinson (1989, p. 423), Hairer, Nørsett & Wanner (1993, p. 134), Kaw & Kalu (2008, §8.4) and Stoer & Bulirsch (2002, p. 476) leave out the factor h in the definition of the stages. Ascher & Petzold (1998, p. 81), Butcher (2008, p. 93) and Iserles (1996, p. 38) use the y values as stages.
4. Süli & Mayers 2003, p. 328
5. Press et al. 2007, p. 907
6. Iserles 1996, p. 38
7. Iserles 1996, p. 39
8. Iserles 1996, p. 39
9. As a counterexample, consider any explicit 2-stage Runge-Kutta scheme with ${\displaystyle b_{1}=b_{2}=1/2}$ and ${\displaystyle c_{1}}$ and ${\displaystyle a_{21}}$ randomly chosen. This method is consistent and (in general) first-order convergent. On the other hand, the 1-stage method with ${\displaystyle b_{1}=1/2}$ is inconsistent and fails to converge, even though it trivially holds that ${\displaystyle \sum _{j=1}^{i-1}a_{ij}=c_{i}{\text{ for }}i=2,\ldots ,s.}$.
10. Butcher 2008, p. 187
11. Butcher 2008
12. Süli & Mayers 2003, p. 352
13. Hairer, Nørsett & Wanner (1993, p. 138) refer to Kutta (1901).
14. Süli & Mayers 2003, p. 327
15. Lambert 1991, p. 278
16. Dormand, J. R.; Prince, P. J. (October 1978). "New Runge-Kutta Algorithms for Numerical Simulation in Dynamical Astronomy". Celestial Mechanics. 18 (3): 223–232. doi:10.1007/BF01230162.
17. Fehlberg, E. (October 1974). Classical seventh-, sixth-, and fifth-order Runge-Kutta-Nyström formulas with stepsize control for general second-order differential equations (報告) (第NASA TR R-432版). Marshall Space Flight Center, AL: National Aeronautics and Space Administration.
18. Süli & Mayers 2003, pp. 349–351
19. Iserles 1996, p. 41; Süli & Mayers 2003
20. Süli & Mayers 2003, p. 353
21. Iserles 1996
22. Iserles 1996, p. 47
23. Hairer & Wanner 1996
24. Hairer & Wanner 1996, p. 40
25. Iserles 1996, p. 60
26. Iserles 1996
27. Iserles 1996, p. 63
28. Lambert 1991, p. 275
29. Lambert 1991, p. 274
30. PDF reporting this derivation