倭算数理研究所

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

ルンゲ=クッタ法を導く (1) : 準備としての高次テイラー法

次回に、イマイチ導出が謎な『ルンゲ=クッタ (Runge-Kutta)法*1を導いてみます。 今回はそれを行うのに必要な準備として、高次テイラーを簡単に導出します*2

目次

  1. 準備としての高次テイラー
  2. 問題設定
  3. 1次の場合
  4. 2次の場合
  5. 3次の場合
  6. 4次の場合

微分方程式

問題設定はオイラー法を導出したときと同じです:

次の微分方程式を解け:

  { \displaystyle\begin{align*}
    \frac{dy}{dx} = f(x,\,y)
\end{align*}}

{ f(x, y) }{ C^\infty } 級の関数としましょう。 初期条件はここでは省略(ホントは大切ですけど)。 この微分方程式をシミュレーションで解く際に必要なのは、{ f(x,\,y),\,h } を用いて { (x_n,\, y_n) } の組から { (x_{n+1},\,y_{n+1}) } を計算することです。 { x_{n+1} }

  { \displaystyle\begin{align*}
    x_{n+1} = x_n + h
\end{align*}}

で与えられるので、メインとなるのは { y_n } から { y_{n+1} } を計算することになります。 今回導く高次テイラー法や次回から導くルンゲ=クッタ法は、この { y_{n+1} } を計算するアルゴリズムが異なるだけです({ x_{n+1} } の表式は同じ)。

高次 Taylor 法

問題の微分方程式の解が

  { \displaystyle\begin{align*}
    y = y(x)
\end{align*}}

であるとします。 このとき

  { \displaystyle\begin{align*}
    \begin{cases}
        y_n = y(x_n) \\
        y_{n+1} = y(x_{n+1}) = y(x_n + h)
    \end{cases}
\end{align*}}

となります。 ここで、{ y_{n+1} }{ x_n } の周りでテイラー展開すると

  { \displaystyle\begin{align*}
    y_{n+1}
        &= y(x_n+h) \\
        &= y(x_n) + h y'(x_n) + \frac{h^2}{2!} y''(x_n) + \frac{h^3}{3!} y'''(x_n) + \cdots \\
        &= y_n + h \sum_{m=0}^\infty \frac{h^m}{(m+1)!} y^{(m+1)}(x_n)
\end{align*}}

さらに問題の微分方程式 { y' = f(x,\,y) } より

  { \displaystyle\begin{align*}
    y^{(m+1)}(x) = \frac{d^m}{dx^m}f(x,\,y)
\end{align*}}

なので、結局次式を得ます:

  { \displaystyle\begin{align*}
    y_{n+1}
        &= y_n + hf(x_n,\,y_n) + \frac{h^2}{2!}\left\{\frac{df}{dx}(x_n,\,y_n)\right\}
            + \frac{h^3}{3!}\left\{\frac{d^2f}{dx^2}(x_n,\,y_n)\right\} + \cdots \\[2mm]
        &= y_n + h\sum_{m=0}^\infty \frac{h^m}{(m+1)!}\left\{\frac{d^m}{dx^m}(x_n,\,y_n)\right\}
\end{align*}}

ここで幾つか注意を。 まず

  { \displaystyle\begin{align*}
    \frac{d^mf}{dx^m}(x_n, \, y_n)
\end{align*}}

{ f(x,\,y) }{ x }{ m }微分した後、{ (x_n,\,y_n) } を代入するという意味です。 また、

  { \displaystyle\begin{align*}
    \frac{d^mf}{dx^m}
\end{align*}}

偏微分ではなく微分です。 つまり、{ x } の引数だけでなく、{ y } を通した { x } への依存性に関しても微分を行います:

  { \displaystyle\begin{align*}
    \frac{d}{dx}f(x,\,y) = \frac{\partial}{\partial x} f(x,\,y) + \left\{\frac{\partial}{\partial y} f(x,\,y)\right\}\frac{dy}{dx}
\end{align*}}

ちなみに、高次テイラー法で { h } の2次以上を切り捨てると オイラー法 (Euler method)になります。

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

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

計算物理学入門

計算物理学入門

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

*1:『ルンゲークッタ 法』とは常微分方程式をシミュレーションで解く際に使う技法。

*2:Javaによるアルゴリズム事典』とこちらのサイトを参考にしました。