喺數值分析 當中， Runge-Kutta法 （德文 ：Runge-Kutta-Verfahren ，英文 ：Runge-Kutta methods ）、又可以叫龍俄-厥他法 （粵拼 ：lung4 ngo4 - kyut3 taa1 faat3 RK法 ，查實係一系列隱式同顯式 迭代方法，其中包括眾所周知嘅喊做Euler法 嘅例程，呢個例程用於幫常微分方程 嘅近似解做時間離散化 。 [1] Carl Runge 跟Wilhelm Kutta 於1900年左右開發出嚟嘅。
微分方程y'= sin（t）^ 2 * y嘅Runge-Kutta法嘅比較（橙色為精確解） 經典RK法 [ 編輯 ] 經典Runge-Kutta法使用嘅斜率 Runge-Kutta系列裏便最廣為人知嘅一種通常着喊做「 RK4」、「經典Runge-Kutta法」抑或直程「 Runge-Kutta法」，「4」係因爲有四階。
指定一個初值問題 ：
d
y
d
t
=
f
(
t
,
y
)
,
y
(
t
0
)
=
y
0
.
{\displaystyle {\frac {dy}{dt}}=f(t,y),\quad y(t_{0})=y_{0}.}
其中
y
{\displaystyle y}
t
{\displaystyle t}
y
{\displaystyle y}
d
y
d
t
{\displaystyle {\frac {dy}{dt}}}
t
{\displaystyle t}
y
{\displaystyle y}
t
0
{\displaystyle t_{0}}
y
{\displaystyle y}
y
0
{\displaystyle y_{0}}
f
{\displaystyle f}
初始條件
t
0
{\displaystyle t_{0}}
y
0
{\displaystyle y_{0}}
而今揀一個𨂾長h > 0並定義埋：
y
n
+
1
=
y
n
+
1
6
h
(
k
1
+
2
k
2
+
2
k
3
+
k
4
)
,
t
n
+
1
=
t
n
+
h
{\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}}}
對於n = 0、1、2、3，...，使用[2]
k
1
=
f
(
t
n
,
y
n
)
,
k
2
=
f
(
t
n
+
h
2
,
y
n
+
h
k
1
2
)
,
k
3
=
f
(
t
n
+
h
2
,
y
n
+
h
k
2
2
)
,
k
4
=
f
(
t
n
+
h
,
y
n
+
h
k
3
)
.
{\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] 呢度嘅
y
n
+
1
{\displaystyle y_{n+1}}
y
(
t
n
+
1
)
{\displaystyle y(t_{n+1})}
y
n
+
1
{\displaystyle y_{n+1}}
y
n
{\displaystyle y_{n}}
加權平均值 嚟決定，其中每個增量都係一個乘積、由間隔大細h 乘埋一個估出嘅斜率得出嘅，個斜率k 由微分方程右側嘅函數 f 指定。
k
1
{\displaystyle k_{1}}
y
{\displaystyle y}
歐拉法 ）;
k
2
{\displaystyle k_{2}}
y
{\displaystyle y}
k
1
{\displaystyle k_{1}}
k
3
{\displaystyle k_{3}}
y
{\displaystyle y}
k
2
{\displaystyle k_{2}}
k
4
{\displaystyle k_{4}}
y
{\displaystyle y}
k
3
{\displaystyle k_{3}}
喺對四個斜率求平均值陣時，對中點度嘅斜率賦予大啲嘅權重。若果
f
{\displaystyle f}
y
{\displaystyle y}
Simpson法則 [4]
種RK4方法係四階方法，噉樣代表局部截尾誤差同階於
O
(
h
5
)
{\displaystyle O(h^{5})}
總累積誤差 同階於
O
(
h
4
)
{\displaystyle O(h^{4})}
喺許多實際應用裏便，隻函數
f
{\displaystyle f}
t
{\displaystyle t}
自治系統 （autonomous system）抑或時唔變系統，特別係喺物理學裏便），而且佢哋嘅增量根本唔會得到計算，亦都唔會傳遞畀函數
f
{\displaystyle f}
t
n
+
1
{\displaystyle t_{n+1}}
顯式RK法 [ 編輯 ] 顯式Runge-Kutta法家族係上述RK4方法嘅概括，之係噉樣得出：
y
n
+
1
=
y
n
+
h
∑
i
=
1
s
b
i
k
i
,
{\displaystyle y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},}
其中[5]
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
)
)
,
⋮
k
s
=
f
(
t
n
+
c
s
h
,
y
n
+
h
(
a
s
1
k
1
+
a
s
2
k
2
+
⋯
+
a
s
,
s
−
1
k
s
−
1
)
)
.
{\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}}}
（注意：以上公式喺某些文本裏便可能具有唔同但等效嘅定義）。 要指定一個特定嘅方法，就要提供整數s （級數）、係數aij （對於 1 ≤ j < i ≤ s ）、bi （對於i = 1，2，...， s ）同埋ci （對於i = 2，3，...， s ）。矩陣[aij ]着喊做Runge-Kutta矩陣 ，而 bi 、ci 着喊做權重 、節點 。 [6] Butcher表 （Butcher tableau ）嘅形式嚟排列（根據返John C. Butcher個名）：
0
{\displaystyle 0}
c
2
{\displaystyle c_{2}}
a
21
{\displaystyle a_{21}}
c
3
{\displaystyle c_{3}}
a
31
{\displaystyle a_{31}}
a
32
{\displaystyle a_{32}}
⋮
{\displaystyle \vdots }
⋮
{\displaystyle \vdots }
⋱
{\displaystyle \ddots }
c
s
{\displaystyle c_{s}}
a
s
1
{\displaystyle a_{s1}}
a
s
2
{\displaystyle a_{s2}}
⋯
{\displaystyle \cdots }
a
s
,
s
−
1
{\displaystyle a_{s,s-1}}
b
1
{\displaystyle b_{1}}
b
2
{\displaystyle b_{2}}
⋯
{\displaystyle \cdots }
b
s
−
1
{\displaystyle b_{s-1}}
b
s
{\displaystyle b_{s}}
泰勒級數 展開表明，Runge-Kutta法相容（consistent）若且唯若有：
∑
i
=
1
s
b
i
=
1.
{\displaystyle \sum _{i=1}^{s}b_{i}=1.}
若果要求呢個方法有返一定階數p ，噉亦都有相應嘅要求，噉樣意味住局部截尾誤差為O(hp +1 )。啲可以從截尾誤差本身嘅定義裏便得出。譬如，若果b 1 + b 2 = 1， b 2 c 2 = 1/2， b 2 a 21 = 1/2，噉兩級方法嘅階數為2。 [7] [8]
∑
j
=
1
i
−
1
a
i
j
=
c
i
for
i
=
2
,
…
,
s
.
{\displaystyle \sum _{j=1}^{i-1}a_{ij}=c_{i}{\text{ for }}i=2,\ldots ,s.}
之不過，僅衹呢個條件唔夠保持相容性。 [9]
一般嚟講，若果
s
{\displaystyle s}
p
{\displaystyle p}
s
≥
p
{\displaystyle s\geq p}
p
≥
5
{\displaystyle p\geq 5}
s
≥
p
+
1
{\displaystyle s\geq p+1}
[10] 清晰 （sharp ）。譬如，所有啲已知嘅8階方法至少有11級，即使有可能存在啲方法少啲級嘅。 （上便個範圍表明到可能係有9級嘅方法；但亦都可能係個範圍根本唔清晰。 ）實質上，確切嘅最細級數
s
{\displaystyle s}
p
{\displaystyle p}
[11]
p
1
2
3
4
5
6
7
8
min
s
1
2
3
4
6
7
9
11
{\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}}}
上便啲證明得嘅範圍就意味住搵唔到階數
p
=
1
,
2
,
…
,
6
{\displaystyle p=1,2,\ldots ,6}
p
=
7
{\displaystyle p=7}
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
種Runge-Kutta法嘅其中一種熹微變體亦都歸功於1901年Kutta，喊做3/8規則。[13]
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
+
h
f
(
t
n
,
y
n
)
{\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n})}
兩級二階方法 [ 編輯 ] 中點法係一個兩級二階方法嘅示例：
y
n
+
1
=
y
n
+
h
f
(
t
n
+
1
2
h
,
y
n
+
1
2
h
f
(
t
n
,
y
n
)
)
.
{\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).}
相應嘅圖表係：
中點法唔係唯一嘅兩級二階Runge-Kutta法。有一系列噉樣嘅方法，由α參數設定並由以下公式畀出[14]
y
n
+
1
=
y
n
+
h
(
(
1
−
1
2
α
)
f
(
t
n
,
y
n
)
+
1
2
α
f
(
t
n
+
α
h
,
y
n
+
α
h
f
(
t
n
,
y
n
)
)
)
.
{\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 )}.}
個Butcher表係：
0
α
{\displaystyle \alpha }
α
{\displaystyle \alpha }
(
1
−
1
2
α
)
{\displaystyle (1-{\tfrac {1}{2\alpha }})}
1
2
α
{\displaystyle {\tfrac {1}{2\alpha }}}
喺噉樣嘅方法族入便，
α
=
1
2
{\displaystyle \alpha ={\tfrac {1}{2}}}
α
=
1
{\displaystyle \alpha =1}
Heun法 。 [4]
譬如，考慮α= 2/3嘅兩級二階Runge-Kutta法，之亦都喊做Ralston方法 。佢由表格畀出：
用到相應嘅方程式
k
1
=
f
(
t
n
,
y
n
)
,
k
2
=
f
(
t
n
+
2
3
h
,
y
n
+
2
3
h
k
1
)
,
y
n
+
1
=
y
n
+
h
(
1
4
k
1
+
3
4
k
2
)
.
{\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}}}
呢個方法用於解決初值問題
d
y
d
t
=
tan
(
y
)
+
1
,
y
0
=
1
,
t
∈
[
1
,
1.1
]
{\displaystyle {\frac {dy}{dt}}=\tan(y)+1,\quad y_{0}=1,\ t\in [1,1.1]}
其中𨂾長h = 0.025，因此呢個方法需要執行四𨂾。
呢個方法進行如下：
t
0
=
1
:
{\displaystyle t_{0}=1\colon }
y
0
=
1
{\displaystyle y_{0}=1}
t
1
=
1.025
:
{\displaystyle t_{1}=1.025\colon }
y
0
=
1
{\displaystyle y_{0}=1}
k
1
=
2.557407725
{\displaystyle k_{1}=2.557407725}
k
2
=
f
(
t
0
+
2
3
h
,
y
0
+
2
3
h
k
1
)
=
2.7138981400
{\displaystyle k_{2}=f(t_{0}+{\tfrac {2}{3}}h,\ y_{0}+{\tfrac {2}{3}}hk_{1})=2.7138981400}
y
1
=
y
0
+
h
(
1
4
k
1
+
3
4
k
2
)
=
1.066869388
_
{\displaystyle y_{1}=y_{0}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.066869388}}}
t
2
=
1.05
:
{\displaystyle t_{2}=1.05\colon }
y
1
=
1.066869388
{\displaystyle y_{1}=1.066869388}
k
1
=
2.813524695
{\displaystyle k_{1}=2.813524695}
k
2
=
f
(
t
1
+
2
3
h
,
y
1
+
2
3
h
k
1
)
{\displaystyle k_{2}=f(t_{1}+{\tfrac {2}{3}}h,\ y_{1}+{\tfrac {2}{3}}hk_{1})}
y
2
=
y
1
+
h
(
1
4
k
1
+
3
4
k
2
)
=
1.141332181
_
{\displaystyle y_{2}=y_{1}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.141332181}}}
t
3
=
1.075
:
{\displaystyle t_{3}=1.075\colon }
y
2
=
1.141332181
{\displaystyle y_{2}=1.141332181}
k
1
=
3.183536647
{\displaystyle k_{1}=3.183536647}
k
2
=
f
(
t
2
+
2
3
h
,
y
2
+
2
3
h
k
1
)
{\displaystyle k_{2}=f(t_{2}+{\tfrac {2}{3}}h,\ y_{2}+{\tfrac {2}{3}}hk_{1})}
y
3
=
y
2
+
h
(
1
4
k
1
+
3
4
k
2
)
=
1.227417567
_
{\displaystyle y_{3}=y_{2}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.227417567}}}
t
4
=
1.1
:
{\displaystyle t_{4}=1.1\colon }
y
3
=
1.227417567
{\displaystyle y_{3}=1.227417567}
k
1
=
3.796866512
{\displaystyle k_{1}=3.796866512}
k
2
=
f
(
t
3
+
2
3
h
,
y
3
+
2
3
h
k
1
)
{\displaystyle k_{2}=f(t_{3}+{\tfrac {2}{3}}h,\ y_{3}+{\tfrac {2}{3}}hk_{1})}
y
4
=
y
3
+
h
(
1
4
k
1
+
3
4
k
2
)
=
1.335079087
_
.
{\displaystyle y_{4}=y_{3}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.335079087}}.}
數值解對應於帶下劃線啲值。
自適應RK法 [ 編輯 ] 自適應方法旨喺產生單個Runge-Kutta𨂾步嘅局部截尾誤差嘅估計。噉樣可以通過兩種方法來完成，一種係有序嘅
p
{\displaystyle p}
p
−
1
{\displaystyle p-1}
喺積分期間，𨂾長會着調整，嚟令估計嘅誤差保持喺用戶定義嘅閾值以下：若果誤差太大，噉攞細啲嘅𨂾長重複執行過一𨂾；若果誤差細得多，就增加𨂾長嚟慳時間。噉樣（幾乎）可以畀出最佳𨂾長，同時慳返計算時間。而且，用戶唔使嘥時間來搵返啱嘅𨂾長。
低階𨂾步由以下式得出：
y
n
+
1
∗
=
y
n
+
h
∑
i
=
1
s
b
i
∗
k
i
,
{\displaystyle y_{n+1}^{*}=y_{n}+h\sum _{i=1}^{s}b_{i}^{*}k_{i},}
其中
k
i
{\displaystyle k_{i}}
e
n
+
1
=
y
n
+
1
−
y
n
+
1
∗
=
h
∑
i
=
1
s
(
b
i
−
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},}
噉樣係
O
(
h
p
)
{\displaystyle O(h^{p})}
b
i
∗
{\displaystyle b_{i}^{*}}
0
c
2
{\displaystyle c_{2}}
a
21
{\displaystyle a_{21}}
c
3
{\displaystyle c_{3}}
a
31
{\displaystyle a_{31}}
a
32
{\displaystyle a_{32}}
⋮
{\displaystyle \vdots }
⋮
{\displaystyle \vdots }
⋱
{\displaystyle \ddots }
c
s
{\displaystyle c_{s}}
a
s
1
{\displaystyle a_{s1}}
a
s
2
{\displaystyle a_{s2}}
⋯
{\displaystyle \cdots }
a
s
,
s
−
1
{\displaystyle a_{s,s-1}}
b
1
{\displaystyle b_{1}}
b
2
{\displaystyle b_{2}}
⋯
{\displaystyle \cdots }
b
s
−
1
{\displaystyle b_{s-1}}
b
s
{\displaystyle b_{s}}
b
1
∗
{\displaystyle b_{1}^{*}}
b
2
∗
{\displaystyle b_{2}^{*}}
⋯
{\displaystyle \cdots }
b
s
−
1
∗
{\displaystyle b_{s-1}^{*}}
b
s
∗
{\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
但係，最簡單嘅自適應Runge-Kutta法涉及Heun法 （階數2）佮Euler法 （階數1）嘅組合。佢嘅擴展後嘅Butcher表係：
其他啲自適應嘅Runge-Kutta法有Bogacki-Shampine法 （3、2階）同埋Cash-Karp法 跟Dormand-Prince法 （5、4階）。
非融合RK法 [ 編輯 ] 若果
c
i
,
i
=
1
,
2
,
…
,
s
{\displaystyle c_{i},\,i=1,2,\ldots ,s}
融合嘅 （nonconfluent ）[15]
Runge-Kutta-Nyström法 [ 編輯 ] Runge-Kutta-Nyström法係特定嘅Runge-Kutta法，之針對二階微分方程進行唨改良[16] [17]
d
2
y
d
t
2
=
f
(
y
,
y
˙
,
t
)
.
{\displaystyle {\frac {d^{2}y}{dt^{2}}}=f(y,{\dot {y}},t).}
隱式RK法 [ 編輯 ] 到目前為止，所有啲提到嘅Runge-Kutta法都係顯式方法 。顯式嘅Runge-Kutta法通常唔適合於攞嚟求解啲剛性方程 （Stiff equations），因為佢哋嘅絕對穩定域好細。特別係佢係有界嘅。 [18] 偏微分方程求解 裏便尤為重要。
顯式Runge-Kutta法個唔穩定性激唨啲隱式方法嘅發展。隱式Runge-Kutta法具有以下形式：
y
n
+
1
=
y
n
+
h
∑
i
=
1
s
b
i
k
i
,
{\displaystyle y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},}
其中
k
i
=
f
(
t
n
+
c
i
h
,
y
n
+
h
∑
j
=
1
s
a
i
j
k
j
)
,
i
=
1
,
…
,
s
.
{\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] 顯式方法嘅區別喺於，喺顯式方法裏便， j 嘅總和僅衹等於i -1。呢個亦都顯示喺張Butcher表裏便，當中顯式方法個係數矩陣
a
i
j
{\displaystyle a_{ij}}
j 度嘅總和上升到s 而且係數矩陣唔係三角形，係噉畀出對應嘅Butcher表形式[12]
c
1
a
11
a
12
…
a
1
s
c
2
a
21
a
22
…
a
2
s
⋮
⋮
⋮
⋱
⋮
c
s
a
s
1
a
s
2
…
a
s
s
b
1
b
2
…
b
s
b
1
∗
b
2
∗
…
b
s
∗
=
c
A
b
T
{\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}}}
請參閱返上方嘅自適應Runge-Kutta法 嚟獲取有關
b
∗
{\displaystyle b^{*}}
呢種差異導致，每一𨂾都必須求解代數方程組。噉樣大大增加唨計算成本。若果使用s 級嘅方法求解m 個分量嘅微分方程，噉個代數方程組將具有ms 個分量。噉樣可以戥隱式線性多步法 （ODE嘅另一種大方法系列）形成對比：隱式s 𨂾線性多步方法需要求解只有m 個分量嘅代數方程組，因此系統嘅大細唔會隨住𨂾數增加而增加。 [20]
隱式Runge-Kutta法嘅最簡單示例係後褪Euler法 ：
y
n
+
1
=
y
n
+
h
f
(
t
n
+
h
,
y
n
+
1
)
.
{\displaystyle y_{n+1}=y_{n}+hf(t_{n}+h,\ y_{n+1}).\,}
張Butcher表好簡單：
1
1
1
{\displaystyle {\begin{array}{c|c}1&1\\\hline &1\\\end{array}}}
呢張Butcher表戥以下公式相對應
k
1
=
f
(
t
n
+
h
,
y
n
+
h
k
1
)
and
y
n
+
1
=
y
n
+
h
k
1
,
{\displaystyle 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表為：
0
0
0
1
1
2
1
2
1
2
1
2
1
0
{\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}}}
梯形規則係一種選點法 （似文中所述）。所有選點法都係隱式Runge-Kutta法，但並非所有隱式Runge-Kutta法都係選點法。 [21]
Gauss-Legendre法 構成唨一系列基於Gauss正交 嘅搭配方法。s 級嘅Gauss-Legendre法嘅階數為2 s （因此，可以構造任意高階嘅方法）。 [22]
1
2
−
1
6
3
1
4
1
4
−
1
6
3
1
2
+
1
6
3
1
4
+
1
6
3
1
4
1
2
1
2
1
2
+
1
2
3
1
2
−
1
2
3
{\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] 戥顯式方法相比，隱式Runge-Kutta法嘅優勢喺於其更大嘅穩定性，尤其係喺應用於剛性方程 嗰陣。考慮線性測試方程y' = λy 。 應用到呢條方程嘅Runge-Kutta法簡化成迭代式
y
n
+
1
=
r
(
h
λ
)
y
n
{\displaystyle y_{n+1}=r(h\lambda )\,y_{n}}
r 由下式得出：
r
(
z
)
=
1
+
z
b
T
(
I
−
z
A
)
−
1
e
=
det
(
I
−
z
A
+
z
e
b
T
)
det
(
I
−
z
A
)
,
{\displaystyle r(z)=1+zb^{T}(I-zA)^{-1}e={\frac {\det(I-zA+zeb^{T})}{\det(I-zA)}},}
[23] 其中e 代表1嘅向量。函數r 喊做穩定性函數 。 [24] s 級， 噉r 係兩個次數為s 嘅多項式個商。 顯式方法有嚴格意義上嘅下三角矩陣A ，噉樣意味住det（ I - zA ）= 1而且穩定性函數係多項式。 [25]
若果|r (z )| < 1，z = h λ ，線性測試方程嘅數值解將衰減到零。噉樣嘅z 嘅集合喊做絕對穩定性域 。特別嚟講，若果啲Re(z ) < 0嘅z 係喺嗮絕對穩定域，呢個方法就着認為係絕對穩定嘅。顯式Runge-Kutta法嘅穩定性函數係多項式，因此顯式Runge-Kutta法永遠唔會係A穩定 嘅。 [25]
若果方法嘅階數為p ，噉穩定性函數滿足
z
→
0
{\displaystyle z\to 0}
r
(
z
)
=
e
z
+
O
(
z
p
+
1
)
{\displaystyle r(z)={\textrm {e}}^{z}+O(z^{p+1})}
Padé近似值 。若且唯若'm≤N≤M' m 次嘅分子同埋n 次嘅分母嘅先係A穩定嘅。 [26]
s 級嘅Gauss-Legendre法嘅階數為2 s ，因此佢個穩定性函數為m = n = s 嘅Padé近似值。因此，呢個方法係A穩定嘅。 [27] 線性多步法 嘅階數唔超得過2。
B穩定性 [ 編輯 ] 微分方程解嘅A穩定性概念同線性自治方程有關
y
′
=
λ
y
{\displaystyle y'=\lambda y}
G穩定性 同埋B穩定性 （Butcher，1975），分別對應多步方法（同埋戥佢相關嘅啲單腿方法（one-leg methods））同埋Runge-Kutta法。 Runge-Kutta法應用於非線性系統
y
′
=
f
(
y
)
{\displaystyle y'=f(y)}
⟨
f
(
y
)
−
f
(
z
)
,
y
−
z
⟩
<
0
{\displaystyle \langle f(y)-f(z),\ y-z\rangle <0}
‖
y
n
+
1
−
z
n
+
1
‖
≤
‖
y
n
−
z
n
‖
{\displaystyle \|y_{n+1}-z_{n+1}\|\leq \|y_{n}-z_{n}\|}
設
B
{\displaystyle B}
M
{\displaystyle M}
Q
{\displaystyle Q}
s
×
s
{\displaystyle s\times s}
B
=
diag
(
b
1
,
b
2
,
…
,
b
s
)
,
M
=
B
A
+
A
T
B
−
b
b
T
,
Q
=
B
A
−
1
+
A
−
T
B
−
A
−
T
b
b
T
A
−
1
.
{\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}.}
若果
B
{\displaystyle B}
M
{\displaystyle M}
代數穩定嘅 [28] [29]
B
{\displaystyle B}
Q
{\displaystyle Q}
四階RK法嘅推導 [ 編輯 ] 一般來講，
s
{\displaystyle s}
y
t
+
h
=
y
t
+
h
⋅
∑
i
=
1
s
a
i
k
i
+
O
(
h
s
+
1
)
,
{\displaystyle 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
⋅
∑
j
=
1
s
β
i
j
f
(
k
j
,
t
n
+
α
i
h
)
{\displaystyle k_{i}=y_{t}+h\cdot \sum _{j=1}^{s}\beta _{ij}f\left(k_{j},\ t_{n}+\alpha _{i}h\right)}
係通過評估
y
t
{\displaystyle y_{t}}
i
{\displaystyle i}
如前所述，喺任一間隔
(
t
,
t
+
h
)
{\displaystyle (t,\ t+h)}
s
=
4
{\displaystyle s=4}
[30]
α
i
β
i
j
α
1
=
0
β
21
=
1
2
α
2
=
1
2
β
32
=
1
2
α
3
=
1
2
β
43
=
1
α
4
=
1
{\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}}}
否則
β
i
j
=
0
{\displaystyle \beta _{ij}=0}
y
t
+
h
1
=
y
t
+
h
f
(
y
t
,
t
)
y
t
+
h
2
=
y
t
+
h
f
(
y
t
+
h
/
2
1
,
t
+
h
2
)
y
t
+
h
3
=
y
t
+
h
f
(
y
t
+
h
/
2
2
,
t
+
h
2
)
{\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}}}
其中
y
t
+
h
/
2
1
=
y
t
+
y
t
+
h
1
2
{\displaystyle y_{t+h/2}^{1}={\dfrac {y_{t}+y_{t+h}^{1}}{2}}}
y
t
+
h
/
2
2
=
y
t
+
y
t
+
h
2
2
{\displaystyle y_{t+h/2}^{2}={\dfrac {y_{t}+y_{t+h}^{2}}{2}}}
k
1
=
f
(
y
t
,
t
)
k
2
=
f
(
y
t
+
h
/
2
1
,
t
+
h
2
)
=
f
(
y
t
+
h
2
k
1
,
t
+
h
2
)
k
3
=
f
(
y
t
+
h
/
2
2
,
t
+
h
2
)
=
f
(
y
t
+
h
2
k
2
,
t
+
h
2
)
k
4
=
f
(
y
t
+
h
3
,
t
+
h
)
=
f
(
y
t
+
h
k
3
,
t
+
h
)
{\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}}}
對於先前嘅關係，可以證明以下等式滿足得
O
(
h
2
)
{\displaystyle {\mathcal {O}}(h^{2})}
k
2
=
f
(
y
t
+
h
/
2
1
,
t
+
h
2
)
=
f
(
y
t
+
h
2
k
1
,
t
+
h
2
)
=
f
(
y
t
,
t
)
+
h
2
d
d
t
f
(
y
t
,
t
)
k
3
=
f
(
y
t
+
h
/
2
2
,
t
+
h
2
)
=
f
(
y
t
+
h
2
f
(
y
t
+
h
2
k
1
,
t
+
h
2
)
,
t
+
h
2
)
=
f
(
y
t
,
t
)
+
h
2
d
d
t
[
f
(
y
t
,
t
)
+
h
2
d
d
t
f
(
y
t
,
t
)
]
k
4
=
f
(
y
t
+
h
3
,
t
+
h
)
=
f
(
y
t
+
h
f
(
y
t
+
h
2
k
2
,
t
+
h
2
)
,
t
+
h
)
=
f
(
y
t
+
h
f
(
y
t
+
h
2
f
(
y
t
+
h
2
f
(
y
t
,
t
)
,
t
+
h
2
)
,
t
+
h
2
)
,
t
+
h
)
=
f
(
y
t
,
t
)
+
h
d
d
t
[
f
(
y
t
,
t
)
+
h
2
d
d
t
[
f
(
y
t
,
t
)
+
h
2
d
d
t
f
(
y
t
,
t
)
]
]
{\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}}}
其中：
d
d
t
f
(
y
t
,
t
)
=
∂
∂
y
f
(
y
t
,
t
)
y
˙
t
+
∂
∂
t
f
(
y
t
,
t
)
=
f
y
(
y
t
,
t
)
y
˙
+
f
t
(
y
t
,
t
)
:=
y
¨
t
{\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}}
係對應時間嘅
f
{\displaystyle f}
若果而今使用啱啱導出嘅內容來表達通用公式，噉會得到：
y
t
+
h
=
y
t
+
h
{
a
⋅
f
(
y
t
,
t
)
+
b
⋅
[
f
(
y
t
,
t
)
+
h
2
d
d
t
f
(
y
t
,
t
)
]
+
+
c
⋅
[
f
(
y
t
,
t
)
+
h
2
d
d
t
[
f
(
y
t
,
t
)
+
h
2
d
d
t
f
(
y
t
,
t
)
]
]
+
+
d
⋅
[
f
(
y
t
,
t
)
+
h
d
d
t
[
f
(
y
t
,
t
)
+
h
2
d
d
t
[
f
(
y
t
,
t
)
+
h
2
d
d
t
f
(
y
t
,
t
)
]
]
]
}
+
O
(
h
5
)
=
y
t
+
a
⋅
h
f
t
+
b
⋅
h
f
t
+
b
⋅
h
2
2
d
f
t
d
t
+
c
⋅
h
f
t
+
c
⋅
h
2
2
d
f
t
d
t
+
+
c
⋅
h
3
4
d
2
f
t
d
t
2
+
d
⋅
h
f
t
+
d
⋅
h
2
d
f
t
d
t
+
d
⋅
h
3
2
d
2
f
t
d
t
2
+
d
⋅
h
4
4
d
3
f
t
d
t
3
+
O
(
h
5
)
{\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}}}
並捉佢戥
y
t
{\displaystyle y_{t}}
y
t
+
h
{\displaystyle y_{t+h}}
泰勒級數 比較，有 ：
y
t
+
h
=
y
t
+
h
y
˙
t
+
h
2
2
y
¨
t
+
h
3
6
y
t
(
3
)
+
h
4
24
y
t
(
4
)
+
O
(
h
5
)
=
=
y
t
+
h
f
(
y
t
,
t
)
+
h
2
2
d
d
t
f
(
y
t
,
t
)
+
h
3
6
d
2
d
t
2
f
(
y
t
,
t
)
+
h
4
24
d
3
d
t
3
f
(
y
t
,
t
)
{\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}}}
獲有一個係數約束系統：
{
a
+
b
+
c
+
d
=
1
1
2
b
+
1
2
c
+
d
=
1
2
1
4
c
+
1
2
d
=
1
6
1
4
d
=
1
24
{\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}}}
個系統解出陣時畀出
a
=
1
6
,
b
=
1
3
,
c
=
1
3
,
d
=
1
6
{\displaystyle a={\frac {1}{6}},b={\frac {1}{3}},c={\frac {1}{3}},d={\frac {1}{6}}}
↑ DEVRIES, Paul L. ; HASBUN, Javier E. A first course in computational physics. Second edition. Jones and Bartlett Publishers: 2011. p. 215.
↑ Press et al. 2007 , p. 908; Süli & Mayers 2003 , p. 328↑ 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.0 4.1 Süli & Mayers 2003 , p. 328↑ Press et al. 2007 , p. 907↑ Iserles 1996 , p. 38↑ Iserles 1996 , p. 39↑ Iserles 1996 , p. 39↑
As a counterexample, consider any explicit 2-stage Runge-Kutta scheme with
b
1
=
b
2
=
1
/
2
{\displaystyle b_{1}=b_{2}=1/2}
c
1
{\displaystyle c_{1}}
a
21
{\displaystyle a_{21}}
b
1
=
1
/
2
{\displaystyle b_{1}=1/2}
∑
j
=
1
i
−
1
a
i
j
=
c
i
for
i
=
2
,
…
,
s
.
{\displaystyle \sum _{j=1}^{i-1}a_{ij}=c_{i}{\text{ for }}i=2,\ldots ,s.}
↑ Butcher 2008 , p. 187↑ Butcher 2008 ↑ 12.0 12.1 Süli & Mayers 2003 , p. 352↑ Hairer, Nørsett & Wanner (1993 , p. 138) refer to Kutta (1901) .↑ Süli & Mayers 2003 , p. 327↑ Lambert 1991 , p. 278↑ 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 . ↑ 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. ↑ Süli & Mayers 2003 , pp. 349–351↑ Iserles 1996 , p. 41; Süli & Mayers 2003 ↑ 20.0 20.1 Süli & Mayers 2003 , p. 353↑ Iserles 1996 ↑ Iserles 1996 , p. 47↑ Hairer & Wanner 1996 ↑ Hairer & Wanner 1996 , p. 40↑ 25.0 25.1 Iserles 1996 , p. 60↑ Iserles 1996 ↑ Iserles 1996 , p. 63↑ Lambert 1991 , p. 275↑ Lambert 1991 , p. 274↑ PDF reporting this derivation
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–453Ascher, 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-8 978-0-89871-412-8 Atkinson, Kendall A. (1989), An Introduction to Numerical Analysis (第2版), New York: John Wiley & Sons , ISBN 978-0-471-50023-0 978-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-7 978-0-470-72335-7 Cellier, F.; Kofman, E. (2006), Continuous System Simulation , Springer Verlag , ISBN 0-387-26102-8 0-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 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-0 978-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-5 978-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-2 978-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-5 0-471-92990-5 Kaw, Autar; Kalu, Egwu (2008), Numerical Methods with Applications 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-8 978-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-3 978-0-387-95452-3 Süli, Endre; Mayers, David (2003), An Introduction to Numerical Analysis , Cambridge University Press , ISBN 0-521-00794-1 0-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） 連出去 [ 編輯 ] 
![{\displaystyle {\frac {dy}{dt}}=\tan(y)+1,\quad y_{0}=1,\ t\in [1,1.1]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4a9f561fa8e63f2b86d6226df93dd25d29f76fc0)
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d7fab64b114fcd6b0ff7ccc1bb7d9754e3eace2d)
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c5d443e3cc50f9e0636ccb9094ff5c4c304b455d)
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d5fde6d4547411755a49dbe4fb9078743a9f1a08)
