python - n체 문제(n-body problem)

이번 학기에 가장 인상적이었던 실습 내용입니다. 

실습이라기엔 사실 기말고사 마지막 문제 내용이었습니다... 

https://gfycat.com/ko/shorttermhighcapybara

n-body problem, n체 문제, 혹은 다체 문제는 질량을 가진 여러 개체들이 공간 안에서 만유인력을 주고받으며 발생하는 운동 궤적을 알아내는 것입니다. 

천체물리학에서 매우 중요한 문제입니다만, 개체가 3개 이상일 때 해석적인 궤적을 얻기는 매우 어렵거나 불가능하므로 일반적으로 수치적인 방법을 이용합니다. 

 

저는 각 질점이 동일한 질량을 가지고 회전운동은 고려하지 않는 상황을 가정하였습니다. 

 

---

import numpy as np
import matplotlib.pyplot as plt

def accel(x,y,z):
    r32 = (x*x + y*y + z*z)**1.5
    return mi*np.array([-x/r32, -y/r32, -z/r32])

def Energy(vx,vy,vz):
    return mi/2*sum(vx**2+vy**2+vz**2)

def W(x,y,z):
    val = 0
    x1 = x*1
    y1 = y*1
    z1 = z*1
    for i in range(n-1):
        x1 = np.roll(x1,1)
        y1 = np.roll(y1,1)
        z1 = np.roll(z1,1)
        r = ( (x1-x)**2 + (y1-y)**2 + (z1-z)**2 )**0.5
        val -= 1/r
    return sum(val*mi**2/2)

def r_dist(x):
    return x**0.5
    
n = 50
mi = 1/n

r1, r2, r3 = np.random.random(n), np.random.random(n), np.random.random(n)
r = r_dist(r1)
theta = np.arccos(2*r2-1)
phi = 2*np.pi*r3

x = r * np.sin(theta) * np.cos(phi)
y = r * r_dist(r1)* np.sin(theta) * np.sin(phi)
z = r * np.cos(theta)
m = np.ones(n)*mi

r4, r5 = np.random.random(n), np.random.random(n)
eta1 = np.arccos(2*r4-1)
eta2 = 2*np.pi*r5
vx = 0.5 * np.sin(eta1) * np.cos(eta2)
vy = 0.5 * np.sin(eta1) * np.sin(eta2)
vz = 0.5 * np.cos(eta1)

plt.figure( figsize=(6,6) )
plt.plot(x, y, c='deepskyblue', marker='o', \
         ms=10, mfc='paleturquoise', ls='None')

plt.xlabel('x', size=15)
plt.ylabel('y', size=15)
plt.title('fig.11', size=25)
plt.tight_layout()

t,tmax = 0,2
h = 0.001
hh = 0.5*h

tsol = [t*1]
xsol,ysol,zsol= [x*1],[y*1],[z*1]
Esol,Wsol = [Energy(vx,vy,vz)],[W(x,y,z)]

count = 0
while t < tmax :
    vhx,vhy,vhz = vx*1, vy*1, vz*1
    for i in range(1,n):# v+0.5
        x1 = np.roll(x, i)
        y1 = np.roll(y, i)
        z1 = np.roll(z, i)
        
        f = accel(x-x1, y-y1, z-z1)
        vhx += hh*f[0]
        vhy += hh*f[1]
        vhz += hh*f[2]
        
    x += h*vhx
    y += h*vhy
    z += h*vhz
    xsol.append(x*1)
    ysol.append(y*1)    
    zsol.append(z*1)
    vx,vy,vz = vhx*1, vhy*1, vhz*1

    for i in range(1,n):
        x1 = np.roll(x, i)
        y1 = np.roll(y, i)
        z1 = np.roll(z, i)
        
        f = accel(x-x1, y-y1, z-z1)
        vx += hh*f[0]
        vy += hh*f[1]
        vz += hh*f[2]
    t += h
    Esol.append( Energy(vx,vy,vz) )
    Wsol.append( W(x,y,z) )
    tsol.append(t)
    count += 1
    if count%1000==0 : print(count//1000)

xsol = np.array(xsol)
ysol = np.array(ysol)
tsol = np.array(tsol)
color = ['r','orange','y','g','b','midnightblue','purple','r','orange','y','g','b','midnightblue','purple']
k=0
plt.figure( figsize=(8,7) )
for i in range(k,k+7):
    plt.plot(xsol.T[i][tsol<0.5], ysol.T[i][tsol<0.5], c=color[i-k], ls='--', lw=1, alpha=0.1)
    plt.plot(xsol.T[i][abs(tsol-0.8)<0.3], ysol.T[i][abs(tsol-0.8)<0.3], c=color[i-k], ls='--', lw=1, alpha=0.3)
    plt.plot(xsol.T[i][abs(tsol-1.3)<0.2], ysol.T[i][abs(tsol-1.3)<0.2], c=color[i-k], ls='--', lw=1, alpha=0.5)
    plt.plot(xsol.T[i][tsol>=1.5], ysol.T[i][tsol>=1.5], c=color[i-k], ls='--', lw=1)
    plt.plot(xsol.T[i][-1], ysol.T[i][-1], c=color[i-k], marker='o', ms=15, mfc='w')

plt.xlim(-1,1)
plt.ylim(-1,1)
plt.tight_layout()

xsol = np.array(xsol)
ysol = np.array(ysol)
zsol = np.array(zsol)
tsol = np.array(tsol)

plt.figure(figsize=(10,10))
ax = plt.axes(projection='3d')
ax.set(facecolor='None')
k=3
for i in range(k,k+7):
    ax.plot3D(xsol.T[i][tsol<0.5], ysol.T[i][tsol<0.5], zsol.T[i][tsol<0.5], \
              c=color[i-k], ls='--', lw=1, alpha=0.1)
    ax.plot3D(xsol.T[i][abs(tsol-0.8)<0.3], ysol.T[i][abs(tsol-0.8)<0.3], zsol.T[i][abs(tsol-0.8)<0.3], \
              c=color[i-k], ls='--', lw=1.1, alpha=0.3)
    ax.plot3D(xsol.T[i][abs(tsol-1.3)<0.2], ysol.T[i][abs(tsol-1.3)<0.2], zsol.T[i][abs(tsol-1.3)<0.2],\
              c=color[i-k], ls='--', lw=1.2, alpha=0.5)
    ax.plot3D(xsol.T[i][tsol>=1.5], ysol.T[i][tsol>=1.5], zsol.T[i][tsol>=1.5], \
              c=color[i-k], ls='--', lw=1.4)
    ax.plot3D(xsol.T[i][-1], ysol.T[i][-1], zsol.T[i][-1], c=color[i-k], marker='o', ms=15, mfc='w')
plt.xlim(-1,1)
plt.ylim(-1,1)
plt.tight_layout()

---

가속도의 반영에는 계산이 빠른 leapfrog 알고리즘을 사용하였습니다. 

질점의 분포는 원점을 중심으로 반지름 1 이내의 영역으로 하였고, 중심으로부터 거리에 따른 밀도가 거리에 반비례하도록 설정하였습니다.  

초기 속도는 크기가 0.5이고 방향은 랜덤으로 설정하였습니다. 

 

묘사를 애니메이션으로 할 수 있으면 좋을텐데, 그래도 일단 투명도를 조절한 선으로 잘 표현할 수 있는 듯해서 좋습니다. 

 

 

n=4, 4체 문제

n=2, 2체 문제

n=7, 7체 문제