数値計算
数学

あれもこれも陽的ルンゲ=クッタ法

はじめに

何を書くか迷っていたのですが、一つ前の方が数値計算の導入の話を書かれていたので、それに乗っかりたいと思います。

オイラー法/ホイン法/ルンゲクッタ法をつかった常微分方程式の数値解析超入門

こちらの記事で紹介して頂いたオイラー法、ホイン法、ルンゲ=クッタ法(RK4)は、より一般化された陽的ルンゲ=クッタ法の枠組みで捉えることができます。この記事では陽的ルンゲ=クッタ法を理解するためのポイントについて書きたいと思います。

1つ目のポイント

陽的ルンゲ=クッタ法の1つ目のポイントです。

  1. テイラー展開の高階導関数を使うことなく
  2. 様々な位置での1階導関数の値に
  3. 係数をかけて足したものを使うことで
  4. 高次の精度を実現する

これだけではわけがわからないので、具体例を用いて説明します。

前進オイラー法

まずは1階導関数が$x$のみに依存する関数$f(x)$で表される場合を考えます。

\frac{dy}{dx} = f(x)

このとき前進オイラー法は以下の式で表されます。

y_{n+1} = y_n + hf(x_n)

$h$は刻み幅です。図にすると下記のようになります。

オイラー法では、$x_n$での傾き$f(x_n)$のみを使います。

ホイン法

同様に、ホイン法は以下の式で表されます。

y_{n+1} = y_n + h\frac{f(x_n) + f(x_{n+1})}{2}  

ホイン法は、$x_n$での傾き$f(x_n)$と$x_{n+1}$での傾き$f(x_{n+1})$に、係数$\frac{1}{2}$をかけて足すことで平均を取っています。

前進オイラー法の精度

前進オイラー法は、$y(x_{n+1})$のテイラー展開を1階導関数で打ち切った式と一致します。

y(x_{n+1}) = y(x_n) + hy'(x_n) + O(h^2)

局所打切り誤差$O(h^2)$が$h$の2乗なので、大域打切り誤差は1次下がった$O(h)$となります。
よってオイラー法は1次精度です。

ホイン法の精度

同様に、$y(x_{n+1})$をテイラー展開します。今回は2階の導関数まで書きます。

y(x_{n+1}) = y(x_n) + hy'(x_n) + \frac{h^2}{2}y''(x_n) + O(h^3)

2階の導関数を消すために、$y'(x_{n+1})$のテイラー展開を使います。

y'(x_{n+1}) = y'(x_n) + hy''(x_n)  + O(h^2)

この式を先程の式に代入して2階の導関数を消します。

\begin{align}
y(x_{n+1}) &= y(x_n) + hy'(x_n) + \frac{h^2}{2}y''(x_n) + O(h^3) \\
                      &= y(x_n) + \frac{h}{2} y'(x_n) + \frac{h}{2} \{ y'(x_n) + hy''(x_n) \}+ O(h^3) \\
                      &= y(x_n) + \frac{h}{2} y'(x_n) + \frac{h}{2} y'(x_{n+1})+ O(h^3) \\
                      &= y(x_n) + h\frac{y'(x_n) + y'(x_{n+1})}{2}  + O(h^3)
\end{align}

なんということでしょう、2次精度になりました。

ホイン法を1つ目のポイントで振り返りましょう。

  1. テイラー展開の2階導関数$y''(x_n)$を使うことなく
  2. $x_n$と$x_{n+1}$での1階導関数$y'(x_n) , y'(x_{n+1})$の値に
  3. 係数$\frac{1}{2}$をかけて足したもの$\frac{y'(x_n) + y'(x_{n+1})}{2}$を使うことで
  4. 2次の精度を実現する

2つ目のポイント

続いて、1階導関数が$y(x)$にも依存する関数$f(x,y(x))$で表される場合を考えます。

\frac{dy}{dx} = f(x,y(x))

この時、ホイン法は先程と同じように、下記のように書けそうです。

y_{n+1} = y_n + h\frac{f(x_n, y_n) + f(x_{n+1}, y_{n+1})}{2}  

しかし$x_{n+1}$での傾き$f(x_{n+1}, y_{n+1})$には、求めたい$y_{n+1}$が入っているため使えません。仕方ないので$x_{n+1}$での$y$の値は、傾き$f(x_n, y_n)$から求めた値$y_n + hf(x_n,y_n)$を使います。

y_{n+1} = y_n + h\frac{f(x_n, y_n) + f(x_{n+1}, y_n+hf(x_n, y_n) )}{2}  

これが1階導関数が$f(x,y(x))$のような関数で表される場合のホイン法です。
図にすると下記のようになります。

以上より2つ目のポイントは、様々な位置における1次導関数の値を求めるのに必要な$y$の値も、既知の1次導関数の値を使って求めるということです。

陽的ルンゲ=クッタ法

それでは以上の2点のポイントを踏まえて、陽的ルンゲ=クッタ法を見ていきましょう。
1階導関数が$y(x)$にも依存する関数$f(x,y(x))$で表される場合、

\frac{dy}{dx} = f(x,y(x))

$s$段陽的ルンゲ=クッタ法は以下の式で表されます。

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

ここで$k$は様々な位置での1階導関数の値です。$s$は段数と呼ばれ、$k$の数です。
$c_i$により$k_i$の$x$の位置($x_n + hc_i$)が決まります。この時の$y$の値($y_n + h \sum_{j=1}^{i-1} a_{ij} k_j$)は、既知の$k_j$の値に、係数$a_{ij}$をかけて足したものを使って求めます。総和の範囲が$j=1,2,...,i-1$となっているため、$k_i$を求めるときの、$k_j$の値は全て求まっています。最後に、求めた$k_i$の値に、係数$b_i$をかけて足すことで、最終的な傾き($\sum_{i=1}^s b_ik_i$)を求めます。

陽的ルンゲ=クッタ法の係数$a, b, c$を変えるだけで、前進オイラー法やホイン法を含めた様々な種類・精度の手法で計算できます。また係数$a, b, c$をわかりやすく表記するためにブッチャー配列が使われます。$s$段陽的ルンゲ=クッタ法のブッチャー配列は下記になります。

\begin{array}
{c|cccccc}
0\\
c_2    & a_{21} \\
c_3    & a_{31} & a_{32} \\
\vdots & \vdots & \vdots & \ddots \\
c_s    & a_{s1} & a_{s2} & \cdots & a_{s,s-1} \\
\hline
& b_1 & b_2 & \cdots & b_{s-1} & b_s
\end{array}

それでは、前進オイラー法、ホイン法、ルンゲ=クッタ法(RK4)が、陽的ルンゲ=クッタ法でどのように表されるか見ていきます。

1段1次ルンゲ=クッタ法 : 前進オイラー法

1段ルンゲ=クッタ法

\begin{align}
y_{n+1} &= y_n + h (b_1k_1) \\
k_1     &= f(x_n, y_n)
\end{align}
\begin{array}
{c|c}
0\\
\hline
& b_1 
\end{array}

1段1次ルンゲ=クッタ法 : 前進オイラー法

\begin{align}
y_{n+1} &= y_n + h (1 \cdot k_i) \\
k_1     &= f(x_n, y_n)
\end{align}
\begin{array}
{c|c}
0\\
\hline
& 1
\end{array}

2段2次ルンゲ=クッタ法 : ホイン法

2段ルンゲ=クッタ法

\begin{align}
y_{n+1} &= y_n + h (b_1k_1+b_2k_2) \\
k_1     &= f(x_n, y_n) \\
k_2     &= f(x_n + c_2h, y_n + h(a_{21}k_1) )
\end{align}
\begin{array}
{c|cc}
0\\
c_2 & a_{21}\\
\hline
& b_1 & b_2
\end{array}

2段2次ルンゲ=クッタ法 : ホイン法

\begin{align}
y_{n+1} &= y_n + h (\frac{1}{2}k_1+\frac{1}{2}k_2) \\
k_1     &= f(x_n, y_n) \\
k_2     &= f(x_n + 1 \cdot h, y_n + h(1 \cdot k_1) )
\end{align}
\begin{array}
{c|cc}
0\\
1 & 1\\
\hline
& 1/2 & 1/2
\end{array}

4段4次ルンゲ=クッタ法 : RK4

4段ルンゲ=クッタ法

\begin{align}
y_{n+1} &= y_n + h (b_1k_1+b_2k_2+b_3k_3+b_4k_4) \\
k_1     &= f(x_n, y_n) \\
k_2     &= f(x_n + c_2h, y_n + h(a_{21}k_1) ) \\
k_3     &= f(x_n + c_3h, y_n + h(a_{31}k_1 + a_{32}k_2) ) \\
k_4     &= f(x_n + c_4h, y_n + h(a_{41}k_1 + a_{42}k_2 + a_{43}k_3) ) \\
\end{align}
\begin{array}
{c|cccc}
0\\
c_2 & a_{21}\\
c_3 & a_{31} & a_{32}\\
c_4 & a_{41} & a_{42} & a_{43} \\
\hline
& b_1 & b_2 & b_3 & b_4
\end{array}

4段4次ルンゲ=クッタ法 : RK4

\begin{align}
y_{n+1} &= y_n + h (\frac{1}{6} k_1 + \frac{1}{3} k_2 + \frac{1}{3} k_3 + \frac{1}{6} k_4) \\
k_1     &= f(x_n, y_n) \\
k_2     &= f(x_n + \frac{1}{2} h, y_n + h(\frac{1}{2}  k_1) ) \\
k_3     &= f(x_n + \frac{1}{2} h, y_n + h(0 \cdot k_1 + \frac{1}{2}  k_2) ) \\
k_4     &= f(x_n + 1 \cdot h, y_n + h(0 \cdot k_1 + 0 \cdot k_2 + 1 \cdot k_3) ) \\
\end{align}
\begin{array}
{c|cccc}
0\\
1/2 & 1/2\\
1/2 & 0 & 1/2\\
1   & 0 & 0 & 1 \\
\hline
& 1/6 & 1/3 & 1/3 & 1/6
\end{array}

計算例

最後に計算結果を示して終わりたいと思います。

微分方程式

\frac{dy}{dx} = f(y) = \lambda y

初期値

y(0) = y_0

厳密解

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

パラメータ

\begin{align}
\lambda &= -1 \\
y_0     &= 1  \\
h       &= 0.5   
\end{align}

計算結果