三体問題のプログラムを改良して高精度化と高速化。
もともと任意の個数の物体について運動を計算できるようにしていたので、ラグランジュポイントのシミュレーションもしてみました。
意外に初期値設定が難しい。円運動をさせるのに、地球と月の重力の両方を考慮した初期速度を与える必要がある。やってみないとわからんもんです。
L1,L2は理論どおり動径方向に不安定。月の裏側のジオン公国はまだましですが、月の表側のサイド5は速攻で地球の重力に惹かれて行きました。ルウム戦役の残骸なんかすぐにどっかにいってしまいそうです。
それはそうと、今日は朝から6時からマウンテンバイクで遠出。芹生峠を越えて黒谷降りて、林道でさらに北に行って、大きな杉を見てきました。
帰りは京北に出て、周山街道で帰宅。2時半ぐらい。結構疲れました。
# coding: utf-8 import numpy as np from scene import * import math class Simulation(Scene): def setup(self): self.disp_interval=200 self.delta=100 self.background_color = 'black' # self.setLagrange() self.setThreeBody() self.mass_=np.repeat(self.mass0.reshape((-1,1)),3,axis=1) self.mass=np.zeros((self.N,self.N-1,3)) for i in range(self.N-1): self.mass[:,i,:]=np.roll(self.mass_,i+1,axis=0) self.cells=[] for i in range(self.N): self.cells.append(SpriteNode('shp:Circle')) for i,cell in enumerate(self.cells): cell.position=(self.pos[i,0]*self.size.w/10,self.pos[i,1]*self.size.h/10) size=np.clip(self.mass0[i]**0.33*15,10,100)#質量の1/3乗に比例した大きさ。あまり小さいと見えないのでクリップ cell.size=(size,size) cell.color=self.colors[i] self.add_child(cell) self.disp=np.zeros((self.N,self.N-1,3)) def setThreeBody(self): self.G=0.1#重力定数 self.N=3 self.mass0=np.ones(self.N) self.pos=np.random.rand(self.N,3)*3 self.vel=np.random.rand(self.N,3)/3 self.colors=["white","blue","red","pink","gray"] def setLagrange(self): self.G=0.05#重力定数、適当 self.N=7 self.mass0=np.ones(self.N) self.mass0[0]=81.3#地球質量 self.mass0[1]=1#月質量 self.mass0[2]=0.0000001 #コロニー質量、適当 self.mass0[3]=0.0000001 #コロニー質量、適当 self.mass0[4]=0.0000001 #コロニー質量、適当 self.mass0[5]=0.0000001 #コロニー質量、適当 self.mass0[6]=0.0000001 #コロニー質量、適当 R=4#月軌道半径、適当 self.pos=np.zeros((self.N,3)) self.vel=np.zeros((self.N,3)) self.pos[:,:]=0 self.pos[1,1]=R self.pos[2,1]=-R self.vel[:,:]=0 self.vel[1,0]=np.sqrt(self.G*self.mass0[0]/R)#円運動になるように設定 self.vel[2,0]=-np.sqrt(self.G*self.mass0[0]/R+self.G*self.mass0[1]/R/2)#円運動になるように設定 t=math.pi*2/3 rotation=np.array([[np.cos(t),np.sin(t),0],[-np.sin(t),np.cos(t),0],[0,0,1]]) self.pos[3]=np.dot(rotation,self.pos[2]) self.vel[3]=np.dot(rotation,self.vel[2]) t=-math.pi*2/3 rotation=np.array([[np.cos(t),np.sin(t),0],[-np.sin(t),np.cos(t),0],[0,0,1]]) self.pos[4]=np.dot(rotation,self.pos[2]) self.vel[4]=np.dot(rotation,self.vel[2]) r=R*(self.mass0[1]/3/self.mass0[0])**(1/3) self.pos[5,1]=R-r self.vel[5,0]=np.sqrt(self.G*self.mass0[0]/(R-r)-0.99*self.G*self.mass0[1]/(r**2/(R-r)))#円運動 self.pos[6,1]=R+r self.vel[6,0]=np.sqrt(self.G*self.mass0[0]/(R+r)+self.G*self.mass0[1]/(r**2/(R+r)))#円運動 self.vel[2,0]=-np.sqrt(self.G*self.mass0[0]/R+self.G*self.mass0[1]/4/R)#円運動になるように設定 self.colors=["skyblue","white","gray","gray","gray","gray","gray"] def update(self): dt=self.dt/self.delta for j in range(self.disp_interval): for i in range(self.N-1):#np.rollを使うとメモリーコピーが生じて遅くなるので、参照で対応 self.disp[0:i+1,i,:]=self.pos[-(i+1)::,:]-self.pos[0:i+1,:] self.disp[i+1::,i,:]=self.pos[0:-(i+1),:]-self.pos[i+1::,:] r2=(np.sum(self.disp**2,axis=2)+1e-8).reshape(self.N,self.N-1,1) f=np.sum(self.disp/r2**1.5*self.mass,axis=1)*self.G self.vel+=f*dt self.pos+=self.vel*dt self.pos=self.pos-np.sum(self.pos*self.mass_,axis=0)/np.sum(self.mass_,axis=0)+5 for i,cell in enumerate(self.cells): cell.position=(self.pos[i,0]*self.size.w/10,self.pos[i,1]*self.size.h/10) if __name__ == '__main__': run(Simulation(), PORTRAIT, show_fps=False)
