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 となっている)。
mac mini (yosemite) 不具合：Google Japanese Inputが原因か？
mac miniの調子がすこぶる悪い。
(1)セーフモードで起動＝起動時に「Shift」を押しておく
>>> 起動に2時間近くかかったりする。
しばらくは普通に動くのだが、
(2)普通に起動
>>> 起動にx時間かかる(x>1)
(3)Apple Diagnostics=起動時に「D」を押しておく
>>> 「異常なし」
＝＝＝
しばらくセーフモードで騙し騙し使っていたのだが、あるとき、
- Cannot start conversion engine. Please restart your computer.
- Failed to get current config values
などのエラーメッセージがでてきたので、これをもとに調べると、どうもGoogle Japanese Input(Google日本語入力)が悪さをしている模様。
>>> Google Japanese Input をアンインストール。
どうなることか。
>>> 治った。[追記。]
python loadtext
data_0= numpy.loadtxt(inputFileName_0, delimiter=' ', comments='#', usecols=range(13))
データの行の末尾が半角スペースになってしまう場合は、読み込むカラムの個数を与える(この場合は13カラム)。
軌道計算アニメーション
ネット上の例にあったもの
0.1sec x 10 = 1secごとにパラパラアニメ
$ convert -delay 10 ./Pic/*.png test.gif
-----
[基本]
自分の場合
$ convert -delay 1 1with_bar/*.png 1snapshot_with_bar.gif
例えば
1with_bar
というフォルダに複数のpngファイルがあったとき、それらをファイルの名前順につなげてくれるようである。
- google chromeなどのブラウザで見ると動いてくれる。
- Macのpreviewで見ると、combineされた(100枚の画像なら100ページの)画像として見れる。
gifはストップ・再生ができないので（ブラウザの更新ボタンで最初からなら再生できる）、snapshotを見る上ではpreviewがいい。
-----
[少し応用]
$ convert -delay 20 -resize 50% omega_m*Gyr/0*/fig*Lz/*png reduced50_delay20.gif
などとすると、もとのpngファイルのサイズを50%にしたもののgifファイルを作れる。gifが重すぎるときに便利。
-delay 20 ::: 20/100 secごと
-resize 50% ::: 画像サイズが50%
Python いろいろなパッケージのインストール
(1) fitsio
$pip install fitsio
(2) esutil [ https://github.com/esheldon/esutil ]
$ pip install esutil
(3) isodist [ ]
$ python setup.py build
$ python setup.py install
(4) periodictable-1.4.1 [ https://pypi.python.org/pypi/periodictable ]
$ sudo python setup.py install
lambda式
Pythonの「ラムダ文」が難しいｗw
import scipy.optimize as op
def lnprob(theta, ell_data, bee_data):
lp = lnprior(theta)
if not np.isfinite(lp):
return -np.inf
return lp + lnlike(theta, ell_data, bee_data)
theta=[5.05, 250., -200., 1.1] #探索するパラメータの初期値。
minus_lnprob = lambda theta,*args: -lnprob(theta,*args)
result = op.minimize(minus_lnprob, theta, args=(ell_data, bee_data), method='TNC', bounds=bnds, options={'maxiter': 10000,'xtol': 1e-6, 'ftol':1e-6})
print(result)
print(result["x"])
試行錯誤の末、上記に辿り着き、なんとかうまく行ったが、理解できないｗｗ
minus_lnprob = lambda theta,*args: -lnprob(theta,*args)
の部分は
def minus_lnprob(theta, ell_data, bee_data):
return -lnprob(theta, ell_data, bee_data)
と同値のようである。
