Julia で楽してピタゴラス三体問題

はじめに

先日、周期的境界条件を使って多体問題のコードを書いた。そこでは、４次のRunge-Kuttaを使っていたため、運動は周期的境界条件を抜きにしても、適当だった。そこで、正確な運動を知るためのシミュレーションを行った。
JuliaにはDifferentialEquationsという非常に便利なパッケージがある。今回はそれを使う。

ピタゴラスの三体問題

質量比と距離がともに$3:4:5$になるように配置した三体問題である。以下が結果。

ソースコード

# 質点の質量
m1, m2, m3 = 3.0, 4.0, 5.0

# 重力定数（この問題のために適当に選択）
G = 1.0

# 微分方程式の定義
function three_body!(du, u, p, t)
    # u = [x1, y1, x2, y2, x3, y3, vx1, vy1, vx2, vy2, vx3, vy3]
    x1, y1, x2, y2, x3, y3, vx1, vy1, vx2, vy2, vx3, vy3 = u
    r12 = ((x1 - x2)^2 + (y1 - y2)^2)^1.5
    r13 = ((x1 - x3)^2 + (y1 - y3)^2)^1.5
    r23 = ((x2 - x3)^2 + (y2 - y3)^2)^1.5

    # dr/dt = v
    du[1] = vx1
    du[2] = vy1
    du[3] = vx2
    du[4] = vy2
    du[5] = vx3
    du[6] = vy3
    
    # dv/dt = F/m
    du[7] = G * m2 * (x2 - x1) / r12 + G * m3 * (x3 - x1) / r13
    du[8] = G * m2 * (y2 - y1) / r12 + G * m3 * (y3 - y1) / r13

    du[9] = G * m1 * (x1 - x2) / r12 + G * m3 * (x3 - x2) / r23
    du[10] = G * m1 * (y1 - y2) / r12 + G * m3 * (y3 - y2) / r23

    du[11] = G * m1 * (x1 - x3) / r13 + G * m2 * (x2 - x3) / r23
    du[12] = G * m1 * (y1 - y3) / r13 + G * m2 * (y2 - y3) / r23
end

# 初期条件の設定
u0 = [1.0, 3.0, -2.0, -1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 
# [x1, y1, x2, y2, x3, y3, vx1, vy1, vx2, vy2, vx3, vy3]

# 時間範囲の設定
tspan = (0.0, 40.0)

# 問題の定義
prob = ODEProblem(three_body!, u0, tspan)

# 解の計算（Runge-Kutta-Fehlberg法）
sol = solve(prob, AutoTsit5(Rosenbrock23()))

結果のプロット

plot(sol, idxs=(1, 2), label="Body 1")
plot!(sol, idxs=(3, 4), label="Body 2")
plot!(sol, idxs=(5, 6), label="Body 3", xlabel="x", ylabel="y", xlims=(-5,5),ylims=(-5,5))
# 結果の保存
savefig("three_body.png")

animation

anim=@animate for i in 1:length(sol.t)
    plot(xlim=(-5, 5), ylim=(-5, 5), legend=:bottomright, xlabel="x", ylabel="y")
    scatter!([sol.u[i][1]], [sol.u[i][2]], color=:blue)
    scatter!([sol.u[i][3]], [sol.u[i][4]], color=:red)
    scatter!([sol.u[i][5]], [sol.u[i][6]], color=:green)
end every 1

gif(anim, "three_body.gif", fps = 15)

解説

ここでは、Runge-Kutta-Fehlberg 法の刻み幅の自動制御を行い、ヤバそうなところ（発散してしまうところ）は時間刻みを小さくしそうでないところはざっくりした刻みにしている。

まとめ

このように、DifferentialEquationsパッケージを使えば、自分で実装するのが面倒なことでも上記のように書くだけでやってくれる。