7 数値積分法

ここまでは、「相互作用をいつ、どうやって計算するか」という話をしてきた。 しかし、実際に数値解を求めるには求まった重力加速度を使って各粒子の軌道 を積分する必要がある。以下、その方法について簡単に述べる。

7.1 The Euler method

これは知らない人はいないと思いたいが、ある微分方程式の初期値問題

$$dx/dt = f(x,t)$$ (2)

$$x(0) = x_0$$ (3)

があるときに、時刻 $t+\Delta t$ での近似解 $x_{i+1}$ を、 時刻 $t$ での近似解$x_i$を使って

$$x_{i+1} = x_i + f(x_i,t)$$ (4)

とする方法である。局所打ち切り誤差は $O(\Delta t^2)$ であり、したがっ て大域誤差は $O(\Delta t)$ということになる。プログラミングは簡単である が、誤差があまりに大きくほとんど実用にはならない。

7.2 The leapfrog method

これはセパラブルなハミルトン系、すなわちハミルトニアンが運動量の関数と 位置の関数に分かれるものに対する特別な方法である。通常の、運動量に対す るハミルトニアンが2次形式である場合には、

$$v_{i+1/2} = v_{i-1/2} + \Delta t a(x_i)$$ (5)

$$x_{i+1} = x_{i} + \Delta t v_{i+1/2}$$ (6)

と書くのが普通である。これでは速度と位置がずれた時間でしか定義されない が、出発用公式として

$$v_{1/2} = v + \Delta t a(x_0)/2$$ (7)

を使い、さらに終了用公式として

$$v_{i} = v_{i-1/2} + \Delta t a(x_i)/2$$ (8)

を使うことで最初と最後を合わせることが出来る。この形は、実は

$$x_{i+1} = x_{i} + \Delta t v_{i} + \Delta t^2 a(x_i)/2$$ (9)

$$v_{i+1} = v_{i} + \Delta t [a(x_i)+a(x_{i+1})]/2$$ (10)

と数学的に等価である（証明してみること）また、以下のようにも書ける

$$v_{i+1/2} = v_i + \Delta t a(x_i)/2$$ (11)

$$x_{i+1} = x_{i} + \Delta t v_{i+1/2} + \Delta t^2 a(x_i)/2$$ (12)

$$v_{i+1} = v_{i+1/2} + \Delta ta(x_{i+1})/2$$ (13)

さらにまた、速度を消去して$x_{i-1}$, $x_{i}$, $x_{i+1}$ の関係式の形で かいてあることもあるかもしれない。なお、すべて違う名前がついていたりす るが、結果は丸め誤差を別にすれば完全に同じである。

この公式は局所誤差が $O(\Delta t^3)$、大域誤差が$O(\Delta t^2)$である。

さて、局所誤差という観点からはこれは決して良い公式というわけではないが、 現実には特に無衝突系の計算ではこの公式は非常に広く使われている。衝突系 の計算で、独立時間刻みを使うような場合にはこの方法は使われないが、それ 以外の非常に多くの問題はこの方法で解かれているといってよい。

これは、別にもっとよい方法を知らないからとかではなく、実はこの方法がい くつかの意味で非常に良い方法であるからである。いくつかの意味とは、例え ば力学平衡の系でエネルギーや角運動量が非常によく保存するということであ る。これらの保存量については多くの場合に誤差がある程度以上増えない。

この、ある程度以上誤差が増えないというのは目覚しい性質である。通常の方 法では、エネルギーの誤差は時間に比例して増えていく。従って、長時間計算 をしようとすればそれだけ正確な計算をする必要がある。ところが、エネルギー の誤差は溜っていかないのならば、かならずしも精度を上げる必要はないとも 考えられる。

もちろん、エネルギーが保存していればそれだけで計算が正しいということに はならない。が、理論的には、いくつかの重要な結果が得られている。まず、

1. leapfrog は symplectic method のもっとも簡単なものの一つである。
2. leapfrog は symmetric method のもっとも簡単なものの一つである。

Symplectic method については、以下のようなことが知られている

1. symplectic method は、すくなくともある種のハミルトニアンに対して 使った場合に、それに近い別のハミルトン系に対する厳密解を与えることがあ る。
2. 周期解を持つハミルトン系に対して使った場合に、どんな量でも誤差が 最悪で時間に比例してしか増えない。
3. 時間刻を変えると上のようなことは成り立たなくなる

これに対し、 symmetric method については以下のようなことが知られている

1. 周期解を持つ時間対称な系に対して使った場合に、どんな量でも誤差が 最悪で時間に比例してしか増えない。
2. 時間刻を変えてもうまくいくようにすることも出来る。

というわけで、 leap frog が、それよりも局所的には誤差が小さいルンゲ・ クッタとかに比べて良い結果を与えるのは当然のことである。

7.3 The Aarseth method

さて、独立時間刻みを使った場合には、話はそう簡単にはいかない。というの は、 symplectic method の良い点も symmetric method の良い点も失われる からである。このため、通常はより高次の方法が用いられる。まあ、 cosmological な計算なんかでは独立時間刻みで2次の公式を使っている例もな いわけではないが、結果を信用できるかどうかには注意したほうがよいかもし れない。

高次の方法としては、伝統的には Aarseth が60年代に開発した方法が用いら れてきた。これは、いわゆる線形多段法 linear multistep method に基づく ものである。具体的な実現についてはここでは省くが、基本的に、過去数ステッ プの加速度から補外多項式を作り、それを積分して位置と速度を更新する。独 立時間刻みを使うに当たっては、他の粒子の位置は補外多項式で計算する。

7.4 The Hermite method

Aarseth method は広く使われてきたが、結構面倒であるし高次にしていくと 計算量が増え、また安定性が急激に悪くなる。もう少し簡単な方法はないかと いうことで考えられたのが Hermite method である。これでは、加速度の他に その時間導関数まで直接計算し、それから補外、補間多項式を作る。これは最 近かなり広く使われるようになってきた。