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

倭算数理研究所

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

ルンゲ=クッタ法を導く (2) : 問題設定

ルンゲ=クッタ (Runge-Kutta) 法を導くシリーズ(目次)。 今回は「{ s } 次のルンゲ=クッタ (Runge-Kutta) 法 」の表式(ここで扱う形)の紹介と、前回の「高次テイラー法」を踏まえたルンゲ=クッタ法の導出のための問題設定をします。 具体的な導出は次回以降(目次)。

{ s } 次のルンゲ=クッタ法

{ s } 次のルンゲ=クッタ法は以下の式を用いて { y_n } から { y_{n+1} } を計算します:

  { \displaystyle\begin{align*}
    &y_{n+1} = y_n + h \sum_{i=1}^s a_i k_i \\[4mm]
    &\begin{cases}
        k_1 = f(x_n,\,y_n) \\[2mm]
        k_i = f(x_n + c_i h,\, y_n + c_i h k_{i-1}) & (2 \le i \le s)
   \end{cases}
\end{align*}}

ここで、{ a_i\; (1 \le i \le s) }{ c_i \; (2 le i \le s) } は定数です*1。 中間値の定理のイメージ。 これらの定数は何でもいいわけではなく、必要な精度で前回やった高次 Taylor 法と等しくなるようにとります。 どのようにこれらの定数を決めればそうなるかを導くのが ルンゲ=クッタ法の導出です。 ただし、これらの定数の値は一意に定まらないので、通常

  { \displaystyle\begin{align*}
    &0 < c_i < 1 & (2 \le i \le s-1) \\
    &c_s = 1
\end{align*}}

という制限もかけます。

ちなみに、高次テイラー法と同じ精度になるなら敢えてルンゲ=クッタ法を使う有難味はどこにあるのかというと、ルンゲ=クッタ法は関数 { f(x,\,y) } だけを用いて評価でき、特にその導関数を使う必要がないというところです。 シミュレーションでは { f(x,\,y) } の計算アルゴリズムが与えられただけではその導関数を計算するのが簡単ではないので、この利点はかなり大きくなります。 高次導関数が必要な 高次 Taylor 法は実際にはほとんど使われてないかと思います。

ちょっとした制限
上記の記述は、{ f(x,\,y) }{ y } に対する増加の表式が、wikipedia:ルンゲ=クッタ法 に掲載されている「{ s } 次の ルンゲ=クッタ法」とは異なっています:

  • { y } の増加の式が { k_i } の線形結合でなくなった
  • { y } の増加の式に { x } の増加の式の係数 { c_i } と同じものを使っている

表式としてはいくらか一般性を失っていますが、通常の(4次の)ルンゲ=クッタ法を導く際には、特に不都合になることはありません。

問題設定
ルンゲ=クッタ法の導出」としてやるべき事は、上記 { y_{n+1} }{ h }{ s } 次まで「高次テイラー法」の { y_{n+1} }

  { \displaystyle\begin{align*}
    y_{n+1} = y_n + h \sum_{m=0}^\infty \frac{h^m}{(m+1)!}\left\{\frac{d^mf}{dx^m}(x_n,\,y_n)\right\}
\end{align*}}

に一致するように { a_i \; (1 \le i \le s),\;c_i\; (2 \le i \le s) } を決定する、です。 「{ s } 次のルンゲ=クッタ法」の { y_{n+1} } の表式を用いて書き換えると

  { \displaystyle\begin{align*}
  &\sum_{i=1}^s a_i k_i
     = \sum_{m=0}^\infty \frac{h^m}{(m+1)!}\left\{\frac{d^mf}{dx^m}(x_n,\,y_n)\right\} \\[4mm]
  &\begin{cases}
    k_1 = f(x_n,\,y_n) \\[2mm]
    k_i = f(x_n + c_i h,\, y_n + c_i h k_{i-1}) & (2 \le i \le s)
  \end{cases}
\end{align*}}


{ h }{ s-1 } 次まで成り立つように { a_i\; (1 \le i \le s),\;c_i\; (2 \le i \le s) } を決定する

となります。 ここで、両辺を { h }{ m }微分して { h = 0 } とおくと(これは大雑把に言って 「{ h } で冪展開したときの係数」を比べてます。 後の補足参照)

{ s } 次のルンゲ=クッタ法を導くためには、{ 0 \le m \le s-1 } を満たす整数 { m } に対して


  { \displaystyle\begin{align*}
    &\sum_{i = 1}^s a_i \left[\frac{d^m k_i}{dh^m}\right]_{h=0} = \frac{1}{m+1}\left\{\frac{d^mf}{dx^m}(x_n,\,y_n)\right\} \\[4mm]
    &\begin{cases}
        k_1 = f(x_n,\,y_n) \\[2mm]
        k_i = f(x_n + c_i h,\, y_n + c_i h k_{i-1}) & (2 \le i \le s)
    \end{cases}
\end{align*}}


が成り立つように { a_i \; (1 \le i \le s),\;c_i\; (2 \le i \le s) } を決定する

右辺の { m } についての和がなくなり、代わりに { s } 個({ m } の取り得る値の数)の条件式になりました。

補足

{ x } の関数 { f(x) } の冪展開が

  { \displaystyle\begin{align*}
    f(x) = \sum_{n=0}^\infty a_n x^n
\end{align*}}

で与えられるとき、両辺を { m }微分すると

  { \displaystyle\begin{align*}
    \frac{d^m f}{dx^m}(x) = \sum_{n=0}^\infty \frac{n!}{(n-m)!}a_n x^{n-m}
\end{align*}}

更に両辺で { x = 0 } とおくと

  { \displaystyle\begin{align*}
    \frac{d^m f}{dx^m}(0) = m! \, a_m
\end{align*}}

よって、{ f(x) }{ x } で冪展開したときの { x^n } の係数は以下で与えられます:

  { \displaystyle\begin{align*}
    a_n = \frac{1}{n!}\frac{d^n f}{dx^n}(0)
\end{align*}}

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

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

計算物理学入門

計算物理学入門

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

*1:Wikipedia の記事とは文字の用い方が異なっています。 Wikipedia の記事からここで使っている表式に移るには、次の置き換えを行います:[tex:{ a_{i,i-1\} \rightarrow c_i,\,b_i \rightarrow a_i