Excelでシンプレクティック公式(8次)を計算してみる。

エネルギーが保存するハミルトン系の微分方程式の数値解法であるシンプレティック法だけれど、高次のものの数表がいまいち見つからない。吉田春夫さんの論文　"Construction of higher order symplectic integrators", Physics Letters A,Volume 150, number 5,6,7(1990)から引っ張ってくることにした。いろいろ書かれている中で、8次の公式のSolution Cを使うことにした。論文ではwiで書いてあるのを使いやすいようにciとdiに書き直す。

i	ci	di
0	0	0
1	0.314515325	0.62903065
2	0.999190057	1.369349464
3	0.152381158	-1.064587148
4	0.299385476	1.6633581
5	-0.007805591	-1.678969283
6	-1.61921866	-1.559468038
7	-0.623838613	0.311790812
8	0.985390848	1.658990885
9	0.985390848	0.311790812
10	-0.623838613	-1.559468038
11	-1.61921866	-1.678969283
12	-0.007805591	1.6633581
13	0.299385476	-1.064587148
14	0.152381158	1.369349464
15	0.999190057	0.62903065
16	0.314515325	0

（間違っていても責任はとりません。。。）ここで、

q_i=q_i-1 + dt c_i p_i-1/m

p_i=p_i-1 - dt d_i ∂U(q_i)/∂q

で、dtすすめるのに16ステップ使う（結構時間かかるな）。とりあえず単振動dp/dt=-qを計算する。

エネルギーもぴったり。結構すごいな。