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

倭算数理研究所

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

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

シミュレーション技法

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

満たすべき条件式

満たすべき条件式は

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

ただし

  { \displaystyle\begin{align*}
    k_1 &= f(x_n,\,y_n) \\
    k_2 &= f(x_n+c_2h,\,y_n+c_2hk_1) \\
    k_3 &= f(x_n+c_3h,\,y_n+c_3hk_2)
\end{align*}}

(1) 式:{ k } の0階微分の条件式

s = 2 の場合と同様にして

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

(2) 式:k の1階微分の条件式

これも s = 2 の場合と同様にして

  { \displaystyle\begin{align*}
    a_2 c_2 +a_3 c_3= \tfrac{1}{2}
\end{align*}}

(3) 式:{ k } の2階微分の条件式

{ k_1 }微分
{ k_1 }{ h } に依存しないので簡単。

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

{ k_i }微分
{ k_i\; (i = 2,\,3) } の場合

  { \displaystyle\begin{align*}
    \left[\frac{d^2k_i}{dh^2}\right]_{h=0}
        = c_i^2\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) + 2c_i\frac{\partial f}{\partial y}\left[\frac{dk_{i-1}}{dh}\right]_{h=0}
\end{align*}}

ここで

  { \displaystyle\begin{align*}
    \left[\frac{dk_1}{dh}\right]_{h=0} &= 0, & \left[\frac{dk_2}{dh}\right]_{h=0} &= c_2\frac{df}{dx}
\end{align*}}

より次式を得ます:

  { \displaystyle\begin{align*}
    \left[\frac{d^2k_2}{dh^2}\right]_{h=0}
        &= c_2^2\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) \\[4mm]
    \left[\frac{d^2k_3}{dh^2}\right]_{h=0}
        &= c_3^2\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) + 2c_3c_2\frac{\partial f}{\partial y}\frac{df}{dx}
\end{align*}}

{ f(x,\,y) } の2階微分
{ f(x,\,y) } の2階微分は次のようになります:

  { \displaystyle\begin{align*}
    \frac{d^2f}{dx^2}
        &= \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}
\end{align*}}

得られる関係式
上記で得られた関係を (3) 式に代入して整理すると

  { \displaystyle\begin{align*}
    &\left(a_2c_2^2 + a_3c_3^2\right)\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\} + 2a_3c_3c_2\left\{\frac{\partial f}{\partial y}\frac{df}{dx}\right\} \\[2mm]
    &\qquad = \frac{1}{3}\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*}}

したがって次式が得られます:

  { \displaystyle\begin{align*}
    a_2c_2^2 + a_3c_3^2 &= \tfrac{1}{3}, &
    2a_3c_3c_2 &= \tfrac{1}{3}
\end{align*}}

{ a,\,c } が満たすべき関係式

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

  { \displaystyle\begin{align*}
    a_1 + a_2 +a_3 &= 1 \\
    a_2 c_2 +a_3 c_3 &= \tfrac{1}{2} \\
    a_2c_2^2 + a_3c_3^2 &= \tfrac{1}{3} \\
    2a_3c_3c_2 &= \tfrac{1}{3}
\end{align*}}

となります。

{ a,\,c } の値

上記の関係式は、文字が5つで式が4つなので、完全には文字の値が決まりません。 なので、文字1つを適当に決めてやる必要があります。 ただ、任意の値でいいわけではないようで、例えば { c_3=1 } を課すと { c_2 } が実数で求まりません。 ちょっと面倒な(そんなに難しくありませんが)解析をしてやると、{ c_2 } が実数として存在するためには

  { \displaystyle\begin{align*}
     c_3 \le \tfrac{3}{4}
\end{align*}}

という条件が満たされている必要があります。 これが充分条件というわけではありませんが、1つの文字の値を決めるのには適度な目安でしょう。

{ c_3=\frac{3}{4} } の場合
この場合、各変数の値は

  { \displaystyle\begin{align*}
    (a_1,\,a_2,\,a_3,\,c_2,\,c_3) = \left(\tfrac{2}{9},\,\tfrac{1}{3},\,\tfrac{4}{9},\,\tfrac{1}{2},\,\tfrac{3}{4}\right)
\end{align*}}

と定まります。 なんとなく、{ a,\,c } それぞれの組の数字が綺麗に並んでる気も。 このとき { y_{n+1} } は以下のようになります:

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

{ c_3 = \frac{2}{3} } の場合
このとき、各変数の値は次の2通りに決まります:

  { \displaystyle\begin{align*}
     (a_1,\,a_2,\,a_3,\,c_2,\,c_3)
        = \left(\tfrac{1}{4},\,0,\,\tfrac{3}{4},\,\tfrac{1}{3},\,\tfrac{2}{3}\right),\,
            \left(\tfrac{1}{4},\,\tfrac{3}{8},\,\tfrac{3}{8},\,\tfrac{2}{3},\,\tfrac{2}{3}\right)
\end{align*}}

このとき、それぞれの場合の { y_{n+1} }

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

もしくは

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

となります。 1組目の { y_{n+1} } には直接 { k_2 } が含まれてませんが、{ k_3 } を計算する際に必要になります。 一応、少しは計算量減りまする。

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

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

計算物理学入門

計算物理学入門

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