Google ColabのJulia言語でDifferentialEquationsパッケージを使ってピタゴラスの三体問題を35段14次ルンゲクッタFeagen法で計算し、Colab上でGIFアニメにしてみた。dtminを小さくしないと途中で止まってしまうことにハマった…

今回はピタゴラスの三体問題。質量mが3,4,5の三体が辺の長さが3,4,5の直角三角形の頂点におかれている場合にどのような動きをするか、というもの。

以前、Python、Visual C#, Unityで計算したがJuliaでもやってみよう。

まずはめちゃくちゃ時間がかかるが

using Pkg

Pkg.add("DifferentialEquations")

としてコードはこんな感じ。使うのはもちろん35段14次ルンゲクッタFeagen法で、今回はアダプティブありでabstol=reltol=1e-30とした。BigFloatも使用。

using DifferentialEquations

using Plots

function threebody!(du, u, p, t)

    m1 = p[1]

    m2 = p[2]

    m3 = p[3]

    qx1 = u[1]

    vx1 = u[2]

    qy1 = u[3]

    vy1 = u[4]

    qx2 = u[5]

    vx2 = u[6]

    qy2 = u[7]

    vy2 = u[8]

    qx3 = u[9]

    vx3 = u[10]

    qy3 = u[11]

    vy3 = u[12]

    Ra = sqrt((qx2 - qx1)^2.0 + (qy2 - qy1)^2.0)^3.0

    Rb = sqrt((qx3 - qx1)^2.0 + (qy3 - qy1)^2.0)^3.0

    du[1] = vx1

    du[2] = m2 * (qx2 - qx1) / Ra + m3 * (qx3 - qx1) / Rb

    du[3] = vy1

    du[4] = m2 * (qy2 - qy1) / Ra + m3 * (qy3 - qy1) / Rb

    Ra = sqrt((qx1 - qx2)^2.0 + (qy1 - qy2)^2.0)^3.0

    Rb = sqrt((qx3 - qx2)^2.0 + (qy3 - qy2)^2.0)^3.0

    du[5] = vx2

    du[6] = m1 * (qx1 - qx2) / Ra + m3 * (qx3 - qx2) / Rb

    du[7] = vy2

    du[8] = m1 * (qy1 - qy2) / Ra + m3 * (qy3 - qy2) / Rb

    Ra = sqrt((qx1 - qx3)^2.0 + (qy1 - qy3)^2.0)^3.0

    Rb = sqrt((qx2 - qx3)^2.0 + (qy2 - qy3)^2.0)^3.0

    du[9] = vx3

    du[10] = m1 * (qx1 - qx3) / Ra + m2 * (qx2 - qx3) / Rb

    du[11] = vy3

    du[12] = m1 * (qy1 - qy3) / Ra + m2 * (qy2 - qy3) / Rb

    nothing

end

tspan = (0, 100)

p = [3.0, 4.0, 5.0]

u0 = [big(1.0), big(0.0), big(3.0), big(0.0), big(-2.0), big(0.0), big(-1.0), big(0.0), big(1.0), big(0.0), big(-1.0), big(0.0)]

prob = ODEProblem(threebody!, u0, tspan, p)

sol = solve(prob, Feagin14(), dtmin = 1.0e-20, abstol = big(1.0e-30), reltol = big(1.0e-30))

これが何回やってもまともな精度がでないまま、途中でステップが小さくなり過ぎだと出て止まる。止まらないように精度を落とすと全然まともな波形が出ない。どのアルゴリズム使っても一緒。

なぜ？と思ったらデフォルトの最小刻み幅設定が～10^（-6)で、それはこの問題に対しては大きすぎていたから。ここにめっちゃくちゃハマった。dtminを設定するとうまく動いた。

Juliaは結果を簡単にGIFアニメにできるのでやってみた。ただ膨大なデータになったので500個飛ばしにしないといつまでたっても処理が終わらなかった。結果はこれ。