はじめに
何を書くか迷っていたのですが、一つ前の方が数値計算の導入の話を書かれていたので、それに乗っかりたいと思います。
オイラー法/ホイン法/ルンゲクッタ法をつかった常微分方程式の数値解析超入門
こちらの記事で紹介して頂いたオイラー法、ホイン法、ルンゲ=クッタ法(RK4)は、より一般化された陽的ルンゲ=クッタ法の枠組みで捉えることができます。この記事では陽的ルンゲ=クッタ法を理解するためのポイントについて書きたいと思います。
1つ目のポイント
陽的ルンゲ=クッタ法の1つ目のポイントです。
- テイラー展開の高階導関数を使うことなく
- 様々な位置での１階導関数の値に
- 係数をかけて足したものを使うことで
- 高次の精度を実現する
これだけではわけがわからないので、具体例を用いて説明します。
前進オイラー法
まずは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})$のテイラー展開を１階導関数で打ち切った式と一致します。
y(x_{n+1}) = y(x_n) + hy'(x_n) + O(h^2)
局所打切り誤差$O(h^2)$が$h$の2乗なので、大域打切り誤差は１次下がった$O(h)$となります。
よってオイラー法は１次精度です。
ホイン法の精度
同様に、$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つ目のポイントで振り返りましょう。
- テイラー展開の2階導関数$y''(x_n)$を使うことなく
- $x_n$と$x_{n+1}$での１階導関数$y'(x_n) , y'(x_{n+1})$の値に
- 係数$\frac{1}{2}$をかけて足したもの$\frac{y'(x_n) + y'(x_{n+1})}{2}$を使うことで
- 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つ目のポイントは、様々な位置における１次導関数の値を求めるのに必要な$y$の値も、既知の１次導関数の値を使って求めるということです。
陽的ルンゲ=クッタ法
それでは以上の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$は様々な位置での１階導関数の値です。$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}
計算結果
