Runge-Kutta法

喺數值分析當中， **Runge-Kutta法**（，）、又可以叫**龍俄-厥他法**（），簡稱**RK法**，查實係一系列隱式同顯式迭代方法，其中包括眾所周知嘅喊做Euler法嘅例程，呢個例程用於幫常微分方程嘅近似解做時間離散化。 啲方法係由德國數學家Carl Runge跟Wilhelm Kutta於1900年左右開發出嚟嘅。

## 經典RK法

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

指定一個初值問題：

:  \frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0. 

其中y係時間t嘅未知函數（標量抑或向量）  ，有待逼近；已知y嘅變化值\frac{dy}{dt}係t跟y本身嘅函數。喺初始時間t_0，個y值係y_0 。畀有函數f跟初始條件t_0 ，y_0。

而今揀一個𨂾長*h* &#x3E; 0並定義埋：

: \begin{align}
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{align}

對於*n* = 0、1、2、3，...，使用

: 
\begin{align}
 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{align}

: *（注意：以上等式喺唔同嘅文本裏便具有唔同但等效嘅定義）。*

呢度嘅y_{n+1}係y(t_{n+1})嘅RK4近似值 ，同埋下一個值（ y_{n+1} ）由當前值（ y_n ）加埋四個增量嘅加權平均值嚟決定，其中每個增量都係一個乘積、由間隔大細*h*乘埋一個估出嘅斜率得出嘅，個斜率*k*由微分方程右側嘅函數 *f* 指定。

* k_1係區間開頭嘅斜率，用到 y  （歐拉法）;
* k_2係區間中點嘅斜率，用到 y 跟 k_1  ;
* k_3又係中點度嘅斜率，之用到 y 跟 k_2  ;
* k_4係區間孻尾嘅斜率，用到 y 跟 k_3  。

喺對四個斜率求平均值陣時，對中點度嘅斜率賦予大啲嘅權重。若果f獨立於y ，即微分方程相當於一個簡單嘅積分，係噉RK4就係。 

種RK4方法係四階方法，噉樣代表局部截尾誤差同階於O(h^5) ，而總累積誤差同階於O(h^4) 。

喺許多實際應用裏便，隻函數f獨立於t （所謂嘅自治系統（autonomous system）抑或時唔變系統，特別係喺物理學裏便），而且佢哋嘅增量根本唔會得到計算，亦都唔會傳遞畀函數f ，只有對t_{n+1}個最終公式用到。

## 顯式RK法
顯式Runge-Kutta法家族係上述RK4方法嘅概括，之係噉樣得出：

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

其中

: 
\begin{align}
 k_1 & = f(t_n, y_n), \\
 k_2 & = f(t_n+c_2h, y_n+h(a_{21}k_1)), \\
 k_3 & = f(t_n+c_3h, y_n+h(a_{31}k_1+a_{32}k_2)), \\
   & \ \ \vdots \\
 k_s & = f(t_n+c_sh, y_n+h(a_{s1}k_1+a_{s2}k_2+\cdots+a_{s,s-1}k_{s-1})).
\end{align}

: *（注意：以上公式喺某些文本裏便可能具有唔同但等效嘅定義）。* 

要指定一個特定嘅方法，就要提供整數*s*（級數）、係數*a<sub>ij</sub>*（對於 1 ≤ *j* i*（對於*i* = 1，2，...， *s* ）同埋*ci* （對於*i* = 2，3，...， *s* ）。矩陣[*aij*]着喊做*'Runge-Kutta矩陣**，而 ''b<sub>i</sub>'' 、''c<sub>i</sub>''着喊做**權重**、**節點**。 啲數據通常係攞一套助記手段喊做**Butcher表**（''Butcher tableau'' ）嘅形式嚟排列（根據返John C. Butcher個名）：

[[泰勒級數]]展開表明，Runge-Kutta法相容（consistent）若且唯若有：

: <math>\sum_{i=1}^{s} b_{i} = 1.</math>

若果要求呢個方法有返一定階數''p'' ，噉亦都有相應嘅要求，噉樣意味住局部截尾誤差為O(''h<sup>p</sup>''<sup>+1</sup>)。啲可以從截尾誤差本身嘅定義裏便得出。譬如，若果''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，噉兩級方法嘅階數為2。 請注意，確定係數嘅常用條件係：

: <math>\sum_{j=1}^{i-1} a_{ij} = c_i \text{ for } i=2, \ldots, s.</math>

之不過，僅衹呢個條件唔夠保持相容性。 

一般嚟講，若果<math>s</math>階龍格-庫塔方法有<math>p</math>階，噉就證明得到級數必須滿足<math>s \ge p</math> ， 而若果<math>p \ge 5</math> ， 就有<math>s \ge p+1</math> 。 但係，仲唔知道喺所有情況下啲界限係唔係都咁**清晰**（''sharp''）。譬如，所有啲已知嘅8階方法至少有11級，即使有可能存在啲方法少啲級嘅。 （上便個範圍表明到可能係有9級嘅方法；但亦都可能係個範圍根本唔清晰。 ）實質上，確切嘅最細級數<math>s</math>仲係一個開放問題（open problem）、對於一個顯式<math>p</math>階Runge-Kutta法嚟講，因爲仲未發現啲方法、滿足得上述啲範圍入便啲等號嘅。一些已知嘅值係： 

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

上便啲證明得嘅範圍就意味住搵唔到階數<math>p=1, 2, \ldots, 6</math>嘅方法，所需級數少過目前搵到嘅同階方法嘅。但係，可以想像到，搵得到一種階數<math>p=7</math>之只有8級嘅方法，之不過如表所示而今已知嘅方法至少要到9級。

=== 例子 ===
RK4法屬於呢個框架。佢個Butcher表係：

種Runge-Kutta法嘅其中一種熹微變體亦都歸功於1901年Kutta，喊做3/8規則。呢個方法嘅主要優點係，幾乎所有啲誤差係數都細過啲流行方法入便嘅誤差係數，但每一𨂾時間需要稍爲多啲嘅FLOP（浮點運算）。佢個Butcher表係：

但係，最簡單嘅Runge-Kutta法梗係（正向）[[歐拉法]]，公式係<math> y_{n+1} = y_n + hf(t_n, y_n) </math> 。呢個係唯一相容（consistent）嘅一級顯式Runge-Kutta法。個相應圖表係：


=== 兩級二階方法 ===
中點法係一個兩級二階方法嘅示例：

: <math> y_{n+1} = y_n + hf\left(t_n+\frac{1}{2}h, y_n+\frac{1}{2}hf(t_n,\ y_n)\right). </math>

相應嘅圖表係：

中點法唔係唯一嘅兩級二階Runge-Kutta法。有一系列噉樣嘅方法，由α參數設定並由以下公式畀出：

: <math> y_{n+1} = y_n + h\bigl( (1-\tfrac1{2\alpha}) f(t_n, y_n) + \tfrac1{2\alpha} f(t_n + \alpha h, y_n + \alpha h f(t_n, y_n)) \bigr). </math>

個Butcher表係：

喺噉樣嘅方法族入便，<math>\alpha=\tfrac12</math>嗰陣畀出中點法，而<math>\alpha=1</math>嗰陣係[[Heun法]]。 

== 使法 ==
譬如，考慮α= 2/3嘅兩級二階Runge-Kutta法，之亦都喊做[[Ralston方法]]。佢由表格畀出：

用到相應嘅方程式

: <math> \begin{align}
k_1 &= f(t_n,\ y_n), \\
k_2 &= f(t_n + \tfrac{2}{3}h,\ y_n + \tfrac{2}{3}h k_1), \\
y_{n+1} &= y_n + h\left(\tfrac{1}{4}k_1+\tfrac{3}{4}k_2\right).
\end{align} </math>

呢個方法用於解決初值問題

: <math> \frac{dy}{dt} = \tan(y)+1,\quad y_0=1,\ t\in [1, 1.1]</math>

其中𨂾長''h'' = 0.025，因此呢個方法需要執行四𨂾。

呢個方法進行如下：

數值解對應於帶下劃線啲值。

== 自適應RK法 ==
自適應方法旨喺產生單個Runge-Kutta𨂾步嘅局部截尾誤差嘅估計。噉樣可以通過兩種方法來完成，一種係有序嘅<math>p</math>還有一個訂單<math>p - 1</math> 。啲方法係交織嘅，即佢哋有共同嘅中間𨂾步。因此，戥採用高階方法嘅一𨂾相比，估計誤差嘅計算成本好細抑或可忽略唔計。

喺積分期間，𨂾長會着調整，嚟令估計嘅誤差保持喺用戶定義嘅閾值以下：若果誤差太大，噉攞細啲嘅𨂾長重複執行過一𨂾；若果誤差細得多，就增加𨂾長嚟慳時間。噉樣（幾乎）可以畀出最佳𨂾長，同時慳返計算時間。而且，用戶唔使嘥時間來搵返啱嘅𨂾長。

低階𨂾步由以下式得出：

: <math> y^*_{n+1} = y_n + h\sum_{i=1}^s b^*_i k_i, </math>

其中<math>k_i</math>戥高階方法相同。然之後個誤差係：

: <math> e_{n+1} = y_{n+1} - y^*_{n+1} = h\sum_{i=1}^s (b_i - b^*_i) k_i, </math>

噉樣係<math>O(h^p)</math> 。幫呢種方法擴展Butcher表，嚟畀出<math>b^*_i</math> ：

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

但係，最簡單嘅自適應Runge-Kutta法涉及[[亨氏法|Heun法]]（階數2）佮[[歐拉法|Euler法]]（階數1）嘅組合。佢嘅擴展後嘅Butcher表係：

其他啲自適應嘅Runge-Kutta法有[[Bogacki-Shampine方法|Bogacki-Shampine法]]（3、2階）同埋[[現金卡普法|Cash-Karp法]]跟[[休眠-王子法|Dormand-Prince法]]（5、4階）。

== 非融合RK法 ==
若果<math>c_i,\,i=1,2,\ldots,s</math>互相之間唔同，噉Runge-Kutta法就着認爲係非''融合嘅''（''nonconfluent''）。

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

: <math> \frac{d^2 y}{d t^2}= f(y,\dot{y},t).</math>

== 隱式RK法 ==
到目前為止，所有啲提到嘅Runge-Kutta法都係[[顯式和隱式方法|顯式方法]]。顯式嘅Runge-Kutta法通常唔適合於攞嚟求解啲[[剛性方程]]（Stiff equations），因為佢哋嘅絕對穩定域好細。特別係佢係有界嘅。 噉樣個問題喺[[偏微分方程嘅數值分析|偏微分方程求解]]裏便尤為重要。

顯式Runge-Kutta法個唔穩定性激唨啲隱式方法嘅發展。隱式Runge-Kutta法具有以下形式：

: <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> 

顯式方法嘅區別喺於，喺顯式方法裏便， ''j''嘅總和僅衹等於''i'' -1。呢個亦都顯示喺張Butcher表裏便，當中顯式方法個係數矩陣<math>a_{ij}</math>係下三角形式。喺隱式方法裏便， ''j''度嘅總和上升到''s''而且係數矩陣唔係三角形，係噉畀出對應嘅Butcher表形式：

: <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\\
    & b^*_1 & b^*_2 & \dots & b^*_s\\
\end{array} =

\begin{array}{c|c}
\mathbf{c}& A\\
\hline
     & \mathbf{b^T} \\
\end{array}
</math>

請參閱返上方嘅[[Runge-Kutta法#自適應Runge-Kutta法|自適應Runge-Kutta法]]嚟獲取有關<math>b^*</math>排嘅介紹。

呢種差異導致，每一𨂾都必須求解代數方程組。噉樣大大增加唨計算成本。若果使用''s''級嘅方法求解''m''個分量嘅微分方程，噉個代數方程組將具有''ms''個分量。噉樣可以戥隱式[[線性多步法]]（ODE嘅另一種大方法系列）形成對比：隱式''s''𨂾線性多步方法需要求解只有''m''個分量嘅代數方程組，因此系統嘅大細唔會隨住𨂾數增加而增加。 

=== 例子 ===
隱式Runge-Kutta法嘅最簡單示例係[[後向歐拉法|後褪Euler法]]：

: <math>y_{n + 1} = y_n + h f(t_n + h,\ y_{n + 1}). \,</math>

張Butcher表好簡單：

: <math>
\begin{array}{c|c}
1 & 1 \\
\hline
 & 1 \\
\end{array}
</math>

呢張Butcher表戥以下公式相對應

: <math> k_1 = f(t_n + h,\ y_n + h k_1) \quad\text{and}\quad y_{n+1} = y_n + h k_1, </math>

可以重新排列嚟獲得上便列出嘅後褪Euler法嘅公式。

隱式Runge-Kutta法嘅另一個示例係[[梯形法則（微分方程）|梯形規則]]。佢個Butcher表為：

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

梯形規則係一種[[選點法]]（似文中所述）。所有選點法都係隱式Runge-Kutta法，但並非所有隱式Runge-Kutta法都係選點法。 

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

: <math>
\begin{array}{c|cc}
\frac12 - \frac16 \sqrt3 & \frac14         & \frac14 - \frac16 \sqrt3 \\
\frac12 + \frac16 \sqrt3 & \frac14 + \frac16 \sqrt3 & \frac14 \\
\hline
             & \frac12         & \frac12 \\
             & \frac12+\frac12 \sqrt3  & \frac12-\frac12 \sqrt3
\end{array}
</math> 

=== 穩定 ===
戥顯式方法相比，隱式Runge-Kutta法嘅優勢喺於其更大嘅穩定性，尤其係喺應用於[[剛性方程]]嗰陣。考慮線性測試方程''y'<nowiki/>'' = λ''y***。*'應用到呢條方程嘅Runge-Kutta法簡化成迭代式 y_{n+1} = r(h\lambda) \, y_n  ，其中*r*由下式得出：

:  r(z) = 1 + z b^T (I-zA)^{-1} e = \frac{\det(I-zA+zeb^T)}{\det(I-zA)},  

其中*e*代表1嘅向量。函數*r*喊做**穩定性函數**。 從公式可以得出，若果呢個方法有*s*級*，*噉*r*係兩個次數為*s*嘅多項式個商*。*顯式方法有嚴格意義上嘅下三角矩陣*A* ，噉樣意味住det（ *I* - *zA* ）= 1而且穩定性函數係多項式。 

若果|*r*(*z*)|  z \to 0 陣時 r(z) = \textrm{e}^z + O(z^{p+1})  。因此，研究畀定次數嘅多項式個商、之逼近指數函數嘅，會好有趣。啲喊做Padé近似值。若且唯若***'m≤N≤M***' + 2嗰陣，一個Padé近似值有*m*次嘅分子同埋*n*次嘅分母嘅先係A穩定嘅。 

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

## B穩定性
微分方程解嘅A穩定性概念同線性自治方程有關y'=\lambda y 。 Dahlquist進行唨一項研究針對啲數值方案嘅穩定性、啲應用於滿足單調性條件嘅非線性系統嘅。相應嘅概念被定義為：**G穩定性**同埋**B穩定性**（Butcher，1975），分別對應多步方法（同埋戥佢相關嘅啲單腿方法（one-leg methods））同埋Runge-Kutta法*。* Runge-Kutta法應用於非線性系統y'=f(y) ，有\langle f(y)-f(z),\ y-z \rangle；若果呢種情況意味住\|y_{n+1}-z_{n+1}\|\leq\|y_{n}-z_{n}\|兩個數值解，噉就喊做B穩定。

設B ， M，Q係三個s\times s矩陣，定義爲：

: B=\operatorname{diag}(b_1,b_2,\ldots,b_s),\, M=BA+A^TB-bb^T,\, Q=BA^{-1}+A^{-T}B-A^{-T}bb^TA^{-1}.

若果B跟M矩陣都係非負定嘅，噉呢個Runge-Kutta法就可以講係**代數穩定嘅** 。 B穩定性嘅充分條件係： B跟Q係非負定嘅。

## 四階RK法嘅推導
一般來講，s階Runge-Kutta法寫得成：

: y_{t + h} = y_t + h \cdot \sum_{i=1}^s a_i k_i +\mathcal{O}(h^{s+1}),

其中：

: k_i =y_t + h \cdot \sum_{j = 1}^s \beta_{ij} f\left(k_j,\ t_n + \alpha_i h \right)

係通過評估y_t喺第i 階嘅導數嚟獲得嘅增量。

如前所述，喺任一間隔(t,\ t +h)嘅起點、中間點、末點進行求值，捉s=4代入通用公式可以得到四階Runge-Kutta方法。因此，揀：

: 
\begin{align}
&\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{align}


否則\beta_{ij} = 0。首先定義以下數量：

: \begin{align}
y^1_{t+h} &= y_t + hf\left(y_t,\ t\right) \\
y^2_{t+h} &= y_t + hf\left(y^1_{t+h/2},\ t+\frac{h}{2}\right) \\
y^3_{t+h} &= y_t + hf\left(y^2_{t+h/2},\ t+\frac{h}{2}\right)
\end{align}

其中y^1_{t+h/2} = \dfrac{y_t + y^1_{t+h}}{2}跟y^2_{t+h/2} = \dfrac{y_t + y^2_{t+h}}{2} 。若果定義：

: \begin{align}
k_1 &= f(y_t,\ t) \\
k_2 &= f\left(y^1_{t+h/2},\ t + \frac{h}{2}\right) = f\left(y_t + \frac{h}{2} k_1,\ t + \frac{h}{2}\right) \\
k_3 &= f\left(y^2_{t+h/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^3_{t+h},\ t + h\right) = f\left(y_t + h k_3,\ t + h\right)
\end{align}

對於先前嘅關係，可以證明以下等式滿足得\mathcal{O}(h^2) ：

: \begin{align}
k_2 &= f\left(y^1_{t+h/2},\ 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^2_{t+h/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^3_{t+h},\ t + h\right) = f\left(y_t + h f\left(y_t + \frac{h}{2} k_2,\ t + \frac{h}{2}\right),\ t + h\right) \\
&= f\left(y_t + h f\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{align}

其中：

: \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

係對應時間嘅f個總導數（total derivative）。

若果而今使用啱啱導出嘅內容來表達通用公式，噉會得到：

: \begin{align}
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 h f_t + b \cdot h f_t + b \cdot \frac{h^2}{2} \frac{df_t}{dt}  + c \cdot h f_t+ c \cdot \frac{h^2}{2} \frac{df_t}{dt} + \\
&{}+ c \cdot \frac{h^3}{4} \frac{d^2f_t}{dt^2} + d \cdot h f_t + d \cdot h^2 \frac{df_t}{dt} + d \cdot \frac{h^3}{2} \frac{d^2f_t}{dt^2} + d \cdot \frac{h^4}{4} \frac{d^3f_t}{dt^3} + \mathcal{O}(h^5)
\end{align}

並捉佢戥y_t附近y_{t+h}嘅泰勒級數比較，有 ：

: \begin{align}
 y_{t+h} &= y_t + h \dot y_t + \frac{h^2}{2} \ddot y_t + \frac{h^3}{6} y^{(3)}_t + \frac{h^4}{24} y^{(4)}_t + \mathcal{O}(h^5) = \\
&= y_t + h f(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{align}

獲有一個係數約束系統：

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

個系統解出陣時畀出a = \frac{1}{6}, b = \frac{1}{3}, c = \frac{1}{3}, d = \frac{1}{6}，似之前噉。

## 睇埋
* 歐拉法
* Runge-Kutta法一覽
* 常微分方程嘅數值方法
* 隨機微分方程嘅Runge-Kutta法
* 一般線性方法
* 李羣積分器

## 註

## 考

*  。
*  。
* 978-0-89871-412-8 。
* 978-0-471-50023-0 。
*  。
*  。
* 978-0-470-72335-7 。
* 0-387-26102-8 。
*  。
*  （請參閱第6章）。
* 978-3-540-56670-0 。
* 978-3-540-60452-5 。
* 978-0-521-55655-2 。
* 0-471-92990-5
*  。
* 978-0-521-88068-8 。另外，第17.2節。 Runge-Kutta嘅自適應𨂾長控制。
* 978-0-387-95452-3 。
* 0-521-00794-1 。
*  。
* 高級離散數學ignou參考書（代碼mcs033）

## 連出去

*  
* 龍格-庫塔四階方法
* Matlab裏便嘅跟踪器組件庫實現`RungeKStep`實現32種嵌入式Runge Kutta算法，喺RungeKNystroemSStep裏便實現24種嵌入式Runge- `RungeKNystroemSStep` 4種通用Runge- `RungeKNystroemGStep`算法。