삼체(The Three-Body Problem) 정말 불가능할까?

삼체의 궤적

다음은 중국의 작가 류츠신이 지은 SF 장편 소설, 삼체(The Three-Body Problem)를 읽으면서 궁금한 사항을 모아 정리한 것이다.  

삼체문제(Three-body problem)란? 

삼체문제는 아이작 뉴턴이 역학을 집대성한 자신의 저서, 프린키피아(Principia)를 통해 처음 언급하였다.

태양과 지구 간에는 상호 중력이 작용한다. 이의 영향으로 지구는 태양의 주위를 도는 공전을 하게 된다. 이와 같은 원리로 달은 지구를 중심으로 공전을 한다. 이와 같이 2개의 물체 간에 중력이 상호 어떤식으로 작용하는지, 그 결과로 어떠한 궤도 움직임을 보일 것인지는 쉽게 예측할 수 다. 이를 2개의 물체간의 문제, 즉 이체문제(Two-body problem)라고한다. 

 

이제 물체를 하나더 추가하여, 3개의 물체를 가정해보자. 3개의 물체 간에 상호 중력이 작동한다면, 그들의 궤도 움직임을 어떠할까? 태양과 지구의 움직임처럼 쉽게 예측 가능할까? 

1980년 수학자 앙리 푸앵카레(Henri Poincaré)는 삼체 문제의 일반해를 구하는 것은 불가능하다는 사실을 증명하였다. 즉, 예측이 사실상 불가능하다는 것이다. 

 

실제로 삼체문제가 어떤 궤도를 만들어내는지 궁금하다면, 아래의 링크를 참고 : https://javalab.org/three_body_problem/

 

삼체문제를 파이썬 코드로 구현하기

삼체문제를 시각화하는 파이썬 코드를 다음과 같다.  이 코드를 실행하면 세 천체가 서로 중력에 의해 영향을 받으며 움직이는 궤적을 보여주는 그래프가 생성됩니다. 각 천체의 초기 위치는 점으로 표시되고, 시간이 지남에 따라 궤적이 그려진다.

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp

# Constants and initial conditions
G = 1  # Gravitational constant for simplicity

# Masses of the three bodies
m1, m2, m3 = 1, 1, 1

# Initial positions and velocities
r1 = np.array([0.97000436, -0.24308753])
v1 = np.array([0.466203685, 0.43236573])
r2 = np.array([-0.97000436, 0.24308753])
v2 = np.array([0.466203685, 0.43236573])
r3 = np.array([0, 0])
v3 = np.array([-2*0.466203685, -2*0.43236573])

# State vector: [r1, r2, r3, v1, v2, v3]
y0 = np.hstack([r1, r2, r3, v1, v2, v3])

def deriv(t, y):
    r1 = y[0:2]
    r2 = y[2:4]
    r3 = y[4:6]
    v1 = y[6:8]
    v2 = y[8:10]
    v3 = y[10:12]
    
    dr1dt = v1
    dr2dt = v2
    dr3dt = v3
    
    dv1dt = G * m2 * (r2 - r1) / np.linalg.norm(r2 - r1)**3 + G * m3 * (r3 - r1) / np.linalg.norm(r3 - r1)**3
    dv2dt = G * m1 * (r1 - r2) / np.linalg.norm(r1 - r2)**3 + G * m3 * (r3 - r2) / np.linalg.norm(r3 - r2)**3
    dv3dt = G * m1 * (r1 - r3) / np.linalg.norm(r1 - r3)**3 + G * m2 * (r2 - r3) / np.linalg.norm(r2 - r3)**3
    
    return np.hstack([dr1dt, dr2dt, dr3dt, dv1dt, dv2dt, dv3dt])

# Integrate the equations of motion
t_span = (0, 10)
t_eval = np.linspace(t_span[0], t_span[1], 1000)
sol = solve_ivp(deriv, t_span, y0, t_eval=t_eval, rtol=1e-9, atol=1e-9)

# Extract the positions
r1_sol = sol.y[0:2, :]
r2_sol = sol.y[2:4, :]
r3_sol = sol.y[4:6, :]

# Plot the trajectories
plt.figure(figsize=(10, 10))
plt.plot(r1_sol[0], r1_sol[1], label='Body 1')
plt.plot(r2_sol[0], r2_sol[1], label='Body 2')
plt.plot(r3_sol[0], r3_sol[1], label='Body 3')

# Initial positions
plt.scatter(*r1, color='blue')
plt.scatter(*r2, color='orange')
plt.scatter(*r3, color='green')

# Legends and labels
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('Three Body Problem')
plt.grid()
plt.show()