三体問題

三体問題がどのくらい予測不能かを見てみたくなったのでプログラムを書いてシミュレーションしてみた。

三次元空間で計算しているけど表示は2次元に射影。
確かに動きの先がよめない。太陽が三つある世界の人は大変だ。

import numpy as np
from scene import *

class Simulation(Scene):
    def setup(self):
        self.counter=0
        self.background_color = 'black'
        self.N=3
        self.cells=[]
        for i in range(self.N):
            self.cells.append(SpriteNode('shp:Circle'))       
        self.mass=1
        self.delta=1
        self.pos=np.random.randn(self.N,3)
        self.vel=np.random.randn(self.N,3)/2
        self.disp=np.zeros((self.N,self.N-1,3))
        colors=["white","#FFCC33","#FF3333"]
        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)
            cell.color=colors[i]
            self.add_child(cell)
    def update(self):
        dt=self.dt/self.delta #self.dtはSceneクラスに最初から定義されている前のアップデートからの時間(秒)
        for i in range(self.N-1):
            self.disp[:,i,:]=np.roll(self.pos,i+1,axis=0)-self.pos
        r2=(np.sum(np.abs(self.disp),axis=2)+1e-8).reshape(self.N,self.N-1,1)
        f=np.sum(self.disp/r2**1.5,axis=1)
        self.vel+=f*dt/self.mass
        self.pos+=self.vel*dt
        if self.counter%100==0:
            self.pos=self.pos-np.average(self.pos,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=True)