經典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})}
嘅RK4近似值 ,同埋下一個值(
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}
,即微分方程相當於一個簡單嘅積分,係噉RK4就係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}}
個最終公式用到。
顯式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}
仲係一個開放問題(open problem)、對於一個顯式
p
{\displaystyle p}
階Runge-Kutta法嚟講,因爲仲未發現啲方法、滿足得上述啲範圍入便啲等號嘅。一些已知嘅值係: [ 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}
之只有8級嘅方法,之不過如表所示而今已知嘅方法至少要到9級。
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] 呢個方法嘅主要優點係,幾乎所有啲誤差係數都細過啲流行方法入便嘅誤差係數,但每一𨂾時間需要稍爲多啲嘅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
+
h
f
(
t
n
,
y
n
)
{\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n})}
。呢個係唯一相容(consistent)嘅一級顯式Runge-Kutta法。個相應圖表係:
中點法係一個兩級二階方法嘅示例:
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}}.}
數值解對應於帶下劃線啲值。
若果
c
i
,
i
=
1
,
2
,
…
,
s
{\displaystyle c_{i},\,i=1,2,\ldots ,s}
互相之間唔同,噉Runge-Kutta法就着認爲係非融合嘅 (nonconfluent )[ 15] 。
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).}
到目前為止,所有啲提到嘅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] 兩級(即有四階)嘅方法會有噉樣嘅Butcher表:
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' + 2嗰陣,一個Padé近似值有m 次嘅分子同埋n 次嘅分母嘅先係A穩定嘅。 [ 26]
s 級嘅Gauss-Legendre法嘅階數為2 s ,因此佢個穩定性函數為m = n = s 嘅Padé近似值。因此,呢個方法係A穩定嘅。 [ 27] 噉樣表明A穩定嘅Runge-Kutta可以有任意高嘅階。相反,A穩定線性多步法 嘅階數唔超得過2。
微分方程解嘅A穩定性概念同線性自治方程有關
y
′
=
λ
y
{\displaystyle y'=\lambda y}
。 Dahlquist進行唨一項研究針對啲數值方案嘅穩定性、啲應用於滿足單調性條件嘅非線性系統嘅。相應嘅概念被定義為: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穩定。
設
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}
矩陣都係非負定嘅,噉呢個Runge-Kutta法就可以講係代數穩定嘅 [ 28] 。 B穩定性[ 29] 嘅充分條件係:
B
{\displaystyle B}
跟
Q
{\displaystyle Q}
係非負定嘅。
一般來講,
s
{\displaystyle s}
階Runge-Kutta法寫得成:
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}
代入通用公式可以得到四階Runge-Kutta方法[ 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}
個總導數(total derivative)。
若果而今使用啱啱導出嘅內容來表達通用公式,噉會得到:
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}
and
c
1
{\displaystyle c_{1}}
and
a
21
{\displaystyle a_{21}}
randomly chosen. This method is consistent and (in general) first-order convergent. On the other hand, the 1-stage method with
b
1
=
1
/
2
{\displaystyle b_{1}=1/2}
is inconsistent and fails to converge, even though it trivially holds that
∑
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–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-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 (請參閱第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-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 (第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-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)