수치 해법

오일러 방법과 국소·전역 오차

초기값 문제 $y' = f(t, y)$, $y(t_0) = y_0$의 해를 시각 $t_n = t_0 + nh$에서의 근삿값 $y_n$으로 구합니다. $h$를 걸음 폭(step size)이라 합니다. 해 $y(t)$를 $t_n$에서 테일러 전개하면(CalculusBook 12장)

마지막 항을 버리고, 모르는 $y(t_n)$ 대신 이미 계산한 $y_n$을 넣은 것이 오일러 방법입니다.

버린 항이 오차의 원천입니다. 그런데 오차에는 두 종류가 있다는 것을 구별해야 합니다. 한 걸음에서 새로 생기는 오차와, 지금까지 쌓인 오차입니다.

국소 오차는 위의 테일러 나머지 그 자체라서 $|\tau_{n+1}| \le \frac{M}{2}h^2$($M = \max|y''|$)입니다. 구간 $[t_0, T]$를 지나려면 $N = (T - t_0)/h$걸음이 필요하니, 국소 오차를 그냥 더하면 $N\cdot\frac{M}{2}h^2 = \frac{M(T - t_0)}{2}h$, 즉 $h$의 1제곱입니다. 실제로는 앞선 오차가 다음 걸음에서 해곡선의 벌어짐에 따라 커지거나 줄어드므로 한 가지 인자가 더 붙습니다.

정리의 메시지는 둘입니다. 첫째, 전역 오차는 $h$에 비례합니다. $h$를 반으로 줄이면 오차도 반이 됩니다. 국소 오차는 $h^2$인데 걸음 수가 $1/h$이므로 하나가 깎입니다. 둘째, $e^{L(t - t_0)}$ 인자 때문에 해곡선이 서로 벌어지는 방정식($\partial f/\partial y > 0$)에서는 오래 적분할수록 오차가 지수적으로 커질 수 있습니다. 해곡선이 모이는 방정식에서는 이 상한이 지나치게 비관적입니다.

아래 시뮬레이터에서 각 걸음의 국소 오차(노란 막대)와 끝점의 전역 오차(청록 막대)를 보세요. 점선은 각 걸음의 출발점을 지나는 정확한 해 $u_n$입니다. $y' = y$에서는 점선들이 위로 벌어지므로 국소 오차가 쌓이며 커지고, $y' = t - y$에서는 점선들이 한데 모이므로 앞선 오차가 저절로 줄어듭니다.

오차 차수: 로그–로그 그래프의 기울기

정리 5.3은 오일러 방법의 오차가 $Ch$ 꼴이라는 것을 말합니다. 더 좋은 방법은 $Ch^2$, $Ch^4$처럼 $h$의 더 높은 거듭제곱으로 줄어듭니다. 이 지수를 방법의 차수라 합니다.

$E \approx Ch^p$의 양변에 로그를 취하면 $\log E \approx \log C + p\log h$. 그래서 오차를 $h$에 대해 로그–로그 눈금으로 그리면 직선이 되고, 그 기울기가 차수 $p$입니다. 두 간격 $h$와 $h/2$의 오차만 있어도 $p \approx \log_2(E_h / E_{h/2})$로 차수를 추정할 수 있습니다. 예제 5.1에서 $\log_2 1.92 = 0.94$, 1에 가깝습니다.

이 기울기는 새 프로그램을 검증하는 가장 강력한 도구입니다. RK4를 구현했는데 기울기가 4가 아니라 1로 나오면 어딘가에 버그가 있습니다(흔히 $k_2$를 계산할 때 $t + h/2$ 대신 $t$를 넣은 경우). 반대로 $h$를 너무 줄이면 오차가 더 이상 줄지 않고 들쭉날쭉해지는 바닥이 나타납니다. 컴퓨터의 실수가 약 16자리만 저장하기 때문에 생기는 반올림 오차로, 걸음 수가 많을수록 쌓입니다(MathBook: 컴퓨터는 어떻게 틀리는가).

개선 오일러(호인)와 중점법

미분방정식을 한 걸음 적분하면 $y(t_{n+1}) = y(t_n) + \int_{t_n}^{t_{n+1}} f(t, y(t))\,dt$입니다. 오일러 방법은 이 적분을 왼쪽 직사각형으로 근사한 것입니다. 사다리꼴 공식을 쓰면 오른쪽 끝의 기울기 $f(t_{n+1}, y_{n+1})$이 필요한데, 그 $y_{n+1}$이 바로 지금 구하려는 값입니다. 해결책은 오일러로 먼저 어림해서(예측) 그 점의 기울기를 재고, 평균을 쓰는(수정) 것입니다.

왜 2차일까요? 참 해의 테일러 전개를 2차 항까지 쓰면 $y(t + h) = y + hf + \frac{h^2}{2}(f_t + f_y f) + O(h^3)$입니다(연쇄 법칙으로 $y'' = f_t + f_y y' = f_t + f_y f$). 한편 호인의 $k_2$를 전개하면 $k_2 = f + h(f_t + f_y f) + O(h^2)$이므로 $y + \frac{h}{2}(k_1 + k_2) = y + hf + \frac{h^2}{2}(f_t + f_y f) + O(h^3)$. 2차 항까지 정확히 일치합니다.

아래 시뮬레이터는 세 방법(그리고 다음 절의 RK4)을 같은 방향장 위에서 나란히 걷게 합니다. 출발점을 끌고 걸음 폭을 바꿔 보세요. 회색 굵은 선은 아주 작은 걸음으로 계산한 기준 해입니다.

룽게–쿠타 4차: 기울기 네 개의 가중 평균

같은 생각을 더 밀고 나가면, 한 걸음 안에서 기울기를 네 번 재서 테일러 전개의 4차 항까지 맞출 수 있습니다. 이것이 공학 계산의 표준 도구인 룽게–쿠타 4차 방법(classical Runge–Kutta, RK4)입니다.

가중치 $1, 2, 2, 1$(합 6)이 낯익습니다. $f$가 $t$에만 의존하면($y' = f(t)$) $k_2 = k_3 = f(t_n + h/2)$이고 RK4는 $\frac{h}{6}\bigl(f(t_n) + 4f(t_n + \frac h2) + f(t_n + h)\bigr)$, 바로 심프슨 공식입니다. 심프슨 공식이 4차인 것처럼 RK4도 4차입니다.

표의 마지막 열이 차수입니다. $h$를 반으로 줄이면 RK4의 오차는 $2^4 = 16$분의 1이 됩니다. 아래 시뮬레이터에서 기울기 네 개가 어디서 재어지는지 단계별로 보세요.

연립 방정식과 고계 방정식에 적용하기

요령은 두 가지입니다. 첫째, 고계 방정식은 1계 연립으로 바꿉니다. 각도 $\theta$와 각속도 $\omega = \theta'$를 함께 상태로 보면 $\theta' = \omega$, $\omega' = -\frac{g}{L}\sin\theta$입니다. 일반적으로 $y^{(m)} = F(t, y, y', \dots, y^{(m-1)})$은 $\mathbf{u} = (y, y', \dots, y^{(m-1)})$로 두면 $\mathbf{u}' = \mathbf{f}(t, \mathbf{u})$ 꼴의 1계 연립이 됩니다(11장에서 다시 씁니다). 둘째, 오일러·호인·RK4의 식은 그대로 벡터에 적용됩니다. $y_n$을 벡터 $\mathbf{u}_n$으로, $f$를 벡터 함수로 읽기만 하면 됩니다.

표는 RK4로 계산한 값입니다. 진폭 30°에서도 단진동 공식과 1.7%밖에 차이 나지 않으니 괘종시계의 작은 흔들림에는 단진동 근사가 충분하고, 180°에 가까워지면 진자가 꼭대기에서 오래 머물러 주기가 한없이 길어집니다. 위상 평면의 그림은 12장에서 다시 봅니다.

적응 간격: 걸음 폭을 스스로 고르기

걸음 폭을 고르려면 지금 걸음의 오차를 알아야 하는데, 참 해를 모르니 오차를 직접 잴 수는 없습니다. 대신 차수가 다른 두 방법으로 같은 걸음을 계산해 그 차이를 오차의 추정값으로 씁니다. 차수가 높은 쪽이 훨씬 정확하므로 차이는 거의 낮은 쪽의 오차입니다.

마지막 식의 뜻은 간단합니다. 국소 오차가 $Ch^{p+1}$에 비례하므로, 오차를 $\mathrm{tol}$로 맞추려면 $h$에 $(\mathrm{tol}/\mathrm{err})^{1/(p+1)}$을 곱하면 됩니다. 0.9는 다음 걸음이 다시 버려지지 않도록 두는 안전 계수입니다. 실무의 표준은 6단 공식 안에 4차와 5차 근사를 함께 담은 RK45 계열이며, ode45와 solve_ivp의 기본값인 도르만드–프린스 쌍(Dormand–Prince)도 그중 하나입니다.

시뮬레이터에서 이심률을 0.9, tol을 10⁻⁶으로 두면 적응 방법은 한 바퀴에 54걸음을 채택하고(24번 버림) 위치 오차 약 $9\times10^{-4}$으로 돌아옵니다. 고정 간격 RK4로 같은 오차를 얻으려면 걸음 수를 두 배씩 늘려 찾았을 때 3200걸음이 필요합니다. 걸음 폭은 근일점 근처의 10⁻³ 수준부터 원일점 근처의 0.5 가까이까지 수백 배 바뀝니다. 적응 간격이 $h$를 고르는 일을 사용자에게서 덜어 줄 뿐 아니라 계산량도 몇 분의 1로 줄인 것입니다.

안정성 영역과 강성, 음적 오일러

이 현상을 이해하는 열쇠는 가장 단순한 방정식, 시험 방정식 $y' = \lambda y$입니다. $\lambda$는 복소수일 수도 있습니다(진동하는 성분). 실수부가 음수면 참 해 $e^{\lambda t}$는 0으로 줄어듭니다. 수치 해도 줄어드는가가 질문입니다. 앞 절의 계산처럼 어떤 룽게–쿠타 방법도 이 방정식에서는 $y_{n+1} = R(z)\,y_n$, $z = h\lambda$ 꼴이 됩니다.

RK4의 영역은 더 넓어서 실수축에서는 $-2.785 \le z \le 0$, 허수축에서는 $|z| \le 2\sqrt2 \approx 2.83$까지 안정합니다. 그래도 크기가 유한하다는 점은 같습니다. 명시적 방법(다음 값을 이미 아는 값으로만 계산하는 방법)의 $R(z)$는 다항식이고, 다항식은 $|z| \to \infty$에서 커지므로 영역은 언제나 유한합니다.

이제 문제 상자의 회로를 이해할 수 있습니다. 다음 방정식은 해가 $\cos t$를 따라가는 느린 성분과 $e^{-\lambda t}$로 사라지는 빠른 과도 성분을 함께 가집니다.

$\lambda = 50$이면 과도 성분은 $t = 0.1$ 무렵 이미 $e^{-5} \approx 0.007$로 사라지고, 그 뒤의 해는 $\cos t$처럼 매끄럽습니다. 정확도만 보면 $h = 0.1$이면 충분합니다. 그러나 오일러 방법의 오차는 매 걸음 $1 - h\lambda = -4$배로 증폭되므로 $h \le 2/50 = 0.04$를 지켜야 합니다. 해의 모양이 아니라 이미 사라진 성분의 빠르기가 걸음 폭을 묶는 것, 이것이 강성입니다.

처방은 오른쪽 끝의 기울기를 쓰는 것입니다. $y_{n+1} = y_n + hf(t_{n+1}, y_{n+1})$은 구하려는 $y_{n+1}$이 양변에 있는 방정식이라 매 걸음 그것을 풀어야 합니다(선형이면 나눗셈 한 번, 비선형이면 뉴턴 반복). 대신 시험 방정식에서 $y_{n+1} = y_n + hλy_{n+1}$, 곧 $R(z) = 1/(1 - z)$이고, $\operatorname{Re}z \lt 0$이면 $|1 - z| > 1$이라 언제나 $|R| \lt 1$입니다.

회로 시뮬레이터 SPICE가 사다리꼴 공식이나 후방 차분 공식(BDF) 같은 음적 방법을 쓰는 것, MATLAB이 강성 문제용으로 ode15s를 따로 두는 것이 이 때문입니다. 화학 반응(빠른 반응과 느린 반응이 섞인 계), 열 방정식을 격자로 나눈 연립 방정식(16장의 $r \le 1/2$ 조건)도 대표적인 강성 문제입니다.

에너지 보존: 심플렉틱 오일러와 궤도 계산

용수철 $x'' = -x$(위상 평면의 원 운동)에서 정확히 계산해 볼 수 있습니다. 오일러 방법은 $x_{n+1} = x_n + hv_n$, $v_{n+1} = v_n - hx_n$이고, 에너지 $x^2 + v^2$를 계산하면

$h = 0.1$로 100걸음($t = 10$, 1.6바퀴)이면 $1.01^{100} = 2.70$배, 에너지가 거의 세 배가 됩니다. 이제 아주 작은 수정을 해 봅시다. 속도를 먼저 갱신하고, 새 속도로 위치를 갱신합니다.

$x'' = -x$에서 한 걸음은 행렬 $\begin{bmatrix}x_{n+1}\\v_{n+1}\end{bmatrix} = \begin{bmatrix}1 - h^2 & h\\ -h & 1\end{bmatrix}\begin{bmatrix}x_n\\v_n\end{bmatrix}$이고 행렬식은 $(1 - h^2) + h^2 = 1$입니다. 오일러 행렬 $\begin{bmatrix}1 & h\\ -h & 1\end{bmatrix}$의 행렬식은 $1 + h^2$입니다. 행렬식이 1이라는 것은 위상 평면의 넓이를 보존한다는 뜻이고(LinearBook: 행렬식), 이것이 2차원에서 "심플렉틱"의 의미입니다. 참 해의 흐름도 넓이를 보존하는데(해밀턴 계의 리우빌 정리), 오일러는 이 성질을 깨고 심플렉틱 오일러는 지킵니다.

RK4는 어떨까요? 정확도가 높아 에너지 오차가 매우 작지만 심플렉틱은 아닙니다. $x'' = -x$에서 $|R(ih)|^2 = 1 - h^6/72 + h^8/576$이라 에너지가 매 걸음 아주 조금씩 줄어듭니다($h = 0.1$이면 걸음당 $1.4\times10^{-8}$). 짧은 계산에서는 RK4가 압도적으로 낫지만, 수억 걸음을 도는 태양계 장기 적분이나 분자 동역학에서는 이 꾸준한 표류가 문제가 되어 심플렉틱 방법(속도 베를레 등)을 씁니다. 게임 물리 엔진이 흔히 반음적 오일러를 쓰는 것도 같은 이유입니다.

이 도구가 쓰이는 곳

교과서의 방정식은 대부분 손으로 풀리지만, 실제 공학의 방정식은 대부분 수치로만 풀립니다. 그래서 어떤 방법이 어떤 문제에 맞는지 고르는 눈이 중요합니다.