一 連 の ル ン ゲ = ク ッ タ 公 式 の 中 で 最 も 広 く 知 ら れ て い る の が 、 古 典 的 ル ン ゲ = ク ッ タ 法 ( R K 4 、 も し く は 単 に 狭 義 の ル ン ゲ = ク ッ タ 法 、 英 : t h e ( c l a s s i c a l ) R u n g e – K u t t a m e t h o d ) な ど と 呼 ば れ る 4 次 の 公 式 で あ る 。
次 の 初 期 値 問 題 を 考 え る 。
y
′
=
f
(
t
,
y
)
,
y
(
t
0
)
=
y
0
.
{\displaystyle y'=f(t,y),\quad y(t_{0})=y_{0}.}
但 し 、 y ( t ) が 近 似 的 に 求 め た い 未 知 関 数 で あ り 、 そ の t に お け る 勾 配 は f ( t , y ) に よ っ て t 及 び y ( t ) の 関 数 と し て 与 え ら れ て い る 。 時 刻 t 0 に お け る 初 期 値 は y 0 で 与 え ら れ て い る 。
今 、 時 刻 t n に お け る 値 y n = y ( t n ) が 既 知 の と き 、 十 分 に 小 さ な ス テ ッ プ 幅 h に 対 し て y n + 1 , t n + 1 を 以 下 の 式 で 与 え る と 、 y n + 1 は y ( t n + 1 ) の 4 次 精 度 の 近 似 に な っ て い る 。 こ の ス テ ッ プ を 逐 次 的 に 繰 り 返 す こ と に よ っ て 、 初 期 値 y 0 か ら 任 意 の 時 刻 t n に お け る 近 似 値 y n が 求 め ら れ る 。
y
n
+
1
=
y
n
+
h
6
(
k
1
+
2
k
2
+
2
k
3
+
k
4
)
,
t
n
+
1
=
t
n
+
h
.
{\displaystyle {\begin{aligned}y_{n+1}&=y_{n}+{h \over 6}(k_{1}+2k_{2}+2k_{3}+k_{4}),\\t_{n+1}&=t_{n}+h.\end{aligned}}}
こ こ で 、
k
1
=
f
(
t
n
,
y
n
)
,
k
2
=
f
(
t
n
+
h
2
,
y
n
+
h
2
k
1
)
,
k
3
=
f
(
t
n
+
h
2
,
y
n
+
h
2
k
2
)
,
k
4
=
f
(
t
n
+
h
,
y
n
+
h
k
3
)
{\displaystyle {\begin{aligned}k_{1}&=f\left(t_{n},y_{n}\right),\\k_{2}&=f\left(t_{n}+{h \over 2},y_{n}+{h \over 2}k_{1}\right),\\k_{3}&=f\left(t_{n}+{h \over 2},y_{n}+{h \over 2}k_{2}\right),\\k_{4}&=f\left(t_{n}+h,y_{n}+hk_{3}\right)\end{aligned}}}
で あ る 。 次 の 値 ( y n + 1 ) は 、 現 在 の 値 ( y n ) に 増 分 を 加 え た も の で あ り 、 増 分 は 勾 配 の 推 定 値 に 間 隔 h を 乗 じ た も の に な っ て い る 。 勾 配 の 推 定 値 は 、 k 1 , . . . , k 4 の 4 つ の 勾 配 の 重 み 付 け 平 均 で 求 め る 。 k 1 , . . . , k 4 の そ れ ぞ れ の 勾 配 は 、 特 定 の ( t , y ) に 対 す る f に よ っ て 与 え ら れ 、 以 下 の よ う に 解 釈 で き る 。
● k 1 は 区 間 の 最 初 t n に お け る 勾 配 で あ る ( オ イ ラ ー 法 で 用 い る 勾 配 に 一 致 す る ) 。
● k 2 は 区 間 の 中 央 t n + h / 2 に お け る 勾 配 の 近 似 値 で あ る ( 中 点 法 で 用 い る 勾 配 ) 。 計 算 に 用 い る 中 央 の y の 値 は 、 初 期 位 置 の 勾 配 k 1 を 用 い て オ イ ラ ー 法 で 推 定 す る 。
● k 3 は 区 間 の 中 央 に お け る 勾 配 の も う 一 つ の 近 似 値 で あ る 。 中 央 の y を k 2 の 値 か ら 推 定 し て 用 い る 。
● k 4 は 区 間 の 最 後 t n + h に お け る 勾 配 の 近 似 値 で あ り 、 k 3 の 値 か ら 推 定 さ れ た 最 後 の 点 の y の 値 を 用 い る 。
重 み 付 き 平 均 で は 、 中 央 の 勾 配 に 対 し て 大 き な 重 み を 用 い る 。 シ ン プ ソ ン 則 を 用 い た 平 均 と 同 等 の 形 に な る [ 3 ] 。
slope
=
k
1
+
2
k
2
+
2
k
3
+
k
4
6
.
{\displaystyle {\mbox{slope}}={\frac {k_{1}+2k_{2}+2k_{3}+k_{4}}{6}}.}
R K 4 は 4 次 の 方 法 で あ る 。 厳 密 解 と R K 4 の テ イ ラ ー 展 開 が 4 次 の 項 ま で 一 致 し 、 1 ス テ ッ プ の 推 定 誤 差 は
O
(
h
5
)
{\displaystyle O(h^{5})}
の オ ー ダ ー に な る 。 目 的 の 時 刻 の y を 求 め る の に 必 要 な ス テ ッ プ 数 は
O
(
h
−
1
)
{\displaystyle O(h^{-1})}
に な る の で 、 全 体 の 推 定 誤 差 は
O
(
h
4
)
{\displaystyle O(h^{4})}
に な る 。
前 述 の R K 4 の 一 般 化 と し て 、 以 下 の 形 式 を 持 つ s 段 の ル ン ゲ ゠ ク ッ タ 法 を 構 成 す る こ と が で き る 。 整 数 s を そ の ル ン ゲ ゠ ク ッ タ 法 の 段 数 ( s t a g e ) と い う 。
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},}
但 し [ 4 ] 、
k
i
=
f
(
t
n
+
h
c
i
,
y
n
+
h
∑
j
=
1
s
a
i
j
k
j
)
,
i
=
1
,
…
,
s
.
{\displaystyle k_{i}=f\left(t_{n}+hc_{i},y_{n}+h\sum _{j=1}^{s}a_{ij}k_{j}\right),\quad i=1,\ldots ,s.}
( 文 献 に よ っ て 、 等 価 で あ る が 上 と 異 な る 定 義 の 仕 方 を し て い る も の が あ る こ と に 注 意 す る ) [ 5 ] 。
具 体 的 な ル ン ゲ ゠ ク ッ タ 公 式 は 、 解 が で き る だ け 高 い 精 度 を 持 つ よ う に 適 切 に 選 ば れ た 係 数 a ij ( ル ン ゲ ゠ ク ッ タ 行 列 ) 、 b i ( 重 み ) 、 c j ( 節 点 ) で 指 定 さ れ る (
i
,
j
≤
s
{\displaystyle i,j\leq s}
) [ 6 ] 。 特 に
i
≤
j
{\displaystyle i\leq j}
に 対 し て
a
i
j
=
0
{\displaystyle a_{ij}=0}
を 満 た す 方 法 が 広 く 用 い ら れ 、 総 称 し て 陽 的 ル ン ゲ ゠ ク ッ タ 法 ( E R K 、 英 : e x p l i c i t R u n g e – K u t t a m e t h o d s ) と 呼 ぶ 。 そ う で な い も の を 陰 的 ル ン ゲ ゠ ク ッ タ 法 ( I R K 、 英 : i m p l i c i t R u n g e – K u t t a m e t h o d s ) と 呼 ぶ 。
近 似 値 y n を y n + 1 か ら 計 算 す る と き に 発 生 す る 誤 差 の 大 き さ が
O
(
h
p
+
1
)
{\displaystyle O(h^{p+1})}
の と き 、 そ の ル ン ゲ ゠ ク ッ タ 公 式 は p 次 精 度 を 持 つ と い い 、 p を 次 数 ( ま た は 位 数 ) と 呼 ぶ 。 p 次 の ル ン ゲ ゠ ク ッ タ 公 式 は 、 誤 差 の 大 き さ の 条 件 に 誤 差 の 表 式 を 代 入 し 、 係 数 の 条 件 を 求 め る こ と に よ っ て 得 ら れ る 。 例 え ば 、 2 段 の 陽 的 方 法 が 2 次 精 度 を 持 つ た め の 係 数 に 対 す る 条 件 は 、
b
1
+
b
2
=
1
{\displaystyle b_{1}+b_{2}=1}
,
b
2
c
2
=
1
/
2
{\displaystyle b_{2}c_{2}=1/2}
, か つ
b
2
a
21
=
1
/
2
{\displaystyle b_{2}a_{21}=1/2}
で あ る [ 7 ] 。 こ の 条 件 を 満 た す 範 囲 内 で 様 々 な 方 法 を 考 え 得 る 。
こ れ ら の 係 数 を 分 か り や す く 表 記 す る 方 法 と し て 、 以 下 の よ う な 形 式 の ブ ッ チ ャ ー 配 列 ( B u t c h e r t a b l e a u ) が 知 ら れ て い る [ 8 ] :
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
=
c
A
b
T
{\displaystyle {\begin{aligned}{\begin{array}{c|ccccc}c_{1}&a_{11}&a_{12}&\cdots &a_{1s}\\c_{2}&a_{21}&a_{22}&\cdots &a_{2s}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{s1}&a_{s2}&\cdots &a_{ss}\\\hline &b_{1}&b_{2}&\cdots &b_{s}\end{array}}&={\begin{array}{c|c}\mathbf {c} &A\\\hline &\mathbf {b} ^{\mathrm {T} }\end{array}}\end{aligned}}}
実 際 に は 、 そ れ ぞ れ の ル ン ゲ ゠ ク ッ タ 法 に つ い て 各 要 素 に 具 体 的 な 値 を 入 れ て 用 い る 。 陽 的 な 方 法 で は ル ン ゲ ゠ ク ッ タ 行 列 の 上 三 角 成 分 は 常 に 0 に な る の で 表 記 を 省 略 す る 。 例 え ば R K 4 は ブ ッ チ ャ ー 配 列 を 用 い て 以 下 の よ う に 表 現 さ れ る 。
0
1
/
2
1
/
2
1
/
2
0
1
/
2
1
0
0
1
1
/
6
1
/
3
1
/
3
1
/
6
{\displaystyle {\begin{aligned}{\begin{array}{c|ccccc}0&\\1/2&1/2\\1/2&0&1/2\\1&0&0&1\\\hline &1/6&1/3&1/3&1/6\end{array}}\end{aligned}}}
ル ン ゲ ゠ ク ッ タ 法 が 一 貫 し て い る た め に は 、 次 の 条 件 が 満 た さ れ て い る 必 要 が あ る 。
∑
j
=
1
s
a
i
j
=
c
i
,
i
=
1
,
…
,
s
.
{\displaystyle \sum _{j=1}^{s}a_{ij}=c_{i},\;i=1,\ldots ,s.}
ル ン ゲ = ク ッ タ 法 は 、 数 値 積 分 に お け る 求 積 法 ( q u a d r a t u r e ) と 深 く 繋 が る 。 時 刻 t n で の 値 か ら t n + 1 = t n + h で の 値 を 求 め る と き の 方 程 式 は 以 下 の よ う に 定 め る 。
y
′
=
f
(
t
,
y
)
,
y
(
t
n
)
=
y
n
.
{\displaystyle y'=f(t,y),\;y(t_{n})=y_{n}.}
求 積 法 は 、 与 え ら れ た 区 間 で の 定 積 分 の 値 を 被 積 分 関 数 の 値 の 線 型 結 合 と し て 近 似 す る 方 法 で あ る 。 簡 単 の た め に 、 区 間 を [ 0 , 1 ] と す る 。 よ っ て 求 積 法 の 公 式 は
∫
0
1
f
(
t
)
d
t
≈
∑
i
=
1
s
b
i
f
(
c
i
)
{\displaystyle \int _{0}^{1}f(t )dt\approx \sum _{i=1}^{s}b_{i}f(c_{i})}
と な る 。 こ こ で 、 b i と c i は 先 に 選 ば れ た 定 数 で あ り 、 前 述 の 重 み と 節 点 に 対 応 す る 。 上 記 式 に 対 し 、 等 式 が す べ て の p − 1 次 以 下 の 多 項 式 に 成 立 す る と き ︵ す な わ ち 、 誤 差 が 0 の と き ︶ 、 そ の 求 積 法 は p 次 精 度 で あ り 、 p を 次 数 と 呼 ぶ 。 節 点 が s 個 の 時 、 最 大 次 数 は 2 s で あ り 、 そ の 方 法 は s 次 ガ ウ ス ・ ル ジ ャ ン ド ル 公 式 と 呼 ば れ る 。
そ し て 上 記 の 方 程 式 を 積 分 形 式 に 変 形 し 、 求 積 法 を 用 い る と 次 の 公 式 と な る [ 6 ] 。
y
(
t
n
+
1
)
=
y
n
+
∫
y
n
y
n
+
1
f
(
t
,
y
(
t
)
)
d
t
=
y
n
+
h
∫
0
1
y
(
t
n
+
h
τ
,
y
(
t
n
+
h
τ
)
)
d
τ
=
y
n
+
h
∑
i
=
1
s
b
i
f
(
t
n
+
c
i
h
,
y
(
t
n
+
c
i
h
)
)
{\displaystyle y(t_{n+1})=y_{n}+\int _{y_{n}}^{y_{n+1}}f(t,y(t ))dt=y_{n}+h\int _{0}^{1}y(t_{n}+h\tau ,y(t_{n}+h\tau ))d\tau =y_{n}+h\sum _{i=1}^{s}b_{i}f(t_{n}+c_{i}h,y(t_{n}+c_{i}h))}
k
i
=
f
(
t
n
+
c
i
h
,
y
(
t
n
+
c
i
h
)
)
{\displaystyle k_{i}=f(t_{n}+c_{i}h,y(t_{n}+c_{i}h))}
と お く 。 k i あ る い は y ( t n + c i h ) を 適 切 に ︵ 線 型 結 合 と し て ︶ 近 似 す る こ と で ル ン ゲ = ク ッ タ 法 の 公 式 と な る 。 そ の 上 、 係 数 を テ イ ラ ー 展 開 よ り 正 し く 選 択 す る と 、 方 法 の 収 束 性 も 求 積 法 の 収 束 性 よ り 保 証 さ れ る 。 し か し 、 局 所 誤 差 の オ ー ダ ー や 上 界 は 、 方 法 に よ っ て 大 き く 異 な る の で 、 方 法 別 に 計 算 し な け れ ば な ら な い 。
陽的ルンゲ゠クッタ法では
a
i
j
=
0
{\displaystyle a_{ij}=0}
(
i
≤
j
{\displaystyle i\leq j}
) が満たされるので、ブッチャー配列は以下の形になる。
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}}
こ の 時 、 各 勾 配 k 1 , . . . , k s を 以 下 の よ う に 逐 次 的 に 求 め る こ と が で き る 。
k
1
=
f
(
t
n
,
y
n
)
,
k
2
=
f
(
t
n
+
c
2
h
,
y
n
+
a
21
h
k
1
)
,
k
3
=
f
(
t
n
+
c
3
h
,
y
n
+
a
31
h
k
1
+
a
32
h
k
2
)
,
⋮
k
s
=
f
(
t
n
+
c
s
h
,
y
n
+
a
s
1
h
k
1
+
a
s
2
h
k
2
+
⋯
+
a
s
,
s
−
1
h
k
s
−
1
)
.
{\displaystyle {\begin{aligned}k_{1}&=f(t_{n},y_{n}),\\k_{2}&=f(t_{n}+c_{2}h,y_{n}+a_{21}hk_{1}),\\k_{3}&=f(t_{n}+c_{3}h,y_{n}+a_{31}hk_{1}+a_{32}hk_{2}),\\&\vdots \\k_{s}&=f(t_{n}+c_{s}h,y_{n}+a_{s1}hk_{1}+a_{s2}hk_{2}+\cdots +a_{s,s-1}hk_{s-1}).\end{aligned}}}
( 文 献 に よ っ て 、 等 価 で あ る が 上 と 異 な る 定 義 の 仕 方 を し て い る も の が あ る こ と に 注 意 す る ) [ 5 ] 。
陽 的 ル ン ゲ ゠ ク ッ タ 法 が 目 的 の 精 度 p に な る 係 数 を 持 つ た め に は 、 十 分 に 大 き な 段 数 s が 必 要 に な る 。 少 な く と も 次 数 以 上 の 段 数 が 必 要 (
s
≥
p
{\displaystyle s\geq p}
) で あ り 、 更 に 5 次 以 上 の 場 合 に は 段 数 は 次 数 よ り も 大 き く 取 ら な け れ ば な ら な い (
s
>
p
{\displaystyle s>p}
) こ と が 知 ら れ て い る 。 各 次 数 p に 対 す る 最 低 段 数 s の 一 般 式 は 未 解 決 問 題 で あ る 。 具 体 的 な 値 が 幾 つ か 知 ら れ て い る :
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}}}
オ イ ラ ー 法 、 ホ イ ン 法 、 古 典 的 ル ン ゲ = ク ッ タ 法 ( R K 4 ) の 相 対 誤 差 の 比 較 。 初 期 値
y
(
0
)
=
0
{\displaystyle y(0)=0}
の も と で の 常 微 分 方 程 式
d
y
d
t
=
cos
(
y
)
{\displaystyle {\frac {dy}{dt}}=\cos(y )}
の 数 値 解 の
t
=
1
{\displaystyle t=1}
で の 値 を ス テ ッ プ サ イ ズ の 関 数 と し て 対 数 プ ロ ッ ト し た 。 各 手 法 が そ れ ぞ れ 1 次 精 度 、 2 次 精 度 、 4 次 精 度 で あ る こ と に 対 応 し て 、 傾 き 1 , 2 , 4 で 誤 差 が 減 少 し て い る [ 13 ] 。
最 も 単 純 な ル ン ゲ ゠ ク ッ タ 法 が 1 段 1 次 の 方 法 で あ り 、 一 意 に 定 ま る 。 ( 前 進 ) オ イ ラ ー 法 と 等 価 に な る 。
y
n
+
1
=
y
n
+
h
f
(
t
n
,
y
n
)
.
{\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}).}
以 下 の ブ ッ チ ャ ー 配 列 で 表 現 さ れ る 。
2段2次の方法は1パラメータの自由度 α を持ち、以下のブッチャー配列で与えられる。
0
α
{\displaystyle \alpha }
α
{\displaystyle \alpha }
(
1
−
1
2
α
)
{\displaystyle (1-{\tfrac {1}{2\alpha }})}
1
2
α
{\displaystyle {\tfrac {1}{2\alpha }}}
対応する公式は
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 )}.}
である。
α
=
1
/
2
{\displaystyle \alpha =1/2}
の場合が 中点法 (または修正オイラー法 , modified Euler method)
y
n
+
1
=
y
n
+
h
f
(
t
n
+
h
2
,
y
n
+
h
2
f
(
t
n
,
y
n
)
)
{\displaystyle y_{n+1}=y_{n}+hf\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {h}{2}}f(t_{n},y_{n})\right)}
に対応し、以下のブッチャー配列で与えられる。
α
=
1
{\displaystyle \alpha =1}
の場合は ホイン法 (英語版 ) (または改良オイラー法 , improved Euler method)として知られ、ブッチャー配列は以下の通りである[3] 。
古典的ルンゲ゠クッタ法 (RK4) は既に述べた通り以下のブッチャー配列で与えられる[8] :
0
1/2
1/2
1/2
0
1/2
1
0
0
1
1/6
1/3
1/3
1/6
修正版としてクッタの3/8公式 (英 : Kutta's 3/8-rule ) が知られている[15] 。
0
1/3
1/3
2/3
-1/3
1
1
1
-1
1
1/8
3/8
3/8
1/8
他に使われている方法は、ルンゲ=クッタ法のリスト に参照。
前 述 の 通 り 、 2 段 の 陽 的 方 法 が 2 次 精 度 を 持 つ た め の 係 数 に 対 す る 条 件 は 、
b
1
+
b
2
=
1
{\displaystyle b_{1}+b_{2}=1}
,
b
2
c
2
=
1
/
2
{\displaystyle b_{2}c_{2}=1/2}
, か つ
b
2
a
21
=
1
/
2
{\displaystyle b_{2}a_{21}=1/2}
で あ る [ 7 ] 。 例 と し て 、 そ れ ら の 条 件 の 導 出 を 見 て み よ う 。
一 般 的 に 2 段 陽 的 ル ン ゲ = ク ッ タ 法 に 対 応 す る 配 列 は 以 下 の 通 り で あ る ( 一 貫 性 を 用 い て
c
1
=
0
{\displaystyle c_{1}=0}
が わ か る ) 。
0
c
2
a
21
b
1
b
2
{\displaystyle {\begin{array}{c|cc}0&&\\c_{2}&a_{21}&\\\hline &b_{1}&b_{2}\\\end{array}}}
公 式 か ら
k
2
=
f
(
t
n
+
c
2
h
,
y
n
+
h
a
21
f
(
t
n
,
y
n
)
)
{\displaystyle k_{2}=f(t_{n}+c_{2}h,y_{n}+ha_{21}f(t_{n},y_{n}))}
の f ( t n , y n ) に 関 す る テ イ ラ ー 展 開 を 考 え る
k
2
=
f
(
t
n
,
y
n
)
+
c
2
h
f
t
(
t
n
,
y
n
)
+
h
a
21
f
(
t
n
,
y
n
)
f
y
(
t
n
,
y
n
)
+
O
(
h
2
)
{\displaystyle k_{2}=f(t_{n},y_{n})+c_{2}hf_{t}(t_{n},y_{n})+ha_{21}f(t_{n},y_{n})f_{y}(t_{n},y_{n})+O(h^{2})}
こ こ で 、 f t は f の t に 関 す る 偏 微 分 で あ る 。 そ れ を
y
n
+
1
=
y
n
+
h
(
b
1
k
1
+
b
2
k
2
)
{\displaystyle y_{n+1}=y_{n}+h(b_{1}k_{1}+b_{2}k_{2})}
に 代 入 し 、
k
1
=
f
(
t
n
,
y
n
)
{\displaystyle k_{1}=f(t_{n},y_{n})}
を 用 い て 整 理 す る と
y
n
+
1
=
y
n
+
h
(
b
1
+
b
2
)
f
(
t
n
,
y
n
)
+
h
2
b
2
(
c
2
f
t
(
t
n
,
y
n
)
+
a
21
f
(
t
n
,
y
n
)
f
y
(
t
n
,
y
n
)
)
+
O
(
h
3
)
{\displaystyle y_{n+1}=y_{n}+h(b_{1}+b_{2})f(t_{n},y_{n})+h^{2}b_{2}(c_{2}f_{t}(t_{n},y_{n})+a_{21}f(t_{n},y_{n})f_{y}(t_{n},y_{n}))+O(h^{3})}
と な る 。 同 時 に 、 方 程 式
y
′
(
t
)
=
f
(
t
,
y
)
{\displaystyle y'(t )=f(t,y)}
の t に 関 す る 微 分 を 取 っ て y ' を f ( t , y ) に 置 き 換 え る と
y
″
(
t
)
=
f
t
(
t
,
y
)
+
f
(
t
,
y
)
f
y
(
t
,
y
)
{\displaystyle y''(t )=f_{t}(t,y)+f(t,y)f_{y}(t,y)}
と な る 。 厳 密 解 を y ( t ) と す る 。 上 記 式 を 用 い て y ( t n + 1 ) の t n に 関 す る テ イ ラ ー 展 開 を 考 え る
y
(
t
n
+
1
)
=
y
(
t
n
)
+
h
f
(
t
n
,
y
n
)
+
1
2
h
2
(
f
t
(
t
n
,
y
n
)
+
f
(
t
n
,
y
n
)
f
y
(
t
n
,
y
n
)
)
+
O
(
h
3
)
{\displaystyle y(t_{n+1})=y(t_{n})+hf(t_{n},y_{n})+{\frac {1}{2}}h^{2}(f_{t}(t_{n},y_{n})+f(t_{n},y_{n})f_{y}(t_{n},y_{n}))+O(h^{3})}
2 次 精 度 を 持 つ 条 件 は 、 局 所 誤 差
e
n
+
1
,
h
=
y
(
t
n
+
1
)
−
y
n
+
1
=
O
(
h
3
)
{\displaystyle e_{n+1,h}=y(t_{n+1})-y_{n+1}=O(h^{3})}
で あ る 。 二 つ の 展 開 式 の 係 数 を 比 較 す る と 、 そ の 条 件 が
b
1
+
b
2
=
1
{\displaystyle b_{1}+b_{2}=1}
,
b
2
c
2
=
1
/
2
{\displaystyle b_{2}c_{2}=1/2}
, か つ
b
2
a
21
=
1
/
2
{\displaystyle b_{2}a_{21}=1/2}
で あ る こ と が わ か る 。
い ま ま で に 述 べ た 方 法 は す べ て 陽 公 式 で あ る 。 陽 的 ル ン ゲ = ク ッ タ 方 法 の 絶 対 安 定 性 領 域 ︵ r e g i o n o f a b s o l u t e s t a b i l i t y ︶ が 小 さ く て 有 界 で あ る た め 、 硬 い 方 程 式 の 解 を 計 算 す る 場 合 、 方 法 は 不 適 切 で あ る 。 こ の 問 題 は 特 に 偏 微 分 方 程 式 の 解 を 計 算 す る と き に 重 要 で あ る 。
陽 的 方 法 の 不 安 定 さ は 陰 的 ル ン ゲ = ク ッ タ 法 の 開 発 の 動 機 と な る 。 陰 的 ル ン ゲ = ク ッ タ 法 は 上 述 の 一 般 的 な 公 式 と 同 じ く 、 以 下 の 形 で 与 え ら れ る
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
+
h
c
i
,
y
n
+
h
∑
j
=
1
s
a
i
j
k
j
)
,
i
=
1
,
…
,
s
.
{\displaystyle k_{i}=f\left(t_{n}+hc_{i},y_{n}+h\sum _{j=1}^{s}a_{ij}k_{j}\right),\quad i=1,\ldots ,s.}
で あ り [ 4 ] 、 対 応 す る ル ン ゲ = ク ッ タ 行 列 は 厳 密 な 下 三 角 行 列 で は な い 。 結 果 と し て 、 一 時 刻 ご と に 代 数 方 程 式 系 を 解 か な け れ ば な ら な く な る 。 応 じ て 、 計 算 コ ス ト も か な り 上 が る 。 s 段 法 を 使 っ て m 個 の 成 分 か ら な る 微 分 方 程 式 系 ︵ す な わ ち
y
→
=
(
y
1
,
…
,
y
m
)
{\displaystyle {\vec {y}}=(y_{1},\ldots ,y_{m})}
の と き ︶ に 適 用 す る 場 合 に 対 応 す る 代 数 方 程 式 の 数 は ms と な る 。 こ れ は 陰 的 線 型 多 段 法 と も 比 較 で き る ‥ s 段 陰 的 線 型 多 段 法 を 使 っ て 同 じ 方 程 式 系 に 適 用 す る 場 合 に 対 応 す る 代 数 方 程 式 の 数 は m で あ り 、 段 数 s に 依 存 し な い [ 20 ] 。
一 番 単 純 な 陰 的 ル ン ゲ = ク ッ タ 法 は 、 後 退 オ イ ラ ー 法
y
n
+
1
=
y
n
+
h
f
(
t
n
+
1
,
y
n
+
1
)
{\displaystyle y_{n+1}=y_{n}+hf(t_{n+1},y_{n+1})}
で あ る 。 対 応 す る ブ ッ チ ャ ー 配 列 も 単 純 に 以 下 の 通 り で あ る 。
1
1
1
{\displaystyle {\begin{array}{c|c}1&1\\\hline &1\\\end{array}}}
他 の 例 と し て 、 台 形 公 式 と 呼 ば れ る 方 法 は 以 下 の 公 式 で 与 え ら れ る 。
y
n
+
1
=
y
n
+
1
2
h
(
f
(
t
n
,
y
n
)
+
f
(
t
n
+
1
,
y
n
+
1
)
)
{\displaystyle y_{n+1}=y_{n}+{\frac {1}{2}}h(f(t_{n},y_{n})+f(t_{n+1},y_{n+1}))}
対 応 す る 配 列 は 以 下 の 通 り で あ る 。
0
0
0
1
1
2
1
2
1
2
1
2
{\displaystyle {\begin{array}{c|cc}0&0&0\\1&{\frac {1}{2}}&{\frac {1}{2}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\\end{array}}}
そ の 対 応 は 、
y
n
+
1
=
y
n
+
1
2
k
1
+
1
2
k
2
{\displaystyle y_{n+1}=y_{n}+{\frac {1}{2}}k_{1}+{\frac {1}{2}}k_{2}}
を 覚 え て て 、 公 式 を 以 下 の よ う に 書 き 換 え る と 明 ら か に な る 。
y
n
+
1
=
y
n
+
1
2
f
(
t
n
,
y
n
)
+
1
2
f
(
t
n
+
h
,
y
n
+
1
2
k
1
+
1
2
k
2
)
{\displaystyle y_{n+1}=y_{n}+{\frac {1}{2}}f(t_{n},y_{n})+{\frac {1}{2}}f\left(t_{n}+h,y_{n}+{\frac {1}{2}}k_{1}+{\frac {1}{2}}k_{2}\right)}
台 形 公 式 は 、 選 点 法 で あ る 。 選 点 法 は す べ て 陰 的 ル ン ゲ = ク ッ タ 法 で あ る け ど 、 陰 的 ル ン ゲ = ク ッ タ 法 が す べ て 選 点 法 で あ る わ け で は な い 。
選 点 法 の 中 で は 、 ガ ウ ス 求 積 に 基 づ い た ガ ウ ス ・ ル ジ ャ ン ド ル 法 ︵ 英 語 版 ︶ が 一 番 次 数 の 高 い 方 法 で あ る 。 s 段 ガ ウ ス ・ ル ジ ャ ン ド ル 法 の 次 数 は 2 s と な る ︵ よ っ て 、 任 意 に 高 い 次 数 を 持 つ 方 法 を 構 造 で き る よ う に な る ︶ [ 22 ] 。 例 と し て 、 2 段 ガ ウ ス ・ ル ジ ャ ン ド ル 法 は 次 の 配 列 で 与 え ら れ る 。
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
{\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}}\\\end{array}}}
[22]
数値分析における安定性 は、それぞれ異なる定義が複数存在する。ルンゲ=クッタ法の安定性を反映できる概念は、主に以下の二つである。
線 型 テ ス ト 方 程 式
y
′
=
λ
y
{\displaystyle y'=\lambda y}
を 考 え る 。 一 つ の ル ン ゲ = ク ッ タ 法 を 使 っ て こ の 方 程 式 に 適 用 す る と
y
n
+
1
=
r
(
h
λ
)
y
n
{\displaystyle y_{n+1}=r(h\lambda )y_{n}}
と な る 。 こ こ で 、
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)}}}
は 安 定 性 関 数 と 呼 ば れ る C 上 の 有 理 関 数 で あ る ︵ e は す べ て の 成 分 が 1 の ベ ク ト ル [ 要 曖 昧 さ 回 避 ] で あ る ︶ 。 s 段 法 の 場 合 、 行 列 式 の 展 開 に よ っ て r ( z ) は 二 つ の s 次 多 項 式 の 商 と な る 。 陽 的 方 法 の 場 合 、 対 応 す る ル ン ゲ = ク ッ タ 行 列 が 狭 義 下 三 角 行 列 で あ る た め 、
det
(
I
−
z
A
)
=
1
{\displaystyle \det(I-zA)=1}
が わ か る 。 し た が っ て 、 r ( z ) は s 次 多 項 式 と な る [ 25 ] 。
ル ン ゲ = ク ッ タ 法 に よ る 上 記 の 方 程 式 の 数 値 解 が ゼ ロ に 減 衰 す る 条 件 が
|
r
(
z
)
|
<
1
{\displaystyle |r(z )|<1}
(
z
=
h
λ
{\displaystyle z=h\lambda }
) で あ る 。 上 記 の 条 件 を 満 た す z の 集 合 は 方 法 に 対 す る 絶 対 安 定 性 領 域 ︵ r e g i o n o f a b s o l u t e s t a b i l i t y ︶ で あ る 。 特 に 、 絶 対 安 定 性 領 域 が 左 複 素 数 平 面 を 含 む と き 、 そ の 方 法 は A - 安 定 ( A - s t a b l e ) と い う 。 陽 的 ル ン ゲ = ク ッ タ 法 は 、 安 定 性 関 数 が 多 項 式 で あ る た め 、 A - 安 定 な 方 法 に な れ な い [ 25 ] 。
も っ と 一 般 的 に 、 方 法 の 次 数 が p の と き に 、 安 定 性 関 数 に 対 し
r
(
z
)
=
e
z
+
O
(
z
p
+
1
)
{\displaystyle r(z )=e^{z}+O(z^{p+1})}
が 成 立 す る 。 ゆ え に 、 指 数 関 数 e z の 、 与 え ら れ た 次 数 の 多 項 式 の 商 か ら な る 有 理 関 数 の 中 で の 最 善 近 似 が 重 要 だ と 考 え ら れ る 。 そ の よ う な 有 理 関 数 は パ デ 近 似 と 呼 ば れ る 。 分 子 の 次 数 が m で 分 母 の 次 数 が n の パ デ 近 似 式 が A - 安 定 性 の 条 件
|
r
(
z
)
|
<
1
{\displaystyle |r(z )|<1}
,
R
e
(
z
)
≤
0
{\displaystyle \mathrm {Re} (z )\leq 0}
[ 26 ] を 満 足 す る た め の 必 要 十 分 条 件 は 、
m
≤
n
≤
m
+
2
{\displaystyle m\leq n\leq m+2}
で あ る 。
s 段 ガ ウ ス ・ ル ジ ャ ン ド ル 法 の 次 数 は 前 述 通 り 2 s で あ る 。 よ っ て 安 定 性 関 数 の 分 子 と 分 母 は 同 じ く s 次 多 項 式 と な る 。 す な わ ち 、
m
=
n
=
s
{\displaystyle m=n=s}
で あ る 。 上 記 の 条 件 が 満 た さ れ る の で 、 ガ ウ ス ・ ル ジ ャ ン ド ル 法 は A - 安 定 で あ る こ と が わ か る 。 故 に 任 意 に 高 い 次 数 を 持 つ 、 A - 安 定 な ル ン ゲ = ク ッ タ 法 が 存 在 す る 。 比 べ て 、 A - 安 定 性 を 持 つ 線 型 多 段 法 の 次 数 は 2 以 下 で あ る [ 29 ] 。
以 上 の こ と か ら 、 陰 的 ル ン ゲ = ク ッ タ 法 は 陽 的 な 方 法 よ り も 優 れ た 安 定 性 を 持 つ こ と も わ か る 。
A - 安 定 性 と い う 概 念 は 線 型 自 励 方 程 式
y
′
=
λ
y
{\displaystyle y'=\lambda y}
を 考 え る 結 果 で あ る 。 ダ ー ル キ ス ト ︵ 英 語 版 ︶ は 、 数 値 的 方 法 を 使 っ て 、 と あ る 単 調 性 条 件 を 満 た す 非 線 型 方 程 式 系 に 適 用 す る と き の 安 定 性 を 主 張 し た 。 対 応 す る 安 定 性 は 、 線 型 多 段 法 の 場 合 に G - 安 定 性 ︵ G - s t a b i l i t y ︶ で 、 ル ン ゲ = ク ッ タ 法 の 場 合 に B - 安 定 性 ( B - s t a b i l i t y ) と 呼 ぶ 。 一 般 的 な テ ス ト 方 程 式
y
′
=
f
(
x
,
y
)
{\displaystyle y'=f(x,y)}
と 単 調 性 条 件
⟨
f
(
x
,
y
)
−
f
(
x
,
z
)
,
y
−
z
⟩
≤
0
{\displaystyle \langle f(x,y)-f(x,z),y-z\rangle \leq 0}
を 満 足 す る f を 考 え る 。 も し す べ て の
h
≥
0
{\displaystyle h\geq 0}
に 対 し 、 不 等 式
‖
y
n
+
1
−
z
n
+
1
‖
≤
‖
y
n
−
z
n
‖
{\displaystyle \|y_{n+1}-z_{n+1}\|\leq \|y_{n}-z_{n}\|}
が 成 立 す る と き 、 そ の ル ン ゲ = ク ッ タ 法 は B - 安 定 と い う 。 こ こ で 、 y n と z n は そ れ ぞ れ の 初 期 値 に 対 す る 数 値 解 で あ る 。
f
(
x
,
y
)
=
λ
y
{\displaystyle f(x,y)=\lambda y}
に 方 法 を 適 用 す る こ と で 、 B - 安 定 性 は A - 安 定 性 よ り 強 い 条 件 で あ る こ と が わ か る 。
(一) ^ P r e s s e t a l . 2 0 0 7 , p . 9 0 8 .
(二) ^ S ü l i & M a y e r s 2 0 0 3 , p . 3 2 8 .
(三) ^ a b S ü l i & M a y e r s 2 0 0 3 , p . 3 2 8 .
(四) ^ a b I s e r l e s 2 0 0 8 , p . 4 1 ; S ü l i & M a y e r s 2 0 0 3 , p p . 3 5 1 – 3 5 2 .
(五) ^ a b A t k i n s o n ( 1 9 8 9 , p . 4 2 3 ) , H a i r e r , N ø r s e t t & W a n n e r ( 1 9 9 3 , p . 1 3 4 ) , K a w & K a l u ( 2 0 0 8 , § 8 . 4 ) a n d S t o e r & B u l i r s c h ( 2 0 0 2 , p . 4 7 6 ) l e a v e o u t t h e f a c t o r h i n t h e d e f i n i t i o n o f t h e s t a g e s . A s c h e r & P e t z o l d ( 1 9 9 8 , p . 8 1 ) , B u t c h e r ( 2 0 0 8 , p . 9 3 ) a n d I s e r l e s ( 2 0 0 8 , p . 3 8 ) u s e t h e y v a l u e s a s s t a g e s .
(六) ^ a b I s e r l e s 2 0 0 8 , p . 3 8 .
(七) ^ a b I s e r l e s 2 0 0 8 , p . 3 9 .
(八) ^ a b S ü l i & M a y e r s 2 0 0 3 , p . 3 5 2 .
(九) ^ I s e r l e s 2 0 0 8 , p . 3 4 .
(十) ^ P r e s s e t a l . 2 0 0 7 , p . 9 0 7 .
(11) ^ B u t c h e r 2 0 0 8 , p . 1 8 7 .
(12) ^ B u t c h e r 2 0 0 8 , p p . 1 8 7 – 1 9 6 .
(13) ^ 齊 藤 宣 一 ﹃ 数 値 解 析 ︵ 共 立 講 座 数 学 探 検 17 ︶ ﹄ 共 立 出 版 、 2 0 1 7 年 、 1 0 7 頁 。 I S B N 9 7 8 4 3 2 0 9 9 2 7 4 0 。
(14) ^ S ü l i & M a y e r s 2 0 0 3 , p . 3 2 7 .
(15) ^ H a i r e r , N ø r s e t t & W a n n e r ( 1 9 9 3 , p . 1 3 8 ) r e f e r t o K u t t a ( 1 9 0 1 ) .
(16) ^ I s e r l e s 2 0 0 8 , p . 1 0 7 .
(17) ^ I s e r l e s 2 0 0 8 , p . 1 0 9 .
(18) ^ H a i r e r , N ø r s e t t & W a n n e r 1 9 9 3 , p . 1 7 7 .
(19) ^ S ü l i & M a y e r s 2 0 0 3 , p p . 3 4 9 – 3 5 1 .
(20) ^ S ü l i & M a y e r s 2 0 0 3 , p . 3 5 3 .
(21) ^ I s e r l e s 2 0 0 8 , p p . 4 3 – 4 4 .
(22) ^ a b I s e r l e s 2 0 0 8 , p . 4 7 .
(23) ^ H a i r e r & W a n n e r 1 9 9 6 , p p . 4 0 – 4 1 .
(24) ^ H a i r e r & W a n n e r 1 9 9 6 , p . 4 0 .
(25) ^ a b I s e r l e s 2 0 0 8 , p . 6 0 .
(26) ^ I s e r l e s ( 2 0 0 8 ) は そ の よ う な 有 理 関 数 を A - a c c e p t a b l e と 呼 ぶ 。
(27) ^ I s e r l e s 2 0 0 8 , p p . 6 2 – 6 3 .
(28) ^ I s e r l e s 2 0 0 8 , p . 6 3 .
(29) ^ こ の 結 果 ︵ 時 々 に s e c o n d D a h l q u i s t b a r r i e r ︶ は 、 D a h l q u i s t ( 1 9 6 3 ) に よ る も の で あ る 。
(30) ^ B u t c h e r 1 9 7 5 .
(31) ^ H a i r e r & W a n n e r 1 9 9 6 , p . 1 8 1 .
Atkinson, Kendall A. (1989), An Introduction to Numerical Analysis (2nd ed.), New York: John Wiley & Sons , ISBN 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 , http://journals.cambridge.org/action/displayAbstract?fromPage=online&aid=4922056 .
Butcher, John C. (1975), “A stability property of implicit Runge-Kutta methods”, BIT 15 : 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 .
Dahlquist, Germund (1963), “A special stability problem for linear multistep methods”, BIT 3 : 27–43, doi :10.1007/BF01963532 , ISSN 0006-3835 .
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 .
Hairer, Ernst; Wanner, Gerhard (1996), Solving ordinary differential equations II: Stiff and differential-algebraic problems (2nd ed.), Berlin, New York: Springer-Verlag , ISBN 978-3-540-60452-5 .
Iserles, Arieh (2008), A First Course in the Numerical Analysis of Differential Equations (Second Edition) , Cambridge University Press , ISBN 978-0-521-73490-5 .
Kaw, Autar; Kalu, Egwu (2008), Numerical Methods with Applications (1st ed.), autarkaw.com, http://numericalmethods.eng.usf.edu/topics/textbook_index.html .
Kutta, Martin Wilhelm (1901), “Beitrag zur näherungsweisen Integration totaler Differentialgleichungen”, Zeitschrift für Mathematik und Physik 46 : 435–453 .
Press, William H.; Flannery, Brian P.; Teukolsky, Saul A. ; Vetterling, William T. (2007), “Section 17.1 Runge-Kutta Method” , Numerical Recipes: The Art of Scientific Computing (3rd ed.), Cambridge University Press , ISBN 978-0-521-88068-8 , http://apps.nrbook.com/empanel/index.html#pg=907 . Also, Section 17.2. Adaptive Stepsize Control for Runge-Kutta .
Stoer, Josef; Bulirsch, Roland (2002), Introduction to Numerical Analysis (3rd ed.), Berlin, New York: Springer-Verlag , ISBN 978-0-387-95452-3 .
Süli, Endre; Mayers, David (2003), An Introduction to Numerical Analysis , Cambridge University Press , ISBN 0-521-00794-1 .