Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その5） ピタゴラスの三体問題を計算する。rtolとatolを設定しないと無茶苦茶になる。

非常に計算が難しいので有名なピタゴラスの三体問題。

質量mが3,4,5の三体が辺の長さが3,4,5の直角三角形の頂点におかれている場合にどのような動きをするか、というもの。(Burrauさんが問題提起して、Szebehelyさんが数値的に解決したということです）。

詳細はこちら。

http://www.ucolick.org/~laugh/oxide/projects/burrau.html

でDOP853で相対＆絶対誤差を非常に小さく(rtolとatol)すると割と簡単に計算できる。

こんな感じ。

ソースはこちら。

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt

def pythagoras(t, x, m1, m2, m3): #odeintのときとt,xの並びが逆
    f=np.zeros(12)

    qx1=x[0]
    vx1=x[1]
    qy1=x[2]
    vy1=x[3]

    qx2=x[4]
    vx2=x[5] 
    qy2=x[6]
    vy2=x[7]

    qx3=x[8]
    vx3=x[9]
    qy3=x[10]
    vy3=x[11]

    Ra=np.sqrt((qx2-qx1)**2 + (qy2-qy1)**2)**3
    Rb=np.sqrt((qx3-qx1)**2 + (qy3-qy1)**2)**3
    f[0]=vx1
    f[1]=m2*(qx2-qx1)/Ra + m3*(qx3-qx1)/Rb
    f[2]=vy1
    f[3]=m2*(qy2-qy1)/Ra + m3*(qy3-qy1)/Rb

    Ra=np.sqrt((qx1-qx2)**2 + (qy1-qy2)**2)**3
    Rb=np.sqrt((qx3-qx2)**2 + (qy3-qy2)**2)**3
    f[4]=vx2
    f[5]=m1*(qx1-qx2)/Ra + m3*(qx3-qx2)/Rb
    f[6]=vy2
    f[7]=m1*(qy1-qy2)/Ra + m3*(qy3-qy2)/Rb

    Ra=np.sqrt((qx1-qx3)**2 + (qy1-qy3)**2)**3
    Rb=np.sqrt((qx2-qx3)**2 + (qy2-qy3)**2)**3
    f[8]=vx3
    f[9]=m1*(qx1-qx3)/Ra + m2*(qx2-qx3)/Rb
    f[10]=vy3
    f[11]=m1*(qy1-qy3)/Ra + m2*(qy2-qy3)/Rb

    return f

t0=0
tmax = 100
N=10000

x0=[1,0,3,0,-2,0,-1,0,1,0,-1,0]
solver=ode(pythagoras)
solver.set_integrator('dop853',atol=1.E-12,rtol=1.E-12)
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆
solver.set_f_params(3,4,5)

t=np.linspace(0, tmax, N)
sol= np.empty((N, 12))
sol[0] = x0

k=1
while solver.successful() and solver.t < tmax:
    solver.integrate(t[k])
    sol[k] = solver.y
    k+= 1


# Plot
fig = plt.figure(figsize=(12,12))
plt.rcParams["font.size"] = 18
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Pythagoras 3 body problem by DOP853")
plt.grid()
plt.plot(sol[:,0],sol[:,2])
plt.plot(sol[:,4],sol[:,6])
plt.plot(sol[:,8],sol[:,10])
plt.xlim(-6,6)
plt.ylim(-6,6)

plt.show()