leapfrog 法による軌道計算

(1) leapfrog法

- 「Kick-Drift-Kick法」をimplementするために、簡単なKepler問題を解く。

- 重力場：質量1e11 Msunの質点

- 初期条件：(x,y)=(10, 0) kpc, (v_x, v_y)=(0, 100) km/s [eccentricity=0.77]

- 積分時間： 1 Gyr

- 積分timestep：(a) 1 Myr [1e3 steps] ; (b) 0.1 Myr [1e4 steps]

- (a) precession顕著。(b) precession無視可能。

- エネルギー誤差：(a) 15 %; (b) 0.5 % (during 1 Gyr)

- 計算時間：KDK (a) 7 ms; (b) 50 ms

参考までに、GSLの微分方程式のルーチンを使うと、20 ms (1 Myr / step)、 153 ms(0.1 Myr / step)。GSLでは、precessionはどちらのtimestepの場合も見られない。

>>> Kepler 問題を解く上では、GSLで何ら問題ない。ただし、多粒子になると、GSLのルーチンだけでは、計算時間が顕著に増えると予想される。「h=0.1 Myr, leapfrog」は計算精度、計算速度のバランス面で優れている。

>>> KDKの積分を行う関数を、「0.1 Myrの積分を10回行う」と記述すると、エネルギー誤差をおさえたまま9 ms 程度で計算できることが判明。つまり、KDKの計算時間は「関数の呼び出し回数」で決まる。 [追記2016/01/01]

(2) KDK法

timestep = h とすると、

v_(n+1/2) = v_n + a(x_n) h/2

x_(n+1) = x_n + v_(n+1/2) h

v_(n+1) = v_(n+1/2) + a(x_(n+1)) h/2

Springelの講義ノートでは、この２番目の式が間違っている(h が h/2 となっている)。