資料室を作ってみた。

morinoseikatu2003-12-09

どうかな?
メモ帳代わりに使おうと思ってんだけど。
フラクタルとカオスに関する本
複雑系のカオス的シナリオ」
 津田一郎
「カオス的脳観」
「解析概論」
笠原晧司 著
 新微分方程式対話 [新版]
 日評数学選書
 1995, 日本評論社
 ISBN 4-535-60119-4

Euler法
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法の右辺を とする代わりに、いわば のように半ステップだけ進めた値を用いるのである(具体的なことは 常微分方程式のプログラム例を参考にしてほしい)。この方法は「改良Euler法」あるいは「中点法」と呼ばれ、 Euler法に比べるとはるかに精度が高い。たとえば、 として、 Euler法と改良Euler法の誤差を比べてみよう。

Euler法
1ステップあたりの誤差は 程度だが、 たとえば まで 計算しようとすると、1000ステップ分の誤差が蓄積することになる。 蓄積した誤差の大きさは 0.0001 の 1000倍なので 0.1 となり、 たとえば に対して 10%の誤差を生じる。
改良Euler法
1ステップあたりの誤差は 程度なので、 同じく までの計算で 1000ステップ分の誤差が蓄積したとしても、 0.1%の誤差を生じるに過ぎない。
なお、Heun法という計算法でも、改良Euler法と同じくらいの精度が得られる。くわしくは参考文献を見ること。


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法は次のように書ける:
反復 {
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 に反比例して増えるから、 結果として累積する誤差は 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