Google ColabのJulia言語でシンプレクティック8次の公式（吉田春夫さんの）を使って三体問題の周期解（8の字解、クラスI.A.1のbutterfly I）を計算。

今回はシンプレクティック8次で計算。

こちらの論文から引っ張ってきた。

"Construction of higher order symplectic integrators", Physics Letters A,Volume 150, number 5,6,7(1990)

コードはこんな感じで。

using Plots

#力の計算（x方向）

function fx(k)

    G = 1.0

    f = 0.0

    for i in 1:n

        if i != k

            R = sqrt((qx[i] - qx[k]) ^ 2 + (qy[i] - qy[k]) ^ 2)

            f = f + G * m[k] * m[i] * (qx[i] - qx[k]) / (R ^ 3)

        end

    end

    return f

end

#力の計算（y方向）

function fy(k)

    G = 1.0

    f = 0.0

    for i in 1:n

        if i != k

            R = sqrt((qx[i] - qx[k]) ^ 2 + (qy[i] - qy[k]) ^ 2)

            f = f + G * m[k] * m[i] * (qy[i] - qy[k]) / (R ^ 3)

        end

    end

    return f

end

# 初期設定

n = 3 #星の個数

m = zeros(Float64, n)

qx = zeros(Float64, n)

qy = zeros(Float64, n)

px = zeros(Float64, n)

py = zeros(Float64, n)

dt = 3.0 / 10000.0

#8の字解

#param = [ 0.97000436,-0.5 * 0.93240737,-0.24308753,-0.5 * 0.86473146, -0.97000436,-0.5 * 0.93240737,0.24308753,-0.5 * 0.86473146,0,0.93240737, 0,0.86473146]

#mparam = [1.0, 1.0, 1.0]

#tmax = 10000 #繰り返し回数

#13個見つかった解のうち、クラスI.A.1のbutterfly I

vx1 = 0.30689

vy1 = 0.12551

param =[-1.0, vx1, 0.0, vy1, 1.0, vx1, 0.0, vy1, 0.0 , -2.0*vx1, 0.0, -2.0*vy1]

mparam = [1.0, 1.0, 1.0]

tmax = 30000 #繰り返し回数

m[1] = mparam[1]

m[2] = mparam[2]

m[3] = mparam[3]

qx[1] = param[1]

qy[1] = param[3]

px[1] = param[2]

py[1] = param[4]

qx[2] = param[5]

qy[2] = param[7]

px[2] = param[6]

py[2] = param[8]

qx[3] = param[9]

qy[3] = param[11]

px[3] = param[10]

py[3] = param[12]

#表示用

x = zeros(Float64, n, tmax)

y = zeros(Float64, n, tmax)

x[1, 1] = qx[1]

y[1, 1] = qy[1]

x[2, 1] = qx[2]

y[2, 1] = qy[2]

x[3, 1] = qx[3]

y[3, 1] = qy[3]

#　8次のシンプレクティック積分公式の係数

w = zeros(Float64, 8)

d = zeros(Float64, 16)

c = zeros(Float64, 16)

w[1] = 0.311790812418427

w[2] = -1.55946803821447

w[3] = -1.6789692825964

w[4] = 1.66335809963315

w[5] = -1.06458714789183

w[6] = 1.36934946416871

w[7] = 0.629030650210433

w[8] = 1.0 - 2.0 * (w[1] + w[2] + w[3] + w[4] + w[5] + w[6] + w[7])

d[1] = w[7]

d[15] = w[7]

d[2] = w[6]

d[14] = w[6]

d[3] = w[5]

d[13] = w[5]

d[4] = w[4]

d[12] = w[4]

d[5] = w[3]

d[11] = w[3]

d[6] = w[2]

d[10] = w[2]

d[7] = w[1]

d[9] = w[1]

d[8] = w[8]

d[16] = 0.0

c[1] = 0.5 * w[7]

c[16] = 0.5 * w[7]

c[2] = 0.5 * (w[7] + w[6])

c[15] = 0.5 * (w[7] + w[6])

c[3] = 0.5 * (w[6] + w[5])

c[14] = 0.5 * (w[6] + w[5])

c[4] = 0.5 * (w[5] + w[4])

c[13] = 0.5 * (w[5] + w[4])

c[5] = 0.5 * (w[4] + w[3])

c[12] = 0.5 * (w[4] + w[3])

c[6] = 0.5 * (w[3] + w[2])

c[11] = 0.5 * (w[3] + w[2])

c[7] = 0.5 * (w[2] + w[1])

c[10] = 0.5 * (w[2] + w[1])

c[8] = 0.5 * (w[1] + w[8])

c[9] = 0.5 * (w[1] + w[8])

for i in 2:tmax

    for k in 1:16

        for j in 1:n

            qx[j] = qx[j] + dt * c[k] * px[j] / m[j]

            qy[j] = qy[j] + dt * c[k] * py[j] / m[j]

        end

        for j in 1:n

            px[j] = px[j] + dt * d[k] * fx(j)

            py[j] = py[j] + dt * d[k] * fy(j)

        end

    end

    for j in 1:n

        x[j, i] = qx[j]

        y[j, i] = qy[j]

    end

end

アニメーションにするには

@gif for i in 1:200:tmax

    plot(x[1,1:i],y[1,1:i], xlim=(-1.2, 1.2), ylim=(-0.5,0.5), label="", xlabel="x", ylabel="y", size=(800,500))

    plot!(x[2,1:i],y[2,1:i], label="")

    plot!(x[3,1:i],y[3,1:i], label="")

end

とする。

まずは有名な8の字解。

次はこちらの論文に乗っていた13個見つかった解のうち、クラスI.A.1のbutterfly I

Three Classes of Newtonian Three-Body Planar Periodic Orbits

他のもやろうと思ったが、刻み幅一定だと、いくら8次でも衝突が多いような解は発散してしまう…

ちょっといまいち。

だがDifferentialEquations.jlのアダプティブな刻み華で、誤差を小さくしたものなら計算できる。

クラスII.C.3a yin-yang IIならばこんな風に面白い動きが再現できる。