仮想計算機構

IT業界と無縁な派遣社員のブログ

PythonによるN体シミュレーション

はじめに

N体シミュレーションを行いました。
初期値等は以下の情報を参考にしています。

上記サイトとの違いは、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

f:id:riverta1992:20210606163149p:plain

t=0.5

f:id:riverta1992:20210606163251p:plain

t=1.0

f:id:riverta1992:20210606163323p:plain

t=2.0

f:id:riverta1992:20210606163720p:plain

アニメーション

下記を参考に png ファイルから動画を作成しました。

結果は以下の通りです。