격자를 두 배로 촘촘히 했더니 오차가 세 배가 됐다 — Jacobi 반복의 증분 수렴 판정이 숨기는 것

증분이 es보다 작아져도 남은 오차는 es/(1-ρ)이다. Jacobi의 1-ρ는 h²로 줄어서, 같은 es로 격자만 줄이면 답이 오히려 나빠진다. 최적 ω의 SOR은 이 배율을 1/h로 낮춘다.

촘촘한 격자가 더 틀린 답을 냈다#

2003년 어느 학부 전산유체 숙제 리포트에 이런 문장이 있다. "격자 간격이 더 좁아서 계산량은 늘었지만, 그만큼의 정확도를 보여주지는 않았다." 2D 정상 열전도를 Jacobi 반복(이웃 네 점의 평균으로 값을 갈아끼우는 반복법)으로 푼 결과다. 촘촘한 격자가 더 정확해야 한다는 상식과 반대다. 이 글은 그 원인이 격자가 아니라 멈추는 규칙에 있음을 보인다. 끝까지 읽으면 "수렴했다"는 판정을 격자와 함께 어떻게 바꿔야 하는지 알게 된다.

2003년 숙제 리포트의 한 줄#

리포트의 문제는 단순하다. 직사각형 판의 아래 벽은 100 °C, 나머지 세 벽은 0 °C다. 내부 온도는 라플라스 방정식 $\nabla^2 T = 0$ 을 따른다. 5점 차분(중심점과 상하좌우 네 점으로 만든 2차 정확도 근사)으로 이산화하면 연립방정식이 된다. 리포트는 이를 Jacobi, 점 Gauss–Seidel, 선 Gauss–Seidel, PSOR, LSOR로 풀었다.

멈추는 규칙은 모든 방법에 같았다. 한 번 반복할 때 바뀐 양이 허용치 $e_s$ 보다 작으면 멈춘다. 리포트는 $e_s = 10^{-4}$ 을 썼다. 그리고 두 가지를 관찰했다. 격자를 줄이면 오차가 오히려 커졌다. SOR의 $\omega$ 는 0.8, 1.0, 1.2 중 1.2가 가장 빨랐다. 두 관찰은 같은 뿌리에서 나온다.

증분은 오차가 아니다#

반복 $k$ 에서 해를 $T^k$ , 이산식의 참해를 $T^*$ 라 하자. 우리가 원하는 것은 오차 $e_k = T^k - T^*$ 다. 하지만 $T^*$ 를 모르니 증분 $d_k = T^{k+1} - T^k$ 를 대신 잰다. 둘의 관계는 반복 행렬의 스펙트럼 반경(가장 큰 고유값의 절댓값) $\rho$ 가 정한다.

e_{k+1} = \rho\, e_k \quad\Rightarrow\quad d_k = e_{k+1} - e_k = (\rho - 1)\, e_k

|e_k| \approx \frac{|d_k|}{1-\rho}

여기서 $\rho$ 는 가장 느리게 줄어드는 오차 모드의 감소율이다. 몇 번 반복하면 그 모드만 남으므로 이 근사가 잘 맞는다. 결국 $|d_k| < e_s$ 에서 멈추면 남은 오차는 약 $e_s/(1-\rho)$ 다. $\rho$ 가 1에 가까우면 이 배율이 커진다.

아래 시뮬레이션에서 직접 조작해보자. 가장 느린 모드 하나만 떼어낸 그림이다.

cells per side n = 40 (h = 0.0250) stop tolerance es = 1e-3

iterating | rho = 0.99692 | k = 0 | increment = 0.00e+0 | error left = 1.00e+2 | predicted es/(1-rho) = 3.24e-1

$n$ 을 10에서 160으로 올려보자. 파란 증분 곡선은 거의 바닥에 붙는데, 주황 오차 곡선은 그대로 남는다. 오른쪽 그래프에서 멈춘 순간 노란 막대가 그 간격, 즉 $1/(1-\rho)$ 배다.

ρ는 격자가 정한다#

Jacobi의 $\rho$ 는 계산할 수 있다. 크기 $L_x \times L_y$ 판에 간격 $h$ 인 균일 격자라면 다음과 같다.

\rho_J = \tfrac12\left[\cos\frac{\pi h}{L_x} + \cos\frac{\pi h}{L_y}\right] \approx 1 - \frac{\pi^2 h^2}{4}\left(\frac{1}{L_x^2} + \frac{1}{L_y^2}\right)

$h$ 는 격자 간격, $L_x, L_y$ 는 판의 가로·세로 길이다. 핵심은 $1-\rho_J \propto h^2$ 이라는 점이다. $h$ 를 반으로 줄이면 $1-\rho_J$ 는 1/4이 된다. 같은 $e_s$ 로 멈추면 남은 반복 오차가 네 배가 된다.

반면 이산화 오차(5점 차분 자체의 오차)는 $h^2$ 로 줄어든다. 그래서 총오차는 두 항의 합이 된다.

\|T^k - T_{\text{exact}}\| \lesssim C h^2 + \frac{e_s}{1-\rho} \approx C h^2 + \frac{16\, e_s}{5\pi^2 h^2}

마지막 식은 $L_x=1$ , $L_y=2$ 일 때다. 앞 항은 줄고 뒤 항은 는다. $h$ 가 어느 값 아래로 내려가면 격자를 줄일수록 답이 나빠진다. 2003년 리포트가 본 것이 이것이다.

코드: 같은 판, 세 격자, 세 반복법#

리포트의 21×41 격자를 정사각 셀로 보고 판을 1×2로 두었다. 격자는 21×41, 41×81, 81×161 세 단계로 균일하게 줄인다. 멈추는 규칙은 $\max|d_k| < 10^{-3}$ 이다. 리포트의 $e_s$ 는 상대 오차로 보이지만, 여기서는 해석을 쉽게 하려고 절대값을 썼다.

참해는 두 개를 준비한다. 하나는 5점 이산식의 정확해 $T_{\text{discrete}}$ 다. 반복 오차는 이것과 비교한다. 다른 하나는 연속 문제의 푸리에 급수 해다. 총오차는 이것과 비교한다. 아래 벽 양 끝 모서리는 온도가 불연속이라 오차가 크므로, 총오차는 $y \ge L_y/8$ 영역에서만 잰다. SOR은 빨강–검정 순서(체스판처럼 두 색으로 나눠 번갈아 갱신)로 벡터화했다.

import numpy as np
 
LX, LY, TB = 1.0, 2.0, 100.0   # 1 x 2 판, 아래 벽 100 °C, 나머지 벽 0
 
def exact_plate(nx, ny, terms=399):
    # 연속 라플라스 문제의 해석해 (푸리에 급수)
    x = np.linspace(0, LX, nx); y = np.linspace(0, LY, ny)
    X, Y = np.meshgrid(x, y, indexing="ij")
    T = np.zeros_like(X)
    for n in range(1, terms + 1, 2):
        k = n * np.pi / LX
        decay = np.exp(-k * Y) * (1 - np.exp(-2 * k * (LY - Y))) / (1 - np.exp(-2 * k * LY))
        T += 4 * TB / (n * np.pi) * np.sin(k * X) * decay
    return T
 
def discrete_plate(nx, ny):
    # 5점 이산식의 정확해 (x는 사인 변환, y는 점화식)
    T = np.zeros((nx, ny)); i = np.arange(1, nx - 1); J = ny - 1; j = np.arange(ny)
    for p in range(1, nx - 1):
        s = np.sin(p * np.pi * i / (nx - 1))
        c = 2 / (nx - 1) * np.sum(TB * s)
        lam = 2 - 2 * np.cos(p * np.pi / (nx - 1))
        r = (2 + lam + np.sqrt((2 + lam) ** 2 - 4)) / 2
        T[1:-1, :] += np.outer(s, c * (r ** (j - J) - r ** (J - j)) / (r ** -J - r ** J))
    T[1:-1, 0] = TB
    return T
 
def jacobi_sweep(T, w, masks):
    Tn = T.copy()
    Tn[1:-1, 1:-1] = 0.25 * (T[2:, 1:-1] + T[:-2, 1:-1] + T[1:-1, 2:] + T[1:-1, :-2])
    return Tn
 
def rb_sor_sweep(T, w, masks):
    T = T.copy()
    for m in masks:                       # 빨강 셀 먼저, 그다음 검정 셀
        gs = np.zeros_like(T)
        gs[1:-1, 1:-1] = 0.25 * (T[2:, 1:-1] + T[:-2, 1:-1] + T[1:-1, 2:] + T[1:-1, :-2])
        T[m] = (1 - w) * T[m] + w * gs[m]
    return T
 
def omega_opt(nx, ny):
    rho_j = 0.5 * (np.cos(np.pi / (nx - 1)) + np.cos(np.pi / (ny - 1)))
    return rho_j, 2 / (1 + np.sqrt(1 - rho_j ** 2))
 
def solve_until(nx, ny, sweep, w, es):
    T = np.zeros((nx, ny)); T[1:-1, 0] = TB
    I, Jg = np.meshgrid(np.arange(nx), np.arange(ny), indexing="ij")
    inner = (I > 0) & (I < nx - 1) & (Jg > 0) & (Jg < ny - 1)
    masks = [inner & ((I + Jg) % 2 == c) for c in (0, 1)]
    for k in range(1, 100000):
        Tn = sweep(T, w, masks)
        d = np.abs(Tn - T).max(); T = Tn
        if d < es:
            return T, k
 
es = 1e-3                                 # max |T^(k+1) - T^k| < es 이면 멈춘다
print(f"stop rule: max|dT| < {es}   error measured on y >= LY/8 (away from hot corners)")
print(f"{'grid':>8} {'method':>10} {'omega':>6} {'iters':>6} {'es/(1-rho)':>11} {'iter_err':>9} {'disc_err':>9} {'total':>8}")
for nx, ny in [(21, 41), (41, 81), (81, 161)]:
    Td, Ta = discrete_plate(nx, ny), exact_plate(nx, ny)
    band = slice((ny - 1) // 8, None)
    disc = np.abs(Td - Ta)[:, band].max()
    rho_j, w_opt = omega_opt(nx, ny)
    for name, sweep, w, rho in [("Jacobi", jacobi_sweep, 1.0, rho_j),
                                ("SOR", rb_sor_sweep, 1.2, None),
                                ("SOR-opt", rb_sor_sweep, w_opt, w_opt - 1)]:
        T, k = solve_until(nx, ny, sweep, w, es)
        bound = f"{es / (1 - rho):11.3f}" if rho is not None else f"{'-':>11}"
        print(f"{nx:>3}x{ny:<4} {name:>10} {w:6.3f} {k:6d} {bound} "
              f"{np.abs(T - Td).max():9.4f} {disc:9.4f} {np.abs(T - Ta)[:, band].max():8.4f}")

stop rule: max|dT| < 0.001   error measured on y >= LY/8 (away from hot corners)
    grid     method  omega  iters  es/(1-rho)  iter_err  disc_err    total
 21x41       Jacobi  1.000    629       0.130    0.1264    0.1691   0.1525
 21x41          SOR  1.200    256           -    0.0422    0.1691   0.1639
 21x41      SOR-opt  1.780     47       0.005    0.0021    0.1691   0.1689
 41x81       Jacobi  1.000   1831       0.519    0.4922    0.0439   0.4850
 41x81          SOR  1.200    791           -    0.1683    0.0439   0.1616
 41x81      SOR-opt  1.883     90       0.009    0.0033    0.0439   0.0435
 81x161      Jacobi  1.000   4815       2.075    1.8097    0.0112   1.8078
 81x161         SOR  1.200   2260           -    0.6486    0.0112   0.6468
 81x161     SOR-opt  1.940    174       0.017    0.0048    0.0112   0.0107

위 코드를 브라우저에서 직접 돌려보자. 왼쪽은 Jacobi, 오른쪽은 SOR이 같은 판에서 동시에 진행한다.

es = 1e-3 SOR omega = 1.780 sweeps/frame = 4

Jacobi : k = 0 | increment = - | iter_err = - | total err (y ≥ LY/8) = - | es/(1-rho) = 0.130
SOR    : k = 0 | increment = - | iter_err = - | total err (y ≥ LY/8) = - | es/(1-rho) = 0.005
discretization error of this grid (y ≥ LY/8) is the floor no solver can beat

격자 버튼을 21×41에서 81×161로 바꾸면 Jacobi의 total err가 0.15에서 1.81까지 커진다. $\omega$ 를 1.2로 내리면 SOR도 같은 함정에 빠진다. show error field를 켜면 남은 오차가 판 한가운데 큰 혹으로 남아 있는 것이 보인다.

1.81이라는 숫자가 말하는 것#

출력 표에서 Jacobi 줄만 보자. 이산화 오차는 0.169 → 0.044 → 0.011로 정확히 1/4씩 준다. 격자는 제 할 일을 했다. 그런데 반복 오차는 0.126 → 0.492 → 1.810으로 네 배씩 는다. 예측값 $e_s/(1-\rho)$ 0.130, 0.519, 2.075와 거의 같다.

총오차는 둘을 합친 결과다. 21×41에서는 이산화 오차가 지배해서 0.15다. 41×81부터는 반복 오차가 지배한다. 0.49, 1.81로 커진다. 가장 촘촘한 격자의 답이 가장 거친 격자보다 12배 나쁘다. 반복 횟수는 629 → 1831 → 4815로 늘었는데 말이다.

21×41에서 총오차 0.1525가 이산화 오차 0.1691보다 작은 것은 우연이다. 반복 오차의 부호가 이산화 오차를 일부 상쇄했다. 이런 상쇄를 정확도로 착각하면 안 된다.

ω=1.2에서 멈추지 말았어야 했다#

리포트는 SOR의 $\omega$ 를 1.2까지만 시험했다. 그런데 이 문제의 최적값은 훨씬 크다.

\omega_{\text{opt}} = \frac{2}{1+\sqrt{1-\rho_J^2}}, \qquad \rho_{\text{SOR}} = \omega_{\text{opt}} - 1 \approx 1 - c\,h

$\rho_{\text{SOR}}$ 은 최적 $\omega$ 일 때 SOR의 스펙트럼 반경이다. $1-\rho$ 가 $h^2$ 이 아니라 $h$ 에 비례한다. 출력에서 최적 $\omega$ 는 1.780, 1.883, 1.940이다. 반복 횟수는 47, 90, 174로 Jacobi보다 13~28배 적다.

더 중요한 것은 멈출 때 남는 오차다. $e_s/(1-\rho)$ 가 0.005, 0.009, 0.017로 $1/h$ 만큼만 는다. 그래서 총오차가 0.169 → 0.044 → 0.011로 이산화 오차를 그대로 따라간다. 격자를 줄인 만큼 답이 좋아진다. $\omega=1.2$ 는 $\rho$ 를 조금만 줄인다. 그래서 반복 오차가 0.042 → 0.168 → 0.649로 여전히 네 배씩 큰다.

최적 $\omega$ 를 구하는 공식은 이 글처럼 직사각형 판에 균일 격자일 때만 정확하다. 실제 문제에서는 리포트 말처럼 시행착오가 필요하다. 그래도 방향은 분명하다. 1보다 조금 큰 값이 아니라 2에 가까운 값을 먼저 시험해야 한다. 정보 전달 순서가 수렴을 바꾸는 다른 예는 LU-SGS 셀 번호 재정렬 글에서 다뤘다.

멈추는 규칙을 격자와 함께 바꾸는 법#

해법은 세 가지다. 첫째, 증분 대신 오차를 추정한다. 연속한 두 증분의 비 $\hat\rho = |d_k|/|d_{k-1}|$ 로 $\rho$ 를 추정한다. 그다음 $|d_k|\,\hat\rho/(1-\hat\rho) < e_s$ 로 판정한다. 코드 한 줄이면 된다.

둘째, 허용치를 격자에 맞춘다. 반복 오차가 이산화 오차보다 충분히 작으면 된다. Jacobi라면 $e_s \propto h^4$ , 최적 SOR이라면 $e_s \propto h^3$ 로 줄여야 한다. 격자를 바꿀 때마다 $e_s$ 를 그대로 두는 것이 이번 함정의 직접 원인이다.

셋째, 판정량을 잔차로 바꾼다. 잔차 $r = b - AT$ 는 오차에 $A$ 를 곱한 것이다. 이것도 격자에 따라 배율이 변하므로, 초기 잔차 대비 상대값을 쓴다. SIMPLE 류 외부 반복에서도 같은 문제가 생긴다. 언더릴랙세이션이 수렴해를 바꾸는 글과 함께 보면, 멈춤 판정과 이완 계수가 결과에 남기는 흔적을 구분할 수 있다.

격자 수렴 연구라면 이 점이 특히 중요하다. 반복 오차가 섞이면 관측 차수가 2차가 아니게 나온다. 검증 케이스가 실제로 재는 것에서 본 차수 오염과 같은 종류다.

내 솔버의 "converged"를 믿기 전에#

로그에 converged가 찍혔다면 두 숫자를 더 출력해보자. 마지막 두 증분의 비 $\hat\rho$ 와 $1/(1-\hat\rho)$ 다. 이 배율이 수백을 넘으면 허용치는 실제 오차보다 그만큼 낙관적이다. 격자를 한 단계 줄였는데 이 배율이 네 배 뛰었다면 Jacobi 계열의 $h^2$ 함정이다. 2003년 리포트의 관찰은 틀리지 않았다. 다만 범인이 격자가 아니라 멈추는 규칙이었다.