5.5.1.0.3 埋め込み型の公式

E. Fehlberg (1970) は次の公式 RKF45 を提案した^5.4。

$$x_{j+1} = x_j + h\left( \frac{16}{135}k_1 +\frac{6656}{... ...28561}{56430}k_4 -\frac{9}{50}k_5 +\frac{2}{55}k_6 \right),$$ (5.4)

ただし

$$\begin{eqnarray*} k_1 &=& f(t_j,x_j) \\ k_2 &=& f\left(t_j+\frac14 h, x_j+\frac... ...65}k_3 +\frac{1859}{4104}k_4-\frac{11}{40}k_5 \right) \right) \end{eqnarray*}$$

これは $6$ 段 $5$ 次の公式であるが、それだけでなく

$$x_{j+1}^\ast=x_j+h \left( \frac{25}{216}k_1+\frac{1408}{2565}k_3+\frac{2197}{4104}k_4-\frac15k_5 \right)$$ (5.5)

という値を作ると、$x_{j+1}^\ast$ は $x(t_{j+1})$ に対して $4$ 次の近似 値となる。この次数の差を利用してステップ幅の自動調節をしよう、というア イディア。

いま (5.4), (5.5) の局所離散化誤差 をそれぞれ $\tau$, $\tau^\ast$ とおくと

$$\tau(t,h) = C(t)h^5+O(h^6)$$ (5.6)

$$\tau^\ast(t,h) = C^\ast(t)h^4+O(h^5)$$

である。あらかじめ許容誤差限界 $\eps_{\rm TOL}$ を決めておき、第 $j$ ステップでは

$$\Vert x_{j+1}^\ast-x_{j+1}\Vert\le \eps_{\rm TOL}$$

が成り立ったとする。これは (5.6) より高次項を無視して

$$\Vert C^\ast(t_j)h^5\Vert\le \eps_{\rm TOL}$$

が成り立つことを意味する。そのとき、次のステップにおけるステップ幅 $\hat h$ は、やはり

$$\Vert C^\ast(t_{j+1})\hat h^5\Vert= \eps_{\rm TOL}$$

となるように選ぶべきであろう。ステップ幅が小さいときにはもっともである と思われる $C^\ast(t_j)=C^\ast(t_{j+1})$ を仮定して

$$\Vert C^\ast(t_{j})\hat h^5\Vert= \eps_{\rm TOL}.$$

ところで近似的に

$$\Vert C^\ast(t_j)\Vert=\frac{\Vert x_{j+1}^\ast-x_{j+1}\Vert}{h^5}$$

とみなせるから、これを代入して

$$\frac{\Vert x_{j+1}^\ast-x_{j+1}\Vert}{h^5} \hat h^5= \eps_{\rm TOL}.$$

すなわち

$$\hat h=h\sqrt[5]{\frac{\eps_{\rm TOL}}{\Vert x_{j+1}^\ast-x_{j+1}\Vert}}.$$ (5.7)

すなわち第 $j+1$ ステップでは、(5.7) で定められ る ステップ幅 $\hat h$ で公式を適用すれば、 $\Vert x_{j+2}-x_{j+2}^\ast\Vert$ も $\eps_{\rm TOL}$ の限界内にあるであろう、 と考える。実際には安全率を見込んで

$$\hat h=\alpha h\sqrt[5]{\frac{\eps_{\rm TOL}}{\Vert x_{j+1}^\ast-x_{j+1}\Vert}}.$$ (5.8)

とする。ここで $\alpha$ は $1$ より小さい正数で普通 $0.8\sim 0.9$ とす る。こうして、次のアルゴリズムが得られる。

許容限界 $\eps_{\rm TOL}$ を定めて、 (5.4), (5.5) により $x_{j+1}$, $x_{j+1}^\ast$ を求めよ。
次のステップ幅を(5.8) で求めて、計算を続行せよ。

このアイディアのキーは、 $s$ 段 $m$ 次公式に、$f$ の値を計算することな く $(m-1)$ 次公式を付随させるところにある。このような Runge-Kutta 型 公式を、$(m-1)$ 次公式が $m$ 次公式に埋め込まれているといい、埋め込み 型 Runge-Kutta 法と呼ぶ。