📂数値解析

ルンゲ-クッタ法で係数を決定する方法

概要¹

次のように与えられた常微分方程式について考える。$y$は$t$に対する関数であり、プライム$(^{\prime})$は$t$に対する微分を意味する。

$$y^{\prime} = f(t,y),\quad t \ge t_{0},\quad y(t_{0}) = y_{0}, \quad t_{n+1} = t_{n} + h$$

明示的ルンゲ・クッタ法とは、与えられた$y_{n} = y(t_{n})$に対して$y_{n+1} = y(t_{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$$

ここでそれぞれの$\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*}$$

RK法では以下の係数は選択して使用するものである。適切な係数を選択する方法について調べる。 $$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} \quad \text{and} \quad \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}$$

$\nu = 2$の場合

まず最も簡単な$\nu = 2$の場合について見てみよう。すると$\xi_{1} = y_{n}$であり$\xi_{2} = y_{n} + ha_{2,1}f(t_{n}, y_{n})$である。従って$y_{n+1}$は次の通りである。

$$\begin{align} y_{n+1} &= y_{n} + h \Big[ b_{1}f(t_{n} + c_{1}h, \xi_{1}) + b_{2}f(t_{n} + c_{2}h, \xi_{2}) \Big] \nonumber \\ &= y_{n} + h \Big[ b_{1}f(t_{n}, y_{n}) + b_{2}f(t_{n} + c_{2}h, y_{n} + ha_{2,1}f(t_{n}, y_{n}) )\Big] \end{align}$$

ここで$c_{1} = 0$である。ベクトル値関数$f$が十分にスムースであると仮定しよう。すると、$f(t_{n} + c_{2}h, y_{n} + ha_{2,1}f(t_{n}, y_{n}) )$の点$(t_{n}, y_{n})$でのテイラー展開は次の通りである。

$$\begin{align*} & f(t_{n} + c_{2}h, y_{n} + ha_{2,1}f(t_{n}, y_{n}) ) \\ &= f(t_{n}, y_{n}) + (c_{2}h, ha_{2,1}f(t_{n}, y_{n}))^{T} \nabla f(t_{n}, y_{n}) + \mathcal{O}(h^{2}) \\ &= f(t_{n}, y_{n}) + h\left[c_{2} \dfrac{\partial f(t_{n}, y_{n})}{\partial t} + a_{2,1}\dfrac{\partial f(t_{n}, y_{n})}{\partial y}f(t_{n}, y_{n}) \right] + \mathcal{O}(h^{2}) \\ \end{align*}$$

$\nabla f$は$f$の勾配を意味する。これを$(1)$に代入すると、

$$\begin{align*} y_{n+1} &= y_{n} + h \bigg[ b_{1}f(t_{n}, y_{n}) + b_{2}f(t_{n}, y_{n}) \\ & \qquad + hb_{2}c_{2} \dfrac{\partial f(t_{n}, y_{n})}{\partial t} + hb_{2}a_{2,1}\dfrac{\partial f(t_{n}, y_{n})}{\partial y}f(t_{n}, y_{n}) + \mathcal{O}(h^{2}) \bigg] \\ &= y_{n} + h (b_{1} + b_{2})f(t_{n}, y_{n}) \\ & \qquad + h^{2}b_{2} \left[ c_{2} \dfrac{\partial f(t_{n}, y_{n})}{\partial t} + a_{2,1}\dfrac{\partial f(t_{n}, y_{n})}{\partial y}f(t_{n}, y_{n}) \right] + \mathcal{O}(h^{3}) \tag{2} \end{align*}$$

ところで、$y^{\prime} = f(t, y)$の全微分は$d (y^{\prime}) = \dfrac{\partial f(t, y)}{\partial t}dt + \dfrac{\partial f(t, y)}{\partial y}dy$であり、$^{\prime} = \dfrac{d}{dt}$なので、$y^{\prime \prime}$は次の通りである。

$$y^{\prime \prime} = \dfrac{\partial f(t, y)}{\partial t} + \dfrac{\partial f(t, y)}{\partial y}\dfrac{dy}{dt} = \dfrac{\partial f(t, y)}{\partial t} + \dfrac{\partial f(t, y)}{\partial y}f(t, y)$$

テイラーの定理により次を得る。

$$\begin{align*} y_{n+1} &= y_{n} + hy^{\prime} + \dfrac{h^{2}}{2}y^{\prime \prime} + \mathcal{O}(h^{3}) \\ &= y_{n} + hf(t_{n},y_{n}) + \dfrac{h^{2}}{2}\left( \dfrac{\partial f(t, y)}{\partial t} + \dfrac{\partial f(t, y)}{\partial y}f(t, y) \right) + \mathcal{O}(h^{3}) \tag{3} \\ \end{align*}$$

$(2)$と$(3)$を比較すると次の関係式を得る。

$$b_{1} + b_{2} = 1, \qquad b_{2}c_{2} = \dfrac{1}{2}, \qquad a_{2,1} = c_{2}$$

これを満足するRKテーブルとしては、以下のものがある。

$$\begin{array}{c|c} \mathbf{c} & A \\ \hline & \mathbf{b}^{t} \end{array} \qquad \begin{array}{c|c} 0 \\ \frac{1}{2} & \frac{1}{2} \\[0.5em] \hline & 0 & 1 \end{array}, \qquad \begin{array}{c|c} 0 \\ \frac{2}{3} & \frac{2}{3} \\[0.5em] \hline \phantom{\dfrac{1}{2}} & \frac{1}{4} & \frac{3}{4} \end{array}, \qquad \begin{array}{c|c} 0 \\ 1 & 1 \\[0.5em] \hline \phantom{\dfrac{1}{2}} & \frac{1}{2} & \frac{1}{2} \end{array}, \qquad$$

$\nu = 3$の場合

まずテイラーの定理により次を得る。

$$\begin{aligned} y_{n+1} &= y_{n} + hy^{\prime} + \dfrac{h^{2}}{2}y^{\prime \prime} + \dfrac{h^{3}}{6}y^{\prime \prime \prime} + \mathcal{O}(h^{4}) \\ &= y_{n} + hf_{n} + \dfrac{h^{2}}{2}\left( \dfrac{\partial f}{\partial t} + \dfrac{\partial f}{\partial y}f \right) + \dfrac{h^{3}}{6}\left( \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{\partial^{2} f}{\partial y^{2}}f^{2} + \dfrac{\partial f}{\partial t}\dfrac{\partial f}{\partial y} \right) + \mathcal{O}(h^{4}) \tag{4} \end{aligned}$$

そして$\sum\limits_{i=1}^{j-1} a_{j,i} = c_{j}$と仮定しよう。すると$\nu = 3$の場合には、

$$c_{1} = 0,\qquad c_{2} = a_{2,1},\qquad c_{3} = a_{3,1} + a_{3,2}$$

$f_{n} = f(t_{n}, y_{n})$と表記すると、

$$\xi_{1} = y_{n}, \quad \xi_{2} = y_{n} + hc_{2} f_{n}, \\[1em] \xi_{3} = y_{n} + ha_{3,1}f_{n} + ha_{3,2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n})$$

従って$y_{n+1}$は次の通りである。

$$\begin{align*} y_{n+1} &= y_{n} + hb_{1}f_{n} + hb_{2}f(t_{n} + c_{2}h, \xi_{2}) + hb_{3}f(t_{n} + c_{3}h, \xi_{3}) \\ &= y_{n} + hb_{1}f_{n} + hb_{2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n}) \\ & \qquad + hb_{3}f(t_{n} + c_{3}h, y_{n} + ha_{3,1}f_{n} + ha_{3,2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n})) \tag{5} \end{align*}$$

$b_{2}$が含まれる項は次の通りである。

$$\begin{align*} & f(t_{n} + c_{2}h, y_{n} + hc_{2}f(t_{n}, y_{n}) ) \\ &= f_{n} + (hc_{2}) \dfrac{\partial f}{\partial t} + (hc_{2}f_{n})\dfrac{\partial f}{\partial y} + \dfrac{(c_{2}h)^{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{(c_{2}h f_{n})^{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}} + \mathcal{O}(h^{3}) \\ \end{align*}$$

これを$(5)$に代入すると、

$$\begin{aligned} y_{n+1} &= y_{n} + h(b_{1} + b_{2})f_{n} + h^{2}b_{2}c_{2} \dfrac{\partial f}{\partial t} + h^{2}b_{2}c_{2} \dfrac{\partial f}{\partial y}f_{n} \\ & \qquad + \dfrac{h^{3}c_{2}^{2}b_{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{h^{3}c_{2}^{2}b_{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}}f_{n}^{2} \\ & \qquad + hb_{3}f(t_{n} + c_{3}h, y_{n} + ha_{3,1}f_{n} + ha_{3,2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n})) + \mathcal{O}(h^{4}) \tag{6} \end{aligned}$$

今度は$b_{3}$が含まれる項を計算してみよう（係数比較に必要ない項は省略する）。

$$\begin{align*} & f(t_{n} + c_{3}h, y_{n} + ha_{3,1}f_{n} + ha_{3,2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n})) \\ &= f_{n} + c_{3}h \dfrac{\partial f}{\partial t} + \left( ha_{3,1}f_{n} + ha_{3,2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n}) \right) \dfrac{\partial f}{\partial y} \\ & \qquad + \dfrac{h^{2}c_{3}^{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{\left( ha_{3,1}f_{n} + ha_{3,2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n}) \right)^{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}} + \mathcal{O}(h^{3}) \\ &= f_{n} + c_{3}h \dfrac{\partial f}{\partial t} + ha_{3,1}f_{n}\dfrac{\partial f}{\partial y} + ha_{3,2} \left[ f_{n} + hc_{2}\dfrac{\partial f}{\partial t} + \mathcal{O}(h^{2}) \right] \dfrac{\partial f}{\partial y} \\ & \qquad + \dfrac{h^{2}c_{3}^{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{(ha_{3,1}f_{n} + ha_{3,2}f_{n})^{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}} + \mathcal{O}(h^{3}) \\ &= f_{n} + c_{3}h\dfrac{\partial f}{\partial t} + hc_{3}f_{n}\dfrac{\partial f}{\partial y} + h^{2}c_{2}a_{3,2}\dfrac{\partial f}{\partial t} \dfrac{\partial f}{\partial y} \\ & \qquad + \dfrac{h^{2}c_{3}^{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{(hc_{3}f_{n})^{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}} + \mathcal{O}(h^{3}) \\ \end{align*}$$

これを$(6)$に代入し整理すると次の通りである。

$$\begin{align*} y_{n+1} &= y_{n} + h(b_{1} + b_{2})f_{n} + h^{2}b_{2}c_{2} \dfrac{\partial f}{\partial t} + h^{2}b_{2}c_{2} \dfrac{\partial f}{\partial y}f_{n} \\ & \qquad + \dfrac{h^{3}c_{2}^{2}b_{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{h^{3}c_{2}^{2}b_{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}}f_{n}^{2} \\ & \qquad + hb_{3} \left[ f_{n} + c_{3}h\dfrac{\partial f}{\partial t} + hc_{3}f_{n}\dfrac{\partial f}{\partial y} + h^{2}c_{2}a_{3,2}\dfrac{\partial f}{\partial t} \dfrac{\partial f}{\partial y} \right. \\ & \qquad + \left. \dfrac{h^{2}c_{3}^{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{(hc_{3}f_{n})^{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}} \right] + \mathcal{O}(h^{4}) \\ &= y_{n} + h(b_{1} + b_{2} + b_{3})f_{n} + h^{2}(b_{2}c_{2} + b_{3}c_{3}) \dfrac{\partial f}{\partial t} + h^{2}(b_{2}c_{2} + b_{3}c_{3}) \dfrac{\partial f}{\partial y}f_{n} \\ & \qquad + \dfrac{h^{3}(b_{2}c_{2}^{2} + b_{3}c_{3}^{2})}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + h^{3}b_{3}c_{2}a_{3,2}\dfrac{\partial f}{\partial t} \dfrac{\partial f}{\partial y} + \dfrac{h^{3}(b_{2}c_{2}^{2} + b_{3}c_{3}^{2})}{2} f_{n}^{2}\dfrac{\partial^{2} f}{\partial y^{2}} + \mathcal{O}(h^{4}) \end{align*}$$

上記の式と$(4)$の係数を比較すると、次の関係式を得る。

$$b_{1} + b_{2} + b_{3} = 1, \quad b_{2}c_{2} + b_{3}c_{3} = \dfrac{1}{2}, \quad b_{2}c_{2}^{2} + b_{3}c_{3}^{2} = \dfrac{1}{3}, \quad b_{3}c_{2}a_{3,2} = \dfrac{1}{6}$$

これを満足するRKテーブルとしては、以下のものがある。

$$\begin{array}{c|c} 0 \\ \frac{1}{2} & \frac{1}{2} \\[0.5em] 1 & -1 & 2 \\[0.5em] \hline \phantom{\dfrac{1}{2}} & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \end{array} \quad (\text{classical RK method})$$ $$\begin{array}{c|c} 0 \\ \frac{2}{3} & \frac{2}{3} \\[0.5em] \frac{2}{3} & 0 & \frac{2}{3} \\[0.5em] \hline \phantom{\dfrac{1}{2}} & \frac{1}{4} & \frac{3}{8} & \frac{3}{8} \end{array} \quad (\text{Nystrom scheme})$$

$\nu = 4$の場合

同様の方法により係数に関する条件を求めることができ、上記で見た規則によれば次の通りであると推測できる。

$$b_{1} + b_{2} + b_{3} + b_{4} = 1, \quad b_{2}c_{2} + b_{3}c_{3} + b_{4}c_{4} = \dfrac{1}{2}, \\[1em] b_{2}c_{2}^{2} + b_{3}c_{3}^{2} + b_{4}c_{4}^{2} = \dfrac{1}{3}, \quad b_{2}c_{2}^{3} + b_{3}c_{3}^{3} + b_{4}c_{4}^{3} = \dfrac{1}{4}$$

RK4で最も有名な係数は次の通りである。以下のテーブルは実際に上記の条件を満たしている。

$$\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次のルンゲ・クッタ法
  • 係数を決定する方法
• 暗示的ルンゲ・クッタ法