はじめに
N体シミュレーションを行いました。
初期値等は以下の情報を参考にしています。
- Create Your Own N-body Simulation (With Python) | by Philip Mocz | The Startup | Medium
- nbody-python/nbody.py at master · pmocz/nbody-python · GitHub
上記サイトとの違いは、6N個の微分方程式を立ててルンゲクッタ法で計算した点です。
実装
環境
実行結果を Sage でプロットしたいので Python を使います。
プログラム
下記プログラム ( main.sage ) を SageShell に load して実行します。
import numpy as np G = 1.0 N = 100 xmax = 2 xmin = -2 h = 0.01 def f_vx(n,m,X,r): i = n-3*N S = 0 for j in range(N): if i==j: continue S += m[j]*(X[j]-X[i])/r[j][i] return S*G def f_vy(n,m,X,r): i = n-4*N S = 0 for j in range(N): if i==j: continue S += m[j]*(X[j+N]-X[i+N])/r[j][i] return S*G def f_vz(n,m,X,r): i = n-5*N S = 0 for j in range(N): if i==j: continue S += m[j]*(X[j+2*N]-X[i+2*N])/r[j][i] return S*G def f(n,m,X,r): if 0<=n and n<3*N: return X[n+3*N] elif 3*N<=n and n<4*N: return f_vx(n,m,X,r) elif 4*N<=n and n<5*N: return f_vy(n,m,X,r) else: return f_vz(n,m,X,r) # 位置の初期化 x = np.random.randn(N) y = np.random.randn(N) z = np.random.randn(N) # 速度の初期化 vx = np.random.randn(N) vy = np.random.randn(N) vz = np.random.randn(N) vx -= np.mean(vx) vy -= np.mean(vy) vz -= np.mean(vz) # 質量 m = [20/N]*N X = np.concatenate([x, y, z, vx, vy, vz], 0) nstep = 1000 g_list = [] r = [[0]*N for i in range(N)] for step in range(nstep): for i in range(N-1): for j in range(i+1,N): r_ij = np.sqrt(0.01+np.power(X[j]-X[i],2)+np.power(X[j+N]-X[i+N],2)+np.power(X[j+2*N]-X[i+2*N],2)) r[i][j] = r_ij**3 r[j][i] = r[i][j] Ka = np.array([h*f(i,m,X,r) for i in range(6*N)]) Kb = np.array([h*f(i,m,X+Ka/2,r) for i in range(6*N)]) Kc = np.array([h*f(i,m,X+Kb/2,r) for i in range(6*N)]) Kd = np.array([h*f(i,m,X+Kc,r) for i in range(6*N)]) new_x = [X[i]+(Ka[i]+2*Kb[i]+2*Kc[i]+Kd[i])/6 for i in range(6*N)] X = np.array(new_x) if step > 20: g_list.pop(0) g = point([(X[i],X[i+N]) for i in range(N)],xmin=xmin,xmax=xmax,ymin=xmin,ymax=xmax,color="blue",size=10) if step>0: g += sum([line([(g_list[i][n][0],g_list[i][n][1]) for i in range(len(g_list))],color="gray",alpha=0.5) for n in range(N)]) g += text("t={:.2f}".format(step*h),(-1.5,1.8),fontsize=20,color="black",bounding_box={'boxstyle':'round', 'fc':'w'}) g.save("image/img_{:06}.png".format(step)) g_list.append([(X[i],X[i+N]) for i in range(N)])
実行結果
今回は時間刻みが 0.01 で for 文を 1000 回まわしているので、実行後の image ディレクトリ配下に1000枚の画像が保存されている状態です。いくつかピックアップしてみます。
t=0
t=0.5
t=1.0
t=2.0