파이썬으로 천체물리학·AI: Rebound를 활용해 N-바디 시뮬레이션 비밀을 풀다

출처: Dev.to

우주는 끊임없는 시계처럼 작동하며, 모든 것을 지배하는 단일하고 보편적인 힘은 중력이다. 이진 별의 밀접한 춤에서부터 은하의 장대한 나선팔까지, 모든 운동은 거대한 몸체들 사이의 쌍대 중력’attraction에 의해 결정된다.

뉴턴 역학의 원리는 아름답게 단순하지만, 많은 상호작용하는 물체들의 운동을 모델링하는 작업(즉, N‑Body Problem)은 단순한 물리 방정식에서 복잡한 수학적·컴퓨팅 난제로 빠르게 전환된다. 이 가이드는 정적인 천체물리 데이터 분석에서 동적·예측 가능한 중력 시스템 모델링으로 초점을 이동한다. 우리는 기술 통계에서 예측 운동학으로 전환하며, REBOUND 시뮬레이션 패키지와 같은 특수화된 고정밀 도구가 필요하다.

이론적 난제: 중력이 왜 어려운가

전문화된 도구가 왜 필요한지를 이해하려면 중력 시스템이 자체적으로 가지고 있는 불안정성을 먼저 알아야 한다.

세 몸체의 혼돈

N=2일 때(예를 들어 지구가 태양을 도는 경우) 시스템은 분석적으로 풀릴 수 있다. 운동 방정식은 케플러 원형 곡선으로 설명되는 정확한 폐곡선 형태의 해를 제공한다.

하지만 제3의 몸체를 도입하면(N=3) 분석적 확실성은 사라진다. 시스템은 초기 조건에 민감해지며, 이는 결정론적 혼돈의 특징이다. 시작 위치나 속도에 미세한 변화가 장기적으로 극단적인 결과 차이를 초래한다. 따라서 복잡한 행성계의 진화를 예측하기 위해서는 수치 통합만 가능하다.

컴퓨팅 스케일링 위기:

두 번째 난제는 원시적인 연산 능력이다. N개의 몸체가 있다면 각 몸이 다른 모든 몸체와 상호작용한다. 고유한 쌍 대의 상호작용 수는 이차적으로 증가한다:

(P = \frac{N(N-1)}{2})

이 O((N^2)) 복잡성은 몸체 수를 두 배로 늘리면 계산 시간이 4배가 된다는 것을 의미한다. 은하 규모 시뮬레이션에서는 Barnes‑Hut과 같은 기술이 계산량을 줄일 수 있지만, 행성계의 장기적 안정성을 모델링하는 데 필요한 정밀도를 희생한다.

불안정 위기: 일반 적분기가 왜 실패하는가

완벽한 힘 계산을 하더라도 일반 수치 적분법(런게‑쿠타 등)은 시스템의 기본 물리 법칙을 준수하지 않아 천체물리학에서는 실패한다. 중력 시스템은 하밀토니계(Hamiltonian)이며, 전체 에너지 보존과 각운동량 보존이 필수적이다.

해결책은 symplectic 적분기이다. 이러한 알고리즘은 하밀토니 역학의 기하학적 구조를 보존하도록 설계되어, 수십억 단계에 걸친 오류가 누적되는 것을 방지한다. 장기적 안정성 연구에는 필수적이며 거부할 수 없다.

REBOUND 소개: 금표 표준

이러한 난관을 고려할 때, 일반 파이썬 라이브러리(scipy.integrate.solve_ivp)는 충분하지 않다. REBOUND(Recursive Bound)라는 강력한 C 기반 시뮬레이션 프레임워크가 robust한 Python 인터페이스를 제공한다.

REBOUND은 중력 역학을 위해 설계되었으며, 다음을 제공한다:

• 전문화된 적분기: WHFast( symplectic, 장거리 안정성용)와 IAS15 (적응형, 고정밀 근접 만남용) 포함.
• 계층적 시간 단계: 몸체마다 다른 시간 단계를 효율적으로 관리한다.
• 데이터 생성: AI 모델을 훈련시켜 시스템 안정성을 예측하는 데 완벽한 도구이다.

Python으로 지구 궤도 시뮬레이션 튜토리얼

코드를 살펴보겠습니다. 가장 간단한 안정적인 N‑체계인 지구가 태양을 도는 시스템을 설정하고, REBOUND를 이용해 1년 동안 통합하여 궤도 안정성을 검증하겠습니다.

선행 조건

라이브러리를 설치합니다:

pip install rebound numpy

Enter fullscreen mode

Exit fullscreen mode

파이썬 코드

아래 스크립트는 시뮬레이션을 초기화하고, 태양과 지구를 추가한 뒤 1년 동안 통합하여 궤도의 안정성을 확인하는 과정을 보여줍니다.

import rebound
import numpy as np
import sys

# 시뮬레이션을 조용히 그리고 효율적으로 실행합니다
rebound.set_status(sys.stdout, color=False)

# --- 1. 시뮬레이션 초기화 및 설정 ---
sim = rebound.Simulation()
sim.units = ('AU', 'Msun', 'year')   # 천체역학 표준 단위
sim.integrator = "IAS15"             # 고정밀 적응형 적분기
sim.dt = 0.001                       # 초기 타임스텝

# --- 2. 입자 추가 ---
# 입자 0: 태양 (질량 = 1.0 태양 질량)
sim.add(m=1.0)

# 입자 1: 지구 (질량 ≈ 3e-6 태양 질량)
# 케플러 원소로 정의됨: 반경축 = 1.0 AU, 이심률 = 0.0167
sim.add(m=3.003e-6, a=1.0, e=0.0167)

# --- 3. 안정화 ---
# 시스템을 중점(COM) 프레임으로 이동하여 드리프트를 방지합니다
sim.move_to_com()

# --- 4. 통합 루프 ---
T_end = 1.0       # 1년 동안 시뮬레이션
N_steps = 100     # 100개 데이터 포인트 기록
times = np.linspace(0, T_end, N_steps)
print(f"시뮬레이션 시작: {T_end} 년(들) ...")
for i, t in enumerate(times):
    sim.integrate(t)  # 포지션을 저장할 수도 있지만, 끝에 안정성을 검증하겠습니다

# --- 5. 결과 검증 및 출력 ---
print("\n--- 시뮬레이션 결과 ---")
print(f"최종 시간 도달: {sim.t:.6f} 년")

# 지구(입자 1)의 최종 상태 확인
final_earth = sim.particles[1]
print("\n지구 최종 상태:")
print(f"위치 (x, y, z): ({final_earth.x:.6f}, {final_earth.y:.6f}, {final_earth.z:.6f}) AU")
print(f"속도 (vx, vy, vz): ({final_earth.vx:.6f}, {final_earth.vy:.6f}, {final_earth.vz:.6f}) AU/년")

# 궤도 파라미터 검증 (초기값과 가깝도록)
print(f"\n최종 반경축 (a): {final_earth.a:.6f} AU")
print(f"최종 이심률 (e): {final_earth.e:.6f}")

# 에너지 보존 검증 (symplectic 적분기에게 중요)
sim.integrate(T_end)  # 정확한 종료 시점을 보장
e_start = -0.000000000000000  # 이 시스템의 이론적 에너지
e_end = sim.calculate_energy()
print(f"\n에너지 보존 검증:")
print(f"종료 시 총 에너지: {e_end:.12f}")

Enter fullscreen mode

Exit fullscreen mode

상세 코드 설명

1. 초기화 및 환경 설정

sim = rebound.Simulation()
sim.units = ('AU', 'Msun', 'year')
sim.integrator = "IAS15"
sim.dt = 0.001

Simulation 객체를 생성하고 단위를 천체역학 표준값(AU, 태양질량, 년)으로 지정한다. 이는 중력 상수 (G) 를 (4\pi^2) 로 정규화하여 부동소수점 오차를 방지한다. IAS15 적분기를 선택해 고정밀·적응형 연산을 수행하고, 초기 타임스텝을 0.001로 설정한다.

2. 입자 추가

sim.add(m=1.0)                         # 태양 (질량 = 1.0 Msun)
sim.add(m=3.003e-6, a=1.0, e=0.0167)   # 지구 (질량 ≈ 3e‑6 Msun, 반경축 1 AU, 이심률 0.0167)

REBOUND는 Cartesian 좌표 혹은 케플러 원소로 입자를 추가할 수 있다. 여기서는 태양을 원점에 두기 위한 암시적 위치와 지구의 궤도 요소를 전달한다. REBOUND가 자동으로 초기 속도를 계산해준다.

3. 안정화

sim.move_to_com()

이 단계는 시스템 전체를 중점(COM) 프레임으로 이동시켜 총 운동량과 드리프트를 없애며, 장기 통합 시 수치 오류가 누적되는 것을 방지한다.

4. 통합 및 검증

T_end = 1.0       # 시뮬레이션 종료 시간 (1년)
N_steps = 100     # 출력 데이터 포인트 수
times = np.linspace(0, T_end, N_steps)

print(f"시뮬레이션 시작: {T_end} 년(들) ...")
for i, t in enumerate(times):
    sim.integrate(t)

시간을 일정한 간격으로 나누어 반복적으로 시스템을 통합한다. 실제 포지션을 저장하지는 않지만, 최종 시점에서의 궤도 안정성을 확인한다.

5. 결과 출력 및 에너지 보존 검증

print("\n--- 시뮬레이션 결과 ---")
print(f"최종 시간 도달: {sim.t:.6f} 년")

final_earth = sim.particles[1]
print("\n지구 최종 상태:")
print(f"위치 (x, y, z): ({final_earth.x:.6f}, {final_earth.y:.6f}, {final_earth.z:.6f}) AU")
print(f"속도 (vx, vy, vz): ({final_earth.vx:.6f}, {final_earth.vy:.6f}, {final_earth.vz:.6f}) AU/년")

print(f"\n최종 반경축 (a): {final_earth.a:.6f} AU")
print(f"최종 이심률 (e): {final_earth.e:.6f}")

sim.integrate(T_end)  # 정확한 종료 시점을 보장
e_start = -0.000000000000000  # 이 시스템의 이론적 에너지
e_end = sim.calculate_energy()
print(f"\n에너지 보존 검증:")
print(f"종료 시 총 에너지: {e_end:.12f}")

최종 위치·속도·궤도 요소(반경축, 이심률)를 출력하고, symplectic 적분기가 보존해야 할 에너지 보존을 확인한다. 이론적 에너지는 거의 0에 가깝고, 시뮬레이션 종료 시 총 에너지도 동일한 범위에 머무르는 것이 정상적인 결과이다.

이 튜토리얼을 통해 REBOUND가 제공하는 고정밀·장기 안정성 특성을 실연적으로 확인할 수 있습니다.