« Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その5) ピタゴラスの三体問題を計算する。rtolとatolを設定しないと無茶苦茶になる。 | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その7) 三体問題の周期解いろいろ(1) 8の字解 »

2020年7月 2日 (木)

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

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

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

 

 

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

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

Pythagoras_20200701222601

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

ソースはこんな感じで。

%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()

« Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その5) ピタゴラスの三体問題を計算する。rtolとatolを設定しないと無茶苦茶になる。 | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その7) 三体問題の周期解いろいろ(1) 8の字解 »

学問・資格」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

« Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その5) ピタゴラスの三体問題を計算する。rtolとatolを設定しないと無茶苦茶になる。 | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その7) 三体問題の周期解いろいろ(1) 8の字解 »

最近の記事

最近のコメント

2021年9月
      1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30    
フォト
無料ブログはココログ