広義Runge–Kutta法の数値安定性

数値解法の安定性とは

(大雑把に言えば)どのくらい刻み幅が大きくなると計算が明らかにおかしくなるか、ということである。

A-stabilityとは

数値解法の安定性を評価する指標の一つである。

次の常微分方程式を考える。
$\displaystyle\frac{\mathrm{d} x(t)}{\mathrm{d} t}=kx(t) \ (k \in \mathbb{C}) \cdots (1)$
これを $x(0)=1$ として解析的に解けば
$\displaystyle x(t)=e^{kt}$

このとき、 $\mathrm{Re}(k)<0$ ならば
$\displaystyle\lim_{t \to \infty} x(t) = 0 \cdots (2)$
である。

ある数値解法で(1)式を解いたとき、刻み幅 $h$ にかかわらず式(2)を満たすならば、その解法はA-stableであるという。

安定領域とは

上の係数 $k$ と、刻み幅 $h$ の積 $hk$ が $\lim_{t \to \infty} x(t) = 0$ を満たす $hk$ の領域である。
広義Runge–Kutta法においてはこの領域は容易に求められる。

広義Runge–Kutta法において次の式が成り立つ。
$\displaystyle x_{n+1}=\phi(hk) \cdot x_{n}$
ここで、 $x_{n}=x(hn)+x_0$ である。 $x_0=x(t_0)$ であり、 $t_0$ は $t$ の初期値である。
また、 $\phi$ は各解法で固有な関数で、安定関数とよばれる。
安定領域は、 $\{z \in \mathbb{C}\ |\ |\phi(z)| < 1\}$ なる領域である。

ある解法がA-stableであるには、安定領域が $\{z \in \mathbb{C}\ |\ \mathbb{Re}(z) < 0\}$ なる領域を含んでいればよい。

各解法の安定領域

今までこのブログで書いてきた各解法(すべて広義Runge–Kutta法である)について、その安定領域を調べる。

陽的Euler法

陽的オイラー法の式は
$\displaystyle x_{n+1} = x_n + f(t_n,x_n)h$

これに、 $f(t_n,x_n) = kx_n$ を代入して、
$\begin{eqnarray} x_{n+1} &=& x_n + kx_n h &=& (1+kh)x_n \end{eqnarray}$

よって、安定関数は
$\displaystyle\phi(z)=1+z$
ゆえに、安定領域は
$\displaystyle|1+z|<1$
(以後の安定関数の計算も全て同様の方法である)

f:id:equal_l2:20150922093556p:plain

(縦軸:虚軸( $\mathrm{Im}$ ) 横軸:実軸( $\mathrm{Re}$ ) 以後の画像も同様)

陰的Euler法

陰的オイラー法の式は
$\displaystyle x_{n+1} = x_n + f(t_{n+1},x_{n+1})h$

安定関数は
$\displaystyle\phi(z)=\frac{1}{1-z}$
ゆえに、安定領域は
$\displaystyle\left|\frac{1}{1-z}\right|<1$
すなわち
$\displaystyle|1-z|>1$

f:id:equal_l2:20150922093612p:plain

この解法はA-stableである。

台形法

台形法の式は
$\displaystyle x_{n+1} = x_n + \frac{h}{2}\left( f(t_n,x_n) + f(t_{n+1},x_{n+1}) \right)$
安定関数は
$\displaystyle\phi(z)=\frac{1+\frac{z}{2}}{1-\frac{z}{2}}$
ゆえに、安定領域は
$\displaystyle\left|\frac{1+\frac{z}{2}}{1-\frac{z}{2}}\right|<1$
変形して、
$\displaystyle\mathrm{Re}(z)<0$

f:id:equal_l2:20150922093624p:plain

この解法はA-stableである。

陽的中点法

陽的中点法の式は
$\displaystyle x_{n+1} = x_{n} + f(t_n+\frac12 h,x_n + \frac12 h f(t_n,x_n))h$
安定関数は
$\displaystyle\phi(z)=\frac12 z^2 + z + 1$

f:id:equal_l2:20150922093634p:plain

陰的中点法

陰的中点法の式は
$\displaystyle x_{n+1} = x_{n} + f(t_n+\frac12 h,\frac12 \left( x_n + x_{n+1} \right))h$
安定関数は
$\displaystyle\phi(z)=\frac{1+\frac{z}{2}}{1-\frac{z}{2}}$
これは台形法と同じ安定関数なので、安定領域は
$\displaystyle\mathrm{Re}(z)<0$
(図省略)

陽的化した陰的Euler法

陽的化した陰的オイラー法の式は
$\displaystyle x_{n+1} = x_n + f(t_{n+1}, x_n + f(t_n,x_n)h)h$
安定関数は
$\displaystyle \phi(z)= z^2 + z + 1$

f:id:equal_l2:20170302060531p:plain

陽的化した台形法 (Heun法)

陽的化した台形法の式は
$\displaystyle x_{n+1} = x_n + \frac{h}{2}\left( f(t_n,x_n) + f(t_{n+1},x_n + f(t_n,x_n)h) \right)$
安定関数は
$\displaystyle\phi(z)=\frac12 z^2 + z + 1$
これは陽的中点法と同じ安定関数である。