適応的時間刻みとイベント検出 | 数値解析 I

1 一歩と二半歩

命題 1.1. $d,p\in\NN$ とし、 $\R^d$ にノルム $\|\cdot\|$ を固定する。 $\Omega\subseteq\R\times\R^d$ を開集合、 $f\colon\Omega\to\R^d$ を連続写像、 $\Phi\colon D\to\R^d$ を $f$ に対する増分関数、 $\Psi_h(t,u)=u+h\Phi(t,u,h)$ をその一歩写像とする。 $t_0<T$ とし、 $J\supseteq[t_0,T]$ を開区間、 $y\colon J\to\R^d$ を、任意の $\tau\in J$ について $(\tau,y(\tau))\in\Omega$ と $y'(\tau)=f(\tau,y(\tau))$ を満たす $C^1$ 級写像とし、 $\delta_y$ を $y$ の局所打切り誤差とする。 $H_0,\rho>0$ 、 $C_1,L_c,\Lambda\ge0$ と写像 $c\colon[t_0,T]\to\R^d$ が次を満たすとする。任意の $s,t\in[t_0,T]$ について $\|c(s)-c(t)\|\le L_c|s-t|$ である。 $t\in[t_0,T)$ と $0<h\le\min\{H_0,T-t\}$ を満たす任意の $t,h$ と、 $\|u-y(t)\|\le\rho$ を満たす任意の $u\in\R^d$ について、 $(t,u,h)\in D$ であり、

\bigl\|\delta_y(t,h)-c(t)h^{p+1}\bigr\|\le C_1h^{p+2},\qquad\bigl\|\Phi(t,u,h)-\Phi(t,y(t),h)\bigr\|\le\Lambda\|u-y(t)\|

が成り立つ。 $M_c:=\max_{s\in[t_0,T]}\|c(s)\|$ 、

C_2:=2^{-p-2}\Bigl(2C_1+L_c+\Lambda\bigl(M_c+\tfrac12C_1H_0\bigr)\Bigr)

と置き、 $H_1\in(0,H_0]$ を $\bigl(M_c+\frac12C_1H_1\bigr)(H_1/2)^{p+1}\le\rho$ を満たす数とする。 $t\in[t_0,T)$ と $0<h\le\min\{H_1,T-t\}$ に対して

Y_1:=\Psi_h(t,y(t)),\qquad Z:=\Psi_{h/2}(t,y(t)),\qquad Y_2:=\Psi_{h/2}\bigl(t+\tfrac h2,Z\bigr)

と置く。

$Y_1,Z,Y_2$ はすべて定まり、 $\bigl\|y(t+h)-Y_2-2^{-p}c(t)h^{p+1}\bigr\|\le C_2h^{p+2}$ である。
$\bigl\|Y_2-Y_1-(1-2^{-p})c(t)h^{p+1}\bigr\|\le(C_1+C_2)h^{p+2}$ である。
$\hat e:=(Y_2-Y_1)/(2^p-1)$ と置くと、 $\bigl\|y(t+h)-Y_2-\hat e\bigr\|\le\frac{2^pC_2+C_1}{2^p-1}h^{p+2}$ である。
$c(t)\ne0$ かつ $C_2h\le2^{-p-1}\|c(t)\|$ ならば、 $y(t+h)\ne Y_2$ であり、 $\frac{\|\hat e-(y(t+h)-Y_2)\|}{\|y(t+h)-Y_2\|}\le\frac{2^{p+1}(2^pC_2+C_1)}{(2^p-1)\|c(t)\|}h$ である。

証明.(1)を示す。 $t':=t+\frac h2$ と置く。 $t'\in[t_0,T)$ であり、 $h\le T-t$ から $\frac h2\le\min\{H_0,T-t'\}$ である。仮定を $(t,h)$ と $(t,\frac h2)$ に $u=y(t)$ として用いると、 $Y_1$ と $Z$ は定まり、 $y(t')-Z=\delta_y(t,\frac h2)$ である。 $h\le H_1$ により

\bigl\|\delta_y(t,\tfrac h2)\bigr\|\le\|c(t)\|(\tfrac h2)^{p+1}+C_1(\tfrac h2)^{p+2}\le\bigl(M_c+\tfrac12C_1H_1\bigr)(\tfrac{H_1}2)^{p+1}\le\rho

であるから、仮定を $(t',\frac h2)$ と $u=Z$ に用いると、 $Y_2$ は定まり、 $\|\Phi(t',Z,\frac h2)-\Phi(t',y(t'),\frac h2)\|\le\Lambda\|\delta_y(t,\frac h2)\|$ である。 $y(t+h)=\Psi_{h/2}(t',y(t'))+\delta_y(t',\frac h2)$ と

\Psi_{h/2}(t',y(t'))-Y_2=y(t')-Z+\tfrac h2\bigl(\Phi(t',y(t'),\tfrac h2)-\Phi(t',Z,\tfrac h2)\bigr)

から、 $R:=\frac h2\bigl(\Phi(t',y(t'),\frac h2)-\Phi(t',Z,\frac h2)\bigr)$ と置くと

y(t+h)-Y_2=\delta_y(t',\tfrac h2)+\delta_y(t,\tfrac h2)+R,\qquad\|R\|\le\tfrac h2\Lambda\bigl\|\delta_y(t,\tfrac h2)\bigr\|

である。 $2c(t)(\frac h2)^{p+1}=2^{-p}c(t)h^{p+1}$ であるから

y(t+h)-Y_2-2^{-p}c(t)h^{p+1}=\bigl(\delta_y(t',\tfrac h2)-c(t')(\tfrac h2)^{p+1}\bigr)+\bigl(c(t')-c(t)\bigr)(\tfrac h2)^{p+1}+\bigl(\delta_y(t,\tfrac h2)-c(t)(\tfrac h2)^{p+1}\bigr)+R

であり、右辺の四つの項のノルムは、仮定と $h\le H_0$ により、それぞれ $C_1(\frac h2)^{p+2}$ 、 $L_c(\frac h2)^{p+2}$ 、 $C_1(\frac h2)^{p+2}$ 、 $\Lambda\bigl(M_c+\frac12C_1H_0\bigr)(\frac h2)^{p+2}$ 以下である。この四つの上界の和は $C_2h^{p+2}$ である。

(2)を示す。 $Y_2-Y_1=\delta_y(t,h)-\bigl(y(t+h)-Y_2\bigr)$ であるから

Y_2-Y_1-(1-2^{-p})c(t)h^{p+1}=\bigl(\delta_y(t,h)-c(t)h^{p+1}\bigr)-\bigl(y(t+h)-Y_2-2^{-p}c(t)h^{p+1}\bigr)

であり、仮定と(1)により右辺のノルムは $(C_1+C_2)h^{p+2}$ 以下である。

(3)を示す。 $D_2:=y(t+h)-Y_2$ と置くと $\hat e=(\delta_y(t,h)-D_2)/(2^p-1)$ であるから

D_2-\hat e=\frac{2^pD_2-\delta_y(t,h)}{2^p-1}=\frac{2^p\bigl(D_2-2^{-p}c(t)h^{p+1}\bigr)-\bigl(\delta_y(t,h)-c(t)h^{p+1}\bigr)}{2^p-1}

であり、(1)と仮定により分子のノルムは $(2^pC_2+C_1)h^{p+2}$ 以下である。

(4)を示す。(1)と $C_2h\le2^{-p-1}\|c(t)\|$ により

\|D_2\|\ge2^{-p}\|c(t)\|h^{p+1}-C_2h^{p+2}\ge2^{-p-1}\|c(t)\|h^{p+1}>0

である。(3)の右辺をこの下界で割ると主張の評価を得る。▨

例 1.2. $d=1$ 、 $\Omega=\R^2$ とし、前進 Euler 法の一歩写像 $\Psi_h(t,u)=u+hf(t,u)$ に命題 1.1の記号 $Y_1,Z,Y_2,\hat e$ を $p=1$ として用いる。このとき $\hat e=Y_2-Y_1$ である。

$f(t,u)=u$ 、 $y(\tau)=e^\tau$ 、 $t=0$ とする。 $Y_1=1+h$ 、 $Z=1+\frac h2$ 、 $Y_2=(1+\frac h2)^2$ であるから $\hat e=h^2/4$ であり、 $y(h)-Y_2-\hat e=e^h-1-h-\frac{h^2}2=\sum_{j\ge3}\frac{h^j}{j!}>0\qquad(h>0)$ である。 $h=0.1$ では $\hat e=0.0025$ 、 $y(h)-Y_2=e^{0.1}-1.1025=0.002670918075\ldots$ である。絶対許容誤差 $0.0026$ と比べると、 $\hat e/0.0026=0.9615\ldots\le1$ である一方、 $y(h)-Y_2>0.0026$ である。
$f(t,u)=t^2$ 、 $y(\tau)=y(0)+\tau^3/3$ とする。 $D=\Omega\times(0,\infty)$ 、 $\Phi(t,u,h)=t^2$ であり、 $\delta_y(t,h)=\bigl((t+h)^3-t^3\bigr)/3-ht^2=th^2+h^3/3$ である。したがって $t_0\le0<T$ を満たす任意の $t_0,T$ と任意の $H_0,\rho>0$ について、 $c(t)=t$ 、 $C_1=\frac13$ 、 $L_c=1$ 、 $\Lambda=0$ が命題 1.1の仮定を満たし、 $C_2=2^{-3}(\frac23+1)=\frac5{24}$ である。 $t=0$ では $c(0)=0$ であり、 $Y_1=Z=y(0)$ 、 $Y_2=y(0)+\frac h2(\frac h2)^2$ から $\hat e=\frac{h^3}8,\qquad y(h)-Y_2=\frac{h^3}3-\frac{h^3}8=\frac{5h^3}{24}$ である。命題 1.1 (1)の評価は等号で成り立つ。 $\bigl(\hat e-(y(h)-Y_2)\bigr)/(y(h)-Y_2)=-\frac25$ は $h$ によらないから、 $c(0)=0$ であるこの場合には命題 1.1 (4)の左辺は $h\to0$ で $0$ に近づかない。
$f(t,u)=t(1-2t)$ 、 $y(\tau)=y(0)+\tau^2/2-2\tau^3/3$ 、 $t=0$ 、 $h=1$ とする。 $f(0,u)=f(\frac12,u)=0$ であるから $Y_1=Z=Y_2=y(0)$ であり、 $\hat e=0$ である。一方 $y(1)-Y_2=\frac12-\frac23=-\frac16$ である。

注意 1.3. $s\in\NN$ とし、 $A\in\R^{s\times s}$ を狭義下三角行列、 $b,\hat b\in\R^s$ とする。連続写像 $f\colon\Omega\to\R^d$ に Butcher 配列 $(A,b)$ と $(A,\hat b)$ の陽的 Runge–Kutta 法を適用すると、二つの方法は段の値 $k_1,\dots,k_s$ を共有し、一歩写像 $\Psi_h$ と $\hat\Psi_h$ の差は $\hat\Psi_h(t,u)-\Psi_h(t,u)=h\sum_{i=1}^s(\hat b_i-b_i)k_i$ である。 $y$ を $f$ の解とし、両者が $(t,y(t),h)$ で定まるとき、局所打切り誤差を $\delta_y$ 、 $\hat\delta_y$ とすると $\hat\Psi_h(t,y(t))-\Psi_h(t,y(t))=\delta_y(t,h)-\hat\delta_y(t,h)$ である。したがって $p\in\NN$ と $\hat C\ge0$ について $\|\hat\delta_y(t,h)\|\le\hat Ch^{p+2}$ ならば、 $\hat\Psi_h(t,y(t))-\Psi_h(t,y(t))$ と $\delta_y(t,h)$ の差のノルムは $\hat Ch^{p+2}$ 以下である。

この推定は一回の試行で $f$ を $s$ 回評価する。 $(A,b)$ の二半歩による推定で $Y_1$ 、 $Z$ 、 $Y_2$ の $s$ 個の段をすべて評価するとき、 $A$ が狭義下三角行列であることから $Y_1$ と $Z$ の第一段はともに $f(t,u)$ であり、この第一段だけを共有すると評価回数は $3s-1$ である。 $s=2$ 、 $a_{21}=1$ 、 $b=(1,0)$ 、 $\hat b=(\frac12,\frac12)$ とすると、 $(A,b)$ の一歩は前進 Euler 法の一歩 $u+hk_1$ に等しく、 $(A,\hat b)$ は Heun 法であり、§E20.27 系 5.2により前進 Euler 法と Heun 法の局所次数はそれぞれ $1$ と $2$ である。このとき差は $\frac h2(k_2-k_1)$ であり、評価回数は $2$ である。

2 刻み制御

定義 2.1. $d,p\in\NN$ 、 $t_0<T$ 、 $u_0\in\R^d$ とする。集合 $\hat D\subseteq\R\times\R^d\times(0,\infty)$ と、 $(t,u,k)\in\hat D$ に対して $\hat\Psi_k(t,u)\in\R^d$ と $\hat e(t,u,k)\in\R^d$ を定める写像を取る。 $a_1,\dots,a_d>0$ 、 $r\ge0$ 、 $\theta\in(0,1)$ 、 $\alpha_{\min}\in(0,1)$ 、 $\alpha_{\max}>1$ 、 $h_{\min}>0$ 、 $h_{\mathrm{init}}>0$ とする。

$u,v,e\in\R^d$ に対して $s_i:=a_i+r\max\{|u_i|,|v_i|\}$ （ $1\le i\le d$ ）と置き、 $E(u,v,e):=\max_{1\le i\le d}\frac{|e_i|}{s_i}$ を尺度付き誤差 (scaled error) という。 $a_i>0$ により $s_i>0$ である。 $a_i$ を成分 $i$ の絶対許容誤差、 $r$ を相対許容誤差という。
$E\ge0$ に対して、 $E=0$ ならば $q(E):=\alpha_{\max}$ 、 $E>0$ ならば $q(E):=\min\Bigl\{\alpha_{\max},\ \max\bigl\{\alpha_{\min},\ \theta E^{-1/(p+1)}\bigr\}\Bigr\}$ と置く。
状態 $(t,u,h)$ を $(t_0,u_0,h_{\mathrm{init}})$ で始め、次の規則を停止するまで繰り返す。 $t=T$ ならば停止する。 $t<T$ かつ $h<h_{\min}$ ならば停止する。それ以外の場合は $k:=\min\{h,T-t\}$ と置き、刻み $k$ の試行 (trial step) を行う。 $(t,u,k)\notin\hat D$ ならば、 $(t,u)$ を変えずに $h$ を $\alpha_{\min}k$ に替える。 $(t,u,k)\in\hat D$ ならば、 $v:=\hat\Psi_k(t,u)$ 、 $E:=E(u,v,\hat e(t,u,k))$ と置き、 $E\le1$ ならば $(t,u,h)$ を $(t+k,v,q(E)k)$ に替え、 $E>1$ ならば $(t,u)$ を変えずに $h$ を $q(E)k$ に替える。 $E\le1$ により $(t,u)$ を替える試行を受理 (acceptance) といい、それ以外の試行を棄却 (rejection) という。この手続きを刻み制御 (step-size control) という。
$\Phi\colon D\to\R^d$ を増分関数とし、 $\Psi_h$ をその一歩写像とする。 $(t,u,k)\in D$ 、 $(t,u,\frac k2)\in D$ 、 $(t+\frac k2,\Psi_{k/2}(t,u),\frac k2)\in D$ を満たす $(t,u,k)$ の全体を $\hat D$ とし、 $\hat\Psi_k(t,u):=\Psi_{k/2}\bigl(t+\tfrac k2,\Psi_{k/2}(t,u)\bigr),\qquad\hat e(t,u,k):=\frac{\hat\Psi_k(t,u)-\Psi_k(t,u)}{2^p-1}$ と置いた刻み制御を、 $\Psi$ の二半歩による刻み制御 (step-size control by step doubling) という。

注意 2.2. $(t,u)$ を固定し、 $\kappa>0$ と区間 $I_0\subseteq(0,\infty)$ について、刻み $k\in I_0$ の試行の尺度付き誤差が $E(k)=\kappa k^{p+1}$ であるとする。 $k\in I_0$ 、 $E(k)>0$ とし、 $\hat k:=\theta E(k)^{-1/(p+1)}k$ が $I_0$ に属するならば、 $E(\hat k)=\theta^{p+1}\kappa k^{p+1}/E(k)=\theta^{p+1}<1$ である。二半歩による刻み制御の推定値は、命題 1.1の仮定の下で、命題 1.1 (3)と命題 1.1 (1)により、解の上で $2^{-p}c(t)k^{p+1}$ との差が $k^{p+2}$ の定数倍以下である。尺度 $s_i$ は候補値 $v$ に依存し、推定値には $k^{p+2}$ の程度の項が加わるから、 $E$ は一般に $\kappa k^{p+1}$ の形でない。 $E=0$ は誤差が $0$ であることを意味しない（例 1.2 (3)）。

$y$ を解とし、受理された状態 $u=y_n$ が時刻 $t_n$ で $y(t_n)$ と異なるとする。 $(t_n,y_n)$ での推定値 $\hat e$ は $y(t_n)$ でなく $y_n$ から計算されるから、この推定値の評価に命題 1.1を用いる場合、同命題の $y$ に当たるのは、時刻 $t_n$ に値 $y_n$ をとり同命題の仮定を満たす解である。刻み $h_n$ の試行で受理された値 $y_{n+1}=\hat\Psi_{h_n}(t_n,y_n)$ について、 $(t_n,y(t_n),h_n)\in\hat D$ であるとき

y(t_{n+1})-y_{n+1}=\bigl(y(t_{n+1})-\hat\Psi_{h_n}(t_n,y(t_n))\bigr)+\bigl(\hat\Psi_{h_n}(t_n,y(t_n))-\hat\Psi_{h_n}(t_n,y_n)\bigr)

であり、受理の判定は、時刻 $t_n$ での誤差 $y(t_n)-y_n$ から生じる第二項を含まない。

命題 2.3.定義 2.1の刻み制御を実数の演算で実行すると、規則の繰返しは有限回で停止し、停止したとき $t=T$ または $h<h_{\min}$ である。受理された試行のうち $t+k<T$ を満たすものの刻み $k$ は $h_{\min}$ 以上であり、受理の回数は $\lceil(T-t_0)/h_{\min}\rceil$ 以下である。

証明. 演習とする（問題 5.1）。▨

系 2.4. $d\in\NN$ とし、 $\R^d$ にノルム $\|\cdot\|$ を固定する。 $\Omega\subseteq\R\times\R^d$ を開集合、 $f\colon\Omega\to\R^d$ を連続写像、 $\hat\Phi\colon\hat D\to\R^d$ を $f$ に対する増分関数、 $\hat\Psi_h$ をその一歩写像とする。 $t_0<T$ とし、 $J\supseteq[t_0,T]$ を開区間、 $y\colon J\to\R^d$ を、任意の $\tau\in J$ について $(\tau,y(\tau))\in\Omega$ と $y'(\tau)=f(\tau,y(\tau))$ を満たす $C^1$ 級写像とする。写像 $\hat e\colon\hat D\to\R^d$ と $\hat\Psi$ による定義 2.1の刻み制御が $t=T$ で停止したとし、受理された試行の後の時刻を $t_1<\dots<t_N=T$ 、状態を $y_1,\dots,y_N$ とし、 $y_0:=u_0$ 、 $h_n:=t_{n+1}-t_n$ と置く。 $\rho\in(0,\infty]$ と $\Lambda\ge0$ を取り、各 $n\in\{0,\dots,N-1\}$ と、 $\|u-y(t_n)\|\le\rho$ を満たす任意の $u\in\R^d$ について

(t_n,u,h_n)\in\hat D,\qquad\|\hat\Psi_{h_n}(t_n,u)-\hat\Psi_{h_n}(t_n,y(t_n))\|\le(1+\Lambda h_n)\|u-y(t_n)\|

が成り立つとする。 $d_n:=y(t_{n+1})-\hat\Psi_{h_n}(t_n,y(t_n))$ 、 $E_n:=\|y_n-y(t_n)\|$ と置き、 $\tau_0,\dots,\tau_{N-1}\ge0$ が $\|d_n\|\le\tau_n$ を満たすとする。

$e^{\Lambda(T-t_0)}\bigl(E_0+\sum_{j=0}^{N-1}\tau_j\bigr)\le\rho$ ならば、 $0\le n\le N$ について $E_n\le e^{\Lambda(t_n-t_0)}\bigl(E_0+\sum_{j=0}^{n-1}\tau_j\bigr)$ である。
$\varepsilon\ge0$ が $0\le j<N$ について $\tau_j\le\varepsilon h_j$ を満たし、 $\varphi_\Lambda$ を§E20.28 補題 1.2 (3)の関数として $e^{\Lambda(T-t_0)}E_0+\varepsilon\varphi_\Lambda(T-t_0)\le\rho$ ならば、 $0\le n\le N$ について $E_n\le e^{\Lambda(t_n-t_0)}E_0+\varepsilon\varphi_\Lambda(t_n-t_0)$ である。

証明. 棄却は $(t,u)$ を変えず、受理は $(t,u)$ を $(t+k,\hat\Psi_k(t,u))$ に替えるから、 $0\le n<N$ について $y_{n+1}=\hat\Psi_{h_n}(t_n,y_n)$ である。§E20.28 定理 1.3の $B_n$ を $r_j=0$ として

B_n=e^{\Lambda(t_n-t_0)}E_0+\sum_{j=0}^{n-1}e^{\Lambda(t_n-t_{j+1})}\|d_j\|

と置く。

(1)を示す。 $0\le t_n-t_{j+1}\le t_n-t_0$ と $\|d_j\|\le\tau_j$ により $B_n\le e^{\Lambda(t_n-t_0)}\bigl(E_0+\sum_{j<n}\tau_j\bigr)$ であり、特に $B_N\le\rho$ である。§E20.28 定理 1.3により $E_n\le B_n$ である。

(2)を示す。§E20.28 補題 1.2 (3)を $a_j:=\|d_j\|\le\varepsilon h_j$ 、 $\alpha:=\varepsilon$ として適用すると $B_n\le e^{\Lambda(t_n-t_0)}E_0+\varepsilon\varphi_\Lambda(t_n-t_0)$ であり、特に $B_N\le\rho$ である。§E20.28 定理 1.3により $E_n\le B_n$ である。▨

注意 2.5. $d=1$ 、 $\Omega=\R^2$ 、 $L>0$ 、 $\varepsilon>0$ 、 $f(t,u)=Lu$ とし、 $\Omega\times(0,\infty)$ 上の増分関数 $\hat\Phi(t,u,h):=\bigl((e^{Lh}-1)u+\varepsilon h\bigr)/h$ の一歩写像 $\hat\Psi_h(t,u)=e^{Lh}u+\varepsilon h$ を考える。解 $y(\tau)=e^{L(\tau-t_0)}$ と任意の格子について、系 2.4の $d_n$ は $d_n=-\varepsilon h_n$ であり、 $\sum_n\|d_n\|=\varepsilon(T-t_0)$ である。 $y_0=y(t_0)$ 、 $y_{n+1}=\hat\Psi_{h_n}(t_n,y_n)$ とすると、 $e_n:=y_n-y(t_n)$ は $e_0=0$ 、 $e_{n+1}=e^{Lh_n}e_n+\varepsilon h_n$ を満たし、

e_N=\varepsilon\sum_{n=0}^{N-1}e^{L(T-t_{n+1})}h_n

である。等間隔 $h_n=h$ では $e_N=\varepsilon h\bigl(e^{L(T-t_0)}-1\bigr)/(e^{Lh}-1)$ であり、 $e^{Lh}-1\le Lhe^{Lh}$ から $e_N\ge\varepsilon e^{-Lh}\bigl(e^{L(T-t_0)}-1\bigr)/L$ である。 $L(T-t_0)=10$ では $e_N$ は局所打切り誤差の和 $\varepsilon(T-t_0)$ の $e^{-Lh}(e^{10}-1)/10$ 倍以上である。部分区間ごとの誤差の和で積分全体の誤差を抑える§E20.21 補題 2.1と異なり、時刻 $t_{n+1}$ で加わった誤差は以後の歩の初期値の誤差となり、 $e^{L(T-t_{n+1})}$ 倍されて $e_N$ に現れる。

3 補間による軌道誤差

命題 3.1. $d\in\NN$ とし、 $u\in\R^d$ に対して $\|u\|_\infty:=\max_{1\le i\le d}|u_i|$ と置く。

$c\in\R$ 、 $h>0$ とし、 $y\colon[c,c+h]\to\R^d$ を $C^4$ 級写像、 $\varepsilon_y,\varepsilon_d\ge0$ とする。 $\tilde y_1,\tilde y_2,\tilde d_1,\tilde d_2\in\R^d$ が $\|\tilde y_1-y(c)\|_\infty,\ \|\tilde y_2-y(c+h)\|_\infty\le\varepsilon_y,\qquad\|\tilde d_1-y'(c)\|_\infty,\ \|\tilde d_2-y'(c+h)\|_\infty\le\varepsilon_d$ を満たすとする。各 $i$ について $\tilde H_i$ をデータ $(c,\tilde y_{1,i},\tilde d_{1,i})$ 、 $(c+h,\tilde y_{2,i},\tilde d_{2,i})$ の Hermite 補間多項式とし、 $\tilde H:=(\tilde H_1,\dots,\tilde H_d)$ 、 $M_4:=\max_i\max_{s\in[c,c+h]}|y_i^{(4)}(s)|$ と置く。このとき $\max_{x\in[c,c+h]}\|\tilde H(x)-y(x)\|_\infty\le\varepsilon_y+\frac{h\varepsilon_d}4+\frac{M_4h^4}{384}$ である。
$\Omega\subseteq\R\times\R^d$ 、 $f\colon\Omega\to\R^d$ 、 $t\in\R$ 、 $u,\tilde u,\tilde d\in\R^d$ 、 $L,\varepsilon_f\ge0$ とし、 $(t,u),(t,\tilde u)\in\Omega$ 、 $\|f(t,\tilde u)-f(t,u)\|_\infty\le L\|\tilde u-u\|_\infty$ 、 $\|\tilde d-f(t,\tilde u)\|_\infty\le\varepsilon_f$ とする。 $y$ が $t$ で微分可能で $y(t)=u$ 、 $y'(t)=f(t,u)$ を満たすならば、 $\|\tilde d-y'(t)\|_\infty\le\varepsilon_f+L\|\tilde u-u\|_\infty$ である。
$t_0<t_1<\dots<t_N$ と $\tilde y_n,\tilde d_n\in\R^d$ （ $0\le n\le N$ ）を取る。各 $n<N$ について $[t_n,t_{n+1}]$ 上の $\tilde y$ を、データ $(t_n,\tilde y_n,\tilde d_n)$ 、 $(t_{n+1},\tilde y_{n+1},\tilde d_{n+1})$ の成分ごとの Hermite 補間多項式で定める。このとき $\tilde y\colon[t_0,t_N]\to\R^d$ は矛盾なく定まる $C^1$ 級写像であり、 $\tilde y(t_n)=\tilde y_n$ 、 $\tilde y'(t_n)=\tilde d_n$ である。

証明.(1)を示す。各 $i$ について $H_i$ を $c,c+h$ における $y_i$ の Hermite 補間多項式とする。§E20.12 命題 7.4 (3)により $[c,c+h]$ 上で $|\tilde H_i-H_i|\le\varepsilon_y+h\varepsilon_d/4$ であり、§E20.12 命題 7.4 (4)を $y_i$ に適用すると $[c,c+h]$ 上で $|y_i-H_i|\le M_4h^4/384$ である。三角不等式により各 $i$ と $x\in[c,c+h]$ で $|\tilde H_i(x)-y_i(x)|$ は右辺以下であり、 $i$ について最大値をとると主張を得る。

(2)を示す。 $y'(t)=f(t,u)$ により $\tilde d-y'(t)=\bigl(\tilde d-f(t,\tilde u)\bigr)+\bigl(f(t,\tilde u)-f(t,u)\bigr)$ であり、二つの項のノルムはそれぞれ $\varepsilon_f$ と $L\|\tilde u-u\|_\infty$ 以下である。

(3)を示す。§E20.12 命題 7.4 (1)により、 $[t_n,t_{n+1}]$ 上の多項式は $t_n$ で値 $\tilde y_n$ と微分係数 $\tilde d_n$ をとり、 $t_{n+1}$ で値 $\tilde y_{n+1}$ と微分係数 $\tilde d_{n+1}$ をとる。したがって $0<n<N$ について、 $t_n$ を共有する二つの区間の多項式は $t_n$ で同じ値 $\tilde y_n$ をとり、 $\tilde y$ は矛盾なく定まって連続である。 $\tilde y$ の $t_n$ での左微分係数と右微分係数はともに $\tilde d_n$ であるから、 $\tilde y$ は $t_n$ で微分可能であり $\tilde y'(t_n)=\tilde d_n$ である。各区間上で $\tilde y'$ は多項式の導関数であって連続であり、 $t_n$ での片側極限はともに $\tilde d_n$ であるから、 $\tilde y'$ は $[t_0,t_N]$ 上で連続である。▨

例 3.2. $c\in\R$ 、 $h>0$ 、 $d=1$ とし、 $w(t):=(t-c)^2(t-c-h)^2$ 、 $f(t,u):=w'(t)=2(t-c)(t-c-h)(2t-2c-h)$ （ $(t,u)\in\R^2$ ）と置く。 $y:=w$ は $y'=f(t,y)$ 、 $y(c)=0$ を満たす。 $y(c)=y(c+h)=y'(c)=y'(c+h)=0$ であるから、§E20.12 命題 7.4 (1)により、この厳密な端点データの Hermite 補間多項式は $\tilde y=0$ である。 $\max_{[c,c+h]}|y-\tilde y|=w(c+\frac h2)=h^4/16$ であり、 $y^{(4)}=24$ から、命題 3.1 (1)の右辺は $\varepsilon_y=\varepsilon_d=0$ のとき $24h^4/384=h^4/16$ に等しい。

$\tilde y$ の欠陥 $\delta(s)=\tilde y'(s)-f(s,\tilde y(s))=-w'(s)$ は $s=c,c+\frac h2,c+h$ で $0$ であるが、軌道誤差は $h^4/16$ である。§E20.28 定理 4.1を $[t_0,T]=[c,c+h]$ 、 $S=[c,c+h]\times\R$ 、 $\ell=0$ として適用すると $|\tilde y(t)-y(t)|\le\int_c^t|w'(s)|\,ds$ である。 $w'$ は $[c,c+\frac h2]$ で非負、 $[c+\frac h2,c+h]$ で非正であるから、右辺は $t=c+\frac h2$ で $w(c+\frac h2)=h^4/16$ となって誤差に等しく、 $t=c+h$ で $h^4/8$ となり、誤差 $0$ より大きい。

4 イベント時刻

定義 4.1. $d\in\NN$ とし、 $I\subseteq\R$ を区間、 $U\subseteq\R^d$ 、 $y\colon I\to U$ を写像、 $g\colon I\times U\to\R$ を関数とする。 $g$ をイベント関数 (event function) といい、 $G(t):=g(t,y(t))$ （ $t\in I$ ）の零点を $y$ のイベント時刻 (event time) という。

補題 4.2. $d\in\NN$ とし、 $\R^d$ にノルム $\|\cdot\|$ を固定する。 $t\in\R$ 、 $U\subseteq\R^d$ とし、 $g$ を $\{t\}\times U$ を含む集合上の実数値関数とする。 $K\ge0$ が任意の $u,v\in U$ について $|g(t,u)-g(t,v)|\le K\|u-v\|$ を満たすとする。 $u,\tilde u\in U$ 、 $\hat g\in\R$ 、 $\varepsilon_y,\varepsilon_g\ge0$ が $\|\tilde u-u\|\le\varepsilon_y$ 、 $|\hat g-g(t,\tilde u)|\le\varepsilon_g$ を満たすならば、 $|\hat g-g(t,u)|\le K\varepsilon_y+\varepsilon_g$ である。

証明. $\hat g-g(t,u)=\bigl(\hat g-g(t,\tilde u)\bigr)+\bigl(g(t,\tilde u)-g(t,u)\bigr)$ であり、二つの項の絶対値はそれぞれ $\varepsilon_g$ と $K\|\tilde u-u\|\le K\varepsilon_y$ 以下である。▨

命題 4.3. $d\in\NN$ とし、 $\R^d$ にノルム $\|\cdot\|$ を固定する。 $I\subseteq\R$ を区間、 $U\subseteq\R^d$ 、 $y\colon I\to U$ と $\tilde y\colon I\to\R^d$ を写像、 $g\colon I\times U\to\R$ を関数とし、 $G(t):=g(t,y(t))$ が $I$ 上で連続かつ $I$ の内部で微分可能であるとする。 $\tau,\hat\tau\in I$ が $G(\tau)=0$ を満たし、 $m>0$ が $\min\{\tau,\hat\tau\}<s<\max\{\tau,\hat\tau\}$ を満たす任意の $s$ について $|G'(s)|\ge m$ を満たすとする。 $K,\varepsilon_y,\varepsilon_g,r\ge0$ と $\hat g\in\R$ が、 $\tilde y(\hat\tau)\in U$ 、任意の $u,v\in U$ についての $|g(\hat\tau,u)-g(\hat\tau,v)|\le K\|u-v\|$ 、 $\|\tilde y(\hat\tau)-y(\hat\tau)\|\le\varepsilon_y$ 、 $|\hat g-g(\hat\tau,\tilde y(\hat\tau))|\le\varepsilon_g$ 、 $|\hat g|\le r$ を満たすならば

|\hat\tau-\tau|\le\frac{K\varepsilon_y+\varepsilon_g+r}m

である。

命題 4.4. $a<b$ とし、 $G\colon[a,b]\to\R$ を連続関数、 $\eta\ge0$ とする。 $\hat g_a,\hat g_b\in\R$ が $|\hat g_a-G(a)|\le\eta$ 、 $|\hat g_b-G(b)|\le\eta$ を満たすとする。

$\hat g_a\hat g_b<0$ かつ $\min\{|\hat g_a|,|\hat g_b|\}>\eta$ ならば、 $G(a)G(b)<0$ であり、 $G$ は $(a,b)$ に零点をもつ。
$G$ が $(a,b)$ で微分可能であり、 $m>0$ が任意の $s\in(a,b)$ について $|G'(s)|\ge m$ を満たすならば、 $G$ の $[a,b]$ における零点は高々一つである。
$\Lambda_G\ge0$ が任意の $s,t\in[a,b]$ について $|G(s)-G(t)|\le\Lambda_G|s-t|$ を満たし、 $|\hat g_a|+|\hat g_b|>\Lambda_G(b-a)+2\eta$ ならば、 $G$ は $[a,b]$ に零点をもたない。

証明.(1)を示す。 $|\hat g_a-G(a)|\le\eta<|\hat g_a|$ であるから $G(a)$ は $\hat g_a$ と同符号であり、同様に $G(b)$ は $\hat g_b$ と同符号である。したがって $G(a)G(b)<0$ である。 $G(a)<0<G(b)$ ならば $G$ に、 $G(a)>0>G(b)$ ならば $-G$ に§D1.12 定理 1.1を適用して、 $(a,b)$ の零点を得る。

(2)を示す。 $\tau_1,\tau_2\in[a,b]$ を $G$ の零点とする。 $\tau_1$ と $\tau_2$ の間の点は $(a,b)$ に属するから、§E20.2 命題 6.1を区間 $[a,b]$ 、零点 $\tau_1$ 、近似点 $\tau_2$ に適用して $|\tau_2-\tau_1|\le|G(\tau_2)|/m=0$ を得る。

(3)を示す。 $s\in[a,b]$ が $G(s)=0$ を満たすとすると、 $|G(a)|=|G(a)-G(s)|\le\Lambda_G(s-a)$ 、 $|G(b)|\le\Lambda_G(b-s)$ であるから

|\hat g_a|+|\hat g_b|\le|G(a)|+|G(b)|+2\eta\le\Lambda_G(b-a)+2\eta

である。仮定はこの不等式の否定であるから、 $G$ は $[a,b]$ に零点をもたない。▨

例 4.5.

$[a,b]=[0,2]$ 、 $G(t)=(t-1)^2$ とする。 $G$ の零点は $1$ だけであり、 $G(0)=G(2)=1$ は同符号である。 $G'(1)=0$ であり、 $|G'(s)|=2|s-1|$ であるから、 $1$ を端点とする開区間上で $|G'|$ の正の下界は存在しない。 $\eta_0\in(0,1)$ に対して $G-\eta_0$ は $G$ との差の絶対値が $\eta_0$ であり、零点 $1\pm\sqrt{\eta_0}$ をもつ。 $\hat\tau:=1+\sqrt{\eta_0}$ は $(G-\eta_0)(\hat\tau)=0$ 、 $G(\hat\tau)=\eta_0$ を満たすが、 $G$ の零点 $1$ との距離は $\sqrt{\eta_0}$ である。 $\eta_0=10^{-8}$ では値の差 $10^{-8}$ に対して時刻の差は $10^{-4}$ である。 $G+\eta_0$ は零点をもたない。
$[a,b]=[0,1]$ 、 $G(t)=(t-\frac15)(t-\frac12)(t-\frac45)$ とする。 $G(0)=-\frac2{25}<0<\frac2{25}=G(1)$ であり、 $G$ は $(0,1)$ に三つの零点をもつ。§E20.4 定義 2.2の二分法は $m_0=\frac12$ で $G(m_0)=0$ となり、§E20.4 定義 2.2 (1)により $\frac12$ を出力して停止する。出力 $\frac12$ は $(0,1)$ の最小の零点 $\frac15$ と異なる。命題 4.4 (2)により、 $(0,1)$ 上で $|G'|\ge m$ を満たす $m>0$ は存在しない。
$[a,b]=[0,1]$ 、 $G(t)=(t-\frac3{10})(t-\frac7{10})$ とする。 $G(0)=G(1)=\frac{21}{100}>0$ であるが、 $G$ は $(0,1)$ に零点 $\frac3{10}$ 、 $\frac7{10}$ をもつ。 $G(s)-G(t)=(s-t)(s+t-1)$ であり、 $s,t\in[0,1]$ では $|s+t-1|\le1$ であるから、 $\Lambda_G=1$ が $[0,1]$ とその部分区間で命題 4.4 (3)の Lipschitz 条件を満たす。 $|G(0)|+|G(1)|=\frac{21}{50}<1$ であるから、 $[0,1]$ では $\eta=0$ の命題 4.4 (3)の仮定は成り立たない。 $[0,\frac15]$ では $|G(0)|+|G(\frac15)|=\frac{21}{100}+\frac1{20}=\frac{26}{100}>\frac15$ であるから、命題 4.4 (3)により $G$ は $[0,\frac15]$ に零点をもたない。

例 4.6. $d=2$ 、 $\Omega=\R\times\R^2$ とし、状態を $(x,v)$ と書いて $f(t,(x,v)):=(v,-1+v^2)$ と置く。初期値 $(x,v)(0)=(1,0)$ の解は $x(t)=1-\log\cosh t$ 、 $v(t)=-\tanh t$ である。 $v<0$ では $-1+v^2=-1-v|v|$ であり、 $x$ と $v$ は、重力加速度と比例定数を $1$ とした速さの二乗に比例する抵抗を受けて落下する物体の高さと速度である。イベント関数 $g(t,(x,v)):=x$ について $G(t)=x(t)$ 、 $G'(t)=-\tanh t$ であり、 $G(0)=1$ と $t>0$ での $G'(t)<0$ により、 $G$ の $[0,\infty)$ の零点は

\tau=\operatorname{arcosh}e=\log\bigl(e+\sqrt{e^2-1}\bigr)=1.657454454153\ldots

ただ一つである。

Heun 法の二半歩による刻み制御を、 $p=2$ 、 $t_0=0$ 、 $T=2$ 、 $u_0=(1,0)$ 、 $a_1=a_2=r=10^{-3}$ 、 $\theta=0.8$ 、 $\alpha_{\min}=0.2$ 、 $\alpha_{\max}=2$ 、 $h_{\min}=10^{-6}$ 、 $h_{\mathrm{init}}=0.5$ として CPython の float で実行した。最初の試行は $k=0.5$ 、 $E=4.3149\ldots$ で棄却されて $h=0.245698\ldots$ となり、以後の8回の試行はすべて受理されて $t=T$ で停止した。8回目の受理の刻みは、終点への切詰めによる $T-t_7$ である。受理後の状態を $(x_n,v_n)$ 、受理の判定に用いた尺度を $s_{n,1},s_{n,2}$ とし、60桁の十進演算で求めた解の値と比べて丸めると次のとおりである。最後の列は $\max\{|x_n-x(t_n)|/s_{n,1},\ |v_n-v(t_n)|/s_{n,2}\}$ である。

$n$	$t_n$	$t_n-t_{n-1}$	$E$	$\lvert x_n-x(t_n)\rvert$	$\lvert v_n-v(t_n)\rvert$	尺度で割った誤差
$1$	$0.245699$	$0.245699$	$0.5243$	$7.283\times10^{-5}$	$6.380\times10^{-4}$	$0.5144$
$2$	$0.489467$	$0.243768$	$0.5676$	$2.312\times10^{-4}$	$1.269\times10^{-3}$	$0.8737$
$3$	$0.724998$	$0.235531$	$0.5278$	$3.327\times10^{-4}$	$1.704\times10^{-3}$	$1.0530$
$4$	$0.958158$	$0.233160$	$0.4626$	$3.396\times10^{-4}$	$1.900\times10^{-3}$	$1.0912$
$5$	$1.199338$	$0.241180$	$0.4137$	$2.753\times10^{-4}$	$1.907\times10^{-3}$	$1.0411$
$6$	$1.458277$	$0.258939$	$0.3754$	$1.634\times10^{-4}$	$1.775\times10^{-3}$	$0.9362$
$7$	$1.745444$	$0.287168$	$0.3419$	$2.251\times10^{-5}$	$1.550\times10^{-3}$	$0.7990$
$8$	$2.000000$	$0.254556$	$0.1461$	$1.770\times10^{-4}$	$1.178\times10^{-3}$	$0.6001$

すべての受理で $E\le1$ であるが、 $n=3,4,5$ では最後の列が $1$ を超え、 $v_n$ の真の誤差は受理の判定に用いた尺度を超えている。

$x_6=0.182000\ldots>0>x_7=-0.082338\ldots$ である。 $\eta:=1.7\times10^{-4}$ と置くと、 $|x_6-x(t_6)|\le\eta$ 、 $|x_7-x(t_7)|\le\eta$ 、 $\min\{|x_6|,|x_7|\}>\eta$ であるから、命題 4.4 (1)により $G$ は $(t_6,t_7)$ に零点をもつ。刻みの端点だけで時刻を $t_7$ とすると、誤差は $t_7-\tau=0.087989\ldots$ である。

$[t_6,t_7]$ 上で $x$ をデータ $(t_6,x_6,v_6)$ 、 $(t_7,x_7,v_7)$ の三次 Hermite 補間多項式 $\tilde x$ で近似する。 $x'=v$ であるから、微分データ $v_n$ の誤差は $|v_n-x'(t_n)|=|v_n-v(t_n)|$ である。 $\sigma:=\operatorname{sech}^2t$ と置くと $x^{(4)}=2\sigma(3\sigma-2)$ であり、 $[t_6,t_7]$ 上で $\sigma<\frac13$ であるから $|x^{(4)}|=2\sigma(2-3\sigma)$ は $\sigma$ について増加し、 $t$ について減少する。したがって $M_4=|x^{(4)}(t_6)|=0.5515\ldots$ である。命題 3.1 (1)を $d=1$ 、 $\varepsilon_y=1.633\ldots\times10^{-4}$ 、 $\varepsilon_d=1.774\ldots\times10^{-3}$ 、 $h=0.287167\ldots$ として適用すると、 $[t_6,t_7]$ 上で

|\tilde x-x|\le\varepsilon_y+\frac{h\varepsilon_d}4+\frac{M_4h^4}{384}=3.005\ldots\times10^{-4}

である。

節点の値は二進浮動小数点数であり有理数であるから、 $\tilde x$ の値を有理数の演算で厳密に計算し、§E20.4 定義 2.2の二分法を $[t_6,t_7]$ 上の $\tilde x$ に60段適用して $\hat\tau:=m_{60}=1.657367744776\ldots$ を得た。 $\tilde x(\hat\tau)$ は厳密に計算した値であり、 $|\tilde x(\hat\tau)|<8.9\times10^{-22}$ 、 $|\tilde x(\hat\tau)-x(\hat\tau)|=8.0628\ldots\times10^{-5}$ である。 $|G'|=\tanh$ は $t>0$ で増加し、 $\hat\tau<\tau$ であるから、 $m:=\tanh\hat\tau=0.929861\ldots$ は $\hat\tau$ と $\tau$ の間で $|G'|\ge m$ を満たす。命題 4.3を $d=1$ 、 $I=[t_6,t_7]$ 、 $U=\R$ 、 $y=x$ 、 $\tilde y=\tilde x$ 、 $g(t,x)=x$ 、 $K=1$ 、 $\hat g=\tilde x(\hat\tau)$ 、 $\varepsilon_g=0$ 、 $\varepsilon_y=|\tilde x(\hat\tau)-x(\hat\tau)|$ 、 $r=|\tilde x(\hat\tau)|$ として適用すると

|\hat\tau-\tau|\le\frac{|\tilde x(\hat\tau)-x(\hat\tau)|+|\tilde x(\hat\tau)|}m=8.67099\ldots\times10^{-5}

であり、実際の誤差は $|\hat\tau-\tau|=8.67093\ldots\times10^{-5}$ である。 $[t_6,t_7]$ 上で $|G'|\ge\tanh t_6=0.8973\ldots$ であるから、同じ命題を $m=\tanh t_6$ として用いると、 $\varepsilon_t>0$ と $[t_6,t_7]$ に属する任意の $\hat\tau$ について、 $|\tilde x(\hat\tau)-x(\hat\tau)|+|\tilde x(\hat\tau)|\le0.89\,\varepsilon_t$ ならば $|\hat\tau-\tau|\le\varepsilon_t$ である。例 4.5 (1)では $|G'|$ の正の下界が存在しないから、この形の十分条件は得られない。上の $\eta$ 、Hermite 補間の評価に用いた $\varepsilon_y,\varepsilon_d$ 、 $|\tilde x(\hat\tau)-x(\hat\tau)|$ は、刻み制御の計算からでなく解の値から得たものである。

5 演習

問題 5.1.命題 2.3の証明を完成させよ。

解答.

$\kappa:=\max\{\alpha_{\min},\theta\}$ と置くと $\kappa<1$ である。棄却では $h$ が $\kappa h$ 以下の値に替わる。実際、 $(t,u,k)\notin\hat D$ の場合の新しい値は $\alpha_{\min}k\le\kappa h$ である。 $E>1$ の場合は $\theta E^{-1/(p+1)}<\theta\le\kappa$ と $\alpha_{\min}\le\kappa$ により $q(E)\le\max\{\alpha_{\min},\theta E^{-1/(p+1)}\}\le\kappa$ であり、新しい値は $q(E)k\le\kappa h$ である。受理では $q(E)\le\alpha_{\max}$ により新しい値は $q(E)k\le\alpha_{\max}h$ である。試行ごとに $k\le T-t$ であるから、状態の $t$ はつねに $T$ 以下であり、受理でだけ増加する。

停止は $t=T$ の規則か $h<h_{\min}$ の規則によってだけ起こるから、停止したとき $t=T$ または $h<h_{\min}$ である。

受理された試行が $t+k<T$ を満たすとする。 $k=\min\{h,T-t\}<T-t$ であるから $k=h$ であり、試行が行われたことから $h\ge h_{\min}$ である。したがって $k\ge h_{\min}$ である。このような受理の回数を $A'$ とすると、これらの受理の後の $t$ はすべて $T$ 未満であり、最後のものの後の $t$ は $t_0+A'h_{\min}$ 以上であるから、 $A'<(T-t_0)/h_{\min}$ 、すなわち $A'\le\lceil(T-t_0)/h_{\min}\rceil-1$ である。 $t+k=T$ を満たす受理の後は $t=T$ であり、次の繰返しで停止するから、そのような受理は高々一回である。したがって受理の回数は $A_{\max}:=\lceil(T-t_0)/h_{\min}\rceil$ 以下である。

$h$ は受理で高々 $\alpha_{\max}$ 倍になり、棄却で減少するから、つねに $h\le H^*:=h_{\mathrm{init}}\alpha_{\max}^{A_{\max}}$ である。 $\kappa^JH^*<h_{\min}$ を満たす $J\in\NN$ を取る。連続する $J$ 回の棄却の後では $h\le\kappa^JH^*<h_{\min}$ であり、棄却は $t$ を変えないから $t<T$ のままであって、次の繰返しで停止する。したがって、停止しない繰返しは、高々 $A_{\max}$ 回の受理と、それらによって分けられた高々 $A_{\max}+1$ 個の区切りの各々での高々 $J$ 回の棄却からなり、その回数は $A_{\max}+(A_{\max}+1)J$ 以下である。これで命題 2.3は示された。▨

§E20.30適応的時間刻みとイベント検出

1 一歩と二半歩

2 刻み制御

3 補間による軌道誤差

4 イベント時刻

5 演習

前提記事