Basic Finite Differential Methods

Reference: Seongjai Kim, Numerical Methods for Partial Differential Equations, Department of Mathematics and Statistics, Mississippi State University, Ch. 2 (Numerical Methods for ODEs)

레벨셋 방법은 결국 인터페이스를 나타내는 함수 $\phi$가 시간에 따라 어떻게 변하는지를 나타내는 편미분방정식(PDE)을 푸는 문제입니다. 그런데 이 PDE를 실제로 풀 때 가장 흔한 접근은, 공간을 grid로 이산화한 뒤 시간 방향으로는 각 격자점에서의 값이 어떻게 변하는지를 나타내는 상미분방정식(ODE)을 푸는 것입니다(이를 method of lines라고 부릅니다). 즉, PDE를 실제로 계산 가능하게 만드는 마지막 단계는 결국 ODE를 얼마나 정확하고 효율적으로 시간 전진(time-stepping)시키느냐의 문제로 귀결됩니다. 이 챕터는 그 시간 전진의 기본 도구인 Euler 방법, Runge-Kutta 방법, multi-step 방법을 다룹니다. 화려한 응용은 없지만, 이후 레벨셋 방법의 수치 스킴을 이해하기 위한 가장 기초적인 언어입니다.

Introduction: 왜 ODE를 수치적으로 풀어야 하는가

이 챕터가 다루는 문제는 초기값 문제(initial value problem, IVP)입니다.

$f$와 $\partial f/\partial y$가 $(x_0,y_0)$ 근방에서 연속이면 이 문제는 그 근방에서 유일한 해를 가진다는 사실이 알려져 있습니다. 문제는 대부분의 $f$에 대해 이 해를 닫힌 형태(closed form)로 구할 수 없다는 점입니다. 그래서 실제로는 해 $y(x)$ 자체를 구하는 대신, 유한한 개수의 점에서 근사값을 계산하는 방식을 씁니다.

이 챕터가 다루는 모든 방법은 step-by-step method라는 공통된 틀을 따릅니다. 구간 $[x_0,T]$를 $h=(T-x_0)/n_t$ 크기의 균일한 간격으로 나누어 mesh point $x_n = x_0+nh$를 만들고, $y_0=y(x_0)$에서 출발해서 $y_1, y_2, \ldots$ 순서로 각 mesh point에서의 근사값 $y_n\approx y(x_n)$을 하나씩 계산해 나갑니다. 즉 한 스텝 앞의 정보를 이용해서 다음 스텝의 값을 구하는 재귀적인 구조입니다. 이 챕터는 이 재귀 규칙을 어떻게 설계하느냐에 따라 Taylor-series 방법, Runge-Kutta 방법, multi-step 방법 세 갈래로 나누어 설명합니다.

Taylor-Series Methods

Euler 방법

가장 단순한 접근은 Taylor 급수를 1차 항까지만 남기고 자르는 것입니다.

여기서 핵심은 $y’(x_0)=f(x_0,y_0)$라는 사실입니다. $x=x_0$에서 2차 이상의 항을 버리면

을 얻습니다. 이 아이디어를 다음 구간에도 그대로 반복 적용하면($y(x_1)$과 $y’(x_1)$을 각각 $y_1$과 $f(x_1,y_1)$로 대체), 일반적으로

이라는 Euler 방법을 얻습니다. 기하학적으로 보면 이는 실제 해 곡선 $y(x)$를 여러 개의 직선 조각으로 이은 다각선으로 근사하는 것과 같습니다. 각 조각의 기울기는 그 구간 시작점에서의 $f$ 값이고, 다음 점은 그 기울기를 따라 $h$만큼 이동한 위치입니다. 즉 매 스텝마다 “지금 위치에서의 접선 방향으로 한 걸음 나아가기”를 반복하는 셈입니다.

Euler 방법의 수렴성

이 방법이 실제로 참값에 수렴하는지가 궁금해집니다. 책에서는 $f$가 두 번째 인자에 대해 Lipschitz 조건을 만족한다는 가정, 즉 어떤 $\lambda>0$이 존재해서

이 성립한다는 가정 하에, 오차 $e_n = y_n-y(x_n)$에 대해 다음의 수렴 정리를 증명합니다.

증명의 구조는 이렇습니다. 참값 $y(x_n)$도 국소적으로는 $y(x_{n+1})=y(x_n)+hf(x_n,y(x_n))+O(h^2)$를 만족하므로, 이를 Euler 방법의 점화식과 빼면 $e_{n+1}=e_n+h[f(x_n,y_n)-f(x_n,y(x_n))]+O(h^2)$가 나옵니다. 여기에 Lipschitz 조건을 적용하면 $\lVert e_{n+1}\rVert \le (1+\lambda h)\lVert e_n\rVert + Ch^2$라는 재귀 부등식을 얻고, 이를 귀납법으로 풀면 위의 수렴 부등식이 나옵니다.

이 결과가 의미하는 바는 명확합니다. $h\to 0$일 때 우변 전체가 0으로 수렴하므로(등호 우변을 展開하면 $O(h)$ 오차가 남습니다) Euler 방법은 수렴하는 방법입니다. 다만 오차의 상한이 $(1+\lambda h)^n$ 형태로, 스텝을 반복할수록 지수적으로 커질 수 있는 구조라는 점도 함께 보여줍니다. 이는 이산 버전의 Gronwall 부등식과 정확히 같은 형태로, “국소적으로 작은 오차가 누적되어도 전체 구간에서는 유계로 통제된다”는 수치해석의 표준적인 논증 패턴입니다.

고차 Taylor 방법

Euler 방법이 1차 정확도에 그치는 이유는 Taylor 급수를 1차 항에서 잘랐기 때문입니다. 더 높은 차수까지 남기면 더 정확한 방법을 얻을 수 있습니다. $y(x)$의 $m$차 Taylor 다항식을 $x_n$ 주변에서 전개하면

여기서 $y^{(k)}(x)=f^{(k-1)}(x,y(x))$라는 사실(연쇄법칙을 반복 적용한 결과)을 사용했습니다. 나머지 항을 버리면 $m$차 Taylor 방법을 얻습니다.

$m=1$이면 앞서 본 Euler 방법이 그대로 나오고, $m=2$이면 $y_{n+1}=y_n+h\left[f(x_n,y_n)+\frac{h}{2}f’(x_n,y_n)\right]$이 됩니다. 책은 $y’=y-x^3+x+1$이라는 구체적인 IVP에 대해 $f’$과 $f’‘$을 연쇄법칙으로 직접 계산해서 $T_3(x,y)$를 구하는 예제를 보여주는데, 이 과정에서 왜 고차 Taylor 방법이 실용적으로 잘 쓰이지 않는지가 드러납니다. 차수 $m$이 늘어날수록 $f$의 고차 도함수를 손으로(또는 기호 계산으로) 미리 구해둬야 하는데, $f$가 복잡한 함수라면 이 도함수 계산 자체가 상당히 번거로운 작업이 됩니다.

Runge-Kutta Methods

앞서 지적한 Taylor 방법의 약점, 즉 $f$의 도함수를 직접 계산해야 한다는 문제를 해결하려는 것이 Runge-Kutta 방법입니다. 핵심 아이디어는, $f$의 도함수 대신 $f$ 자체를 구간 $[x_n,x_{n+1}]$ 안의 여러 점에서 평가한 값들의 가중 평균으로 같은 차수의 정확도를 달성하는 것입니다.

여기서 $w_j\ge 0$이고 $\sum w_j=1$이며, $K_j$는 $f(x,y)$를 재귀적으로 평가한 “기울기” 값들입니다. 이 가중합이 $T_m(x_n,y_n)+O(h^m)$을 근사하도록 $w_j$와 $K_j$ 안의 파라미터들을 결정하는 것이 설계의 목표입니다. 즉 Runge-Kutta 방법은 $m$차 Taylor 방법과 같은 차수의 정확도를 내면서도, $f$의 도함수 대신 $f$ 자체의 여러 번의 함수 평가(function evaluation)만으로 이를 달성하려는 시도라고 이해하면 됩니다.

2차 Runge-Kutta 방법의 유도

이 설계 원리가 실제로 어떻게 작동하는지 2차의 경우로 확인해보겠습니다. 다음 형태를 가정합니다.

목표는 $w_1K_1+w_2K_2 = T_2(x_n,y_n)+O(h^2) = f(x_n,y_n)+\frac{h}{2}f’(x_n,y_n)+O(h^2)$가 되도록 $w_1,w_2,\alpha,\beta$를 정하는 것입니다. 좌변(참값의 Taylor 전개)은 $y’=f$, $y’‘=f_x+f_yf$라는 사실을 이용해서

로 전개됩니다. 우변($K_2$ 안의 $f$를 다시 Taylor 전개)은

이 됩니다. 두 전개식의 $h$, $h^2$ 계수를 각각 맞추면 다음 세 조건이 나옵니다.

여기서 흥미로운 점은, 미지수가 4개($w_1,w_2,\alpha,\beta$)인데 조건은 3개뿐이라는 것입니다. 즉 해가 유일하지 않고, 조건을 만족하는 파라미터 조합이 여러 개 존재합니다. 이는 곧 “2차 정확도를 갖는 Runge-Kutta 방법”이 하나가 아니라 한 묶음의 방법군(family)이라는 뜻이고, 실제로 아래에서 보듯 서로 다른 선택이 서로 다른 이름의 방법으로 이어집니다.

대표적인 2차 방법: Heun의 방법과 수정 Euler 방법

$w_1=w_2=\frac12$, $\alpha=\beta=1$을 선택하면

이 되는데, 이것이 Heun의 방법(2차 Runge-Kutta, RK2)입니다. 직관적으로 보면, $K_1$은 구간 시작점에서의 기울기이고 $K_2$는 Euler 방법으로 한 걸음 미리 나가본 지점에서의 기울기입니다. 이 둘을 평균 내어 사용하는 것은 사다리꼴 공식으로 적분을 근사하는 것과 같은 발상입니다.

반대로 $w_1=0$, $w_2=1$, $\alpha=\beta=\frac12$을 선택하면

이 되는데, 이는 구간의 중점에서 기울기 하나만을 평가해서 사용하는 수정 Euler 방법(modified Euler method)입니다. 두 방법 모두 2차 정확도를 갖지만 $f$를 평가하는 지점과 방식이 다릅니다.

4차 Runge-Kutta 방법 (RK4)

실무에서 가장 널리 쓰이는 것은 4차 Runge-Kutta 방법입니다.

구조를 보면 $K_1$은 시작점, $K_2$와 $K_3$은 중점에서의 기울기를 두 번(한 번은 $K_1$을 이용해서, 한 번은 $K_2$를 이용해서 더 정교하게) 평가하고, $K_4$는 $K_3$를 이용해 예측한 끝점에서의 기울기입니다. 이 네 기울기를 $1:2:2:1$ 비율로 가중 평균하는 것이 Simpson 공식과 형태적으로 닮아 있다는 점도 눈여겨볼 만합니다. RK4의 국소절단오차(local truncation error)는 $\frac{h^5}{5!}y^{(5)}(\xi_n)$이고, 이를 전체 구간에 누적하면 전역오차는 $\frac{(T-x_0)h^4}{5!}y^{(5)}(\xi)$ 수준이 됩니다. 즉 스텝 크기 $h$를 절반으로 줄이면 오차는 대략 $2^4=16$배 줄어드는, 4차 정확도의 방법입니다.

Adaptive 방법

정확도를 높이는 가장 단순한 방법은 스텝 크기 $h$를 줄이는 것이지만, 이는 곧 계산량 증가로 직결됩니다. 그런데 실제 문제에서는 해가 완만하게 변하는 구간에서는 큰 스텝으로도 충분하고, 해가 급격히 변하는 구간에서만 작은 스텝이 필요한 경우가 많습니다. adaptive 방법은 매 스텝마다 국소절단오차를 추정해서 스텝 크기를 그때그때 조절하는 방법입니다. 대표적인 예로 책에서 언급하는 Runge-Kutta-Fehlberg 방법(RKF45)은 RK5 계산 결과를 이용해 RK4의 국소절단오차를 추정하고, 이 추정치를 바탕으로 다음 스텝의 $h$를 늘리거나 줄입니다.

Experiments: 정확도 비교

책은 이 세 가지 one-step 방법(Euler, Heun, RK4)의 정확도를 스프링-질량계의 강제 진동 문제로 실제 비교합니다. 대상 문제는 다음 2차 ODE입니다.

이 문제는 $y_1=y$, $y_2=-y_1’/\omega$로 치환하면 1차 시스템으로 바뀌고, 이 시스템에 대해 각 방법을 적용합니다. 실험에 사용된 파라미터는 다음과 같습니다.

이 값들에서 운동의 기본 주기는 $T=2\pi q/\omega = 1$입니다. 정확도는 $t=1$에서의 $\ell^2$-오차 $\lvert y_{n_t}^h - y(1)\rvert$로 측정하고, 여러 스텝 크기 $h=1/n_t$에 대해 오차를 비교합니다.

괄호 안의 값은 $\alpha := \log(E(2h)/E(h))/\log 2$로 계산한 수렴 차수입니다. 표를 보면 Euler는 약 1차, Heun은 정확히 2차, RK4는 정확히 4차의 수렴 차수를 보여서, 앞서 이론적으로 유도한 각 방법의 정확도 차수가 실제 계산에서도 그대로 재현됩니다. 특히 눈에 띄는 것은 오차의 절대적인 크기 차이입니다. $1/h=100$에서 이미 RK4의 오차($2.61\times10^{-5}$)는 Euler의 오차($1.19$)보다 다섯 자리 이상 작습니다. RK4는 한 스텝에 $f$를 4번 평가해야 해서 Euler(1번)보다 스텝당 비용은 크지만, 목표 정확도에 도달하는 데 필요한 전체 스텝 수가 압도적으로 적기 때문에 실용적으로는 훨씬 효율적입니다. 책이 짚듯이, 이것이 RK4가 널리 쓰이는 이유입니다.

Multi-step Methods

Euler와 Runge-Kutta는 모두 single-step(1스텝) 방법입니다. 즉 $y_{n+1}$을 구하는 데 오직 직전 값 $y_n$(과 그 근방에서의 $f$ 평가)만 사용합니다. 반면 multi-step 방법은 $y_{n+1}$을 구할 때 과거 여러 스텝의 값 $y_n, y_{n-1}, \ldots, y_{n+1-m}$과 그 지점들에서의 $f$ 값을 모두 사용합니다.

여기서 $b_0=0$이면 explicit(open) 방법, $b_0\neq 0$이면 implicit(closed) 방법이라고 부릅니다. $b_0\neq 0$인 경우 우변에 미지값 $y_{n+1}$이 $f(x_{n+1},y_{n+1})$ 형태로 다시 등장하므로, $y_{n+1}$을 구하려면 (일반적으로 비선형인) 방정식을 별도로 풀어야 한다는 뜻입니다.

4차 정확도를 갖는 대표적인 multi-step 방법들은 다음과 같습니다($y_i’=f(x_i,y_i)$로 표기).

• Adams-Bashforth 방법 (explicit): $y_{n+1}=y_n+\frac{h}{24}(55y_n’-59y_{n-1}’+37y_{n-2}’-9y_{n-3}’)$
• Adams-Moulton 방법 (implicit): $y_{n+1}=y_n+\frac{h}{24}(9y_{n+1}’+19y_n’-5y_{n-1}’+y_{n-2}’)$
• Adams-Bashforth-Moulton 방법 (predictor-corrector): 먼저 Adams-Bashforth로 예측값 $y_{n+1}^*$를 구하고, $y_{n+1}’^*=f(x_{n+1},y_{n+1}^*)$를 계산한 뒤, 이를 Adams-Moulton 공식의 $y_{n+1}’$ 자리에 대입해서 $y_{n+1}$을 보정합니다.

predictor-corrector 방식이 흥미로운 이유는, implicit 방법이 갖는 비선형 방정식을 직접 풀어야 하는 부담 없이도 implicit 방법에 가까운 정확도와 안정성을 얻으려는 절충안이기 때문입니다. explicit 방법으로 일단 근사값을 구한 뒤(predict), 그 값을 implicit 공식에 대입해서 한 번 더 다듬는(correct) 구조입니다.

multi-step 방법의 실용적인 매력은 계산 비용에 있습니다. RK4가 한 스텝마다 $f$를 4번 새로 평가해야 하는 반면, multi-step 방법은 과거 스텝에서 이미 계산해 둔 $f$ 값들을 재사용하므로 매 스텝마다 새로 필요한 $f$ 평가는 단 1번뿐입니다. 다만 이 방법은 시작 시점에 $y_1,y_2,y_3$처럼 여러 개의 초기값이 필요한데, 이는 아직 계산되지 않은 시점이므로 스스로 시작(self-start)할 수 없습니다. 그래서 책은 이 초기 몇 개의 값을 RK4로 미리 계산해두고 그다음부터 multi-step 방법으로 넘어가는 방식을 제안합니다. 종합적으로 책은 “RK 방법이 충분히 정확하고 구현이 쉬워서, 실무에서는 multi-step 방법이 자주 쓰이지 않는다”고 평가하면서도, ABM이 특정한 경우 강한 안정성(strong stability)을 보이는 경우가 있다고 덧붙입니다.

고차 미분방정식과 연립방정식으로의 변환

지금까지 다룬 Euler, Runge-Kutta, multi-step 방법은 모두 1차 ODE $y’=f(x,y)$를 전제로 설계되었습니다. 그런데 실제 문제는 앞서 스프링-질량계 예제처럼 2차 이상의 미분방정식으로 주어지는 경우가 많습니다. 다행히 이 경우는 새로운 방법을 만들 필요 없이, $u=y’$로 치환해서 1차 연립방정식으로 바꿔주기만 하면 됩니다.

우변에 더 이상 도함수가 등장하지 않는 1차 시스템이 되었으므로, 앞서 다룬 Euler·Runge-Kutta·multi-step 방법을 벡터 함수 버전으로 그대로 적용할 수 있습니다. 실제로 앞선 정확도 비교 실험에서 스프링-질량계의 2차 ODE를 $y_1=y$, $y_2=-y_1’/\omega$의 1차 시스템으로 바꾸어 풀었던 것이 바로 이 절차의 구체적인 사례입니다.

Comment

이 챕터를 다 읽고 나서 기억해야 할 것은, 정확도 차수를 올리는 방법에는 근본적으로 두 가지 축이 있다는 사실입니다. 하나는 $f$의 고차 도함수를 직접 계산해서 Taylor 급수를 더 많은 항까지 쓰는 것(Taylor 방법)이고, 다른 하나는 도함수 대신 $f$를 여러 지점에서 반복 평가해서 그 가중 평균으로 같은 효과를 내는 것(Runge-Kutta 방법)입니다. 그리고 세 번째 선택지로, 매 스텝마다 새로 계산하는 대신 과거에 이미 계산해 둔 값을 재활용하는 것(multi-step 방법)이 있습니다. 이 세 갈래는 서로 경쟁 관계가 아니라 “계산 비용을 어디서 지불할 것인가”에 대한 서로 다른 답이라고 보는 편이 정확합니다.

특히 2차 Runge-Kutta 방법을 유도하는 과정(두 전개식의 $h$, $h^2$ 계수를 맞추는 것)은 이 책의 다른 곳에서도 반복해서 등장할 법한 표준적인 기법입니다. 유한차분 스텐실의 계수를 정하는 것도 본질적으로 같은 방식(Taylor 전개 후 계수 매칭)이므로, 이 절의 유도 과정을 제대로 따라가 두면 이후 PDE의 유한차분 근사를 이해할 때도 같은 논리를 재사용할 수 있습니다.

다만 이 챕터는 안정성(stability)에 대해서는 깊게 다루지 않습니다. Euler 방법의 수렴 정리는 오차가 유계로 통제된다는 것을 보여주지만, 스텝 크기 $h$가 특정 문제에 비해 너무 크면 방법 자체가 발산하는 안정성 문제는 다른 주제입니다. 특히 레벨셋 방법처럼 PDE를 공간 이산화해서 매우 많은(격자점 수만큼의) 연립 ODE로 만드는 상황에서는, explicit 방법(Euler, RK)의 안정성 조건이 스텝 크기에 엄격한 제약(흔히 CFL 조건이라 부르는)을 거는 경우가 많습니다. 이 챕터는 그 논의의 출발점이 되는 기본 도구들을 정리한 것이고, 안정성과 강성(stiffness) 문제는 이후 다른 장에서 다뤄질 주제로 보입니다. 지금 단계에서 챙겨야 할 것은 “Euler는 1차, Heun/수정 Euler는 2차, RK4는 4차”라는 정확도 차수의 감각과, 왜 RK4가 실무 기본값으로 자리 잡았는지에 대한 이유(도함수 계산 없이 높은 정확도, 그리고 정확도 대비 계산 비용의 이점)입니다.