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

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

[オイラー法/ホイン法/ルンゲクッタ法をつかった常微分方程式の数値解析超入門](https://qiita.com/Dr_ASA/items/d44560cf1367290aee6d)

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

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

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

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

# 前進オイラー法
まずは1階導関数が$x$のみに依存する関数$f(x)$で表される場合を考えます。

```math
\frac{dy}{dx} = f(x)
```
このとき前進オイラー法は以下の式で表されます。

```math
y_{n+1} = y_n + hf(x_n)
```
$h$は刻み幅です。図にすると下記のようになります。

<img src="https://qiita-image-store.s3.amazonaws.com/0/35426/ad730452-3ad6-75ef-f8f7-121a3245efc0.png" width=60%>

オイラー法では、$x_n$での傾き$f(x_n)$のみを使います。
# ホイン法
同様に、ホイン法は以下の式で表されます。

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

<img src="https://qiita-image-store.s3.amazonaws.com/0/35426/134b5a60-5022-5a43-8ff0-aa01583d3711.png" width=60%>

ホイン法は、$x_n$での傾き$f(x_n)$と$x_{n+1}$での傾き$f(x_{n+1})$に、係数$\frac{1}{2}$をかけて足すことで平均を取っています。
# 前進オイラー法の精度
前進オイラー法は、$y(x_{n+1})$のテイラー展開を１階導関数で打ち切った式と一致します。

```math
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階の導関数まで書きます。

```math
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})$のテイラー展開を使います。

```math
y'(x_{n+1}) = y'(x_n) + hy''(x_n)  + O(h^2)
```
この式を先程の式に代入して2階の導関数を消します。

```math
\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}$での１階導関数$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))$で表される場合を考えます。

```math
\frac{dy}{dx} = f(x,y(x))
```
この時、ホイン法は先程と同じように、下記のように書けそうです。

```math
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)$を使います。

```math
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))$のような関数で表される場合のホイン法です。
図にすると下記のようになります。
<img src="https://qiita-image-store.s3.amazonaws.com/0/35426/a34483c1-4d76-92c2-210b-b6ee9cbeef59.png" width=80%>

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

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

```math
\frac{dy}{dx} = f(x,y(x))
```
$s$段陽的ルンゲ=クッタ法は以下の式で表されます。

```math
\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$段陽的ルンゲ=クッタ法のブッチャー配列は下記になります。

```math
\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段ルンゲ=クッタ法**

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

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

```math
\begin{align}
y_{n+1} &= y_n + h (1 \cdot k_i) \\
k_1     &= f(x_n, y_n)
\end{align}
```
```math
\begin{array}
{c|c}
0\\
\hline
& 1
\end{array}
```
# 2段2次ルンゲ=クッタ法 : ホイン法
**2段ルンゲ=クッタ法**

```math
\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}
```
```math
\begin{array}
{c|cc}
0\\
c_2 & a_{21}\\
\hline
& b_1 & b_2
\end{array}
```

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

```math
\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}
```
```math
\begin{array}
{c|cc}
0\\
1 & 1\\
\hline
& 1/2 & 1/2
\end{array}
```
# 4段4次ルンゲ=クッタ法 : RK4
**4段ルンゲ=クッタ法**

```math
\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}
```
```math
\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**

```math
\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}
```
```math
\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}
```
# 計算例
最後に計算結果を示して終わりたいと思います。

**微分方程式**

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

**初期値**

```math
y(0) = y_0
```

**厳密解**

```math
y(x) = y_0 e^{\lambda x}
```
**パラメータ**

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

**計算結果**
<img src="https://qiita-image-store.s3.amazonaws.com/0/35426/25e485c3-702a-7f39-5e77-00bf26cdb1f9.png" width=70%>