Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う（その6） ピタゴラスの三体問題を計算してFuncAminationでGIFアニメにしてみる。

昨日は最終的な軌跡を描いただけでしたが、せっかくなんでアニメにしたい。

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

いろいろ調べて、なんとなくこんな感じでできるんじゃないかということでやってみた（おそらくめっちゃ無駄なことをやっているが、それはおいおい修正するとして、、、）

これ。クリックで動きます。

これならいろいろアニメにして遊べそう。他の三体問題もやってみる。

ソースはこんな感じで。

%matplotlib notebook
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter

INTERVAL = 0.0005 # [s] interval time
FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval between frames
FPS = 1000/FRAME_INTERVAL # frames per second

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

def update(i):
    plt.cla() #現在描写されているグラフを消去
    plt.plot(sol[0:i,0],sol[0:i,2],c='Red',lw=1)
    plt.plot(sol[0:i,4],sol[0:i,6],c='Blue',lw=1)
    plt.plot(sol[0:i,8],sol[0:i,10],c='Green',lw=1)
    plt.scatter(sol[i,0],sol[i,2],c='Red')
    plt.scatter(sol[i,4],sol[i,6],c='Blue')
    plt.scatter(sol[i,8],sol[i,10],c='Green')
    plt.xlim(-8,8)
    plt.ylim(-8,8)
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.title("Pythagoras 3 body problem by DOP853")
    plt.grid()

# Plot
fig = plt.figure(figsize=(6,6))
ani = FuncAnimation(fig, update, frames=range(0,len(t),20), interval=FRAME_INTERVAL, blit=True)
plt.show()
ani.save("pythagoras.gif", writer = 'pillow',fps=FPS)

plt.show()