10.1 常微分方程式の初期値問題の基本数値解法

10.1.1 数値解法の分類

変数 \(t\) を離散化して, \(t_0, t_1, \dots, t_n, \dots\) と表す. これらの点の間隔をステップ幅あるいはステップ幅 (step size) といい, \(h_n\) で表す.
\[
h_n = t_{n+1} – t_n
\] ステップ幅は必ずしも一定でなくてよいが, \(n\) によらず一定の場合には単に \(h\) と表すことがある. 未知関数 \(y(t)\) の \(t_n\) における値 \(y(t_n)\) の近似値を \(y_n\) と表すことにする.

差分法 (微分を差分で近似する方法) により数値解法を導くことを考える. 代表的な差分の求め方には, 次のようなものがある.

• 前進差分: \((y_{n+1} – y_n)/h\)
• 後退差分: \((y_n – y_{n-1})/h\)
• 中心差分: \((y_{n+1} – y_{n-1})/2h\)

まず, 前進差分を使って微分を近似する.
\[
y_n’ = f(t_n, y_n) \simeq (y_{n+1} – {y_n})/h
\] これより次の近似式が得られる.
\[
y_{n+1} = y_n + hf(t_n, y_n)
\] この式を使い, 初期値 \(y_0\) から始めて, \(y_1, y_2, \sim\) と次々に各点における近似解を求めることができることがわかる. これをオイラー法 (Euler method) という.

後退差分を使った場合には, 同様にして次式により近似解を求めることができる.
\[
y_{n+1} = y_n + hf(t_{n+1}, y_{n+1})
\] これを後退オイラー法 (backward Euler method) という.

後退オイラー法では右辺の \(f()\) の中に \(y_{n+1}\) が含まれており, \(y_{n+1}\) を求めるために方程式を解かなければならない (一般的には非線形方程式になる). このような解法を陰解法 (陰的公式, implicit formula) という. これに対して, オイラー法の場合には方程式を解くことなく直接 \(y_{n+1}\) を計算することができた. そのような解法を陽解法 (陽的公式, explicit formula) という.

次に, 中心差分を使った場合には, 次式により近似解を求めることができる.
\[
y_{n+1} = y_{n-1} + 2hf(t_n, y_n)
\] これは中点則とよばれる.

中点則では, \(y_{n+1}\) の計算には \(y_n\) の他に過去の値 \(y_{n-1}\) も必要で, 2 つの値を使うので 2 段法とよばれ, 多段法 (多段階法, multistep method) とよばれる解法のひとつになる. これに対して, オイラー法と後退オイラー法は \(y_n\) だけしか使わないので, 1 段法 (1 段階法, one-step method) とよばれる解法に分類される.

まとめると, 常微分方程式の初期値問題の解法を大きく分類すると 1 段法と多段法があり, それぞれに陽解法と陰解法がある.

10.1.2 数値解法の次数

関数 \(y(t)\) (真の解) を \(t_n\) においてテイラー展開すると次のようになる.
\[
\begin{align}
y(t_{n+1}) & = y(t_n) + hy'(t_n) + (1/2)h^2y”(t_n) + (1/6)h^3y”'(t_n) + \dots + (1/p)h^py^(p)(t_n) + \dots \\
& = y(t_n) + hf(t_n, y(t_n)) + (1/2)h^2f'(t_n, y(t_n))+ (1/6)h^3f”(t_n, y(t_n)) + \dots + (1/p)h^pf^{p-1}(t_n, y(t_n)) + \dots \\
\end{align}
\] 数値解法による近似解 \(y_{n+1}\) の右辺が上のテイラー展開 \(y(t_{n+1})\) の p 次の項まで一致するとき, すなわち次式が成り立つときにその数値解法は p 次の方法であるという.
\[
y_{n+1} – y(t_{n+1}) = O(h^{p+1})
\] オイラー法はテイラー展開の 1 次 \((h)\) の項までで打ち切ったものに等しいから明らかに 1 次の方法である. また, 後退オイラー法は 1 次, 中点則は 2 次の方法であることを確かめることができる.

数値実験 (1)

最も簡単な微分方程式の初期値問題の例を考える.
\[
y’ = y, \space y(0) = 1
\] この問題は解析解がわかっていて \(y = e^t\) である. 上の 3 つの公式を用いてステップ幅を \(0.1, 0.05, 0.01, \dots\) というように小さくしていきながら順に計算する. \(t = 1\) において得られた近似解の真値 \(y = e\) からの誤差と必要としたステップ数の値を対数プロットすると下図のようになった. 横軸は相対誤差 (右に行くほど精度がよい), 縦軸はステップ数 (= 関数評価 (\(f()\) の計算) 回数) である.

オイラー法と後退オイラー法では, 1,000 ステップ (\(h = 10^{-3}\)) でおよそ 3 桁半の精度, 1,000,000 ステップ (\(h = 10^{-6}\)) でおよそ 6 桁半の精度となり, 1 次の方法であることが確かめられた. 一方, 中点則では, 1,000 ステップ (\(h = 10^{-3}\)) でおよそ 6 桁半の精度, 1,000,000 ステップ (\(h = 10^{-6}\)) でおよそ 12 桁半の精度があり, 2 次の方法となっていることがわかる. 同じ計算量であれば中点則はオイラー法の倍の桁数を正しく求めることができる.

10.1.3 1 段法

10.1.3.1 テイラー法

関数 y(t) の t_n におけるテイラー展開を再度次に示す.
\[
y(t_{n+1}) = y(t_n) + hf(t_n, y(t_n)) + (1/2)h^2f'(t_n, y(t_n)) + (1/6)h^3f”(t_n, y(t_n)) + \dots + (1/p)h^pf^{p-1}(t_n, y(t_n)) + \dots
\] オイラー法はテイラー展開を 1 次の項までで打ち切って近似したものに相当した. そこで, より高次の項まで使った近似式を使えばオイラー法より精度のよい, すなわち, より次数の高い解法を得ることができると考えられる. これをテイラー法という.

上式で \(f'(t, y(t))\) と \(f^{\prime\prime}(t, y(t))\) は具体的には次のように計算する.
\[
\begin{align}
& f'(t, y(t)) = f_t(t, y(t)) + f_y(t, y(t))f(t, y(t)) \\
& f^{\prime\prime}(t, y(t)) = f_{tt}(t, y(t)) + 2f_{ty}(t, y(t))f(t, y(t)) + f_{yy}(t, y(t))f^2(t, y(t)) + f_t(t, y(t))f_y(t, y(t)) + f_y^2(t, y(t))f(t, y(t)) \\
\end{align}
\] 2 次の項を計算するのでさえ結構複雑なことがわかる. そのため, \(f()\) の高階の微分が簡単に計算できる関数でない限りこの方法を実際に適用するのは難しい.

数値実験 (2)

数値解法の次数の実験と同じ最も簡単な微分方程式の初期値問題の例を使って 1 次 (オイラー法), 2 次, 3 次のテイラー法による計算を実際に行って比較する.
\[
y’ = y, \space y(0) = 1
\]

横軸は相対誤差 (右に行くほど精度がよい), 縦軸はステップ数の対数プロットである. それぞれ 1 次 ～ 3 次の傾きを示している. ただし, 1 ステップの計算には上の \(f'(t, y(t))\) や \(f^{\prime\prime}(t, y(t))\) の計算が含まれるので, 1 ステップあたりの計算量は, 2 次のテイラー法では \(f()\) の計算が 3 回相当, 3 次のテイラー法では 6 回相当ということになる. そのため計算量を指標とすると一般的にはオイラー法との差はやや縮まる. ただし, 上の例の方程式では微分係数が 0 か 1 と簡単なため 1 ステップあたりの計算量の差は少ない.

10.1.3.2 2 次のルンゲ・クッタ法

ルンゲ・クッタ法 (Runge-Kutta method) は, \(f()\) の微分を計算することなく高次のテイラー展開に相当する精度を得る 1 段法である. いくつかの \(t\) と \(y\) について \(f(t, y)\) を計算しそれらの重み付き平均を作ってテイラー展開の p 次の項まで一致させることにより p 次の公式を作る.

2次のルンゲ・クッタ法の一般形は次式のとおりである. \(a, b_1, b_2, c\) はパラメータである.
\[
\begin{align}
& k_1 = f(t_n, y_n) \\
& k_2 = f(t_n + ch, y_n + h ak_1) \\
& y_n+1 = y_n + h(b_1k_1 + b_2k_2) \\
\end{align}
\] \(k_2\) の \((t_n, y_n)\) におけるテイラー展開は次のようになる.
\[
k_2 = f_n + c h f_{t,n} + a k_1 f_{y,n} + \dots
\] これを上の式に代入して, 次式が得られる.
\[
y_{n+1} = y_n + (b_1 + b_2) h f_n + b_2 h^2(c f_{t,n} + a f_{y,n}f_n + \dots )
\] この式とテイラー展開の式を比較して, \(h\) と \(h^2\) の係数が等しいとすると次式が得られる.
\[
\begin{align}
& b_1 + b_2 = 1 \\
& b_2c = 1/2 \\
& b_2a = 1/2 \\
\end{align}
\] 未知数が 4 個で式が 3 本なので自由度が 1 ある. よく使われるものとしては, \(c = 1\) としたものをホイン法 (Heun method) とよぶ. また, \(c = 1/2\) としたものを修正オイラー法とよぶことがある (2 次のルンゲ・クッタ法とはよばれないが, これらの総称ととらえてよい).

ホイン法の計算式は次のとおりである.
\[
\begin{align}
& k_1 = f(t_n, y_n) \\
& k_2 = f(t_n + h, y_n + h k_1) \\
& y_{n+1} = y_n + h(k_1 + k_2)/2 \\
\end{align}
\] 修正オイラー法の計算式は次のとおりである.
\[
\begin{align}
& k_1 = f(t_n, y_n) \\
& k_2 = f(t_n + h/2, y_n + h k_1/2) \\
& y_{n+1} = y_n + h k_2 \\
\end{align}
\]

数値実験 (3)

再び最も簡単な微分方程式の初期値問題の例でホイン法と修正オイラー法の精度を確認する.
\[
y’ = y, \space y(0) = 1
\]

横軸は相対誤差 (右に行くほど精度がよい), 縦軸はステップ数の対数プロットである. ホイン法と修正オイラー法の計算結果はほぼ同じである. 2 次の傾きを確認でき, 最初の方で示した中点則と似たグラフになった. ただし, 計算式からわかるようにホイン法と修正オイラー法では 1 ステップあたり 2 回の関数評価 (\(f()\) の計算) を行っているのに対して, 中点則では 1 ステップあたり 1 回の関数評価しか行わなくてよい.

10.1.3.3 ルンゲ・クッタ法の一般形

上と同様に 3 次, 4 次 ・・・ と高次のルンゲ・クッタ法の公式を作ることができ, 多くの公式が提案されている.

ルンゲ・クッタ法は一般的に次のように表される.
\[
\begin{align}
& y_{n+1} = y_n + h \sum_{i=1}^s b_ik_i \\
& k_i = f(t_n + c_ih, y_n + h \sum_{j=1}^s a_{ij}k_j) \space (i = 1, 2, \dots, s) \\
\end{align}
\] s は段数(stage)とよばれ, 1 ステップ進むのに必要な \(f(t, y)\) の計算回数を表す (1 段法や多段法の”段”とは異なる). \(a_{ij}\), \(b_i\) および \(c_i\) はパラメータで, \(b_i\) および \(c_i\) は次式を満たす.
\[
\begin{align}
& c_i = \sum_{j=1}^s a_{ij} \space (i = 1, 2, \dots, s) \\
& \sum_{i=1}^s b_i = 1 \\
\end{align}
\]

ルンゲ・クッタ法のパラメータを表すにはブッチャー配列と呼ばれる次のような表記法がある.
\[
\begin{array}{c|cccc}
c_1 & a_11 & a_12 & \dots & a_1s \\
c_2 & a_21 & a_22 & \dots & a_2s \\
c_3 & a_31 & a_32 & \dots & a_3s \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
c_s & a_s1 & a_s2 & \dots & a_ss \\
\hline
& b_1 & b_2 & \dots & b_s \\
\end{array}
\]

\(a_{ij}\) の値によって 3 つの場合に分けられる.

(i) \(a_{ij} = 0 \space (i \le j)\) の場合 (右上三角部分と対角部分が 0):
陽解法. 単に「ルンゲ・クッタ法」あるいは「陽的 (explicit) ルンゲ・クッタ法」という. \(c_1 = 0\) である. ブッチャー配列では \(a_{ij}\) の右上三角部分と対角部分を書かない (空白にする).

(ii) \(a_{ij} = 0 \space (i < j)\) の場合 (右上三角部分が 0):
半陰解法. 「半陰的 (semi-implicit) ルンゲ・クッタ法」という.

(iii) \(a_{ij} \ne 0 \space (i < j)\) の場合 (右上三角部分に 0 でないものがある):
陰解法. 「陰的 (implicit) ルンゲ・クッタ法」という.

以下 (i) について説明する. (ii), (iii) については別途「10.3 スティフな方程式」で触れる.

上で示したホイン法と修正オイラー法をブッチャー配列で表すと次のようになる.

ホイン法

\[
\begin{array}{c|cc}
0 \\
1 & 1 \\
\hline
& 1/2 & 1/2 \\
\end{array}
\]

修正オイラー法

\[
\begin{array}{c|cc}
0 \\
1/2 & 1/2 \\
\hline
& 0 & 1 \\
\end{array}
\] (陽的) ルンゲ・クッタ法ではできるだけ高次の公式になるようにパラメータを決めるが, s 段ならば s 次になるわけではなく, s 段の公式が達成しうる最大の次数 p は次のようになることがわかっている.

段数 s	1	2	3	4	5	6	7	8	9	10	11
最大次数 p	1	2	3	4	4	5	6	6	7	7	8

10.1.3.4 4 次のルンゲ・クッタ法

4 段 4 次の公式を作ろうとすると, パラメータが 12 個で自由度が 4 なので様々な公式を作ることが可能である. その中で最良のものを選ぶための多くの研究が行われた. 中でもクッタが選んだ以下の 2 組の係数が有名である.

(1) 古典的ルンゲ・クッタ法

単に「ルンゲ・クッタ法」というとこの公式を指す.
\[
\begin{array}{c|cccc}
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}
\] 計算式は次のとおりである.
\[
\begin{align}
& 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) \\
& y_{n+1} = y_n + h(k_1 + 2 k_2 + 2 k_3 + k_4)/6 \\
\end{align}
\]

(2) クッタの 3/8 公式

\[
\begin{array}{c|cccc}
0 \\
1/3 & 1/3 \\
2/3 & -1/3 & 1/2 \\
1 & 1 & -1 & 1 \\
\hline
& 1/8 & 3/8 & 3/8 & 1/8 \\
\end{array}
\] 計算式は次のとおりである.
\[
\begin{align}
& k_1 = f(t_n, y_n) \\
& k_2 = f(t_n + h/2, y_n + hk_1/3) \\
& k_3 = f(t_n + h/2, y_n + h(-k_1/3 + k_2)) \\
& k_4 = f(, t_n + h, y_n + h(k_1 – k_2 + k_3)) \\
& y_{n+1} = y_n + h(k_1 + 3k_2 + 3k_3 + k_4)/8 \\
\end{align}
\] ルンゲ・クッタ法一族において最もよく使われているのは 4 次の古典的ルンゲ・クッタ法である. いくつかある 4 段 4 次の公式間での精度の差はほとんどないが, この公式は係数 \(a_{ij}\) に 0 が多く, すべての係数が正であるなどの理由により他のものより使いやすい. また, 5 次の公式は 6 段以上でないと作れないため, これ以上高次のものを求めると複雑になるためこれが普及したのかもしれない.

数値実験 (4)

最も簡単な微分方程式の初期値問題の例で 4 次の古典的ルンゲ・クッタ法および3/8公式の精度を確認する.
\[
y’ = y, \space y(0) = 1
\]

横軸は相対誤差 (右に行くほど精度がよい), 縦軸はステップ数の対数プロットである. 古典的ルンゲ・クッタ法および 3/8 公式の計算結果はほぼ同じである. 両方とも 4 次の傾きをしていることが確認できる. ただし, 1 ステップあたり 4 回の関数評価 (\(f()\) の計算) が必要である.

10.1.4 多段法

1 段法は \(y_{n+1}\) の計算に現在の値 \(y_n\) のみを使う解法であったが, 過去の値 (それまでに計算済の値) \(y_{n-1}, y_{n-2}, \dots\) も使った方がよい結果が得られそうに思われる. このアイデアを実現したのが多段法である.

1 階の常微分方程式の積分形を考える.
\[
y(t_{n+1}) = y(t_n) + \int_{t_n}^{t_{n+1}} f(t, y(t)) dt
\] ここで, 被積分関数 \(f()\) を補間多項式で近似して右辺の積分を求める方法を一般にアダムス法 (Adams method) という.

10.1.4.1 アダムス・バシュフォース法

\(y(t_{n+1})\) と \(y(t_n)\) は近似値 \(y_{n+1}\) と \(y_n\) を使い, 右辺の積分は \((t_i, f_i) \space (i = n – k +1, \dots, n)\) を通る補間多項式 \(p(t)\) で \(f(t, y(t))\) を置き換えると, 次のように k 段の陽解法の公式が得られ, アダムス・バシュフォース法 (Adams-Bashforth method) とよばれる. ただし, \(f_i = f(t_i, y_i)\) を表し, \(y_i \space (i = n – k + 1, \dots, n)\) が計算済なので \(f_i\) はすぐに求めることができる. この公式は, k 次の公式になる.
\[
\begin{align}
y_{n+1} & = y_n + \int_{t_n}^{t_{n+1}} p(t) dt \\
& = y_n + h \sum_{j=0}^{k-1} \gamma_j \nabla^j f_n \\
\end{align}
\] ただし,
\[
\begin{align}
& \nabla^0 f_n = f_n, \nabla^{j+1} f_n = \nabla^j f_n – \nabla^j f_{n-1} \space (後退差分) \\
& \gamma_j = (-1)^j \int_0^1 C(-s, j) ds \\
\end{align}
\] である. ここで, \(C(n, k)\) は二項係数を表す. なお, \(\gamma_j\) は直接計算せずに次の漸化式を用いて求めることができる.
\[
\begin{align}
& \gamma_0 = 1 \\
& \gamma_j + (1/2)\gamma_{j-1} + (1/3)\gamma_{j-2} + \dots + (1/(j + 1))\gamma_0 = 1 \\
\end{align}
\] \(j = 0 \sim 7\) のときの値は次のとおりである.

\(j\)	0	1	2	3	4	5	6	7
\(\gamma_j\)	1	1/2	5/12	3/8	251/720	95/288	19087/60480	5257/17280

この方法を図に示すと次のようになる. 区間 \([t_n, t_{n+1}]\) の積分を計算するために, 区間 \([t_{n-k+1}, t_n]\) の補間データを使用した補間関数 \(p(t)\) を使う. 図で, 黒丸は補間点, 色つきの部分は積分範囲を表す.

\(k = 1 \sim 5\) について整理すると次式のようになる.
\[
\begin{align}
& k = 1 : y_{n+1} = y_n + hf_n \\
& k = 2 : y_{n+1} = y_n + h/2(3f_n – f_{n-1}) \\
& k = 3 : y_{n+1} = y_n + h/12(23f_n – 16f_{n-1} + 5f_{n-2}) \\
& k = 4 : y_{n+1} = y_n + h/24(55f_n – 59f_{n-1} + 37f_{n-2} – 9f_{n-3}) \\
& k = 5 : y_{n+1} = y_n + h/720(1901f_n – 2774f_{n-1} + 2616f_{n-2} – 1274f_{n-3} + 251f_{n-4}) \\
\end{align}
\] k = 1 の場合は 1 段法とみることができオイラー法に一致する.

10.1.4.2 アダムス・モルトン法

アダムス・バシュフォース法では補間関数を使って外挿により積分を求めたが, 精度を上げるために \(t_{n+1}\) を加えて区間 \([t_{n-k+1}, t_{n+1}]\) の補間データを使う, すなわち, \((t_i, f_i) \space (i = n – k + 1, \dots, n, n + 1)\) を通る補間多項式 \(p^*(t)\) を使い, 外挿することを避けたのがアダムス・モルトン法 (Adams-Moulton method) である. こちらの方が精度がよさそうであるが, その代わり右辺に \(y_{n+1}\) が現れ, 陰解法になる. 式は次のとおりである. この公式は, k + 1 次の公式になる.
\[
\begin{align}
y_{n+1} & = y_n + \int_{t_n}^{t_{n+1}} p^*(x) dx \\
& = y_n + h \sum_{j=0}^k \gamma^*_j \nabla^j f_{n+1}
\end{align}
\] ただし,
\[
\begin{align}
& \nabla^0 f_n = f_n, \space \nabla^{j+1} f_n = \nabla^j f_n – \nabla^j f_{n-1} \space (後退差分) \\
& \gamma^*_j = (-1)^j \int_0^1 C(-s+1, j) ds \\
\end{align}
\] である. ここで, \(C(n, k)\) は二項係数である. なお, \(\gamma_j^{*}\) は直接計算せずに次の漸化式を用いて求めることができる.
\[
\begin{align}
& \gamma^*_0 = 1 \\
& \gamma^*_j + (1/2)\gamma^*_{j-1} + (1/3)\gamma^*_{j-2} + \dots + (1/(j + 1))\gamma^*_0 = 0 \\
\end{align}
\] \(j = 0 \sim 7\) のときの値は次のとおりである.

\(j\)	0	1	2	3	4	5	6	7
\(\gamma^*_j\)	1	-1/2	-1/12	-1/24	-19/720	-3/160	-863/60480	-275/24192

この方法を図に示すと次のようになる. 下図で, 黒丸は補間点, 色つきの部分は積分範囲を表す. 積分範囲は補間区間に含まれており, 外挿されていない.

\(k = 0 \sim 5\) について整理すると次式のようになる.
\[
\begin{align}
& k = 0 : y_{n+1} = y_n + hf_{n+1} \\
& k = 1 : y_{n+1} = y_n + h/2(f_{n+1} + f_n) \\
& k = 2 : y_{n+1} = y_n + h/12(5f_{n+1} + 8f_n – f_{n-1}) \\
& k = 3 : y_{n+1} = y_n + h/24(9f_{n+1} + 19f_n – 5f_{n-1} + f_{n-2}) \\
& k = 4 : y_{n+1} = y_n + h/720(251f_{n+1} + 646f_n – 264f_{n-1} + 106f_{n-2} – 19f_{n-3}) \\
& k = 5 : y_{n+1} = y_n + h/1440 (475f_{n+1} + 1427f_n – 798f_{n-1} + 482f_{n-2} – 173f_{n-3} + 27f_{n-4}) \\
\end{align}
\] k = 0 の場合は 1 段法とみることができ後退オイラー法に一致する. k = 1 の場合は台形則とよばれる公式に一致する.

数値実験 (5)

最も簡単な微分方程式の初期値問題の例で k = 4 のアダムス・バシュフォース法 (陽解法) とアダムス・モルトン法 (陰解法) 単体(注)で計算を行い精度を確認する. (注: 実際には後で説明するように組み合わせて使用し, 単体では用いない)
\[
y’ = y, \space y(0) = 1
\]

横軸は相対誤差 (右に行くほど精度がよい), 縦軸はステップ数の対数プロットである. 4 次のルンゲ・クッタ法の計算結果も示した. 多段法では陰解法の方が精度がよいことが確認できる. 多段法のほうが同じ次数では少し精度が悪いように感じられるが, ステップ数で表示している縦軸を関数評価回数 (計算量) で表示し直すとほとんど差がなくなる.

多段法では, はじめに出発値 \(y_1, y_2, \dots, y_{k-1}\) の値を求めておく必要があるが, そのために 1 段解法などを補助的に使わなければならない. この精度が悪いと以降の計算値に影響するため, 使用する多段法の公式の精度に十分見合うだけの精度を持った方法を使用する必要がある. ここでは最初の 4 ステップ目まではルンゲ・クッタ法 (4 次) により計算を行った.

10.1.4.3 予測子・修正子法

多段法においては, 陽解法は計算が容易であるが満足のいく精度が得られない. 一方, 陰解法は精度はよいが一般的には非線形方程式を毎回解かなければならないため計算が難しい.

そこで広く使われている方法が予測子・修正子法 (predictor-corrector method) である. まず予測子 (predictor) とよばれる陽解法により \(y_{n+1}\) の第 1 近似を求める. 次にその値を使って修正子 (corrector) とよばれる陰解法を 1 回～数回適用しこれを改良する方法である.

予測子と修正子の組み合わせにはバリエーションがあるが, 予測子にアダムス・バシュフォース法, 修正子にアダムス・モルトン法を使った予測子・修正子法はアダムス・バシュフォース・モルトン法 (あるいは単にアダムス法) と呼ばれ, 広く使われている.

予測子・修正子法の計算手順は, (1) 予測子を適用して \(y_{n+1}\) を求める (P), (2) それを使って関数値 \(f_{n+1}\) を計算する (E), (3) 修正子を適用して \(y_{n+1}\) を更新する (C), (4) 更新した \(y_{n+1}\) を使って \(f_{n+1}\) を更新する (E) となる. これを PECE と表すことにすると, 最後の E を省略して PEC とする方法, 修正子の適用を 2 回繰り返して PECECE とする方法など, バリエーションがある.

数値実験 (6)

最も簡単な微分方程式の初期値問題の例で k = 4 の予測子・修正子法 (アダムス・バシュフォース・モルトン法) の精度を確認する.
\[
y’ = y, \space y(0) = 1
\] ここでは PEC と PECE の 2 種類について試してみた. 出発値を求めるためには 4 次のルンゲ・クッタ法を使用した.

横軸は相対誤差 (右に行くほど精度がよい), 縦軸は今回は関数評価回数 (ステップ数ではない) の対数プロットとした. ステップあたりの関数評価回数が PEC と PECE で異なるため, 総合的に評価するように縦軸のプロットを変えた.

その PEC と PECE であるが, この例は簡単な関数であるためかほとんど計算結果に差が出なかった. そのため, 2 回目の E の計算のぶんだけ PECE の関数評価回数が増えただけという結果になった. なお, 別途難しい関数で試してみると PECE のほうがよくなることがあった.

10.1.5 使われなくなった解法

よく知られているが現在では使われなくなった解法についてここで触れておく.

10.1.5.1 ルンゲ・クッタ・ジル法

ルンゲ・クッタ・ジル法 (Runge-Kutta-Gill method) は係数が次のように表される 4 次のルンゲ・クッタ法である.
\[
\begin{array}{c|cccc}
0 \\
1/2 & 1/2 \\
1/2 & (√2-1)/2 & (2-√2)/2 \\
1 & 0 & (-√2)/2 & 1+(√2)/2 \\
\hline
& 1/6 & (2-√2)/6 & (2+√2)/6 & 1/6 \\
\end{array}
\] 通常の 4 次のルンゲ・クッタ法と性能 (収束速度) はほぼ同じであるが Gill の提案には次の特長がある.

(i) 通常の 4 次のルンゲ・クッタ法では計算のために 4 個の変数 (連立方程式の場合には配列) が必要なのに対し, この方法では係数を工夫して 3 個の変数 (配列) だけで計算できるようになっている. [メモリの節約]

(ii) 丸め誤差を減らすための手順 (下記参照) が含まれている (これは他のルンゲ・クッタ法にも適用可能である). [精度の向上]

メモリが少なく, おそらく単精度計算が主流であった, 初期の計算機におけるプログラムに適した方法だったと思われる. かつて広く使われていたが, メモリが豊富で倍精度計算が普通の現代では必要性が薄れてきた. さらに, より高次の公式の実用化やステップ幅の自動調節機能の普及もあり出番が少なくなったと思われる.

ここで, (ii)について少し説明する.

常微分方程式の基本的な数値解法では各ステップにおいて関数値の更新差分 (ここでは k と表す) を計算し, 次のようにして関数値を更新するものとする.
\[
y_{n+1} = y_n + k
\] この計算を計算機の浮動小数データ (ワード長は有限) で行うと下図のようになる.

ここで, k をワード長いっぱい正しく求めたとしても, この更新に使用されるのは色のついた一部だけで残りの部分は捨てられてしまう (丸め誤差という). そこで, この捨てられる部分を別の場所に累積しておき, 色のついた部分に頭が出るほどに貯まるたびに k に加えて関数値の更新に使用することにより丸め誤差を減らす方法が考案された. この方法がルンゲ・クッタ・ジル法に使われている.

数値実験 (7)

最も簡単な微分方程式の初期値問題の例でルンゲ・クッタ・ジル法の精度を確認する.
\[
y’ = y, \space y(0) = 1
\]

横軸は相対誤差 (右に行くほど精度がよい), 縦軸は関数評価回数を対数プロットとした. 比較のため普通のルンゲ・クッタ法もプロットしている. 計算は 64 ビット倍精度演算で行われている.

相対誤差が \(10^{-14}\) 付近までは差がないが, 普通のルンゲ・クッタ法では \(10^{-15}\) あたりで頭打ちになり, その先は刻み幅を小さくするとむしろ誤差が大きくなることもある. これは, 64 ビット倍精度では 1 回の演算の精度の限界は \(2^{-16}\) なので, この場合は総合的には \(10^{-15}\) くらいが最高精度だと考えられ, それ以上計算を繰り返すと丸め誤差が積み重なっていきかえって精度が落ちてしまうためと考えられる.

一方で, ルンゲ・クッタ・ジル法では最後のほうの誤差は \(1.6 \times 10^{-16}\) または 0 となった (図で線がとんでいるところは完全に 0 のためプロットされていない). つまり, 最終 1 ビットが違うかどうかのめいっぱいの精度で計算されたことを示している.

10.1.5.2 ミルン法

アダムス法と似ているが, 積分範囲が異なる次の積分形を考える.
\[
y(t_{n+1}) = y(t_{n-1}) + \int_{t_{n-1}}^{t_{n+1}} f(t, y(t)) dt
\] アダムス法と同様に \((t_i, f_i) \space (i = n – k + 1, \dots, n)\) を通る補間多項式 \(p(t)\) で積分を近似すると, 次の陽解法の公式が得られる.
\[
y_{n+1} = y_{n-1} + h \sum_{j=0}^{k-1} \kappa_j \nabla^j f_n
\] ただし,
\[
\kappa_j = (-1)^j \int_{-1}^1 C(-s, j) ds
\] また, \((t_i, f_i) \space (i = n – k +1, \dots, n+1)\) を通る補間多項式 \(p^*(x)\) で積分を近似すると, 次の陰解法の公式が得られる.
\[
y_{n+1} = y_{n-1} + h \sum_{j=0}^k \kappa^*_j \nabla^j f_{n+1}
\] ただし,
\[
\kappa^*_j = (-1)^j \int_{-1}^1 C(-s+1, j) ds
\] 予測子・修正子法の一種にミルン法 (Milne method) とよばれる方法があり, 上の陰解法の式において k = 2 とした次式がその修正子として使われている.
\[
y_{n+1} = y_{n-1} + h/3 (f_{n+1} + 4f_n + f_{n-1})
\] ミルン法の予測子には, 積分範囲を \([t_{n-3}, t_{n+1}]\) に広げた積分形より導かれる次式が使われている.
\[
y_{n+1} = y_{n-3} + (4h/3) (2f_n – f_{n-1} + 2f_{n-2})
\] ミルン法では予測子で求めた近似に対して補正子 (modifier) を適用して加速することが提案されている. すなわち, 予測子で求めた近似を \(y_{n+1}^p\) で表すと, 次のような補正を行う.
\[
y_{n+1} = y_{n+1}^p + (28/29) (y_n – y_n^p)
\] これを M と表すとミルン法の予測子・修正子法は PMECE… となる.

ミルン法は 4 次の公式で, コンピュータが普及する以前の手計算の時代によく使われていた. しかし, コンピュータで使われるようになると安定性の面で問題があることがわかり, 現在では使われなくなった.

なお, 安定性の問題を解決するために, ミルン法の予測子と組み合わせて使う安定な修正子として次のハミング (Hamming) の修正子が提案されている.
\[
y_{n+1} = 9/8 y_n – 1/8 y_{n-2} + h/8 (3f_{n+1} + 6f_n – 3f_{n-1})
\] なお, 補正子は次のようになる.
\[
y_{n+1} = y_{n+1}^p + (112/121) (y_n – y_n^p)
\]

数値実験 (8)

最も簡単な微分方程式の初期値問題の例でミルン法の精度を確認する.
\[
y’ = y, \space y(0) = 1
\]

横軸は相対誤差 (右に行くほど精度がよい), 縦軸は関数評価回数を対数プロットとした. 比較のためルンゲ・クッタ法とアダムス・バシュフォース・モルトン法もプロットしている. なお, 出発値を求めるためには 4 次のルンゲ・クッタ法を使用した.

ミルン法は公式が比較的簡単な割に, 図のように性能のよい解法であることがわかる.

数値実験 (9)

次の初期値問題をミルン法により解く.
\[
y’ = -y + 1, \space y(0) = 2
\] この方程式の解は \(y(t) = e^{-t} + 1\) である. ステップ幅を 0.01 および 0.001 として計算した結果, 次のようになった (見やすさのためデータ点を間引きして表示した. 実際はもっと細かく振動する).

ステップ幅が十分小さいにもかかわらず計算が進んでステップ数が多くなると振動を始め発散した. これは, 解法の差分方程式の一般解が真の解の近似 (有効成分) の他に必要のない成分を含み, 条件によってそれが有効成分よりも大きくなったときに起きる現象である.

この不安定性の対策としてハミングの修正子があるので, それを使って同じ計算を行うと次のようになった.

ハミングの修正子を使うことにより不安定性が収まることが確認された.

数値実験 (1) では中点則の性能がよかったが, 中点則にもこの種の不安定性があることがわかっている. なお, 1 段法ではこの種の不安定性は生じない.