読者です 読者をやめる 読者になる 読者になる

倭算数理研究所

科学・数学・学習関連の記事を、「倭マン日記」とは別に書いていくのだ!

ルンゲ=クッタ法を導く (6) : 4次の場合

シミュレーション技法

ルンゲ=クッタ (Runge-Kutta) 法を導くシリーズ(目次)。 今回は4次のルンゲ=クッタ法。 「問題設定」の { s=4 } の場合。

満たすべき関係式

満たすべき関係式は

  { \displaystyle\begin{align*}
  \sum_{i=1}^4 a_i\left[k_i\right]_{h=0} &= f(x_n,\,y_n) & \cdots (1) \\
  \sum_{i=1}^4 a_i\left[\frac{dk_i}{dh}\right]_{h=0} &= \frac{1}{2}\frac{df}{dx}(x_n,\,y_n) & \cdots (2) \\
  \sum_{i=1}^4 a_i\left[\frac{d^2k_i}{dh^2}\right]_{h=0} &= \frac{1}{3}\frac{d^2f}{dx^2}(x_n,\,y_n) & \cdots (3) \\
  \sum_{i=1}^4 a_i\left[\frac{d^3k_i}{dh^3}\right]_{h=0} &= \frac{1}{4}\frac{d^3f}{dx^3}(x_n,\,y_n) & \cdots (4)
\end{align*}}

ここで

  { \displaystyle\begin{align*}
    k_1 &= f(x_n,\,y_n) \\
    k_i &= f(x_n+c_ih,\,y_n+c_ihk_{i-1}) & (i = 2,\,3,\,4)
\end{align*}
}

です。

(1) ~ (3) 式

s = 3 の場合と同様にして、(1) 式より

  { \displaystyle\begin{align*}
    a_1 + a_2 + a_3 + a_4 = 1
\end{align*}}

(2) 式より

  { \displaystyle\begin{align*}
    a_2c_2 + a_3c_3 + a_4c_4= \tfrac{1}{2}
\end{align*}}

(3) 式より

  { \displaystyle\begin{align*}
    a_2c_2^2 + a_3c_3^2 + a_4c_4^2 &= \tfrac{1}{3} \\
    2\left(a_3c_3c_2 + a_4c_4c_3\right) &= \tfrac{1}{3}
\end{align*}}

(4) 式: { k } の3階微分の関係式

{ k_1 }微分

  { \displaystyle\begin{align*}
    \left[\frac{d^3k_1}{dh^3}\right]_{h=0} = 0
\end{align*}}

{ k_i }微分
{ i \ge 2 } の場合。 長くなるので、所々引数を省略しています:

  { \displaystyle\begin{align*}
    \left[\frac{d^3k_i}{dh^3}\right]_{h=0}
        &= c_i^3\left(\frac{\partial^3 f}{\partial x^3} + 3\frac{\partial^3 f}{\partial x^2 \partial y}f
            + 3\frac{\partial^3 f}{\partial x \partial y^2}f^2 + \frac{\partial^3 f}{\partial y^3}f^3 \right) \\
        &\qquad + 6c_i^2\left(\frac{\partial^2 f}{\partial x \partial y}
            + \frac{\partial^2 f}{\partial y^2}f\right)\left[\frac{dk_{i-1}}{dh}\right]_{h=0}
            + 3c_i\frac{\partial f}{\partial y}\left[\frac{d^2k_{i-1}}{dh^2}\right]_{h=0}
\end{align*}}

よって

  { \displaystyle\begin{align*}
    \left[\frac{d^3k_2}{dh^3}\right]_{h=0}
        = c_2^3\left(\frac{\partial^3 f}{\partial x^3} + 3\frac{\partial^3 f}{\partial x^2 \partial y}f
            + 3\frac{\partial^3 f}{\partial x \partial y^2}f^2 + \frac{\partial^3 f}{\partial y^3}f^3\right)
\end{align*}}

  { \displaystyle\begin{align*}
  \left[\frac{d^3k_3}{dh^3}\right]_{h=0}
    &= c_3^3\left(\frac{\partial^3 f}{\partial x^3} + 3\frac{\partial^3 f}{\partial x^2 \partial y}f
      + 3\frac{\partial^3 f}{\partial x \partial y^2}f^2 + \frac{\partial^3 f}{\partial y^3}f^3 \right) \\
    &\qquad + 6c_3^2c_2\left(\frac{\partial^2 f}{\partial x \partial y} 
      + \frac{\partial^2 f}{\partial y^2}f\right)\frac{df}{dx} \\
    &\qquad + 3c_3c_2^2\frac{\partial f}{\partial y}\left(\frac{\partial^2 f}{\partial x^2}
      + 2\frac{\partial^2 f}{\partial x \partial y}f + \frac{\partial^2 f}{\partial y^2}f^2\right)
\end{align*}}

  { \displaystyle\begin{align*}
  \left[\frac{d^3k_4}{dh^3}\right]_{h=0}
    &= c_4^3\left(\frac{\partial^3 f}{\partial x^3} + 3\frac{\partial^3 f}{\partial x^2 \partial y}f
      + 3\frac{\partial^3 f}{\partial x \partial y^2}f^2 + \frac{\partial^3 f}{\partial y^3}f^3 \right) \\
    &\qquad + 6c_4^2c_3\left(\frac{\partial^2 f}{\partial x \partial y}
      + \frac{\partial^2 f}{\partial y^2}f\right)\frac{df}{dx} \\
    &\qquad + 3c_4c_3^2\frac{\partial f}{\partial y}\left(\frac{\partial^2 f}{\partial x^2}
      + 2\frac{\partial^2 f}{\partial x \partial y}f+ \frac{\partial^2 f}{\partial y^2}f^2\right)\\
    &\qquad + 6c_4c_3c_2\left(\frac{\partial f}{\partial y}\right)^2\frac{df}{dx}
\end{align*}}

{ f(x,\,y) } の3階微分

  { \displaystyle\begin{align*}
  \frac{d^3f}{dx^3}
    &= \frac{\partial^3 f}{\partial x^3} + 3\frac{\partial^3 f}{\partial x^2 \partial y}f
      + 3\frac{\partial^3 f}{\partial x \partial y^2}f^2
      + \frac{\partial^3 f}{\partial y^3}f^3 \\ &\qquad + 3\left(\frac{\partial^2 f}{\partial x \partial y}
      + \frac{\partial^2 f}{\partial y^2}f\right)\frac{df}{dx} + \frac{\partial f}{\partial y}\frac{d^2f}{dx^2} \\
    &= \frac{\partial^3 f}{\partial x^3} + 3\frac{\partial^3 f}{\partial x^2 \partial y}f
      + 3\frac{\partial^3 f}{\partial x \partial y^2}f^2
      + \frac{\partial^3 f}{\partial y^3}f^3 \\ &\qquad + 3\left(\frac{\partial^2 f}{\partial x \partial y}
      + \frac{\partial^2 f}{\partial y^2}f\right)\frac{df}{dx} \\
    &\qquad + \frac{\partial f}{\partial y}\left(\frac{\partial^2 f}{\partial x^2}
      + 2\frac{\partial^2 f}{\partial x \partial y}f + \frac{\partial^2 f}{\partial y^2}f^2
      + \frac{\partial f}{\partial y}\frac{df}{dx}\right)
\end{align*}}

得られる関係式
上記で得られた関係を (4) 式に代入すると次式が得られます:

  { \displaystyle
\begin{align*}
    a_2c_2^3 + a_3c_3^3 + a_4c_4^3 &= \tfrac{1}{4} \\
    2\left(a_3c_3^2c_2 + a_4c_4^2c_3\right) &= \tfrac{1}{4} \\
    3\left(a_3c_3c_2^2 + a_4c_4c_3^2\right) &= \tfrac{1}{4} \\
    6a_4c_4c_3c_2 &= \tfrac{1}{4}
\end{align*}}

{ a,\,c } の値

{ a,\,c } が満たすべき関係式をまとめると

  { \displaystyle\begin{align*}
    a_1 + a_2 + a_3 + a_4 &= 1 \\[2mm]
    a_2c_2 + a_3c_3 + a_4c_4 &= \tfrac{1}{2} \\[2mm]
    a_2c_2^2 + a_3c_3^2 + a_4c_4^2 &= \tfrac{1}{3} \\[2mm]
    2\left(a_3c_3c_2 + a_4c_4c_3\right) &= \tfrac{1}{3} \\[2mm]
    a_2c_2^3 + a_3c_3^3 + a_4c_4^3 &= \tfrac{1}{4} \\[2mm]
    2\left(a_3c_3^2c_2 + a_4c_4^2c_3\right) &= \tfrac{1}{4} \\[2mm]
    3\left(a_3c_3c_2^2 + a_4c_4c_3^2\right) &= \tfrac{1}{4} \\[2mm]
    6a_4c_4c_3c_2 &= \tfrac{1}{4}
\end{align*}}

これらを解いて { a,\,c } の値を求めなければいけませんが、かなり面倒なので、上記の関係式を満たす特別な { a,\,c } の値を与えるだけでヨシとしましょう*1

  { \displaystyle\begin{align*}
    a_1 &= a_4 = \tfrac{1}{6}, &
    a_2 &= a_3 = \tfrac{1}{3} \\[2mm]
    c_2 &= c_3 = \tfrac{1}{2}, &
    c_4 &= 1
\end{align*}}

よって、この場合の { y_{n+1} } の表式は以下のようになります:

4次のルンゲ=クッタ法

  { \displaystyle\begin{align*}
    &y_{n+1} = y_n + \tfrac{1}{6}h\left(k_1 + 2k_2 + 2k_3 + k_4\right) \\[4mm]
    &\begin{cases}
        k_1 = f(x_n,\,y_n) \\[2mm]
        k_2 = f(x_n + \tfrac{1}{2}h,\,y_n + \tfrac{1}{2}hk_1) \\[2mm]
        k_3 = f(x_n + \tfrac{1}{2}h,\,y_n + \tfrac{1}{2}hk_2) \\[2mm]
        k_4 = f(x_n + h,\,y_n + hk_3)
    \end{cases}
\end{align*}}

まぁ、4次のルンゲ=クッタ法の値をそのまま持ってきただけですが。 これらの値は、きちんと上記の条件式を満たします。 以上で 4次のルンゲ=クッタ法の導出完了ってことで。

{ n } 個の条件式も頑張れば導けそうな気がするので、興味ある人はやってみてくださいな。

Javaによるアルゴリズム事典

Javaによるアルゴリズム事典

計算物理学入門

計算物理学入門

  • 作者: ハーベイゴールド,ジャントボチニク,Harvey Gould,Jan Tobochnik,鈴木増雄,石川正勝,溜渕継博,宮島佐介
  • 出版社/メーカー: ピアソンエデュケーション
  • 発売日: 2000/12
  • メディア: 単行本
  • 購入: 1人 クリック: 28回
  • この商品を含むブログ (45件) を見る

*1: {c_4=1,\,c_2=c_3 } くらいの条件を課せば、これらの値を得ることが出来ます。