[RK45 오류 해결] N-Body 시뮬레이션 RK45 구현: 궤적 왜곡 및 적응형 스텝 사이즈 문제 원인과 해결 방법

안녕하세요, 기술 블로그 독자 여러분. N-body 시뮬레이션과 같은 복잡한 물리 시스템을 구현할 때, 수치 해석 기법의 정확성은 매우 중요합니다. 오늘은 스택오버플로우에서 제기된 RK45 (Runge-Kutta-Fehlberg) 메서드를 이용한 N-body 시뮬레이션 중 궤적이 잘못 나타나는 문제에 대해 심층적으로 분석하고, 그 해결 방법을 제시하고자 합니다.

특히 “Golden Pancake Drama”와 같은 알려진 궤적을 재현하려 할 때, 예상과 달리 궤적이 “shifted/squashed” 되거나 닫힌 루프가 형성되지 않는 경험을 하셨다면 이 글이 큰 도움이 될 것입니다. 본 문서는 구글 검색 상위 노출과 개발자 사용자 경험(UX)을 최적화하여 작성되었습니다.

1. 에러 발생 상황

질문자는 Python과 NumPy, Matplotlib을 사용하여 N-body 시뮬레이션(특히 3체 문제)을 위해 RK45 적분기를 직접 구현하는 과정에서 문제를 겪고 있습니다. 알려진 주기 궤도(예: “Golden Pancake Drama”, “Broucke A2”)를 시뮬레이션 했을 때, 궤적의 전반적인 형태는 비슷하지만 실제 값은 잘못되어 플롯이 ‘ shifted/squashed’(이동되거나 찌그러진) 보이거나, 닫힌 루프가 제대로 형성되지 않는 현상이 발생했습니다.

아래 이미지는 예상되는 궤적(RK4 구현)과 RK45 구현으로 출력된 궤적을 비교한 것입니다. 출력된 궤적이 불안정하고 예상 경로에서 벗어난 것을 명확히 볼 수 있습니다.

예상되는 궤적 (RK4 구현):

실제 출력된 궤적 (RK45 구현):

질문자는 rk45_step 함수에서 dt(시간 스텝) 변수가 업데이트된 후 최종적으로 반환되기 전에 이미 사용되는 것이 문제일 수 있다고 추측했습니다. 이는 문제의 일부분에 해당하지만, 더 근본적인 원인이 존재합니다.

2. 명확한 발생 원인

제공된 코드 분석 결과, 궤적 왜곡의 주요 원인은 두 가지입니다.

2.1. RK45 중간 스텝 계산 시 다른 입자들의 상태 참조 오류 (가장 치명적인 원인)

RK45와 같은 Runge-Kutta 계열 적분기는 다음 시간 스텝을 정확히 예측하기 위해 여러 ‘중간 스텝(k1, k2, …, k6)’을 계산합니다. 이때 각 중간 스텝에서는 시스템의 임시적인(trial) 상태를 사용하여 미분값(속도와 가속도)을 평가해야 합니다. 그러나 질문자의 코드에서는 다음과 같은 문제가 발생했습니다.

• arr_calc(inp_arr, index) 함수는 특정 입자(index)의 현재 위치/속도(inp_arr)를 인자로 받습니다.
• 하지만 이 함수 내에서 다른 입자들의 위치(m, x2, y2)를 가져올 때, 전역 변수 state_arr를 직접 참조합니다.
• 이것은 k2, k3 등 중간 스텝을 계산할 때, 대상 입자의 위치는 임시 중간 위치를 사용하지만, 다른 입자들의 위치는 항상 초기(이전) 시간 스텝의 위치를 사용한다는 의미가 됩니다.
• N-body 시뮬레이션에서는 모든 입자가 상호작용하기 때문에, 각 중간 스텝에서 모든 입자의 임시적인 위치를 사용하여 힘과 가속도를 다시 계산해야 합니다. 이러한 불일치로 인해 힘 계산이 부정확해지고 궤적 왜곡이 발생합니다.

2.2. 적응형 스텝 사이즈 조절 로직의 불일치

rk45_step 함수는 dt 값을 적응적으로 조절하는 로직을 포함하고 있습니다. 하지만 구현 방식에 문제가 있습니다.

• k_list (각 입자에 대한 RK45 스텝 결과)는 특정 dt를 기반으로 계산됩니다.
• 이후 에러(r_val)를 평가하고, 이 에러에 따라 dt는 delta 값으로 조정됩니다 (dt *= delta).
• 만약 에러가 허용치 이내이면, k_list를 현재 상태에 적용하지만, 이 k_list는 dt가 조정되기 전의 값으로 계산된 것입니다. 즉, k_list와 실제로 적용된 (혹은 다음 스텝에 사용될) dt 사이에 불일치가 발생합니다.
• RK45의 적응형 스텝 사이즈 로직은 에러가 허용치 이내일 경우 스텝을 수락하고, 에러가 크면 스텝을 거부하고 더 작은 dt로 처음부터 다시 계산해야 합니다. 현재 코드는 스텝을 거부하는 대신 dt를 변경한 후 이전 dt로 계산된 k_list를 적용하는 방식으로 동작하여 정확도를 떨어뜨립니다.

3. 해결 방법 및 코드 예시

위에서 언급된 두 가지 주요 문제를 해결하기 위해 코드를 수정해야 합니다. 핵심은 미분 함수(가속도 계산)가 전체 시스템의 ‘현재 임시 상태’를 정확히 참조하도록 하고, 적응형 스텝 사이즈 로직을 Runge-Kutta-Fehlberg 알고리즘에 맞게 재구성하는 것입니다.

3.1. get_derivatives 함수 재정의: 전체 시스템 상태 반영

가장 먼저, 시스템의 모든 입자에 대한 미분값(속도 및 가속도)을 계산하는 새로운 함수 get_derivatives를 정의합니다. 이 함수는 특정 시점의 전체 시스템 상태 배열(current_full_state_arr)을 인자로 받아 모든 입자들의 미분값을 반환합니다. 이렇게 함으로써 RK45의 중간 스텝 계산 시 다른 입자들의 위치를 정확하게 반영할 수 있습니다.

import numpy as np
from matplotlib import pyplot as plt

# (기존의 softening, big_G, max_step, min_step, err_tolerance 등 전역 변수는 동일)
softening = 1e-6
big_G = 1
max_step = 5e-5
min_step = 1e-12
err_tolerance = 1.5e-11

# 초기 상태 및 시뮬레이션 파라미터 (원 코드와 동일)
state_arr = np.array([
    [1, -1.000000000000, 0.000000000000, 0.151837333844, 0.139577137808],
    [0.6123560888279949, 0.000000000000, 0.000000000000, -0.495911893796, -0.455869192302],
    [1, 1.000000000000, 0.000000000000, 0.151837333844, 0.139577137808]
])
plot_title = 'GOLDEN PANCAKE DRAMA'
delta_t = 1e-4
duration_t = 8.8457131110

def get_derivatives(current_full_state_arr):
    # current_full_state_arr: (N, 5) 배열, 각 입자의 [m, x, y, vx, vy] 포함
    num_bodies = len(current_full_state_arr)
    
    # 각 입자의 [vx, vy, ax, ay]를 저장할 배열
    derivatives = np.zeros((num_bodies, 4)) 

    for i in range(num_bodies):
        m1, x1, y1, vx1, vy1 = current_full_state_arr[i]
        
        tot_acc_x = 0
        tot_acc_y = 0

        for j in range(num_bodies):
            if i != j:
                m2, x2, y2, _, _ = current_full_state_arr[j] # 다른 입자의 질량, 위치만 필요
                
                rx = x1 - x2
                ry = y1 - y2
                r_mag_sqr = rx**2 + ry**2 + softening**2
                
                # 매우 근접한 거리에 대한 처리 (softening이 있지만 0에 가까워지는 경우 방지)
                if r_mag_sqr == 0:
                    continue # 또는 적절한 오류 처리
                
                force_magnitude = big_G * m2 / (r_mag_sqr**1.5)
                tot_acc_x -= force_magnitude * rx
                tot_acc_y -= force_magnitude * ry
        
        derivatives[i] = np.array([vx1, vy1, tot_acc_x, tot_acc_y])
    return derivatives

3.2. rk45_step 함수 재구조화: 정확한 적응형 스텝 구현

이제 rk45_step 함수는 get_derivatives를 사용하여 전체 시스템의 상태를 기반으로 각 k 값을 계산합니다. 또한, 적응형 스텝 로직을 Runge-Kutta-Fehlberg의 원칙에 따라 스텝 수락/거부 방식으로 수정합니다.

def rk45_step(current_state_arr, dt_guess, err_tol, max_h, min_h):
    # current_state_arr: (N, 5) 배열, 스텝 시작 시 시스템 상태
    # dt_guess: 현재 시도하는 시간 스텝
    
    while True:
        # 1. k1 계산 (현재 상태에서의 미분값)
        k1_derivs = get_derivatives(current_state_arr) # (N, 4) 배열
        k1 = k1_derivs * dt_guess # k1은 dt_guess 동안의 변화량 (N, 4)

        # 2. k2 계산 (k1을 이용한 중간 상태에서의 미분값)
        trial_state_k2 = np.copy(current_state_arr)
        trial_state_k2[:, 1:] += k1 / 4 # 질량(0번째 열) 제외, 위치/속도(1번째부터)만 업데이트
        k2_derivs = get_derivatives(trial_state_k2)
        k2 = k2_derivs * dt_guess

        # 3. k3 계산
        trial_state_k3 = np.copy(current_state_arr)
        trial_state_k3[:, 1:] += (k1*3/32 + k2*9/32)
        k3_derivs = get_derivatives(trial_state_k3)
        k3 = k3_derivs * dt_guess

        # 4. k4 계산
        trial_state_k4 = np.copy(current_state_arr)
        trial_state_k4[:, 1:] += (k1*1932/2197 - k2*7200/2197 + k3*7296/2197)
        k4_derivs = get_derivatives(trial_state_k4)
        k4 = k4_derivs * dt_guess

        # 5. k5 계산
        trial_state_k5 = np.copy(current_state_arr)
        trial_state_k5[:, 1:] += (k1*439/216 - k2*8 + k3*3680/513 - k4*845/4104)
        k5_derivs = get_derivatives(trial_state_k5)
        k5 = k5_derivs * dt_guess

        # 6. k6 계산
        trial_state_k6 = np.copy(current_state_arr)
        trial_state_k6[:, 1:] += (-k1*8/27 + k2*2 - k3*3544/2565 + k4*1859/4104 - k5*11/40)
        k6_derivs = get_derivatives(trial_state_k6)
        k6 = k6_derivs * dt_guess

        # 5차 해법을 위한 최종 스텝 변화량 (k_sol) 계산
        k_sol = k1*25/216 + k3*1408/2565 + k4*2197/4104 - k5*1/5

        # 4차와 5차 해법의 차이를 이용한 오차 추정 (error_vec)
        error_vec = k1*1/360 - k3*128/4275 - k4*2197/75240 + k5*1/50 + k6*2/55
        
        # 오차의 크기 계산: 각 입자별 오차 벡터의 크기를 구하고, 그 중 최대값을 사용
        # 원 코드의 arr_mag(r)는 스칼라였지만, 여기서는 각 입자(행)의 오차 벡터 크기를 구함.
        # dt_guess로 나눠 '미분값의 차이'로 정규화하는 것이 일반적임.
        r_val = np.max(np.linalg.norm(error_vec / dt_guess, axis=1)) 

        # 적응형 스텝 사이즈 조절 로직
        # r_val이 0인 경우 (오차가 없는 경우) 무한히 dt를 키우는 것을 방지
        delta = (err_tol / (2 * r_val))**(1/4) if r_val > 0 else 2.0 
        delta = np.clip(delta, 0.1, 10.0) # dt 변화폭 제한 (과도한 증감 방지)

        if r_val <= err_tol:
            # 스텝 수락: k_sol (dt_guess로 계산된)을 사용하여 상태 업데이트
            next_state_arr = np.copy(current_state_arr)
            next_state_arr[:, 1:] += k_sol # 질량 제외하고 위치/속도만 업데이트
            
            # 다음 스텝을 위한 새로운 dt 제안
            dt_next = dt_guess * delta
            dt_next = np.clip(dt_next, min_h, max_h) # 최소/최대 스텝 사이즈 제한
            return next_state_arr, dt_next # 업데이트된 상태와 다음 dt 반환
        else:
            # 스텝 거부: 더 작은 dt_guess로 현재 스텝 다시 시도
            dt_guess = dt_guess * delta # dt_guess를 업데이트하여 while 루프의 다음 반복에서 사용
            dt_guess = np.clip(dt_guess, min_h, max_h)
            # dt_guess가 min_h에 도달했지만 여전히 에러가 큰 경우 (문제 발생 가능성)
            if dt_guess == min_h and r_val > err_tol:
                print(f"Warning: Minimum step size reached ({min_h}) but error tolerance not met (r_val={r_val}).")
            # 루프는 새로운 dt_guess로 현재 스텝을 다시 시도

3.3. 메인 시뮬레이션 루프 (기존과 동일)

메인 시뮬레이션 루프는 rk45_step 함수의 반환값을 받아 state_arr와 delta_t를 업데이트합니다.

x1_plot = []
y1_plot = []
x2_plot = []
y2_plot = []
x3_plot = []
y3_plot = []
t_plot = []

time_val = 0
n = 0
frame_interval = 50

while time_val < duration_t:
    # 수정된 rk45_step 함수 호출
    # rk45_step은 이제 전체 state_arr를 받아 업데이트된 state_arr와 dt_next를 반환
    state_arr, delta_t = rk45_step(state_arr, delta_t, err_tolerance, max_step, min_step)

    # 시간 업데이트는 스텝 수락 후에 이루어져야 함
    time_val += delta_t 
    n += 1

    if n % frame_interval == 0:
        t_plot.append(time_val)

        x1_plot.append(state_arr[0][1])
        y1_plot.append(state_arr[0][2])

        x2_plot.append(state_arr[1][1])
        y2_plot.append(state_arr[1][2])

        x3_plot.append(state_arr[2][1])
        y3_plot.append(state_arr[2][2])

        print(f'TRAJECTORY CALCULATION: {round(100*time_val / duration_t, 3):.3f}%, {round(time_val, 15):.12f}, dt = 1e{np.log10(delta_t):.2f}s')

f=plt.figure(layout='constrained', figsize=(5,5))
plt.plot(x1_plot, y1_plot, c='blue', linewidth=0.25, label='m1')
plt.plot(x2_plot, y2_plot, c='red', linewidth=0.25, label='m2')
plt.plot(x3_plot, y3_plot, c='green', linewidth=0.25, label='m3')
f.legend(loc='outside lower center', fontsize=9, ncols=len(state_arr)+2)
plt.title(plot_title) # plot_title 추가
plt.show()

위와 같이 코드를 수정하면, RK45의 적응형 스텝 로직이 정확하게 동작하고, N-body 시뮬레이션에서 각 입자의 상호작용이 올바르게 계산되어 예상되는 궤적에 훨씬 더 근접한 결과를 얻을 수 있을 것입니다.

4. 향후 예방을 위한 팁

이러한 복잡한 수치 적분기 구현 시 발생하는 문제를 예방하기 위한 몇 가지 팁입니다.

1. 시스템 상태의 정확한 추적: N-body 시뮬레이션과 같이 여러 요소가 상호작용하는 시스템에서는 모든 중간 계산 단계에서 각 요소의 '현재' 또는 '임시' 상태가 정확히 참조되는지 면밀히 확인해야 합니다. 특히 전역 변수를 사용하는 경우, 예상치 못한 부작용을 일으킬 수 있으므로 주의가 필요합니다.
2. 적응형 스텝 로직의 일관성: Runge-Kutta-Fehlberg와 같은 적응형 스텝 사이즈 메서드는 '스텝 수락/거부' 로직이 핵심입니다. 스텝이 거부되면 해당 스텝에서 계산된 모든 중간 결과(k 값)도 함께 거부하고, 새로운 스텝 크기로 처음부터 다시 계산해야 합니다. dt와 k 값의 적용 시점을 명확히 분리하여 혼동을 방지하세요.
3. 모듈화 및 테스트: get_derivatives와 같이 핵심 계산 로직을 독립적인 함수로 모듈화하면 테스트와 디버깅이 훨씬 용이해집니다. 각 함수가 예상대로 작동하는지 단위 테스트를 통해 검증하는 습관을 들이세요.
4. 레퍼런스 코드 참고: SciPy의 solve_ivp와 같은 검증된 라이브러리들은 N-body 시뮬레이션을 위한 고품질 적분기를 제공합니다. 직접 구현하는 것이 학습 목표일지라도, 이러한 라이브러리들의 소스 코드를 참고하거나 구현 방식을 학습하는 것은 매우 큰 도움이 됩니다.
5. 단계별 검증: 복잡한 시스템을 한 번에 구현하기보다는, 먼저 간단한 2체 문제나 더 짧은 시뮬레이션 기간으로 구현을 검증하여 문제를 빠르게 식별하는 것이 좋습니다.

N-body 시뮬레이션은 매우 흥미롭지만, 그 정확성을 보장하기 위해서는 수치 해석 기법에 대한 깊이 있는 이해가 필요합니다. 이 글이 여러분의 시뮬레이션 문제를 해결하고, 더 정확한 물리 시스템을 구현하는 데 도움이 되기를 바랍니다.

---