« Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その4) ローレンツ96モデル 36変数の方程式  | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その6) ピタゴラスの三体問題を計算してFuncAminationでGIFアニメにしてみる。 »

2020年7月 1日 (水)

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)すると割と簡単に計算できる。

こんな感じ。

Pythagoras

 

ソースはこちら。

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

« Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その4) ローレンツ96モデル 36変数の方程式  | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その6) ピタゴラスの三体問題を計算してFuncAminationでGIFアニメにしてみる。 »

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

コメント

コメントを書く

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

« Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その4) ローレンツ96モデル 36変数の方程式  | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その6) ピタゴラスの三体問題を計算してFuncAminationでGIFアニメにしてみる。 »

最近の記事

2021年5月
            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 31          
フォト
無料ブログはココログ