격자를 8배 조여도 2.4674는 그대로였다 — 1807년 푸리에 급수가 남긴 두 상수

격자를 세분하면 L2 오차는 줄지만 톱니 모드의 감쇠는 언제나 π²/4배 느리다. 세분으로 못 고치는 오차가 따로 있다.

격자를 8배 조였더니 한 숫자만 꿈쩍하지 않았다#

1차원 열전도를 격자 64개로 풀었다. 그다음 128, 256, 512로 올렸다. L2 오차는 예상대로 착실히 줄었다.

그런데 표에 같이 찍어 둔 숫자 하나가 움직이지 않았다. 2.4674. 네 번 모두 소수점 넷째 자리까지 같았다.

격자 간격 $\Delta x$ 를 8분의 1로 줄였는데 왜 어떤 오차는 그대로인가. 이 글은 그 숫자의 정체를 따라간다. 답은 $\pi^2/4$ 이고, 출발점은 1807년 파리다.

1807년, 심사위원 셋이 걸고넘어진 것은 수렴이었다#

푸리에(Joseph Fourier)는 1807년 고체 열전도에 관한 논문을 파리 아카데미에 제출했다. 고체 안의 온도 분포가 사인과 코사인의 합으로 적힌다는 내용이었다.

심사위원은 Lagrange, Laplace, Legendre 셋이었다. 논문은 반려됐다. 이유는 물리가 아니라 수학이었다. 급수가 수렴하는지 증명되지 않았다는 것이다.

Lagrange가 직접 증명을 시도했다. 성공하지 못한 채 1813년 세상을 떠났다. 증명은 한참 뒤 푸리에의 제자 Dirichlet이 채웠다.

심사위원들이 까다로웠던 게 아니다. 그들이 본 것은 실제로 거기 있었다. 불연속이 있는 함수에 사인 급수를 씌우면 항을 아무리 늘려도 사라지지 않는 것이 남는다.

폭 0.5 사각파에서 슬라이더 N은 남기는 항 수다. 올리면 봉우리는 점프 쪽으로 붙고, 높이는 거의 안 움직인다.

terms N = 41Fejer averagezoom near jump

max overshoot = 0.000%  · peak at x = 0.0000 (jump at x = 0.25)  · distance to jump = -0.2500

N을 키우면 봉우리가 불연속에 바짝 붙는다. 높이는 8.95% 근처에 머문다. Fejer average를 켜면 봉우리는 사라지고 계단이 퍼진다.

사각파를 사인으로 쌓으면 8.9490%가 남는다#

주기 구간에 폭 0.5짜리 사각파를 놓자. 푸리에 계수는 손으로 적힌다.

f(x) = \frac{1}{2} + \sum_{m=1}^{\infty} \frac{2}{\pi m}\sin\!\left(\frac{m\pi}{2}\right)\cos\big(2\pi m (x - 0.5)\big)

$m$ 은 항 번호, $x$ 는 위치다. 홀수 $m$ 만 살아남는다. 이 급수를 $N$ 항에서 자르면 불연속점 바로 옆에서 부분합이 1을 넘어선다.

이것이 Gibbs 현상이다. 극한값은 8.9490%로, 항수와 무관한 상수다. 정밀도를 사서 없앨 수 있는 오차가 아니다.

열방정식이 사인 급수를 좋아하는 이유#

푸리에가 사인을 고른 것은 취향이 아니었다. 열방정식

\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}

에서 $\alpha$ 는 열확산계수다. 여기에 $T = e^{ikx}$ 를 넣으면 2계도함수가 $-k^2 e^{ikx}$ 가 된다. 즉 $e^{ikx}$ 는 라플라시안의 고유함수이고 고유값은 $-k^2$ 다.

그래서 각 모드는 서로 간섭하지 않고 혼자 $e^{-\alpha k^2 t}$ 로 줄어든다. 편미분방정식이 모드마다 하나씩의 상미분방정식으로 분해된다.

지수 안의 $k^2$ 가 감쇠 속도를 정한다. 파장이 절반이면 감쇠가 4배 빠르다. 짧은 파동일수록 먼저 지워진다. 확산이 하는 일이 이 한 자리에 몰려 있다.

격자는 $-k^2$ 대신 무엇을 돌려주는가#

이제 같은 계산을 격자 위에서 한다. 2계도함수를 2차 중심차분으로 바꾸고 $T_j = e^{ikj\Delta x}$ 를 넣는다.

\frac{T_{j+1} - 2T_j + T_{j-1}}{\Delta x^2} = -\frac{4}{\Delta x^2}\sin^2\!\left(\frac{k\Delta x}{2}\right) e^{ikj\Delta x}

$j$ 는 격자점 번호다. 이산 연산자도 $e^{ikj\Delta x}$ 를 고유함수로 갖는다. 격자에서도 푸리에 분해가 그대로 통한다는 뜻이다.

달라진 것은 고유값이다. $-k^2$ 가 아니라 $-k'^2$ 이고

k'^2 = \frac{4}{\Delta x^2}\sin^2\!\left(\frac{k\Delta x}{2}\right)

이다. 이 $k'$ 를 유효파수(modified wavenumber)라고 부른다. 스킴이 "내가 생각하는 $k$ 는 이것"이라고 고백하는 값이다.

$\theta = k\Delta x$ 가 작으면 유효 곡선은 $\theta^2$ 에 붙는다. 오른쪽으로 갈수록 벌어진다. 같은 파동인데 격자가 내놓은 쪽이 늦게 죽는다.

나이퀴스트에서 $\pi^2/4$ 가 나오는 자리#

격자가 표현할 수 있는 가장 짧은 파장은 $2\Delta x$ 다. 한 칸 올라가고 한 칸 내려가는 톱니 모양이고, 이때 $\theta = k\Delta x = \pi$ 다.

여기에 위 식을 넣으면

\frac{k'^2}{k^2} = \frac{4\sin^2(\pi/2)}{\pi^2} = \frac{4}{\pi^2} = 0.4053

이다. 역수가 $\pi^2/4 = 2.4674$ . 앞에서 꿈쩍하지 않던 그 숫자다.

$\Delta x$ 가 식에서 약분됐다는 점이 중요하다. 비율은 격자 간격에 의존하지 않는다. 격자를 아무리 조여도 "그 격자의 톱니"는 언제나 2.4674배 늦게 죽는다.

차수를 올려도 사정은 비슷하다. 실제 스텐실로 확인한다.

import numpy as np
 
STENCILS = {
    2: [1.0, -2.0, 1.0],
    4: [-1/12, 4/3, -5/2, 4/3, -1/12],
    6: [1/90, -3/20, 3/2, -49/18, 3/2, -3/20, 1/90],
}
 
def modified_k2(order, theta):
    """중심차분 2계도함수의 유효파수 제곱 (dx=1 기준). 정확한 값은 theta**2."""
    c = STENCILS[order]
    h = len(c) // 2
    sym = sum(c[h + j] * np.cos(j * theta) for j in range(-h, h + 1))
    return -sym
 
theta = np.array([np.pi/8, np.pi/4, np.pi/2, 3*np.pi/4, np.pi])
print("theta/pi   exact k^2    order2    order4    order6    (ratio)")
for th in theta:
    ex = th**2
    row = [modified_k2(o, th) for o in (2, 4, 6)]
    print(f"{th/np.pi:6.3f}  {ex:10.5f}  " + "  ".join(f"{v:7.4f}" for v in row)
          + "   " + " ".join(f"{v/ex:6.4f}" for v in row))
 
print()
for o in (2, 4, 6):
    r = modified_k2(o, np.pi) / np.pi**2
    print(f"order {o}: Nyquist ratio = {r:.4f}   decay is {1/r:.3f}x slower")

theta/pi   exact k^2    order2    order4    order6    (ratio)
 0.125     0.15421   0.1522   0.1542   0.1542   0.9872 0.9997 1.0000
 0.250     0.61685   0.5858   0.6144   0.6166   0.9496 0.9960 0.9996
 0.500     2.46740   2.0000   2.3333   2.4222   0.8106 0.9457 0.9817
 0.750     5.55165   3.4142   4.3856   4.8278   0.6150 0.7900 0.8696
 1.000     9.86960   4.0000   5.3333   6.0444   0.4053 0.5404 0.6124
 
order 2: Nyquist ratio = 0.4053   decay is 2.467x slower
order 4: Nyquist ratio = 0.5404   decay is 1.851x slower
order 6: Nyquist ratio = 0.6124   decay is 1.633x slower

표의 $\theta$ 와 차분 차수를 다이얼로 그대로 돌려 보자. 차수를 올리면 저파수 쪽이 $\theta^2$ 에 붙고, 톱니 쪽 비는 조금만 움직인다.

theta/pi = 1.00

k'^2 / k^2 = 0.4053  · grid decays 2.467x slower  · a*t/dx^2 = 0.000

$\theta$ 를 $\pi$ 까지 밀면 격자 복사본이 정확한 해보다 늦게 줄어든다. 6차로 올려도 나이퀴스트 비는 0.6124, 감쇠는 1.633배 느리다.

한 칸에 4개 이상 들어가는 파동( $\theta/\pi \le 0.25$ )은 차수를 올리면 거의 완벽해진다. 톱니는 그렇지 않다. 6차로 올려도 0.6124다.

같은 사각파를 격자 넷에 걸어봤다#

이제 실제로 풀어본다. 시간 적분 오차가 섞이지 않게 시간 방향은 해석적으로 적분하고, 공간 이산화 차이만 남긴다.

import numpy as np
 
L, ALPHA, TEND = 1.0, 1.0, 2.0e-3
 
def tophat_coeffs(m):
    """주기 [0,1)에서 0.25 <= x < 0.75 사각파의 복소 푸리에 계수."""
    m = np.asarray(m, dtype=float)
    out = np.zeros(m.shape, dtype=complex)
    out[m == 0] = 0.5
    mm = m[m != 0]
    out[m != 0] = (np.exp(-2j*np.pi*mm*0.75) - np.exp(-2j*np.pi*mm*0.25)) / (-2j*np.pi*mm)
    return out
 
def gibbs_overshoot(nterm):
    """부분합 2N+1항의 최대 초과량(%). 튀는 지점 근처만 조밀하게 본다."""
    x = np.linspace(0.25, 0.25 + 4.0/nterm, 40001)
    m = np.arange(-nterm, nterm + 1)
    c = tophat_coeffs(m)
    s = (c[None, :] * np.exp(2j*np.pi*m[None, :]*x[:, None])).sum(axis=1).real
    return (s.max() - 1.0) * 100.0
 
def heat_error(ngrid, smooth):
    """열방정식을 시간에 대해 정확히 적분. 공간만 2차 중심차분 대 연속 정확해."""
    dx = L / ngrid
    x = np.arange(ngrid) * dx
    m = np.fft.fftfreq(ngrid, d=1.0/ngrid)
    k = 2*np.pi*m
    kmod2 = (4.0/dx**2) * np.sin(k*dx/2.0)**2          # 격자가 내놓는 k^2
    mref = np.arange(-8192, 8193)
    if smooth:
        t0 = np.cos(2*np.pi*x) + 0.3*np.cos(6*np.pi*x)
        ref = (np.cos(2*np.pi*x)*np.exp(-ALPHA*(2*np.pi)**2*TEND)
               + 0.3*np.cos(6*np.pi*x)*np.exp(-ALPHA*(6*np.pi)**2*TEND))
    else:
        t0 = ((x >= 0.25) & (x < 0.75)).astype(float)
        c = tophat_coeffs(mref)*np.exp(-ALPHA*(2*np.pi*mref)**2*TEND)
        ref = (c[None, :]*np.exp(2j*np.pi*mref[None, :]*x[:, None])).sum(axis=1).real
    num = np.fft.ifft(np.fft.fft(t0)*np.exp(-ALPHA*kmod2*TEND)).real
    nyq = ngrid//2
    return np.sqrt(np.mean((num-ref)**2)), 1.0/(ALPHA*k[nyq]**2), 1.0/(ALPHA*kmod2[nyq])
 
print("N_terms   Gibbs overshoot[%]")
for n in (16, 32, 64, 128, 256, 512):
    print(f"{n:6d}    {gibbs_overshoot(n):9.5f}")
 
for smooth in (True, False):
    print(f"\ninitial condition: {'smooth cosine' if smooth else 'square wave (discontinuous)'}")
    print("grid      L2 error     order     sawtooth decay time exact/grid    ratio")
    prev = None
    for ng in (64, 128, 256, 512):
        l2, t_ex, t_gr = heat_error(ng, smooth)
        order = "  -  " if prev is None else f"{np.log2(prev/l2):5.2f}"
        print(f"{ng:4d}   {l2:.4e}    {order}     {t_ex:.3e} / {t_gr:.3e}   {t_gr/t_ex:.4f}")
        prev = l2

N_terms   Gibbs overshoot[%]
    16      9.01421
    32      8.96527
    64      8.95306
   128      8.95000
   256      8.94924
   512      8.94905
 
initial condition: smooth cosine
grid      L2 error     order     sawtooth decay time exact/grid    ratio
  64   5.3683e-04      -       2.474e-05 / 6.104e-05   2.4674
 128   1.3424e-04     2.00     6.184e-06 / 1.526e-05   2.4674
 256   3.3562e-05     2.00     1.546e-06 / 3.815e-06   2.4674
 512   8.3907e-06     2.00     3.865e-07 / 9.537e-07   2.4674
 
initial condition: square wave (discontinuous)
grid      L2 error     order     sawtooth decay time exact/grid    ratio
  64   2.3383e-02      -       2.474e-05 / 6.104e-05   2.4674
 128   1.1673e-02     1.00     6.184e-06 / 1.526e-05   2.4674
 256   5.8342e-03     1.00     1.546e-06 / 3.815e-06   2.4674
 512   2.9168e-03     1.00     3.865e-07 / 9.537e-07   2.4674

세분하면 줄어드는 오차, 세분해도 안 줄어드는 오차#

매끈한 초기조건에서 수렴차수는 정확히 2.00이다. 스킴은 약속대로 2차다.

같은 스킴에 사각파를 주면 1.00으로 떨어진다. 불연속이 있으면 푸리에 계수가 $1/m$ 로만 줄어들고, 이 느린 꼬리가 앨리어싱으로 되돌아온다. 스킴의 차수가 아니라 해의 매끄러움이 수렴차수를 정한다. 계단 옆의 진동을 죽이는 선택은 flux limiter로 2차 TVD를 세우는 일과 맞닿아 있다.

오른쪽 끝 비율은 여덟 줄 모두 2.4674다. L2 오차는 8배 줄었는데 이 값은 소수점 아래까지 같다.

두 현상을 섞으면 안 된다. L2가 줄어드는 것은 톱니 모드의 진폭이 작아지기 때문이지, 톱니를 다루는 실력이 늘어서가 아니다. 격자를 세분하면 문제를 더 잘 푸는 게 아니라 문제를 옮긴다.

격자 64에서 매끈한 코사인의 L2는 5.3683e-04다. 512로 여덟 배 조이면 8.3907e-06이 된다. 나누면 약 64배다. 로그를 취하면 차수 2.00이 그대로 나온다. 사각파는 같은 구간에서 2.3383e-02에서 2.9168e-03으로 떨어진다. 여덟 배, 차수 1.00이다. 톱니의 감쇠 시간은 격자 64에서 정확해 2.474e-05, 격자 값 6.104e-05다. 512에서는 3.865e-07과 9.537e-07이다. 비는 2.4674에서 소수점 넷째 자리까지 안 움직인다. 나이퀴스트의 정확한 $k^2$ 는 9.86960이고 2차 스텐실은 4.0000만 돌려준다. $4/\pi^2 = 0.4053$ 이다. $\Delta x$ 가 약분된 식은 격자를 조여도 같은 숫자를 찍을 수밖에 없다.

Gibbs 쪽도 같은 종류다. 항 16개에서 초과는 9.01421%다. 512개로 올리면 8.94905%다. 차이는 0.065%p 남짓이다. $\theta/\pi = 0.125$ 줄의 2차는 0.9872라서, 한 파장에 격자 16개가 들어가면 이미 1% 안쪽이다. 6차는 나이퀴스트에서 6.0444를 내서 비는 0.6124, 감쇠는 1.633배 느리다. 차수를 올려 사는 것은 저파수 자릿수이지 톱니의 운명이 아니다.

격자 64의 사각파 L2 2.3383e-02를 매끈한 쪽 5.3683e-04로 나누면 약 44배다. 같은 격자, 같은 스킴인데 초기조건만 다르다. 512에서는 2.9168e-03 대 8.3907e-06, 약 348배다. 세분할수록 불연속 쪽이 상대적으로 커진다. 차수 1과 차수 2가 격자 비의 한 제곱만큼 벌어지기 때문이다. 512/64 = 8이니 오차 비는 약 8배 더 벌어진다. 44 곱하기 8은 352이고, 표의 348과 맞다.

6차로 올려도 0.6124인 이유, 그리고 필터#

차수를 올리는 것은 저파수에서만 남는 장사다. 표의 $\theta/\pi = 0.125$ 줄에서 6차는 1.0000이지만, 나이퀴스트에서는 0.6124다. 어떤 유한 스텐실도 톱니를 정확히 다루지 못한다. $2\Delta x$ 파동은 격자 위에서 미분 자체가 모호하기 때문이다.

그래서 실무는 두 갈래로 간다. 하나는 톱니가 생기지 않게 막는 것이다. 엇갈린 격자나 Rhie–Chow 보간이 그 일을 한다. 오차가 큰 칸만 조여도 그 칸의 톱니 비는 $\pi^2/4$ 로 남는다. 쿼드트리가 격자를 스스로 키우는 방식도 그 비를 지우지는 못한다.

다른 하나는 생긴 톱니를 명시적으로 지우는 것이다. 스펙트럴 필터, 인공점성, WENO의 가중치가 모두 여기 속한다. 위 viz의 Fejer average가 가장 단순한 형태다. Gibbs 봉우리는 사라지지만 계단이 뭉개진다는 대가도 같이 보인다.

어느 쪽이든 공짜는 아니다. 고를 수 있는 것은 오차를 없애는가가 아니라 어디에 둘 것인가다.

Dirichlet이 1829년에 채워 넣은 칸#

Lagrange의 반대는 형식논리가 아니었다. 그는 급수가 불연속을 만났을 때 무슨 일이 벌어지는지를 의심했고, 그 자리에 실제로 8.9490%가 있었다.

Dirichlet은 1829년 조건부 수렴을 증명하며 그 칸을 채웠다. 각 점에서 수렴한다, 다만 균등하지는 않다. Gibbs 봉우리는 반례가 아니라 그 "균등하지 않다"의 모습이었다.

오늘 격자 수렴 시험을 돌릴 때 우리는 같은 질문을 되풀이한다. 이 오차는 세분하면 줄어드는 종류인가, 아니면 $\pi^2/4$ 처럼 격자 간격이 약분되어 사라지지 않는 종류인가. 수렴차수 표 한 줄은 앞의 것만 답해 준다.

수렴차수가 2.00으로 찍혔다고 안심한 적이 있다면, 유효파수를 한 번 그려 볼 값어치가 있다. 스킴이 톱니에 대해 무엇을 고백하는지는 거기에만 적혀 있다.