Runge-Kutta法

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

微分方程y'= sin（t）^ 2 * y嘅Runge-Kutta法嘅比較（橙色為精確解）

經典RK法

經典Runge-Kutta法使用嘅斜率

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 > 0並定義埋：

${\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}}$

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

${\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}}$

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

呢度嘅$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就係Simpson法則。

種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{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}}$

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

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

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

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

$\sum _{i=1}^{s}b_{i}=1.$

若果要求呢個方法有返一定階數p ，噉亦都有相應嘅要求，噉樣意味住局部截尾誤差為O(h^p+1)。啲可以從截尾誤差本身嘅定義裏便得出。譬如，若果b₁ + b₂ = 1， b₂c₂ = 1/2， b₂a₂₁ = 1/2，噉兩級方法嘅階數為2。 請注意，確定係數嘅常用條件係：

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

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

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

${\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}}$

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

例子

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

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

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

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

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

0
	1

兩級二階方法

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

$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

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

$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 )}.$

個Butcher表係：

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

喺噉樣嘅方法族入便，$\alpha ={\tfrac {1}{2}}$嗰陣畀出中點法，而$\alpha =1$嗰陣係Heun法。

使法

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

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

用到相應嘅方程式

${\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}}$

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

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

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

呢個方法進行如下：

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

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

自適應RK法

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

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

低階𨂾步由以下式得出：

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

其中$k_{i}$戥高階方法相同。然之後個誤差係：

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

噉樣係$O(h^{p})$ 。幫呢種方法擴展Butcher表，嚟畀出$b_{i}^{*}$ ：

	0
	$c_{2}$	$a_{21}$
	$c_{3}$	$a_{31}$	$a_{32}$
	$\vdots$	$\vdots$		$\ddots$
	$c_{s}$	$a_{s1}$	$a_{s2}$	$\cdots$	$a_{s,s-1}$
		$b_{1}$	$b_{2}$	$\cdots$	$b_{s-1}$	$b_{s}$
		$b_{1}^{*}$	$b_{2}^{*}$	$\cdots$	$b_{s-1}^{*}$	$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

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

	0
	1	1
		1/2	1/2
		1	0

其他啲自適應嘅Runge-Kutta法有Bogacki-Shampine法（3、2階）同埋Cash-Karp法跟Dormand-Prince法（5、4階）。

非融合RK法

若果$c_{i},\,i=1,2,\ldots ,s$互相之間唔同，噉Runge-Kutta法就着認爲係非融合嘅（nonconfluent）。

Runge-Kutta-Nyström法

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

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

隱式RK法

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

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

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

其中

$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.$

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

${\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}}$

請參閱返上方嘅自適應Runge-Kutta法嚟獲取有關$b^{*}$排嘅介紹。

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

例子

隱式Runge-Kutta法嘅最簡單示例係後褪Euler法：

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

張Butcher表好簡單：

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

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

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

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

隱式Runge-Kutta法嘅另一個示例係梯形規則。佢個Butcher表為：

${\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}}$

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

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

${\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}}$

穩定

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

$r(z)=1+zb^{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)| < 1，z = hλ ，線性測試方程嘅數值解將衰減到零。噉樣嘅z嘅集合喊做絕對穩定性域。特別嚟講，若果啲Re(z) < 0嘅z係喺嗮絕對穩定域，呢個方法就着認為係絕對穩定嘅。顯式Runge-Kutta法嘅穩定性函數係多項式，因此顯式Runge-Kutta法永遠唔會係A穩定嘅。

若果方法嘅階數為p ，噉穩定性函數滿足$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 <0$；若果呢種情況意味住$\|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^{T}B-bb^{T},\,Q=BA^{-1}+A^{-T}B-A^{-T}bb^{T}A^{-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{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}}$

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

${\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}}$

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

${\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}}$

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

${\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}}$

其中：

${\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{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}}$

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

${\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}}$

獲有一個係數約束系統：

${\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法
• 一般線性方法
• 李羣積分器

考

• Runge, Carl David Tolmé (1895), "Über die numerische Auflösung von Differentialgleichungen", Mathematische Annalen, Springer, 46 (2): 167–178, doi:10.1007/BF01446807 。
• Kutta, Martin (1901), "Beitrag zur näherungsweisen Integration totaler Differentialgleichungen", Zeitschrift für Mathematik und Physik, 46: 435–453 。
• Ascher, Uri M.; Petzold, Linda R. (1998), Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, Philadelphia: Society for Industrial and Applied Mathematics, ISBN 978-0-89871-412-8978-0-89871-412-8 。
• Atkinson, Kendall A. (1989), An Introduction to Numerical Analysis (第2版), New York: John Wiley & Sons, ISBN 978-0-471-50023-0978-0-471-50023-0 。
• Butcher, John C. (May 1963), "Coefficients for the study of Runge-Kutta integration processes", Journal of the Australian Mathematical Society, 3 (2): 185–201, doi:10.1017/S1446788700027932 。
• Butcher, John C. (1975), "A stability property of implicit Runge-Kutta methods", BIT, 15 (4): 358–361, doi:10.1007/bf01931672 。
• Butcher, John C. (2008), Numerical Methods for Ordinary Differential Equations, New York: John Wiley & Sons, ISBN 978-0-470-72335-7978-0-470-72335-7 。
• Cellier, F.; Kofman, E. (2006), Continuous System Simulation, Springer Verlag, ISBN 0-387-26102-80-387-26102-8 。
• Dahlquist, Germund (1963), "A special stability problem for linear multistep methods", BIT, 3: 27–43, doi:10.1007/BF01963532, ISSN 0006-3835 。
• Forsythe, George E.; Malcolm, Michael A.; Moler, Cleve B. (1977), Computer Methods for Mathematical Computations, Prentice-Hall （請參閱第6章）。
• Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard (1993), Solving ordinary differential equations I: Nonstiff problems, Berlin, New York: Springer-Verlag, ISBN 978-3-540-56670-0978-3-540-56670-0 。
• Hairer, Ernst; Wanner, Gerhard (1996), Solving ordinary differential equations II: Stiff and differential-algebraic problems (第2版), Berlin, New York: Springer-Verlag, ISBN 978-3-540-60452-5978-3-540-60452-5 。
• Iserles, Arieh (1996), A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, ISBN 978-0-521-55655-2978-0-521-55655-2 。
• Lambert, J.D (1991), Numerical Methods for Ordinary Differential Systems. The Initial Value Problem, John Wiley & Sons, ISBN 0-471-92990-50-471-92990-5
• Kaw, Autar; Kalu, Egwu (2008), Numerical Methods with Applications (第1版), autarkaw.com 。
• Press, William H.; Teukolsky, Saul A.; Vetterling, William T.; Flannery, Brian P. (2007), "Section 17.1 Runge-Kutta Method", Numerical Recipes: The Art of Scientific Computing (第3版), Cambridge University Press, ISBN 978-0-521-88068-8978-0-521-88068-8 。另外，第17.2節。 Runge-Kutta嘅自適應𨂾長控制。
• Stoer, Josef; Bulirsch, Roland (2002), Introduction to Numerical Analysis (第3版), Berlin, New York: Springer-Verlag, ISBN 978-0-387-95452-3978-0-387-95452-3 。
• Süli, Endre; Mayers, David (2003), An Introduction to Numerical Analysis, Cambridge University Press, ISBN 0-521-00794-10-521-00794-1 。
• Tan, Delin; Chen, Zheng (2012), "On A General Formula of Fourth Order Runge-Kutta Method" (PDF), Journal of Mathematical Science & Mathematics Education, 7 (2): 1–10 。
• 高級離散數學ignou參考書（代碼mcs033）

連出去

•  Hazewinkel, Michiel, 編 (2001), "Runge-Kutta method", Encyclopaedia of Mathematics, Kluwer Academic Publishers, ISBN 978-1556080104
• 龍格-庫塔四階方法
• Matlab裏便嘅跟踪器組件庫實現RungeKStep實現32種嵌入式Runge Kutta算法，喺RungeKNystroemSStep裏便實現24種嵌入式Runge- RungeKNystroemSStep 4種通用Runge- RungeKNystroemGStep算法。