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

양해법의 시간 간격은 격자 간격의 제곱에 잠긴다. 1차원에서 그 자물쇠의 비용은 Thomas 스윕 한 번이다.

밤 11시에 켜 둔 계산이 아침에 NaN이었다. 격자를 두 배로 잘게 쪼개고, 시간 간격은 그대로 뒀다. CFL(Courant–Friedrichs–Lewy) 한 줄을 뺀 대가가 클러스터 한 시간이다.

격자 간격이 절반이면 양해 확산의 시간 간격은 4분의 1이다. 점 1024개짜리도, 터지지 말라는 이유만으로 수십만 스텝을 밟는다. 빠져나오는 줄은 짧다. 우변을 미래 시각으로 옮기고, 그 삼대각을 Thomas로 푼다.

r 슬라이더를 0.05에서 5까지 밀면, 같은 가우시안 펄스가 양해와 음해에서 어디부터 갈라지는지 한 화면에 겹친다.

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

0.5 아래에서는 두 곡선이 거의 포개진다. 0.5를 넘는 순간 주황(양해)이 톱니로 떨리다가 1을 넘기면 화면 밖으로 나가고, 시안(음해)만 선이 남는다.

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$ 가 강제된다.

$1024 \to 2048$ 이면 $\Delta t$ 는 $1/4$ 이다. 한 번 더 쪼개 $4096$ 이면, 1024 기준 $\Delta t$ 는 $1/16$ 이고 스텝은 16배다. 점을 8배로 늘리면 같은 시각까지 스텝이 64배가 된다. 청구서는 셀 개수가 아니라 이 제곱이다.

$r = 0.6$ 이면 최단 파장 인자는 $|1 - 4r| = 1.4$ 다. 60스텝은 그 인자의 60제곱이 먹는 자리다. 출력의 $1.33\times 10^{12}$ 에는 초기 펄스가 그 모드에 실린 몫까지 들어 있다. $r = 0.4$ 에서는 인자가 $0.6$ 이라 양해 피크가 $4.21\times 10^{-1}$ 에 머문다.

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

고치는 줄은 하나다. 우변의 $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}

같은 분석으로 증폭 인자는 $1/(1 + 4r)$ 이고, 양수인 $r$ 에서는 항상 1보다 작다. $\Delta t$ 를 키워도 발산하지 않는다. 무조건 안정(unconditionally stable)이다.

대가는 정확도다. 너무 큰 $\Delta t$ 는 안정하지만 정확하지 않다. 그래도 발산하지 않는다는 보장은, 계산을 밤새 두는 사람한테는 작지 않다.

같은 $r = 0.4$ , 60스텝의 피크는 양해 $4.21\times 10^{-1}$ , 음해 $0.4456$ 이다. 안정한 구간인데도 소수점 세 자리에서 갈린다. 시각은 같고, 차이는 절단오차다.

삼대각 행렬과 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번의 연산이고, 메모리는 세 배열만 있으면 된다. 행렬을 통째로 저장할 필요가 없다.

경계 대각만 $1+r$ 이고 내부는 $1+2r$ 이다. 제로 플럭스라 유령 셀을 경계값과 같게 두면 바깥 이웃이 빠진다. 대각의 $2r$ 에서 $r$ 하나와, 그 이웃으로 가던 $-r$ 이 같이 사라진다. 그래서 코드의 b[0] = 1 + r 이다. 내부 식을 경계에 쓰면 총량이 샌다.

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}$ 로 폭발한다. 같은 시간 스텝에서 음해법은 매끄럽게 퍼져나가고, 피크는 $0.3873$ 이다. $r=2$ 에서는 양해가 $2.78\times 10^{38}$ , 음해 피크는 $0.2517$ 까지 내려앉는다.

출력 세 줄은 같은 시각이 아니다#

$\Delta t = r\,\Delta x^2 / D$ 이고 $N = 80$ , $D = 1$ 이면 60스텝의 물리 시간은 $r$ 에 비례한다. $r = 2$ 는 $r = 0.4$ 보다 다섯 배 먼 시각이다. 피크 $0.4456$ 과 $0.2517$ 을 붙여 음해가 더 죽었다고 하면 거짓이다. 시각을 맞추려면 $r$ 에 스텝 수를 곱한 값이 같아야 하고, $r = 2$ 의 12스텝이 $r = 0.4$ 의 60스텝과 같다. 출력은 그 맞춤을 하지 않았다.

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

무조건 안정하다고 음해법이 만능은 아니다. 비용이 세 갈래로 붙는다.

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

Crank–Nicolson의 최단 파장 인자는 $(1-2r)/(1+2r)$ 이다. $r = 2$ 이면 $-0.6$ 이다. 크기는 1보다 작아 안 터지지만 부호가 음이라 위상이 뒤집히고, $r$ 을 더 올리면 $-1$ 에 붙는다.

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

3. 다차원 — 2D/3D에서는 더 이상 삼대각이 아니다. 5점/7점 스텐실은 측면 밴드(side band)를 만든다. 이때는 ADI(Alternating Direction Implicit)로 1D 삼대각을 두세 번 쪼개거나, Jacobi부터 BiCGSTAB까지 이어지는 반복법으로 넘어간다.

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

큰 r을 정답이라고 적으면 실패한다#

$r = 2$ 의 음해 피크 $0.2517$ 은 정답에 가까워진 기록이 아니다. 물리 시간이 다섯 배이고, 시간 오차는 $\Delta t$ 에 비례한다. 곡선을 남기기 전에 그 스텝의 물리 시간을 같이 적지 않으면 비교가 거짓이 된다.

격자가 2차원이 되면 이 삼대각 코드는 그대로 쓸 수 없다. 행렬을 만들지 않고 스윕만 하는 편이 싸 보일 때, 양해가 터진 뒤에서 행렬 없이 음함수를 푸는 LU-SGS는 그 안정 보장을 어디까지 가져갈까. 스윕 방향이 틀리면 그 보장도 같이 깨지는지, 아직 열려 있다.