4.7.2 RKF45 による刻み幅の自動調節 (書き直し版、工事中)

$t=a$ から $t=b$ まで、誤差の大きさの見積が、 計算者が指定した値 $\eps$ より小さくなるように解くことを目標にする。 このとき $\eps$ のことを許容誤差限界と呼ぶ。 $[a,b]$ を $a=t_0<t_1<\cdots<t_N=b$ と分割して解くことにする。 $t=t_j$ における値 $x_j$ まで定まったとする。 $t_{j+1}$ における 5次公式, 4 次公式による近似値をそれぞれ $x_{j+1}$ , $x_{j+1}^*$ とすると、

$$x_{j+1}-x(t_{j+1})=c h^6+O(h^7),\quad x_{j+1}^*-x(t_{j+1})=c' h^5+O(h^6)$$

となる。ただし $h:=t_{j+1}-t_j$ とおいた。これから

$$x_{j+1}-x_{j+1}^*=-c' h^5+O(h^6).$$

通常は、$x_{j+1}$ は $x_{j+1}^*$ よりも格段に精度が良いと期待出来るので、 この式の値が $x_{j+1}^*$ の誤差の良い評価となっていると考えられる:

$$x_{j+1}-x_{j+1}^* =\left(x_{j+1}-x(t_{j+1})\right) +\left(x(t_{j+1})-x_{j+1}^*\right) \kinji x(t_{j+1})-x_{j+1}^*.$$

さて、

$$\delta_{j+1}:=\left\Vert x_{j+1}-x_{j+1}^*\right\Vert$$

と置いておく。

目標は $[a,b]$ で解いたときの許容誤差限界を $\eps$ とすることであるから、 $[t_{j},t_{j+1}]$ で解いたときの許容誤差限界は

$$\eps\cdot \frac{h}{H},\quad H:=b-a$$

とするのが妥当であろう。それゆえ

$$\delta_{j+1}\le \frac{\eps h}{H}$$

であれば良いが、 そうでない場合は、$h$ が大きすぎて十分な精度が得られていないと考え、 もっと $h$ を小さくすることにする。 $h$ の代りに、 $\overline h\ll h$ なる $\overline h$ を用いて

$$\overline t_{j+1}=t_j+\overline h$$

とおき、これに対応した $\overline x_{j+1}$ , $\overline x_{j+1}^*$ を求め、 $\overline\delta_{j+1}:= \left\Vert\overline x_{j+1}-\overline x_{j+1}^*\right\Vert$ とおく。

$$\delta_{j+1}\sim \Vert c'\Vert\left\vert h\right\vert^5,\quad \overline\delta_{j+1}\sim \Vert c'\Vert\left\vert\overline h\right\vert^5$$

となると期待できるから、

$$\overline\delta_{j+1}\sim \left( \frac{\left\vert\overline h\right\vert}{\left\vert h\right\vert} \right)^5 \delta_{j+1}.$$

今度は $\overline \delta_{j+1}\le \eps \overline{h}/H$ となって欲しいわけだが、 この不等式の左辺に推定値を代入した不等式

$$\left( \frac{\left\vert\overline h\right\vert}{\left\vert h\right\vert} \right)^5 \delta_{j+1} \le \frac{\eps \overline h}{H}$$

を $\overline h$ について解くと

$$\overline h \le h \left(\frac{\eps h}{H\delta_{j+1}}\right)^{1/4}$$

を得る。安全のために

$$\overline h:= 0.9 h \left(\frac{\eps h}{H \delta_{j+1}}\right)^{1/4}$$

で $\overline h$ を定めることにする。