#はじめに
この記事は数値計算 Advent Calendar 2018の１日目の記事です.

ルンゲ-クッタ法によって求めた近似解の安定性を調べたいと思います.

オイラー法 (1次のルンゲ-クッタ法)

まずは1次のルンゲ-クッタ法である前進オイラー法(Forward Euler, FE)と後退オイラー法(Backward Euler, BE)の安定性を考えてみたいと思います. 常微分方程式の初期値問題

\begin{align}
\frac{dy}{dx}  &= f(x,y) \\
y(x_0) &= y_0
\end{align}

に対して, 刻み幅を$h$とすると前進オイラー法は,

y_1 = y_0 + hf(x_0,y_0)

となります.
また後退オイラー法は

y_1 = y_0 + hf(x_1,y_1)

です.
#テスト方程式
テスト方程式と呼ばれる以下の常微分方程式に対して, 計算スキームを適用したときの数値解の安定性を考えます.

\begin{align}
\frac{dy}{dx}  &= \lambda y \ \ (\lambda \in \mathbb{C})\\
y(0) &= 1
\end{align}

厳密解は

y(x)  = e^{\lambda x}

です. よって$\lambda$の実部が負のとき,

\text{Re} (\lambda) <0

$x\to \infty$で厳密解は$0$に収束します.

\lim_{x \to \infty} y(x)= 0

#前進オイラー法
テスト方程式に対して前進オイラー法を適用します.

\begin{align}
y_1 &= y_0 + h (\lambda y_0) \\
&= (1 + h\lambda) y_0 
\end{align}

ここで

z = h\lambda \ \ ( z \in \mathbb{C})

とおくと

y_1 = (1 + z) y_0

となります. このときの

R(z) = 1 + z

を安定性関数といいます. また$R(z)$の絶対値が1以下の集合

S = \{ z \in \mathbb{C} : \left|R(z)\right| \leq 1 \}

を安定性領域といいます.
前進オイラー法の安定性領域をプロットすると下記のようになります.

#後退オイラー法
テスト方程式に対して後退オイラー法を適用します.

y_1 = y_0 + h(\lambda y_1)

これを解くと,

y_1 = \frac{1}{1-h\lambda}y_0

となります. よって後退オイラー法の安定性関数は,

R(z) = \frac{1}{1-z}

です.
後退オイラー法の安定性領域をプロットすると下記のようになります.

#計算例
それでは刻み幅$h$を$0.1$に固定して, $\lambda$の値を変えたときに, 厳密解と数値解の振る舞いを見ていきたいと思います.

厳密解 : 収束, 前進オイラー法 : 収束, 後退オイラー法 : 収束
厳密解 : 収束, 前進オイラー法 : 発散, 後退オイラー法 : 収束
厳密解 : 発散, 前進オイラー法 : 発散, 後退オイラー法 : 収束
厳密解 : 発散, 前進オイラー法 : 発散, 後退オイラー法 : 発散

の4パターンをみてみます.
###1. 厳密解 : 収束, 前進オイラー法 : 収束, 後退オイラー法 : 収束
$\lambda = -0.3 + 2i$のときの計算結果は下記になります.

左上が前進オイラー法の安定性領域, バツ印が$z = h\lambda$の値です. 右上が後退オイラー法の安定性領域です. また左下が前進オイラー法の計算結果$y_{{\rm FE}}$(青線)と厳密解$y_{{\rm exact}}$(赤点線), 右下が後退オイラー法の計算結果$y_{{\rm BE}}$(青線)です.

左上と右上のバツ印より前進オイラー法, 後退オイラー法ともに, バツ印が安定性領域内にあります. また$\text{Re} (\lambda) <0$より厳密解は0に収束します(左下, 右下の赤点線). 同様に, 前進オイラー法, 後退オイラー法ともに, 誤差はありますが0に収束します(左下, 右下の青線).

###2. 厳密解 : 発散, 前進オイラー法 : 発散, 後退オイラー法 : 収束
$\lambda = -0.1+ 2i$のときの計算結果は下記になります.

左上の図のバツ印が前進オイラー法の安定性領域の外にでました. すると$\text{Re} (\lambda) <0$より厳密解は0に収束するにも関わらず(左下, 赤点線), 前進オイラー法は, 発散してしまいます(左下, 青線).

###3. 厳密解 : 発散, 前進オイラー法 : 発散, 後退オイラー法 : 収束
$\lambda = 0.1+ 2i$のときの計算結果は下記になります.

今度は$\text{Re} (\lambda) >0$より厳密解は発散するにも関わらず(右下, 赤点線), 後退オイラー法は, 収束してしまいます(右下, 青線).

###4. 厳密解 : 発散, 前進オイラー法 : 発散, 後退オイラー法 : 発散
$\lambda = 0.3+ 2i$のときの計算結果は下記になります.

右上の図のバツ印も後退オイラー法の安定性領域の外にでましたので, すべてのケースで発散するようになりました.

このようにパラメータによっては, 厳密解が収束するにも関わらず, 数値解が発散したり, 逆に厳密解が発散するにも関わらず, 数値解が収束してしまうことがわかりました. 前進オイラー法は, 陽的ルンゲ-クッタ法と呼ばれる手法の中で, 一番精度の低いものであり, 後退オイラー法は, 陰的ルンゲ-クッタ法と呼ばれる手法の中で, 一番精度の低いものです. それではより高次の手法ではどうなるのか見ていきたいと思います.

陽的ルンゲ-クッタ法

常微分方程式に対する$s$段の陽的ルンゲ-クッタ法は以下の式で表されます。

\begin{align}
k_i     &= f\left(x_0 + hc_i, y_0 + h \sum_{j=1}^{i-1} a_{ij} k_j \right) \\
y_{1} &= y_0 + h \sum_{i=1}^s b_ik_i 
\end{align}

ここで$a_{ij},b_i,c_i$はパラメータです.

2段2次のルンゲ-クッタ法(ホイン法)

2段2次のルンゲ-クッタ法(ホイン法)は以下の式となります.

\begin{align}
k_1     &= f(x_0, y_0) \\
k_2     &= f(x_0 + h, y_0 + hk_1 ) \\
y_{1}   &= y_0 + h \frac{k_1+k_2}{2} 
\end{align}

これをテスト方程式に適用すると,

\begin{align}
y_{1}   &= y_0 + (h\lambda) y_0 + \frac{(h\lambda)^2}{2}  y_0
\end{align}

となります. これは厳密解$y(x)=y_0e^{\lambda x}$の$x_0$周りのテイラー展開

\begin{align}
y(x_0+h)&= y_0 + (h\lambda) y_0 + \frac{(h\lambda)^2}{2}  y_0 + \cdots
\end{align}

と2次まで一致します. また$z=h\lambda$とおくと,

\begin{align}
y_{1}   &= \left(1 + z  + \frac{z^2}{2}  \right) y_0
\end{align}

となり安定性関数は,

\begin{align}
R_2(z) = 1 + z  + \frac{z^2}{2} 
\end{align}

です.

s段p次のルンゲ-クッタ法

段数$s$と次数$p$が等しい陽的ルンゲ-クッタ法の安定性関数は,

\begin{align}
R(z) = 1 + z  + \frac{z^2}{2!}  + \cdots + \frac{z^s}{s!} 
\end{align}

となります. これは$z=\lambda x$とした厳密解$e^z$のテイラー展開と$p$次まで一致します.

\begin{align}
e^z = 1 + z  + \frac{z^2}{2!}  + \cdots  + \frac{z^p}{p!} +\cdots
\end{align}

よって3段3次の安定性関数は,

\begin{align}
R_3(z) = 1 + z  + \frac{z^2}{2}  + \frac{z^3}{6} 
\end{align}

また4段4次(古典的ルンゲ-クッタ法)の安定性関数は,

\begin{align}
R_4(z) = 1 + z  + \frac{z^2}{2}  + \frac{z^3}{6} + \frac{z^4}{24} 
\end{align}

です. まとめると,

\begin{align}
R_1(z) &= 1 + z \\
R_2(z) &= 1 + z  + \frac{z^2}{2}  \\
R_3(z) &= 1 + z  + \frac{z^2}{2}  + \frac{z^3}{6}  \\
R_4(z) &= 1 + z  + \frac{z^2}{2}  + \frac{z^3}{6} + \frac{z^4}{24} 
\end{align}

となります(下添字は次数). 4次までの陽的ルンゲ-クッタ法の安定性領域をプロットすると下記のようになります.

上図より次数が上がる程に安定性領域が広がっていくことがわかります. しかし右下の4次であっても, 厳密解が収束する複素平面の左半平面で, 不安定な領域があることがみてとれます.

陰的ルンゲ-クッタ法

常微分方程式に対する$s$段の陰的ルンゲ-クッタ法は, 以下の式で表されます。

\begin{align}
k_i     &= f\left(x_0 + hc_i, y_0 + h \sum_{j=1}^{s} a_{ij} k_j \right) \\
y_{1} &= y_0 + h \sum_{i=1}^s b_ik_i 
\end{align}

陰的ルンゲ-クッタ法は, 常微分方程式の右辺の$f(x,y)$の形によっては, $k_i$を求める非線形連立方程式を解かないといけません. しかしテスト方程式は線形なので, 連立一次方程式を解き, さらにクラメールの公式を使うことで, 下記の安定性関数を求めることができます(参考文献[1], P42, 命題3.2).

$$
R(z) = \frac{\det(I-zA+z\mathbb{1}b^T)}{\det(I-zA)}\tag{1}
$$
ここで, $I$は単位行列で, $b^T=(b_1,\dots,b_s)$, $A=(a_{ij})_{i,j=1}^s$, $\mathbb{1} = (1,\dots,1)^T$です.

陰的ルンゲ-クッタ法の安定性領域をプロットすると下記のようになります.

左上は1段1次の後退オイラー法, 右上1段2次の陰的中点法, 左下は2段3次のラダウⅠA法, 右下は2段4次のガウス-ルジャンドル法です. 高次の手法でも, 後退オイラー法と同じく, 複素平面の左半平面に安定性領域が広がっています. また先程みたように, 後退オイラー法では, 厳密解が発散する右半平面にも安定性領域が広がっていますが(左上), 右側の２つは右半平面に安定性領域が広がっていないことが確認できます.

陰的ルンゲ-クッタ法の安定性関数

先程求めたように, 1段1次の後退オイラー法の安定性関数は,

R_{{\rm BE}}(z) = \frac{1}{1-z}

です. また式(1)より1段2次の陰的中点法(Implicit Midpoint, IM)の安定性関数は,

R_{{\rm IM}}(z) = \frac{1+\frac{z}{2}}{1-\frac{z}{2}}

となります. この安定性関数をそれぞれ原点周りでテイラー展開すると,

\begin{align}
R_{{\rm BE}}(z) &= 1 + z + z^2 + \cdots \\
R_{{\rm IM}}(z) &= 1 + z + \frac{z^2}{2} + \frac{z^3}{4} + \cdots
\end{align}

となります. これは厳密解$e^z$の原点周りのテイラー展開と,後退オイラー法では$1$次まで, 陰的中点法では$2$次まで一致しています. すなわち数値計算手法の次数$p$まで一致します.

式(1)の形からもわかるように, 陰的ルンゲ-クッタ法では, 安定性関数が有理関数(多項式わる多項式)になります. 実は, 上記でみた安定性関数は, 厳密解$e^z$のパデ近似というものになっています. パデ近似は, 関数を有理関数で近似する手法であり, 分母と分子の多項式の次数を足した次数まで, 元の関数のテイラー展開と近似する有理関数のテイラー展開が一致するように, 近似関数を構成します.

$\lambda = -1$のときの厳密解$y(x)=e^{-x}$と, テイラー展開を2次で打ち切った近似関数$y_{{\rm Taylor}}(x)$と, 分母1次・分子1次のパデ近似関数$y_{{\rm Pade}}(x)$のプロットが下記になります. パデ近似関数$y_{{\rm Pade}}(x)$は, 陰的中点法の安定性関数$R_{{\rm IM}}(\lambda x)$と一致します. また$y_{{\rm Taylor}}(x)$は, 2段2次の陽的ルンゲ-クッタ法の安定性関数$R_{{\rm 2}}(\lambda x)$と一致します.

上図からわかるように原点付近では, 近似解は厳密解をよく近似していますが, $y_{{\rm Taylor}}(x)$は$x\to \infty$で発散してしまいます.一方で $y_{{\rm Pade}}(x)$は,
$$
\lim_{x\to \infty} y_{{\rm Pade}}(x) = -1
$$
と収束します. このような性質が陽的ルンゲ-クッタ法に比べて, 陰的ルンゲ-クッタ法の安定性が高い要因になっていると思われます.

まとめ

p次の陽的ルンゲ-クッタ法の安定性関数は, 厳密解のテイラー展開とp次まで一致しました. またp次の陰的ルンゲ-クッタ法では, 安定性関数が有理関数(パデ近似関数)となりましたが, この有理関数をテイラー展開すると, 厳密解のテイラー展開とp次まで一致しました.

ソースコード

後日公開

参考文献

[1] Hairer, E., G. Wanner, "常微分方程式の数値解法 II." (2002).