明示的なルンゲ＝クッタ法
概要
微分方程式ソルバー、ルンゲ＝クッタ法を紹介する。よく使われる 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} $$