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

倭算数理研究所

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

2.2 オイラーのアルゴリズム

ニュートン力学以来、物理法則の多くは微分方程式の形で記述されてきました。 したがって、多くの場合の物理シミュレーションが「微分方程式に従って物理系の状態を遷移させていき、各ステップで物理量を計算して出力する」というのものになっています*1。 このようなシミュレーションに関しては、微分方程式が物理系から引きはがされ、微分方程式が与えられたものとして一般的なシミュレーション手法がいろいろ考案されています。

今回は、それらの手法の中で最も基本となるオイラー法 (Euler's method) を見ていきます。

オイラー

微分方程式が規定する物理系をシミュレートする場合の問題設定は、だいたい以下のようになります:

初期条件 { (x,\,y) = (x_0,\,y_0) } の下で次の微分方程式を解け:

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

物理系の状態は { (x,\,y) } の値の組で指定されるのですが、シミュレーションでは物理系の状態遷移は

  { \displaystyle
\begin{align*}
     (x_0,\,y_0) \rightarrow (x_1,\,y_1) \rightarrow (x_2,\,y_2) \rightarrow \cdots \rightarrow (x_n,\,y_n) \rightarrow (x_{n+1},\, y_{n+1}) \rightarrow \cdots
\end{align*}
}

のようにステップ・バイ・ステップで求めていきます。 初期条件は出発点となる状態を指定し、微分方程式はある状態から次の状態への遷移(のアルゴリズム)を規定します。

{ x } の遷移
{ x } の遷移は通常決まった値だけ増加させるだけです。 ここではその増加幅を { h } としましょう。 このとき

{ x } の遷移:

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

となります。

{ y } の遷移
{ y } の遷移は

から計算します。 微分方程式が規定する物理系に適用されるシミュレーション手法は、(オイラー法に限らず)多くの場合 { x } の増加幅が一定で、{ y } の増加幅をどのように計算するかが変わっているだけです。

オイラーでは、微分を「差分」に置き換えます:

  { \displaystyle
\begin{align*}
    \frac{dy}{dx} \longrightarrow \frac{y_{n+1} - y_n}{h}
\end{align*}
}

この置き換えを微分方程式に施して、それを { y_{n+1} } について解くと以下の式を得ます:

{ y } の増加:

  { \displaystyle
\begin{align*}
    y_{n+1} = y_n + h f(x_n,\,y_n)
\end{align*}
}

結果
結局、オイラー法では以下のような状態遷移を行います:

オイラー法の状態遷移:

  { \displaystyle
\begin{align*}
    \begin{cases}
        x_{n+1} &= x_n + h \\[2mm]
        y_{n+1} &= y_n + h f(x_n,\,y_n)
    \end{cases}
\end{align*}
}

計算物理学入門

計算物理学入門

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

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

*1:まさに、Groops の対象とするようなシミュレーション :-)