음해법과 Thomas 알고리듬 — 1D 확산을 무조건 안정하게 푸는 한 줄의 트릭

양해법이 폭발하는 한계 너머, 삼대각 행렬 한 번 푸는 비용으로 안정성을 산다

밤 11시에 시뮬레이션을 돌려놓고 아침에 출근한 적이 있다. 출근해서 결과 파일을 열었더니 화면이 NaN으로 가득했다. 격자를 두 배로 촘촘하게 만든 직후였고, 시간 스텝은 그대로 두었다. CFL(Courant–Friedrichs–Lewy) 조건을 깜빡한 대가는 한 시간 분량의 클러스터 비용이었다.

확산 방정식의 양해법은 격자가 절반이 되면 시간 스텝을 4분의 1로 줄여야 한다. 격자 1024개짜리 문제를 풀면, 양해법은 폭발하지 않으려는 그 이유 하나로 수십만 스텝을 밟는다. 이 글은 그 함정을 한 줄로 빠져나오는 음해법(implicit method)과, 음해법의 비용을 거의 없애주는 Thomas 알고리듬(tridiagonal solver)을 다룬다. 끝에서 직접 dt를 키워보고, 양해법이 발산하는 동안 음해법은 잠잠한 모습을 본다.

CFL이 묶어버린 미세 격자#

확산 방정식은 다음과 같다.

\frac{\partial q}{\partial t} = D \frac{\partial^2 q}{\partial x^2}

여기서 $q$ 는 농도(또는 온도), $D$ 는 확산 계수다. 가장 단순한 양해법(forward Euler + 중심차분)은

q_i^{n+1} = q_i^{n} + r\,(q_{i+1}^{n} - 2 q_i^{n} + q_{i-1}^{n}),\quad r = \frac{D\,\Delta t}{\Delta x^2}

이 식이 안정하려면 $r \le 1/2$ 여야 한다. von Neumann 분석을 한 줄로 풀면, 가장 짧은 파장(파장 = $2\Delta x$ )의 증폭 인자가 $|1 - 4r|$ 이고 이것이 $\le 1$ 이려면 $r \le 1/2$ 가 강제된다.

문제는 이 조건이 $\Delta x$ 의 제곱에 묶인다는 점이다. 격자를 $1024 \to 2048$ 로 만들면 $\Delta t$ 는 $1/4$ 이 된다. 이건 자릿수의 횡포다.

미래값을 우변에 — 음해법의 단순한 트릭#

해결책은 어이없을 정도로 단순하다. 우변의 $q^n$ 을 $q^{n+1}$ 로 바꿔 쓴다.

q_i^{n+1} = q_i^{n} + r\,(q_{i+1}^{n+1} - 2 q_i^{n+1} + q_{i-1}^{n+1})

미래값이 우변에 있으니 한 셀씩 차례로 풀 수가 없다. 모든 셀이 동시에 결합되어 행렬 방정식이 된다.

A\,\mathbf{q}^{n+1} = \mathbf{q}^{n}

대신 보상이 크다. 같은 von Neumann 분석으로 증폭 인자는 $1/(1 + 4r)$ 이고, 모든 $r$ 에서 1보다 작다. 즉 $\Delta t$ 를 아무리 키워도 발산하지 않는다 — 무조건 안정(unconditionally stable).

대가는 정확도다. 너무 큰 $\Delta t$ 는 안정하지만 정확하지 않다. 그래도 "발산하지 않는다"는 보장은 코드를 돌리는 사람에게 정신 건강을 준다.

삼대각 행렬과 Thomas의 단번 풀이#

행렬 $A$ 를 적어보자. 1차원에서는 셀 $i$ 가 $i-1$ 과 $i+1$ 에만 연결되므로

A = \begin{pmatrix} 1+r & -r & & & \\ -r & 1+2r & -r & & \\ & -r & 1+2r & -r & \\ & & \ddots & \ddots & \ddots \\ & & & -r & 1+r \end{pmatrix}

이 행렬은 **삼대각(tridiagonal)**이다. 대각선 세 줄만 0이 아니다. 일반적인 가우스 소거는 $\mathcal{O}(N^3)$ 이지만, 삼대각이면 $\mathcal{O}(N)$ 로 끝난다. 이 알고리듬을 Thomas algorithm이라고 부른다 — 1949년 Llewellyn Thomas가 IBM 메모로 정리했다.

Thomas의 핵심은 두 단계다.

Forward elimination — 대각 아래 항을 차례로 제거하면서 새로운 대각 $b'$ 와 우변 $d'$ 를 계산한다.
Backward substitution — 마지막 행부터 거꾸로 대입한다.

각 단계가 셀당 4~5번의 연산이고, 메모리는 세 배열만 있으면 된다. 행렬을 통째로 저장할 필요가 없다.

NumPy로 본 양해/음해 비교#

import numpy as np
 
def explicit_diffusion_step(q, r):
    """Forward Euler + central difference."""
    next_q = q.copy()
    next_q[1:-1] = q[1:-1] + r * (q[2:] - 2 * q[1:-1] + q[:-2])
    return next_q
 
def thomas_solver(a, b, c, d):
    """Solve tridiagonal system A x = d in O(N).
    a: sub-diagonal (a[0] is unused)
    b: main diagonal
    c: super-diagonal (c[-1] is unused)
    d: right-hand side
    """
    n = len(d)
    cp = np.zeros(n)
    dp = np.zeros(n)
    cp[0] = c[0] / b[0]
    dp[0] = d[0] / b[0]
    for i in range(1, n):
        m = b[i] - a[i] * cp[i - 1]
        cp[i] = c[i] / m
        dp[i] = (d[i] - a[i] * dp[i - 1]) / m
    x = np.zeros(n)
    x[-1] = dp[-1]
    for i in range(n - 2, -1, -1):
        x[i] = dp[i] - cp[i] * x[i + 1]
    return x
 
def implicit_diffusion_step(q, r):
    """Backward Euler — solve (I - r·L) q^{n+1} = q^n via Thomas."""
    n = len(q)
    a = np.full(n, -r)
    b = np.full(n, 1 + 2 * r)
    c = np.full(n, -r)
    b[0] = 1 + r        # zero-flux boundary
    b[-1] = 1 + r
    return thomas_solver(a, b, c, q.copy())
 
# 가우시안 펄스로 비교
N = 80
x = np.linspace(0, 1, N)
q0 = np.exp(-((x - 0.5) ** 2) / 0.005)
D, dx = 1.0, 1.0 / N
r_values = [0.4, 0.6, 2.0]
for r in r_values:
    dt = r * dx**2 / D
    qe, qi = q0.copy(), q0.copy()
    for _ in range(60):
        qe = explicit_diffusion_step(qe, r)
        qi = implicit_diffusion_step(qi, r)
    print(f"r={r:.2f} | explicit max={np.max(np.abs(qe)):.2e} | implicit max={np.max(np.abs(qi)):.4f}")
# r=0.40 | explicit max=4.21e-01 | implicit max=0.4456
# r=0.60 | explicit max=1.33e+12 | implicit max=0.3873
# r=2.00 | explicit max=2.78e+38 | implicit max=0.2517

r=0.6에서 양해법은 60스텝 만에 $10^{12}$ 로 폭발한다. 같은 시간 스텝에서 음해법은 매끄럽게 퍼져나간다.

직접 dt를 키워보자#

아래 시뮬레이션에서 직접 조작해보자. 슬라이더로 $r = D\,\Delta t/\Delta x^2$ 을 0.05부터 5까지 바꿔본다.

CFL number r = D·dt / dx²: 0.45Explicit: stable (r ≤ 0.5)Implicit: unconditionally stable

$r$ 을 0.5 아래로 두면 두 곡선이 거의 겹친다. 0.5를 넘기는 순간 주황색(양해법)이 톱니로 떨리기 시작하고, 1을 넘기면 화면 밖으로 발사된다. 시안색(음해법)은 어느 $r$ 에서도 무너지지 않는다. 단, $r$ 이 너무 크면 해상도가 떨어져 펄스가 한 스텝에 너무 많이 퍼진다는 점도 같이 보인다.

음해법이 왜 항상 답이 아닌가#

무조건 안정하다고 음해법이 만능은 아니다. 세 가지 비용이 따라붙는다.

1. 정확도 vs 안정성 — 안정성은 진폭이 폭발하지 않는다는 보장일 뿐, 정확하다는 보장이 아니다. 백워드 오일러는 $\mathcal{O}(\Delta t)$ 정확도다. 정확도가 필요하면 Crank–Nicolson( $\mathcal{O}(\Delta t^2)$ , 무조건 안정)을 쓰는 편이 낫다.

2. 비선형성 — 확산 계수가 해에 의존하면( $D = D(q)$ ) 행렬이 매 스텝마다 바뀌고 비선형 솔버(Picard, Newton)가 필요하다.

3. 다차원 — 2D/3D에서는 더 이상 삼대각이 아니다. 5점/7점 스텐실은 측면 밴드(side band)를 만든다. 이때는 ADI(Alternating Direction Implicit)로 1D 삼대각 두/세 번으로 쪼개거나, 켤레기울기(CG)·BiCGSTAB·멀티그리드 같은 반복 솔버를 쓴다.

세 번째는 압축성에서 비압축성으로 넘어갈 때 그대로 다시 등장한다. 압력 Poisson 방정식 $\nabla^2 P = -\rho \nabla\!\cdot(\mathbf{u}\!\cdot\!\nabla\mathbf{u})$ 가 똑같은 구조의 큰 sparse 행렬이기 때문이다. 음해 확산을 풀어본 사람은 압력 풀이 코드도 같은 도구로 쓴다.

기억할 점#

양해 확산법의 $\Delta t$ 제약은 $\Delta x^2$ 에 묶여 있어 미세 격자에서 치명적이다.
음해법은 미래값을 우변으로 옮겨 무조건 안정성을 산다 — 대가는 행렬 풀이 한 번.
1D 삼대각이면 Thomas 알고리듬으로 $\mathcal{O}(N)$ , 격자 100만 개도 한 스텝에 끝난다.