10.2 ステップ幅の自動調節, 密出力, 2階常微分方程式

10.2.1 誤差の推定

求めた数値解の誤差を推定する方法について説明する. 計算結果が得られればそれで終わりというのではなく, その結果の誤差が推定できていれば, 特異点などのために異常な解が得られていないかなどのチェックができ計算結果の利用価値は上がる. また, 現在主流となっているステップ幅の自動調節機能付きルーチンの実装には誤差推定が必須である.

10.2.1.1 1 段法の誤差

初期値問題の数値解の誤差の原因には次の2つがある.

• 離散化誤差 (打ち切り誤差ともいう)
  演算がすべて無限精度で行われたとしたときに, 解法に起因して生じる誤差
  • 局所的離散化誤差
    それまでに得られている値が正確であるとして, 次の 1 ステップで生じる誤差
  • 大域的離散化誤差
    初期値から計算される解と真の解との差
• 丸め誤差
  数値を浮動小数で表し有限桁数で演算することにより生じる誤差

いま, 次の初期値問題を考える.
\[
dy/dt = f(t, y), \space t = t_0 において y = y_0
\] \(t_n\) における計算値を \(y_n\) とし, \(u_n(t_n) = y_n\) とする. すなわち, 計算値 \(y_n\) から決定される解曲線を \(u_n(t\) )とする. \(y_{n+1}\) は \(y_n, t_n, f\) より計算される推定値である. 例えば, オイラー法を使うならば, ステップ幅を \(h\) として, \(y_{n+1} = y_n + hf(t_n, y_n)\) となる.

そうすると, 局所的離散化誤差 \(d_n\) は次のように表すことができる.
\[
d_n = y_{n+1} – u_n(t_{n+1})
\] 局所的離散化誤差は \(p\) 次の解法において \(O(h^{p+1})\) である. すなわち, 次式を満たす数 \(C\) が存在する.
\[
|d_n| ≦ Ch^{p+1}
\] 大域的離散化誤差 \(e_n\) は次のように表される.
\[
e_n = y_n – y(t_n)
\] 以上の関係を下に図示する. これは, 解曲線が少しずれたときに時間が進むにつれて離れていく関数の例である. すなわち, \(e_n > \sum_{i=0}^{n-1} d_i\) となり, ステップが進むにつれて誤差が拡大していく.

下に他の例を示す. これは, 解曲線が少しずれたときに時間が進むにつれて接近していく関数の例である. すなわち, \(e_n < \sum_{i=0}^{n-1} d_i\) となり, ステップが進むにつれて誤差は縮小していく.

2 つの例をあげたが, 関数がどちらの傾向を示すかは \(\partial f/\partial y\) の符号により決まり, 正であれば前者, 負であれば後者になる. もし, \(f\) が \(y\) に依存しない場合, すなわち, \(\partial f/\partial y = 0\) ならば \(e_n = \sum_{i=0}^{n-1} d_i\) となる. また, \(\partial f/\partial y\) の符号が区間内で変化するような複雑な関数もある. \(\partial f/\partial y = 0\) に近ければ, 大域的離散化誤差はおおむね \(O(h^p)\) といえる.

ここまでは丸め誤差を無視して考えてきたが, 有限桁で計算する限り必ず丸め誤差 \(\epsilon_n\) が発生する. 例えば, オイラー法ならば, \(y_{n+1} = y_n + hf(y_n, t_n) + \epsilon_n\) となる.

ここで, 丸め誤差 \(\epsilon_n\) のノルムはある数 \(\epsilon\) 以下であるとする.
\[
\|\epsilon_n\| \le \epsilon
\] いま, 最終点を \(t_N\) とすると, ステップ数 \(N = (t_N – t_0)/h\) である. 最終的な大域的離散化誤差を \(e_N \simeq \sum_{i=0}^{N-1} d_i\) と近似し, 丸め誤差を \(N\epsilon\) と近似すると, 最終的な全誤差は次のように近似することができる.
\[
全誤差 \simeq NCh^{p+1} + N\epsilon = (t_N – t_0)(Ch^p + \epsilon/h)
\]

数値実験 (10)

次の初期値問題をオイラー法により解く.
\[
y’ = y, \space y(0) = 1
\] オイラー法の場合 \(p = 1\) であるから, 全誤差は次のように近似できる.
\[
全誤差 \simeq (t_N – t_0)(Ch + \epsilon/h)
\] \(h\) が大きいときは第 1 項 (離散化誤差) が優勢で, \(h\) が小さいときは第 2 項 (丸め誤差) が優勢になると考えられる. すなわち, ステップ幅を小さくしていくと誤差は小さくなっていくが, ある値 \(((\epsilon/C)^{-2})\) から下では丸め誤差のためにかえって誤差が大きくなると予想される.

ステップ幅を変化させて, \(t = 1\) における相対誤差をプロットした結果を下に示す. 横軸はステップ幅 (右にいくほど小さい), 縦軸は相対誤差である. なお, 丸め誤差の影響を見やすくするために計算は 32 ビット単精度演算 (10 進の精度はおおよそ \(10^{-7}\)) で行った.

図中の青色のプロットのように, この場合 \(10^{-5}\) 付近を境に丸め誤差が支配的になることがわかる.

ルンゲ・クッタ・ジル法の項目で丸め誤差を減らす方法を紹介したが, それをオイラー法に適用して計算した結果をオレンジ色で表示した. 丸め誤差がなければ \(h\) に比例して誤差が減っていくことがわかる.

10.2.1.2 リチャードソンの補外

与えられた点 \(t_n\) における p 次の数値解法の大域的離散化誤差は \(O(h^p)\) となるので, 次のように表すことができる.
\[
y_n – y(t_n) = ch^p + O(h^q), \space q > p
\] ステップ幅 \(h\) で計算した \(y_n\) を \(y_n^{(h)}\) と書くことにして, ステップ幅 \(2h\) でも計算すると次のようになる.
\[
\begin{align}
& y_n^{(h)} = y(t_n) + ch^p + O(h^q) \\
& y_n^{(2h)} = y(t_n) + c(2h)^p + O(h^q) \\
\end{align}
\] これらの差をとると,
\[
y_n^{(2h)} – y_n^{(h)} = (2^p – 1)ch^p + O(h^q)
\] となり, \(y_n^{(h)}\) の誤差は \(O(h^q)\) を無視すると次のように近似できる.
\[
y_n^{(h)} の誤差 \simeq ch^p = \frac{y_n^{(2h)} – y_n^{(h)}}{2^p – 1}
\] これは, ロンバーグ積分で使われるリチャードソンの補外 (Richardson extrapolation) と同じである.

数値実験 (11)

簡単な微分方程式の初期値問題の例をオイラー法, ホイン法, ルンゲ・クッタ法で計算する.
\[
y’ = y, \space y(0) = 1
\] ステップ幅 \(h = 1\) から始めて, \(h = 1/2, 1/4, 1/8 \dots, 1/1024\) というようにステップ幅を半分にして計算を繰り返した. このとき, リチャードソンの補外によりそれぞれの計算結果の誤差を推定して, 実際の誤差と共にプロットする.

この例では, リチャードソンの補外による誤差の推定値は全体として少なめに見積もられているが, 比較的近い値が得られている.

10.2.2 補外法

10.2.2.1 補外法のアルゴリズム

\(H = t – t_0\) とする. 正の整数列 \(n_1 < n_2 < n_3 < \dots\) について \(h_i = H/n_i\) とすると, 整数列に対応するステップ幅 \(h_1 > h_2 > h_3 > \dots\) が得られる.

\(p\) 次の数値解法においてステップ幅 \(h_j\) を使って \(n_j\) ステップ計算して求めた \(t\) における解を \(T_{j1}\) とする. \(T_{11} \sim T_{k1} \space (k \le j)\) の値を使って連立方程式を解けば \(h\) の \(k\) 次補間多項式を作ることができる. そこで \(h = 0\) として補外値を求めればよい近似を得ることができる.

実際の計算においては, 補間多項式そのものは不要で \(h = 0\) における値のみあればよいので, 次のネヴェルのアルゴリズムとよばれる方法を使うことができる.
\[
T_{jk} = T_{j(k-1)} + (T_{j(k-1)} – T_{(j-1)(k-1)})/(n_j/n_{j-k+1} – 1)
\] この漸化式を使って求めた \(T_{jk}\) は解の良い近似となっている. この方法は \(p + k – 1\) 次の数値解法となっており, 補外法 (Extrapolation method) とよばれる.

補外法の計算結果を表す次の表を補外表という.
\[
\begin{array}{c|cccc}
T_{11} \\
T_{21} & T_{22} \\
T_{31} & T_{32} & T_{33} \\
T_{41} & T_{42} & T_{43} & T_{44} \\
\dots & \dots & \dots & \dots & \dots \\
数値解 & 補外値 \dots \\
\end{array}
\] \(T_{21}\) については \(T_{22} (k = 2)\) まで, \(T_{31}\) は \(T_{33} ((k = 3)\) まで, \(\dots\), \(T_{j1}\) は \(T_{jj} (k = j)\) まで補外できる.

上において, \(k = 2, n_1 = 1, n_2 = 2\) の場合はリチャードソンの補外に一致する.

整数列 \(n_i\) のとり方であるが, 次のようにいくつかの方法が提案されている.

ロンバーク列: \(1, 2, 4, 8, 16, 32, 64, 128, 256, 512, \dots\)
ブリアシュ列: \(1, 2, 3, 4, 6, 8, 12, 16, 24, 32, \dots\)
調和数列: \(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, \dots\)

数値実験 (12)

最も簡単な微分方程式の初期値問題の例で補外法の精度を確認する.
\[
y’ = y, \space y(0) = 1
\] 数値解法としてはオイラー法 (p = 1), 整数列としてはロンバーク列 1, 2, …, 128 を使った.

横軸は相対誤差 (右に行くほど精度がよい), 縦軸は関数評価回数の対数プロットとした. 補外に必要な計算量は無視した. 一番左のラインがオイラー法による数値解 \(T_{j1} (k = 1)\)である. 一番右の赤いラインは最高次数で求めた補外値 \(T_{jj} (k = j)\)である.

数値実験 (13)

同じ微分方程式の初期値問題の例で補外法のパラメータ (k の値および整数列の種類) を変化させたときの挙動を確認する.
\[
y’ = y, \space y(0) = 1
\] 数値解法としてはオイラー法 (p = 1) を使った. ロンバーク列で k = 1 (オイラー法), 4, j の場合, ブリアシュ列で k = j, 調和数列で k = j の場合を比較する.

横軸は相対誤差 (右に行くほど精度がよい), 縦軸は関数評価回数の対数プロットとした. 補外に必要な計算量は無視した. ロンバーク列で k = 4 とした (補外表の横幅を 4 に制限した) 場合には途中からこの (対数) プロットで直線となり次数が増えなくなったのがわかる. ブリアシュ列と調和数列はこの例では最初はよかったが途中で失速して誤差が減らなくなった.

10.2.2.2 GBS アルゴリズム

補間多項式が偶数のべき乗のみで表される場合, ネヴェルのアルゴリズムは次のようになり収束が速くなる.
\[
T_{jk} = T_{j(k-1)} + \frac{T_{j(k-1)} – T_{(j-1)(k-1)}}{(n_j/n_{j-k+1})^2 – 1}
\] GBS (Gragg-Bulirsch-Stoer) アルゴリズムはこれを利用したもので, 次の手順で計算を行う.
\[
\begin{align}
& y_1 = y_0 + h_jf(t_0, y_0) \space [オイラー法] \\
& y_{i+1} = y_{i-1} + 2h_jf(t_i, y_i) \space (i = 1, 2, \dots, n_j) \space [中点則] \\
& T_{j1} = (1/4)(y_{n_{j-1}} + 2y_{n_j} + y_{n_{j+1}}) \\
\end{align}
\] このようにすると \(T_{j1}\) は偶数のべき乗のみで表され, 上式を適用して効率よく補外値を求めることができる. なお, 中点則には不安定性の問題があることを 10.1.5.2 で説明したが, この場合は問題は生じない.

整数列 \(n_i\) は偶数になるようにする.

ロンバーク列: \(2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, \dots\)
ブリアシュ列: \(2, 4, 6, 8, 12, 16, 24, 32, 48, 64, \dots\)
調和数列: \(2, 4, 6, 8, 10, 12, 14, 16, 18, 20, \dots\)

数値実験 (14)

最も簡単な微分方程式の初期値問題の例で GBS アルゴリズムによる補外法と通常の補外法の精度を比較する.
\[
y’ = y, \space y(0) = 1
\]

横軸は相対誤差 (右に行くほど精度がよい), 縦軸は関数評価回数の対数プロットである. GBS アルゴリズムの方が明らかに速く誤差が減少している. なお, 図で線が切れているところは誤差が 0 になったためプロットされていない区間であり, 精度いっぱいで解が求められていることを示している.

10.2.3 ステップ幅の自動調節

10.2.3.1 埋め込み型ルンゲ・クッタ法

誤差の推定を行うのにリチャードソンの補外を使うためにはステップ幅を変えて 2 回計算を行わなければならず, 計算量が大きくなるのが欠点である. しかし, ルンゲ・クッタ法自身に誤差推定機能を付加した公式を作ることができれば, 1 回の計算で誤差推定値も得ることができ有用である. 特に, あらかじめ要求精度を設定しておき, それに必要なステップ幅に自動調節するプログラムを作るのに必須である.

ルンゲ・クッタ法は次のように表された. ここで, この公式は p 次であるとする.
\[
\begin{align}
& k_i = f(y_n + h \sum_{j=1}^s a_{ij}k_j, t_n + c_i h) \space (i = 1, 2, \dots, s) \\
& y_{n+1} = y_n + h \sum_{i=1}^s b_ik_i \\
\end{align}
\] 同じ \(a_{ij}\) と \(c_i\) を使い (すなわち, 同じ \(k_i\) を使い), \(b_i\) を \(b^*_i\) に変えた時の近似値 \(y^*_{n+1}\) が \(p^*\) 次 (例えば, \(p^* = p – 1, p^* = p + 1\) など) になっているように定数を決めることができれば, \(k_i\) を計算し直すことなく次式により \(p^*\) 次のときの近似値を求めることができる.
\[
y^*_{n+1} = y_n + h \sum_{i=1}^s b^*_i k_i
\] このとき, 次数が低い方の局所離散化誤差を次のように求めることができる.
\[
y_{n+1} – y^*_{n+1} = h \sum_{i=1}^s (b_i – b^*_i)k_i
\] このような公式は実際に作ることができて, 埋め込み型ルンゲ・クッタ法 (embedded Runge-Kutta method) とよばれる. 主な公式を以下に示す.

10.2.3.1.1 ルンゲ・クッタ・マーソン法

埋め込み型ルンゲ・クッタ法で最初に提案されたのはマーソンの公式 (ルンゲ・クッタ・マーソン法 (Runge-Kutta-Merson method)) とされる. 係数は次のとおりである.
\[
\begin{array}{c|ccccc}
0 \\
1/3 & 1/3 \\
1/3 & 1/6 & 1/6 \\
1/2 & 1/8 & 0 & 3/8 \\
1 & 1/2 & 0 & -3/2 & 2 \\
\hline
y & 1/6 & 0 & 0 & 2/3 & 1/6 \\
y^* & 1/2 & 0 & -3/2 & 2 & 0 \\
\end{array}
\] 通常のルンゲ・クッタ法のブッチャー配列に対して \(y^*\) のための係数が 1 行追加されている.

計算式は次のとおりである.
\[
\begin{align}
& k_1 = f(t_n, y_n) \\
& k_2 = f(t_n + h/3, y_n + hk_1/3) \\
& k_3 = f(t_n + h/3, y_n + h(k_1/6 + k_2/6)) \\
& k_4 = f(t_n + h/2, y_n + h(k_1/8 + 3k_3/8)) \\
& k_5 = f(t_n + h, y_n + h(k_1/2 – 3k_3/2 + 2k_4)) \\
& y_{n+1} = y_n + h(k_1 + 4k_4 + k_5)/6 \space (4次) \\
& y^*_{n+1} = y_n + h(k_1 – 3k_3 + 4k_4)/2 \space (誤差評価用 4次) \\
\end{align}
\] この公式による推定誤差は, \(f(t, y)\) が \(t\) と \(y\) に関して線形ならば \(1/5 (y_{n+1} – y^*_{n+1})\) となる.

1 ステップあたり 5 回の関数評価が必要である. すなわち, 誤差評価を行うために, 通常の 4 次のルンゲ・クッタ法に比べ 1 ステップあたり 1 回余計に関数評価を行うという犠牲を払っている.

10.2.3.1.2 5(4) 次 ルンゲ・クッタ・フェールベルグ法

5(4) 次 ルンゲ・クッタ・フェールベルグ法 (Runge-Kutta-Fehlberg method) は埋め込み型ルンゲ・クッタ法の中でも有名な公式である. RKF45 というプログラムが広く使われている. 係数は次のとおりである.
\[
\begin{array}{c|cccccc}
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 \\
\hline
y & 16/135 & 0 & 6656/12825 & 28561/56430 & -9/50 & 2/55 \\
y^* & 25/216 & 0 & 1408/2565 & 2197/4104 & -1/5 & 0 \\
\end{array}
\] この公式では, 近似値 \(y\) の計算を 5 次で行い, 誤差評価用の値 \(y^*\) の計算を 4 次で行う. これを, 5(4) 次の公式と書く. 1 ステップあたり 6 回の関数評価が必要である.

10.2.3.1.3 5(4) 次 ドルマン・プリンス法

5(4) 次 ドルマン・プリンス法 (Dormand-Prince method) の係数は次のとおりである.
\[
\begin{array}{c|ccccccc}
0 \\
1/5 & 1/5 \\
3/10 & 3/40 & 9/40 \\
4/5 & 44/45 & -56/15 & 32/9 \\
8/9 & 19372/6561 & -25360/2187 & 64448/6561 & -212/729 \\
1 & 9017/3168 & -355/33 & 46732/5247 & 49/176 & -5103/18656 \\
1 & 35/384 & 0 & 500/1113 & 125/192 & -2187/6784 & 11/84 \\
\hline
y & 35/384 & 0 & 500/1113 & 125/192 & -2187/6784 & 11/84 & 0 \\
y^* & 5179/57600 & 0 & 7571/16695 & 393/640 & -92097/339200 & 187/2100 & 1/40 \\
\end{array}
\] 5(4) 次の公式で, 1 ステップあたり 7 回の関数評価が必要である. この公式の特長は, 7 番目 (\(k_7\)) の係数 \(a_{7j}\) が近似値 \(y\) を計算するための係数 \(b_i\) に等しいことである. これは FSAL (First Same As Last) とよばれるアイデアで, \(k_7\) が次の \(k_1 = f(y_{n+1}, t_{n+1})\) に等しくなるよう巧妙に作られている.

10.2.3.1.4 6(5) 次 ルンゲ・クッタ・ヴァーナー法

6(5) 次 ルンゲ・クッタ・ヴァーナー法 (Runge-Kutta-Verner method) の係数は次のとおりである.
\[
\begin{array}{c|cccccccc}
0 \\
1/6 & 1/6 \\
4/15 & 4/75 & 16/75 \\
2/3 & 5/6 & -8/3 & 5/2 \\
5/6 & -165/64 & 55/6 & -425/64 & 85/96 \\
1 & 12/5 & -8 & 4015/612 & -11/36 & 88/255 \\
1/15 & -8263/15000 & 124/75 & -643/680 & -81/250 & 2484/10625 & 0 \\
1 & 3501/1720 & -300/43 & 297275/52632 & -319/2322 & 24068/84065 & 0 & 3850/26703 \\
\hline
y & 3/40 & 0 & 875/2244 & 23/72 & 264/1955 & 0 & 125/11592 & 43/616 \\
y^* & 13/160 & 0 & 2375/5984 & 5/16 & 12/85 & 3/44 & 0 & 0 \\
\end{array}
\] 6(5) 次の公式で, 1 ステップあたり 8 回の関数評価が必要である. 古くから広く使われているサブルーチン DVERK に採用されている公式である.

10.2.3.2 多段法の誤差

多段法における局所誤差は次のように求められる.

解法	次数	局所誤差
アダムス-バシュフォース法	\(k\)	\(\gamma_k h^{k+1} y^{k+1}\)
アダムス-ムルトン法	\(k + 1\)	\(\gamma^*_{k+1} h^{k+2} y^{k+2}\)
後退微分公式 (BDF 法)	\(k\)	\(-1/(k+1) h^{k+1} y^{k+1}\)

例えば, 4 次のアダムス法では, 予測子 (右肩に P) と修正子 (右肩に C) のそれぞれの局所誤差は具体的には次のようになる.
\[
\begin{align}
& y(t_{n+1}) – y^P_{n+1} = \gamma_4h^5y^5 + O(h^6) = 251/720 h^5y^5 + O(h^6) \\
& y(t_{n+1}) – y^C_{n+1} = \gamma^*_4h^5y^5 + O(h^6) = -19/720 h^5y^5 + O(h^6) \\
\end{align}
\] このように多段法では予測子と修正子で局所誤差が異なるので, それを利用して誤差の推定を行うことができる.

両辺の引き算をして高次の項を無視すると次の近似式が得られる.
\[
y^C_{n+1} – y^P_{n+1} \simeq 3/8 h^5y^5
\] これを上の修正子の式に代入して高次の項を無視すると次のような誤差評価式が得られる.
\[
y(t_{n+1}) – y^C_{n+1} \simeq -19/27(y^C_{n+1} – y^P_{n+1})
\]

10.2.3.3 ステップ幅の自動調節機能の実装

一定のステップ幅で計算するのではなく, 目標誤差を満たすように最適なステップ幅に自動調節しながら計算するプログラムが現在の主流である. これは, 目標誤差に対して計算量を最小化することと, 計算の安定性を損なわないようにステップ幅を維持することに効果がある.

上に説明した誤差推定機能を持つ公式を使えばそのようなプログラムを実装することができる. その処理手順はおおまかには次のとおりである.

• ステップ幅の初期値と局所離散化誤差の許容限界を与え計算を開始する
• ステップごとに推定誤差をチェックし, 許容内であれば次ステップに進み計算を続ける (ステップ成功)
• 推定誤差が大きすぎるならば, ステップ幅を小さくしてやり直す (ステップ失敗)
• 推定誤差が小さすぎるならば, ステップ数が不要に増えてしまわないよう次ステップからステップ幅を大きくする

このような処理を繰り返し, 目標精度を満たしつつ最小のステップ数で計算を完了するようにする.

ただし, 実際に実用に供するプログラムを実装するのは簡単ではなく多くの細かい工夫が必要である.

多段法の場合には, ステップ幅を変更すると以前の値の間隔が新しいステップ幅と合わなくなり, 公式がそのままでは使えなくなる. そこで, 不等間隔のステップ点に対する公式に変更した可変ステップ幅アダムス法や可変ステップ幅 BDF 法などを使う必要がある.

数値実験 (15)

ルンゲ・クッタ・マーソン法を使い簡単なステップ幅自動調節機能付きプログラム (RKM) を作成して動作を確認する.

RKM の計算アルゴリズムは次のようにごく簡単なものとした.

• マーソンの公式で 1 ステップ計算後, 推定誤差が許容限界より大きければステップ幅を 1/2 にしてやり直す
• 推定誤差が許容限界の 1/32 より小さいときには次ステップからステップ幅を 2 倍にする
• それ以外ならば 1 ステップ進める

2 種類の微分方程式の初期値問題例を使ってステップ幅の自動調節機能の動作を確認する. 比較のために, 4 次のルンゲ・クッタ法 (RK4) および 5(4) 次のルンゲ・クッタ・フェールベルグ法 (RKF45: XLPack の Derkfa を使用) による計算結果も示す.

横軸は相対誤差 (右に行くほど精度がよい. n = 2 の場合はそれぞれの値の平均), 縦軸は関数評価回数の対数プロットである. RKM と RKF45 では要求精度 (Tol) を少しづつ変えて計算を行い, 得られた結果の実際の誤差と関数評価回数をプロットした. RK4 ではステップ幅を変えていき, そのときの誤差と関数評価回数をプロットした. なお, RKM ではステップ幅の初期値を 50*Tol とした.

(例 1)

変化が少ない関数の例で, 解析解は \(y = e^t\) である. \(t = 1\) における値 (\(e\)) を求め誤差を対数プロットした.
\[
y’ = y, \space y(0) = 1
\]

(例 2)

文献[5]のブラセレータとよばれる変化の激しい方程式の例で \(n = 2\) である. \(t = 16\) における値を求め誤差を対数プロットした.
\[
\begin{align}
& y_1′ = 1 + y_1^2 y_2 – 4y_1 \\
& y_2′ = 3y_1 – y_1^2 y_2 \\
& y_1(0) = 1.5, y_2(0) = 3 \\
\end{align}
\]

RKM と RK4 はどちらも 4 次であるが, 例 1 では RK4 の方が速く, 例 2 では RKM の方が速い. これは, 変化が少ない関数においてはステップ幅の調節を行う必要がないが, RKM では調節のために無駄なステップを費やしているためと思われる.

次に, 上の 2 つの例を計算したときの RKM の動作を詳細に見る.

(例 1)

\(Tol = 10^{-7}\) としたときの動作. 上段は計算値, 下段はステップ幅である.

関数が滑らかなので \(h = 0.082\) になった後はステップ幅は変更されなかった. ステップ幅の初期値 (50*Tol) が小さすぎたようなので工夫が必要であろう. なお, 最後の 1 ステップでステップ幅が変化しているのは誤差の問題ではなく, 終点をちょうど 1 に合わせるために調整されたものである.

(例 2)

\(Tol = 10^{-4}\) としたときの動作. 上段は計算値, 下段はステップ幅である.

急な変化をするところではステップ幅を小さく, ゆるやかな変化をするところではステップ幅を大きく調節しているのがわかる. 成功ステップが 97 回で計算を終えているが, 失敗ステップが 25 回あった. すなわち, 2 割ほど無駄ステップを踏んでいるがそれでも RK4 よりは効率がよかったことになる.

この例で RK4 の方が遅かったのは, 全体の精度を確保するためには急な変化をするところで必要なステップ幅に合わせる必要があり, ゆるやかな変化をするところで必要以上にステップ幅が小さくなり全体としては無駄な計算が多かったためと考えられる.

10.2.4 密出力

ステップ幅を自動調節するプログラムにおいて出力点の間隔が狭いときに効率が悪くなることがある. そのような場合に無駄な計算を避ける方法について述べる.

数値実験 (16)

次の初期値問題の解を求める.
\[
y’ = y, \space y(0) = 1
\] 解析解は \(y = e^t\) である.

ステップ幅を自動調節するプログラム Dopri5a (5(4)次 ドルマン・プリンス法) と Derkfa (5(4) 次 ルンゲ・クッタ・フェールベルグ法) を使って, 要求精度 RTol = ATol = \(10^{-7}\) として, t = 0 ～ 1 で 0.05 ごとに 0.0 ～ 0.05, 0.05 ~ 0.1, …, 0.95 ～ 1.0 というように終点を変えながら繰り返し解を求めた. このときの誤差をプロットすると次のようになった. 赤の実線と左側の目盛は関数形を表す. 点線と右側の目盛は各解法による誤差を表す.

Dopri5a も Derkfa も 0.05 幅の各区間では 1 ステップで解を求め, 誤差は \(10^{-10}\) 付近と要求精度よりかなりよい精度で解を求めた. これは, 自動調節されるべき要求精度に見合ったステップ幅よりも 0.05 という幅が小さく, 終点に合わせるためにステップ幅が常に 0.05 に調節されたためである. すなわち, 要求精度に比べて公式の精度が良すぎるのである. 比較のために 4 次のルンゲ・クッタ法 (RK4) でステップ幅を 0.05 として計算してみると, 誤差はちょうど \(10^{-7}\) 付近になり, RK4 の精度で十分であることがわかる. 関数評価回数は, Dopri5a が 128 回, Derkfa が 124 回で, RK4 は 80 回となった. Dopri5a と Derkfa はこの場合には高精度なゆえに無駄な計算をしていることになる.

そこで, ステップ幅は必要精度の観点から最適なものを選び, ステップに一致しない出力点においては計算を行わずに補間値を使うことにより無駄な計算を避ける方法が考えられる. これを密出力 (dense output) とよぶ.

10.2.4.1 多段法の場合

アダムス・ムルトン法を例にとる. アダムス・ムルトン法の説明で使用した図を再掲する.

\(t_n\) の次のステップ \(t_{n+1} = t_n + h\) における値を求める際には図のように補間多項式 \(p^*(t)\) を求め, これを区間 \([t_n, t_{n+1}]\) で積分することにより値を求めた.
\[
y_{n+1} = y_n + \int_{t_n}^{t_{n+1}} p^*(x) dx
\] 同様に, \(t_n + \theta h (0 \le \theta \le 1)\) における値 \(y(t_n + \theta h)\) は同じ補間式を区間 \([t_n, t_n + \theta h]\) で積分すれば求めることができる. 新たに関数評価を行う必要がないし精度が落ちることもない.

数値実験 (17)

上の例題において 4 次アダムス・ムルトン法で h = 0.1 として t = 0.4, 0.5, …, 1.0 における y の値を求める. t = 0.1, 0.2, 0.3 における値 (出発値) は 4 次のルンゲ・クッタ法で求めた. t = 0.35, 0.45, …, 0.95 における値は上で説明した補間により求める.

青い点がアダムス・ムルトン法による計算値で \(y = e^t\) のグラフになる. オレンジ色の点は補間により求めた点を示す. 緑は計算値および補間値の相対誤差を右側の対数目盛りにより表す. 期待通りに, 誤差が増えることなく補間値が求められている.

10.2.4.2 連続ルンゲ・クッタ法

ルンゲ・クッタ法は 1 段法であり過去の値は使用しないためすぐに使える補間式はない. そこで, 負担の少ない, すなわちこのための追加の関数評価がないまたはごく少数回であるような, 密出力のための補間公式が研究されてきた. 例えば, \(t_n + \theta h (0 \le \theta \le 1)\) における値 \(y(t_n + \theta h)\) を求める次の形の式である.
\[
\begin{align}
& y(t_n + \theta h) = y_n + h \sum_{i=1}^{s^*} b_i(\theta)k_i \\
k& _i = f(t_n + c_ih, y_n + h \sum_{j=1}^{i-1} a_{ij}k_j) \space (i = 1, 2, \dots, s^*) \\
\end{align}
\] ルンゲ・クッタ法の公式と同じ形をしているが, 係数 \(b_i\) は定数ではなく \(\theta\) に依存する. このようなルンゲ・クッタ法の自然な拡張を連続ルンゲ・クッタ法 (Continuous Runge-Kutta (CRK) method) という.

4 次のルンゲ・クッタ法については, \(s^* = s = 4\) に対して密出力公式の次数 \(p^* = 3\) となる次の解が得られている.
\[
\begin{align}
& b_1(\theta) = \theta – (3/2)\theta^2 + (2/3)\theta^3 \\
& b_2(\theta) = b_3(\theta) = \theta^2 – (2/3)\theta^3 \\
& b_4(\theta) = -(1/2)\theta^2 + (2/3)\theta^3 \\
\end{align}
\] \(k_i\) は 4 次ルンゲ・クッタ法の公式と同じである. すなわち, 新たな関数評価を必要としない. なお, 4 次ルンゲ・クッタ法では関数評価を追加しなければ \(p^* = 4\) にはできない.

通常は \(p^* = p – 1\) であれば十分であるが, 導関数 \(y'(t_n + \theta h)\) が正確な必要がある場合には \(p^* = p\) になるようにしなければならない.

数値実験 (18)

上の例題において 4 次ルンゲ・クッタ法で h = 0.1 として t = 0.1, 0.2, …, 1.0 における y の値を求める. t = 0.05, 0.15, …, 0.95 における値は上で示した 3 次の密出力公式を使って求める.

青い点がルンゲ・クッタ法による計算値で \(y = e^t\) のグラフとなる. オレンジ色の点は補間により求めた点を示す. 緑は計算値および補間値の相対誤差を右側の対数目盛りにより表す. 補間値の場合は少し誤差が増える.

10.2.4.3 ステップ幅の自動調節を行う場合

上の 2 つの数値実験でみた補間を使う密出力は, ステップ幅の自動調節を行うプログラムに適用すると効率を上げることができる. すなわち, ステップ幅は必要精度の面から最適なものを選び, ステップに一致しない出力点においては補間値を使うことにより無駄な関数計算を避けることができる.

多段法の場合, 例えば, XLPack に収録されている Deabm (アダムス・バシュフォース・ムルトン法) は, 特に指定しない限り補間値を使うように作られている.

ルンゲ・クッタ法の場合, 例えば, XLPack の Derkfa には 5 次の補間公式を使った密出力機能が組み込まれており, Dopri5a には 4 次の補間公式を使った密出力機能が組み込まれている. これらのプログラムでは密出力のための追加関数評価を必要としない補間公式が採用されている. XLPack に収録されている他の多くのプログラムにも密出力機能が組み込まれている.

数値実験 (19)

数値実験 (16) の例について, Derkfa, Dopri5a および Deabm を使って要求精度 RTol = ATol = \(10^{-7}\) を指定して計算を行い, 0.05 ごとに結果の出力を行う. その際に密出力をオンにするとどう変わるか確認する.

それぞれの図の上段は関数評価回数のプロットである. 青は通常動作, オレンジ色は密出力を行った場合である. どのプログラムでも密出力を行った場合には関数評価回数がおおよそ半分に減っている.

下段は計算結果の相対誤差を表す. 通常動作では, 要求精度 \(10^{-7}\) を満たすステップ幅が 0.05 より大きくなるため, 結果的にステップはその終わりが出力を行う点に調整されて 0.05 にカットされてしまい, 必要以上の精度 \(10^{-9} \sim 10^{-11}\) が出ている. 密出力を行った場合には要求精度に対する最適ステップ幅に調整されているため, 要求どおり \(10^{-7}\) 付近になっている.

密出力においては補間値をどれだけ細かな間隔で求めても関数評価回数は増えない (補間計算程度の計算量しか増えない) ので, 精密なグラフを描くために細かく値を出力したい場合などに有効である. 逆に, 必要なステップ幅が出力点の間隔よりも小さくなる場合には密出力の効果は小さくなるので, 単純に繰り返し計算を行ってもよい.

なお, Deabm の場合には指定しない限り密出力を行うような実装になっている.

10.2.5 2 階常微分方程式

実際に応用で現れる方程式は 2 階常微分方程式が多い.
\[
\begin{align}
& y^{\prime\prime} = f(t, y, y’) \\
& y(t_0) = y_0, \space y'(t_0) = y’_0 \space (初期条件) \\
\end{align}
\] これは「はじめに」の項で説明したように 1 階の連立常微分方程式に変換して解くことができる.

しかし, 次のように \(f()\) が \(y’\) に依存しないような特別な形をした場合には 1 階の連立常微分方程式に変換することなく直接解くことができる方法がある.
\[
\begin{align}
& y^{\prime\prime} = f(t, y) \\
& y(t_0) = y_0, \space y'(t_0) = y’_0 \space (初期条件) \\
\end{align}
\]

10.2.5.1 ニュストレム法

ルンゲ・クッタ法に似ているが係数の決め方が異なるニュストレム法 (Nystrom method) が知られており, 次のように計算を行う.
\[
\begin{align}
& k’_i = f(t_n + c_ih, y_n + c_ihy’_n + h^2 \sum_{j=1}^s \bar{a}_{ij}k’_j) \space (i = 1, 2, \dots, s) \\
& y_{n+1} = y_n + hy’_n + h^2 \sum_{i=1}^s \bar{b}_ik’_i \\
& y’_{n+1} = y’_n + h \sum_{i=1}^s b_i k’_i \\
\end{align}
\] 係数は以下のとおりである.

4 次ニュストレム法

\[
\begin{array}{c|ccc}
0 & & & \bar{a}_{ij} \\
1/2 & 1/8 \\
2/3 & 0 & 1/2 \\
\hline
\bar{b}_i & 1/6 & 1/3 & 0 \\
b_i & 1/6 & 4/6 & 1/6 \\
\end{array}
\]

5 次ニュストレム法

\[
\begin{array}{c|cccc}
0 & & & & \bar{a}_{ij} \\
1/5 & 1/50 \\
2/3 & -1/27 & 7/27 \\
1 & 3/10 & -2/35 & 9/35 \\
\hline
\bar{b}_i & 14/336 & 100/336 & 54/336 & 0 \\
b_i & 14/336 & 125/336 & 162/336 & 35/336 \\
\end{array}
\] ニュストレム法はルンゲ・クッタ法に比べて段数が少なくてすむ (4次: 4段 → 3段, 5次: 6段 → 4段).

10.2.5.2 補外法

1 階の連立常微分方程式に変換して GBS アルゴリズムを適用する.
\[
\begin{align}
& y_1 = y_0 + h_jy’_0 \\
& y’_1 = y’_0 + h_jf(t_0, y_0) \\
\\
& y_{i+1} = y_{i-1} + 2h_jy’_i \space (i = 1, 2, \dots, n_j) \\
& y’_{i+1} = y’_{i-1} + 2h_jf(t_i, y_i) \space (i = 1, 2, \dots, n_j) \\
\\
& T_{j1} = (1/4)(y_{n_{j-1}} + 2y_{n_j} + y_{n_{j+1}}) \\
& T’_{j1} = (1/4)(y’_{n_{j-1}} + 2y’_{n_j} + y’_{n_{j+1}}) \\
\end{align}
\] この場合, 次のようになるので半分だけ (奇数添字のときだけ) 関数評価をすればよくなる.
\[
\begin{align}
& T_{j1} = y_{n_j} \\
& T’_{j1} = (1/2)(y’_{n_{j-1}} + y’_{n_{j+1}}) \\
\end{align}
\]

数値実験 (20)

次の 2 階常微分方程式の初期値問題の例を 4 次, 5 次ニュストレム法および補外法 (GBS, ロンバーク列) により計算する.
\[
\begin{align}
& y^{\prime\prime} = y \\
& y(0) = 1, y'(0) = 1 \\
\end{align}
\]

横軸は相対誤差 (右に行くほど精度がよい), 縦軸はステップ数の対数プロットである. 比較のため 4 次のルンゲ・クッタ法の計算結果も示した.

4 次ニュストレム法は段数は少ないが 4 次のルンゲ・クッタ法と同程度の精度であることが確認できる. 関数評価回数が少ない分だけプロットは下にきている.

10.2.6 遅延微分方程式

遅延微分方程式 (Delay differential equations (DDE)) とは次のような微分方程式をいう.
\[
y'(t) = f(t, y(t), y(t – \tau))
\] \(\tau\) は遅延時間を表す正定数である. 複数の遅延時間が含まれるなどより一般的なものも考えられる.

これをルンゲ・クッタ法で解くためにはステップ幅を \(\tau = kh\) (k は整数) となるように選んで, 計算した \(y(t – \tau)\) の値をとっておくようにすればよい. しかし, ステップ幅を自動調節するプログラムではうまくいかない. その場合, 多段法か密出力機能を持ったルンゲ・クッタ法で過去の計算値をとっておくようにすれば計算できる.

XLPack に収録されている Retarda はドルマン・プリンス法プログラム Dopri5a を遅延微分方程式用に修正したものである. Dopri5a では 1 ステップ前までの過去の値を求められたが, Retarda ではずっと前までの値を得ることができる.

遅延微分方程式の計算例は「10.4 実用プログラムのベンチマーク」を参照せよ.