Euler法とその改良

2000-11-17 更新 / 2008-06-22 微更新


宿題の答えを探して当ページに来られた方へ

もし「○○をルンゲクッタで解け」とかいう宿題を出されていて、 その答えを探してここに来たのなら、以下の内容を読む前に、 宿題を出した先生の意図は何であるかを ちょっと考えてみてください。 もしかすると、先生が本当に求めているのは、 数値計算結果をレポートとして提出することではなく、 あなたが数値計算プログラムを自力で書けるようになること なのではないでしょうか?

そのへんをよくわきまえた上で、先生の意図に応えられるような形で 以下の内容を参考にされたら良いと思います。

参考: 「教えて! goo」に宿題を頼む前に考えてほしいこと


常微分方程式の数値解法でまず始めに習うのが Euler法である。 ただし、Euler法をそのまま用いるのは実用的でなく、 実際には より洗練された形を用いることが多い。 ここでは、改良Euler法(中点法)と、 その発展形であるRunge-Kutta法について説明する。

[「常微分方程式」に戻る] [「微分方程式」に戻る]


Euler法

常微分方程式

  dx/dt = v   (1)


 (x[j+1] - x[j])/dt = v[j]

と近似して解くのがEuler法である。 ただし v=v(x) または v=v(x,t) であり、また x[j] は 時刻 t = t[j] = t[0] + j dt における xの値を示す。 分母を払って
  x[j+1] = x[j] + v[j] dt   (2)

を得る。

例として、 運動方程式 ma = f

dx/dt = v,  dv/dt = a = f/m



として計算することを考えよう。 Euler法でプログラムを書くと
 反復 {
    a[j] = (x[j]に対応する加速度)
    x[j+1] = x[j] + v[j]*dt
    v[j+1] = v[j] + a[j]*dt
 }
のような感じになる。
※ 本当は配列を使わずに書けるが、 ここでは分かりやすいように配列を用いている。

改良Euler法

実際にEuler法で数値計算をしてみると、 そのままでは誤差が大きくて使いものにならない場合が多い。 たとえば、惑星の軌道が楕円軌道になるはずの場合でも、 もとの位置に戻ってこないで少しずれたりする。 もちろん dtをものすごく小さくすれば良いのだが、 こんどはステップ数が多くなって大変である。

そこで、次のように改良してみる:

 反復 {
    a[j] = (x[j]に対応する加速度)
    xh = x[j] + v[j]*(dt/2)
    vh = v[j] + a[j]*(dt/2)
    ah = (xhに対応する加速度)
    x[j+1] = x[j] + vh*dt
    v[j+1] = v[j] + ah*dt
 }
つまり、Euler法の右辺を x[j] とする代わりに、 いわば x[j+1/2]のように 半ステップだけ進めた値を用いるのである(具体的なことは 常微分方程式のプログラム例を 参考にしてほしい)。 この方法は「改良Euler法」 あるいは「中点法」「Heunの矩形公式」などと呼ばれ、 Euler法に比べるとはるかに精度が高い。 たとえば、 dt=0.01 として、 Euler法と改良Euler法の誤差を比べてみよう。
Euler法
1ステップあたりの誤差は dt^2 = 0.0001 程度だが、 たとえば t=10まで 計算しようとすると、1000ステップ分の誤差が蓄積することになる。 蓄積した誤差の大きさは 0.0001 の 1000倍なので 0.1 となり、 たとえば x=1.0 に対して 10%の誤差を生じる。
改良Euler法
1ステップあたりの誤差は dt^3 = 10^{-6} 程度なので、 同じく t=10までの計算で 1000ステップ分の誤差が蓄積したとしても、 0.1%の誤差を生じるに過ぎない。

なお「Heunの台形公式」という計算法(単に「Heun法」とも言う)でも、 改良Euler法と同じくらいの精度が得られる。 くわしくは参考文献を見ること。 また間違いプログラムの例(PDF)も参考になるかもしれない。

Runge-Kutta法

改良Euler法の基本的な考え方は、 中間的なステップを設けることによって 数値計算の精度を高められるという点にある。 この考えを発展させたのが Runge-Kutta法である。

例として、さっきの方程式(1)の数値解法を考える。 改良Euler法は、少し冗長になるが、次のようにも書ける:

 反復 {
    v1dt = v(x[j]          ) * dt
    v2dt = v(x[j] + v1dt/2 ) * dt
    x[j+1] = x[j] + v2dt
 }
また、Heun法(Heunの台形公式)は次のように書ける:
 反復 {
    v1dt = v(x[j]        ) * dt
    v2dt = v(x[j] + v1dt ) * dt
    x[j+1] = x[j] + (v1dt + v2dt)/2
 }
どちらの場合も、 もしいきなり x[j] = x[j] + v1dt とすれば Euler法と同じであるが、 そうせずに補助的に v2dt を求めることで精度を高めている。 これにより、累積誤差を dtの2乗に比例して小さくできる[※]。 改良Euler法やHeun法を総称して「2次Runge-Kutta法」という。
※ 1ステップあたりの誤差は dt^3程度であるが、 結果をだすのに必要な全ステップ数は dt に反比例して増えるから、 結果として累積する誤差は dtの2乗に比例する。

単にRunge-Kutta法と言った場合、次のような 4次Runge-Kutta法を指す場合が多い:

 反復 {
    v1dt = v(x[j]          ) * dt
    v2dt = v(x[j] + v1dt/2 ) * dt
    v3dt = v(x[j] + v2dt/2 ) * dt
    v4dt = v(x[j] + v3dt   ) * dt
    x[j+1] = x[j] + (v1dt + 2*v2dt + 2*v3dt + v4dt)/6
 }
4次Runge-Kutta法では、 一定時間のあいだに累積する誤差は dtの4乗に比例する。 2次と4次のどちらが良いかは、場合にもよるので一概には言えないが、 学部の授業における数値計算としては、 まずは2次の方法をきちんと習得することを勧めたい。

参考文献

W. H. Press ほか: Numerical Recipes in C [C言語による数値計算のレシピ]
技術評論社 1993
第15章: 常微分方程式の数値解法
山本 哲朗: 数値解析入門 (サイエンスライブラリ 現代数学への入門 14)
サイエンス社 1976
8-3. Runge-Kutta法 (pp.107-109)
一松 信: 数値解析 (新数学講座13)
朝倉書店 1982
Sec.29: Runge-Kutta型の公式 (pp.110-115)
三島 信彦: マイコン物理 (物理学 One Point 13)
共立出版 1981

追記

本格的な常微分の計算(MDとか天体計算とか)をしたい人は、 牧野先生「システム数理IV」の講義ノートを読んで勉強するのが良いでしょう。


[このページの冒頭に戻る] [常微分方程式] [微分方程式]