미분방정식의 수치 해법

일관성 + 안정성 = 수렴

ODEBook 5장의 정리 5.3은 오일러 방법에 대해 "국소 오차 $O(h^2)$가 쌓여 전역 오차 $O(h)$"라는 것을 보였습니다. 그 증명에서 몰래 쓴 성질이 하나 있습니다. 앞선 오차가 다음 걸음에서 많아야 $1 + hL$배로 커진다는 것, 곧 오차가 걸음마다 폭발하지 않는다는 성질입니다. 수렴하려면 이것과 "한 걸음이 맞다"는 성질이 둘 다 필요합니다.

다단계 방법은 사정이 다릅니다. $k$단계 선형 다단계 방법(linear multistep method)은

꼴이고, 다음 값이 $k$개의 이전 값으로 정해지는 차분 방정식입니다. $h = 0$으로 두면 $\sum\alpha_j y_{n+j} = 0$, 해는 $\rho$의 근 $\zeta_i$에 대해 $\zeta_i^n$ 꼴입니다. 참 해에 해당하는 근은 $\zeta = 1$(일관성 조건 $\rho(1) = 0$)이고 나머지는 계산이 만들어 낸 기생 해입니다. 반올림 오차나 출발값 오차가 이 기생 해를 조금이라도 깨우면, $|\zeta_i| > 1$일 때 그것이 $|\zeta_i|^n$으로 자랍니다. $n = (T - t_0)/h$이므로 걸음을 줄일수록 더 많이 자랍니다.

이 "일관성 + 안정성 ⇔ 수렴"은 수치해석 전체를 꿰는 원리입니다. 16장에서 편미분방정식의 격자 방법에 대해 같은 모양의 랙스 동등성 정리를 다시 만납니다(16장). 문제 상자의 공식은 $\rho(\zeta) = \zeta^2 + 4\zeta - 5 = (\zeta - 1)(\zeta + 5)$이라 기생 근 $-5$가 있습니다. 2단계 명시적 공식이 3차라는 것 자체가 제1 장벽(명시적이면 차수 ≤ 2)을 넘었으니, 영 안정일 수가 없었던 것입니다.

룽게–쿠타 계열과 부처 표

$s$단 명시적 룽게–쿠타 방법은 단계 기울기 $k_i$를 앞선 기울기의 조합으로 차례로 계산하고, 마지막에 가중 평균으로 걷습니다.

차수 조건은 방법의 한 걸음과 참 해의 테일러 전개를 $h$의 거듭제곱별로 맞춘 것입니다. 스칼라 자율 방정식만 보면 조건이 실제보다 적게 보이므로, 부처는 연립 방정식의 전개에 나타나는 항을 뿌리 있는 나무(rooted tree)로 분류했습니다. 차수 $q$의 나무 하나가 조건 하나이고, 나무의 개수는 $q = 1, 2, 3, 4, 5$에서 1, 1, 2, 4, 9개입니다.

아래 시뮬레이터에서 부처 표의 칸을 직접 고쳐 보세요. 분수(1/6)로 입력할 수 있고, 차수 조건이 하나씩 켜지고 꺼집니다. 그와 동시에 $y' = y\cos t$를 여러 걸음 폭으로 풀어 로그–로그 기울기를 측정합니다. 조건으로 예측한 차수와 측정한 기울기가 일치하는지가 핵심입니다.

임베디드 룽게–쿠타와 오차 제어

ODEBook 5장 6절은 내장 쌍의 원리를 소개했습니다. 같은 단계 기울기에 두 벌의 가중치 $\mathbf b$(p차)와 $\hat{\mathbf b}$($\hat p$차)를 주면 추가 비용 없이 두 근사의 차이 $\mathbf{e} = h\sum(b_i - \hat b_i)k_i$가 낮은 쪽의 국소 오차 추정이 됩니다. 실제 구현에는 세 가지 결정이 더 들어갑니다.

1. 어느 쪽으로 걸을까. 추정한 것은 낮은 차수 쪽 오차지만, 채택된 걸음은 높은 차수 쪽 값으로 진행합니다. 이것을 국소 외삽(local extrapolation)이라 합니다. 도르만드–프린스 5(4)는 5차로, 보가키–샴핀 3(2)는 3차로 걷습니다.
2. 오차를 무엇과 비교할까. 성분마다 크기가 다르므로 척도 $\mathrm{sc}_i = \mathrm{atol} + \mathrm{rtol}\cdot\max(|y_i|, |y_{n+1,i}|)$로 나눈 노름 $\mathrm{err} = \max_i |e_i| / \mathrm{sc}_i$(또는 RMS)를 쓰고, $\mathrm{err} \le 1$이면 채택합니다.
3. 다음 h는. 국소 오차가 $Ch^{q+1}$($q = \min(p, \hat p)$)에 비례하므로 $h_{\text{new}} = h\cdot\min\bigl(f_{\max},\ \max(f_{\min},\ 0.9\,\mathrm{err}^{-1/(q+1)})\bigr)$. 0.9는 안전 계수, $f_{\min} = 0.2$, $f_{\max} = 5$는 한 번에 너무 크게 바뀌지 않게 하는 제한입니다.

function adaptive_rk(f, t0, y0, T, rtol, atol, tableau)
  t ← t0;  y ← y0;  h ← 초기 추정        // 예: 0.01·|y0| / |f(t0, y0)|
  while t < T
    h ← min(h, T − t)
    k[1..s] ← 단계 기울기 (부처 표 A, c)
    y_hi ← y + h·Σ bᵢ kᵢ                  // 높은 차수: 채택되면 이 값으로 진행
    e    ← h·Σ (bᵢ − b̂ᵢ) kᵢ               // 오차 추정 (추가 f 계산 없음)
    err  ← max over i of |eᵢ| / (atol + rtol·max(|yᵢ|, |y_hiᵢ|))
    if err ≤ 1 then  t ← t + h;  y ← y_hi    // 채택
    end if                                // 아니면 버리고 같은 t에서 다시
    h ← h · min(5, max(0.2, 0.9 · err^(−1/(q+1))))
  end while

시뮬레이터는 세 쌍(호인–오일러 2(1), 보가키–샴핀 3(2), 도르만드–프린스 5(4))을 같은 문제에 돌립니다. 아래쪽 그래프는 tol을 $10^{-3}$에서 $10^{-8}$까지 바꾸며 얻은 전역 오차 대 f 계산 횟수, 이른바 작업–정밀도 그림입니다. 반데르폴($\mu = 5$, $t \in [0, 20]$)에서 tol = $10^{-6}$이면 f 계산 횟수가 호인–오일러 36,944회, 보가키–샴핀 4,308회, 도르만드–프린스 1,960회이고(FSAL 절약은 세지 않음), 전역 오차는 각각 $9.1\times10^{-7}$, $2.9\times10^{-6}$, $3.1\times10^{-7}$입니다.

다단계 방법: 애덤스–배시포스, 애덤스–몰턴, 예측–수정

$y(t_{n+1}) = y(t_n) + \int_{t_n}^{t_{n+1}} f(t, y(t))\,dt$에서 피적분 함수를 이미 아는 점들을 지나는 보간 다항식(5장)으로 바꾸고 정확히 적분합니다. 보간점으로 $t_n, t_{n-1}, \dots, t_{n-k+1}$을 쓰면 구간 밖으로 외삽하는 명시적 공식, $t_{n+1}$까지 넣으면 음적 공식이 됩니다.

AB2의 계수는 $f_{n-1}, f_n$을 지나는 직선 $p(t) = f_n + (t - t_n)(f_n - f_{n-1})/h$를 $[t_n, t_n + h]$에서 적분하면 $h f_n + \tfrac h2(f_n - f_{n-1}) = \tfrac h2(3f_n - f_{n-1})$으로 바로 나옵니다. 다단계 방법에는 실무상 두 가지 일이 따라붙습니다. 첫째, 출발: 처음 $k - 1$개의 값은 RK4 같은 일단계 방법으로 만듭니다. 둘째, 걸음 폭 바꾸기: 계수가 등간격을 가정하므로 $h$를 바꾸려면 보간을 다시 하거나 계수를 다시 계산해야 합니다.

음적인 AM은 $y_{n+1}$을 풀어야 합니다. 강성이 없는 문제에서는 뉴턴 반복까지 필요 없고, AB로 예측한 값을 AM의 우변에 넣어 한 번 수정하는 것으로 충분합니다.

// PECE: AB(k) 예측 → f 계산 → AM 수정 → f 계산
for n = k−1, k, …
  y_p ← yₙ + h·Σ βⱼᴬᴮ fₙ₋ⱼ              // P: 예측 (외삽)
  f_p ← f(tₙ₊₁, y_p)                    // E: 계산
  yₙ₊₁ ← yₙ + h·(β₀ᴬᴹ f_p + Σ βⱼᴬᴹ fₙ₊₁₋ⱼ)  // C: 수정
  fₙ₊₁ ← f(tₙ₊₁, yₙ₊₁)                  // E: 다음 걸음용
  est ← |yₙ₊₁ − y_p| · |C_AM| / |C_AB − C_AM|   // 밀른의 오차 추정
end for

마지막 줄이 다단계 방법의 공짜 오차 추정입니다. 예측자와 수정자의 국소 오차가 각각 $C_{AB}h^{p+1}y^{(p+1)}$, $C_{AM}h^{p+1}y^{(p+1)}$이면(같은 차수의 쌍을 쓸 때) 두 값의 차이에서 $y^{(p+1)}$을 소거해 수정자의 오차를 추정할 수 있습니다. 이것을 밀른 장치(Milne's device)라 하며 임베디드 RK의 오차 추정과 같은 역할을 합니다.

표의 출발값은 정확한 해로 주었습니다. 같은 4차라도 PECE 4는 AB4보다 13배 정확하고, f 계산은 두 배입니다. 오차 상수의 비 $\tfrac{251}{19} \approx 13$이 그대로 드러납니다.

안정성 영역과 A-안정

시험 방정식 $y' = \lambda y$에 선형 다단계 방법을 쓰면 $\sum(\alpha_j - z\beta_j)y_{n+j} = 0$($z = h\lambda$)입니다. 이 차분 방정식의 해는 안정성 다항식(stability polynomial) $\pi(\zeta; z) = \rho(\zeta) - z\sigma(\zeta)$의 근 $\zeta_i(z)$의 거듭제곱입니다. 룽게–쿠타의 $|R(z)| \le 1$에 해당하는 조건이 "모든 근이 단위원 안"입니다.

사다리꼴 공식은 A-안정이지만 $R(z) = \frac{1 + z/2}{1 - z/2} \to -1$($z \to -\infty$)이라 L-안정은 아닙니다. 아주 빠른 성분이 지워지지 않고 부호만 바꾸며 남아 회로 시뮬레이터에서 "링잉"이라 불리는 가짜 진동을 만듭니다. BDF1(음적 오일러)과 BDF2는 L-안정입니다. 3차 이상의 BDF는 제2 장벽 때문에 A-안정일 수 없지만 큰 각도의 A(α)-안정을 가집니다.

표의 수치는 시뮬레이터와 같은 방법(격자 위에서 안정성 다항식의 근을 직접 계산하고, BDF의 α는 경계 궤적에서 최소 각을 찾음)으로 계산했습니다. BDF7 이상은 영 안정성마저 잃으므로 실용 BDF 코드는 차수 5 또는 6까지만 씁니다.

강성과 BDF: 음적 방법의 뉴턴 반복

원인은 야코비안 $J = \partial f/\partial y$의 고유값입니다. $t = 40$ 부근에서 $J$의 고유값은 $0$(질량 보존 $y_1 + y_2 + y_3 = 1$), $-0.0214$, $-3393$입니다. 해의 변화는 $-0.02$짜리 느린 모드가 정하지만, 명시적 방법의 걸음은 $h|\lambda_{\max}| \lesssim 3.3$(DP5의 실수축 안정 경계)에 묶여 $h \approx 10^{-3}$을 넘지 못합니다. $40/10^{-3} = 4\times10^4$걸음, 관찰한 걸음 수와 맞습니다. 강성비 $3393/0.0214 \approx 1.6\times10^5$가 걸음 수의 비가 됩니다.

음적 공식은 매 걸음 비선형 방정식 $G(y_{n+1}) = y_{n+1} - \gamma h f(t_{n+1}, y_{n+1}) - \psi_n = 0$($\psi_n$은 과거 값의 조합, BDF2에서 $\gamma = \tfrac23$)을 풉니다. 강성 문제에서는 단순 고정점 반복 $y \leftarrow \psi + \gamma hf(y)$가 $\gamma h\|J\| \lt 1$일 때만 수렴하므로, 바로 그 걸음 제한을 다시 불러옵니다. 그래서 뉴턴 반복(3장, 4장)을 씁니다.

function bdf_step(f, J, t, ψ, γh, y_guess)
  y ← y_guess                         // 과거 값의 외삽(예측자)
  M ← I − γh·J(t, y)                   // 반복 행렬: 걸음마다 한 번만
  LU ← factor(M)                       // 12장의 LU 분해, 여러 반복에서 재사용
  for it = 1, 2, …, maxit
    G ← y − γh·f(t, y) − ψ
    Δ ← solve(LU, −G)                  // 전진·후진 대입만: O(m²)
    y ← y + Δ
    if ‖Δ‖ ≤ 작은 수·‖y‖ then return y
  end for
  return 실패 → h를 줄여 다시

위 구현은 반복 행렬을 걸음마다 한 번만 분해하는 단순화 뉴턴입니다. 수렴은 1차로 떨어지지만 $m$차원 계에서 분해 $O(m^3)$를 반복마다 하지 않아도 되어 훨씬 싸고, 실제 코드는 야코비안을 여러 걸음 동안 재사용하기도 합니다. 큰 계에서 이 분해가 비용의 대부분이라 희소 LU(12장)나 반복 해법(13장)이 강성 ODE 풀이기의 심장이 됩니다.

심플렉틱 적분: 스퇴르머–베를레와 장기 에너지

ODEBook 5장 8절에서 심플렉틱 오일러가 용수철의 수정된 에너지를 정확히 지킨다는 것을 보았습니다. 실무에서 가장 많이 쓰는 심플렉틱 방법은 그것을 대칭으로 만든 2차 방법입니다.

RK4가 에너지를 잃는 이유는 반대로 설명됩니다. RK4의 한 걸음은 넓이를 아주 조금 줄이는 사상이라(선형 진동자에서 걸음당 넓이 비가 $1 - h^6/72 + h^8/576$), 그 작은 수축이 매 걸음 같은 방향으로 쌓입니다. 베를레의 오차는 크지만 한 방향으로 쌓이지 않습니다. 그래서 분자 동역학(수십억 걸음), 행성계의 장기 안정성 연구, 입자 가속기의 궤도 추적, 게임 물리 엔진이 베를레나 그 고차 합성(요시다 4차 등)을 씁니다.

베를레의 오차는 $h$를 반으로 줄일 때마다 4분의 1(2차)이고, RK4는 32분의 1입니다. RK4의 표류는 $t \cdot h^4$에 비례하는데, 이 계수가 작아서 $t = 1000$ 정도에서는 RK4가 여전히 낫습니다. 승부는 적분 시간과 h로 정해지고, 시간이 충분히 길면 반드시 베를레가 이깁니다.

경계값 문제: 슈팅 방법과 유한 차분

첫 번째 방법은 모르는 기울기 $s = y'(0)$을 맞힐 때까지 쏘아 보는 것입니다. 초기값 문제 $y'' = -\lambda e^y$, $y(0) = 0$, $y'(0) = s$를 RK4로 풀어 끝값 $y(1; s)$를 얻고, 빗나간 정도 $F(s) = y(1; s) - 0$의 근을 3장의 근 찾기로 구합니다.

슈팅은 간단하지만 약점이 있습니다. 초기값 문제가 불안정하면(해가 $e^{Lx}$로 벌어지면) $F(s)$가 $s$에 극도로 민감해져 반올림 오차만으로도 맞출 수 없습니다. 구간을 여러 조각으로 나눠 조각마다 쏘고 이음매를 맞추는 다중 슈팅이 처방입니다. 두 번째 방법은 구간 전체를 한꺼번에 푸는 것입니다.