📂数値解析

明示的なルンゲ＝クッタ法

概要

微分方程式ソルバー、ルンゲ＝クッタ法を紹介する。よく使われる 4次ルンゲ＝クッタ法 RK4についての詳しい説明は、別記事にて発表されている。

ビルドアップ

次のように与えられる常微分方程式を考える。$y$は$t$に対する函数であり、$^{\prime}$は$t$に対する微分を意味する。 $$y^{\prime} = f(t,y),\quad t \ge t_{0},\quad y(t_{0}) = y_{0}$$ これを$t_{n}$から$t_{n+1} = t_{n} + h$まで積分すると（混乱を避けるため積分変数を$t$の代わりに$\tau$にしよう）、 $$\int_{t_{n}}^{t_{n+1}}y^{\prime}d\tau = \int_{t_{n}}^{t_{n+1}} f(\tau,y)d\tau$$ 左辺を整理して$\tau \equiv t_{n} + h\tau$に置き換えると次の式を得る。 $$\big[y(\tau)\big]_{t_{n}}^{t_{n+1}} = \int_{t_{n}}^{t_{n+1}} f(\tau,y)d\tau$$ $$\Downarrow \\ y(t_{n+1}) = y(t_{n}) + \int_{t_{n}}^{t_{n+1}} f(\tau,y)d\tau = y(t_{n}) + h\int_{0}^{1} f\big(t_{n} + h\tau,y(t_{n} + h\tau)\big)d\tau$$

今、後置きの積分を和として近似して次の式を得る。

$$y(t_{n+1}) = y(t_{n}) + h\sum_{j=1}^{\nu} b_{j} f\big(t_{n} + c_{j}h,y(t_{n} + c_{j}h)\big),\quad n=0,1,2,\cdots \tag{1}$$

これは、$y_{n} = y(t_{n})$を知っているとき、$y_{n+1} = y(t_{n+1})$を知ることができる近似式である。もちろん、右辺には未知数である$y(t_{n} + c_{j}h)$が含まれているので、これも近似する必要がある。式が複雑なので、次のように$j = 1,2,\dots,\nu$について表記しよう。

$$\xi_{j} = y(t_{n} + c_{j}h)$$

まず$c_{1} = 0$としておこう。すると$\xi_{1} = y_{n}$なので、既知の値である。ここで多段階法を適用して$\xi_{j}$たちを次のように逐次的に近似するとしよう。

$$\begin{align*} \xi_{2} &= y(t_{n} + c_{2}h) \\ &= y_{n} + h a_{2,1}f(t_{n} + c_{1}h, y_{n} + c_{1}h) \\ &= y_{n} + h a_{2,1}f(t_{n}, y_{n}) & (\because c_{1} = 0)\\ &= y_{n} + h a_{2,1}f(t_{n}, \xi_{1}) & (\because \xi_{1}=y_{n}) \\[2em] \xi_{3} &= y(t_{n} + c_{3}h) \\ &= y_{n} + h a_{3,1}f(t_{n} + c_{1}h, y_{n} + c_{1}h) + h a_{3,2}f(t_{n} + c_{2}h, y_{n} + c_{2}h) \\ &= y_{n} + h a_{3,1}f(t_{n}, \xi_{1}) + h a_{3,2}f(t_{n} + c_{2}h, \xi_{2}) \\ &\vdots \\ \xi_{\nu} &= y_{n} + h\sum_{j=1}^{n-1}a_{\nu,j} f(t_{n} + c_{j}h, \xi_{j}) \end{align*}$$

これで$(1)$の右辺の値$y_{n}$、$f(t_{n}, \xi_{1})$、$f(t_{n} + c_{2}h, \xi_{2})$、$\dots$、$f(t_{n} + c_{j-1}h, \xi_{j-1})$全てを求めたことになる。明示的ルンゲ＝クッタ^{explicit Runge-Kutta, ERK}メソッドとは、$y_{n+1}$の値を$y_{n}$、$f(t_{n}, \xi_{1})$、$f(t_{n} + c_{2}h, \xi_{2})$、$\dots$、$f(t_{n} + c_{j-1}h, \xi_{j-1})$の線形結合として近似して得る方法を言う。

定義

明示的ルンゲ＝クッタ法とは、与えられた$y_{n}$に対して$y_{n+1}$を次のように近似する方法である。

$$y_{n+1} = y_{n} + h\sum_{j=1}^{\nu} b_{j} f(t_{n} + c_{j}h, \xi_{j}),\quad n=0,1,2,\cdots \tag{2}$$

ここで、各$\xi_{j}$は、

$$\begin{align*} \xi_{1} &= y_{n} \\ \xi_{2} &= y_{n} + ha_{2,1}f(t_{n}, \xi_{1}) \\ \xi_{3} &= y_{n} + ha_{3,1}f(t_{n}, \xi_{1}) + ha_{3,2}f(t_{n} + c_{2}h, \xi_{2}) \\ \vdots & \\ \xi_{\nu} &= y_{n} + h\sum_{i=1}^{\nu-1}a_{\nu, i}f(t_{n} + c_{i}h, \xi_{i}) \end{align*}$$

この時$\nu \times \nu$ 行列 $A = [a_{j,i}]$をルンゲ＝クッタ行列^{RK matrix}と言う。式の中で空いている成分は$0$である。

$$A = [a_{j,i}] = \begin{bmatrix} \\ a_{2,1} & \\ a_{3,1} & a_{3,2} \\ \vdots & \\ a_{\nu,1} & a_{\nu,2} & \cdots & a_{\nu,\nu-1} & \end{bmatrix}$$

また、以下の$\mathbf{b}$と$\mathbf{c}$をそれぞれRKウェイト^{RK weights}、RKノード^{RK nodes}と言う。

$$\mathbf{b} = \begin{bmatrix} b_{1} \\ b_{2} \\ \vdots \\ b_{\nu} \end{bmatrix} \quad \text{and} \quad \mathbf{c} = \begin{bmatrix} c_{1} \\ c_{2} \\ \vdots \\ c_{\nu} \end{bmatrix}$$

また、$(2)$が**$\nu$ ステージ**^{$\nu$ stages}を持つと言い、このような方法を$\nu$ステージ（または$\nu$次）ルンゲ＝クッタ法と言う。

説明

定義により、RK行列は明らかに下三角行列である。

上記の方法に従って、$\xi_{j}$は$\xi_{1}$から$\xi_{\nu}$まで逐次的に計算することができる。結果として$(2)$を計算でき、$y_{n+1}$が得られる。これを繰り返して、次のステップでの$y$の値$y_{n+2}$、$y_{n+3}$、$\dots$を求めることができる。

RK法では、係数$a_{j,i}$、$b_{j}$、$c_{j}$は求めるものではなく、選んで使用するものである。係数は$\begin{array}{c|c} \mathbf{c} & A \\ \hline & \mathbf{b}^{t} \end{array}$のような形で表記され、これをRKテーブル^{RK tableaux}と言う。よく使われるRK4のRKテーブルは次のようである。

$$\begin{array}{c|cccc} 0 & \\[0.5em] \frac{1}{2} & \frac{1}{2} \\[0.5em] \frac{1}{2} & 0 & \frac{1}{2} \\[0.5em] 1 & 0 & 0 & 1 \\[0.5em] \hline {\displaystyle {\color{white}\dfrac{0}{0}}} & \frac{1}{6} & \frac{1}{3} & \frac{1}{3} & \frac{1}{6} \end{array}$$

参考

• 4次ルンゲ＝クッタ法
• 暗黙のルンゲ＝クッタ法