꺾인 현은 멀쩡했고 잘린 현은 9%가 튀었다 — 1747년 파동방정식 논쟁의 수치적 결말

해를 무엇으로 표현하느냐가 오차의 종류를 정한다. 진행파는 모양을 옮기고, 모드 합은 모서리에서 진동한다.

1747년, 현 하나를 두고 세 사람이 갈라졌다#

달랑베르는 1747년 진동하는 현의 방정식을 세웠다. 인류 최초의 편미분방정식이었다. 그런데 다툼이 붙은 것은 방정식이 아니라 해의 생김새였다. 오일러와 다니엘 베르누이가 각각 다른 답을 들고 나왔고, 세 사람은 30년 넘게 물러서지 않았다.

결론부터 적으면 셋 다 옳았다. 답이 갈라지는 것은 급수를 유한한 항에서 자를 때다. 이 글은 그 논쟁의 두 표현을 같은 코드에 나란히 세우고, 초기조건의 매끄러움이 오차의 종류를 어떻게 바꾸는지 확인한다. 스펙트럴 기법이나 고차 스킴에서 만나는 9%짜리 진동의 출처가 바로 여기다.

달랑베르는 아무것도 전개하지 않았다#

장력 $T$ , 선밀도 $\rho$ 인 현의 미소 구간에 뉴턴 제2법칙을 적용하면 이렇게 된다.

u_{tt} = c^2 u_{xx}, \qquad c = \sqrt{T/\rho}

$u$ 는 현의 가로 변위, $c$ 는 파동 속도다. 달랑베르는 좌표를 $\xi = x - ct$ , $\eta = x + ct$ 로 바꾸면 방정식이 $u_{\xi\eta} = 0$ 이 된다는 것을 알아챘다. 두 번 적분하면 해가 바로 나온다.

u(x,t) = \tfrac{1}{2}\left[\,F(x - ct) + F(x + ct)\,\right]

$F$ 는 초기 변위 $f(x)$ 를 양 끝점에 대해 홀함수로, 주기 $2L$ 로 연장한 함수다. 고정단 조건은 따로 넣지 않는다. 홀함수 연장 자체가 벽에서 부호를 뒤집어 주기 때문이다. 전개도, 계수도, 주파수도 없다. 초기 모양을 절반으로 갈라 양쪽으로 옮기는 것이 전부다.

아래 시뮬레이션에서 직접 조작해보자.

speed0.35

Turn halves off and you see only the string. Turn it on and the same picture is two rigid copies at half height sliding in opposite directions — no frequencies anywhere. Switch to jump: the corners stay razor sharp forever, because transport never smooths anything. The clock reads t·c/L = 0.00; at 2 the shape is exactly back.

halves 버튼으로 두 진행파를 껐다 켜 보면, 흰 선 하나가 사실은 반 높이 복사본 두 개의 합이라는 것이 드러난다. jump 모양에서 모서리가 끝까지 날카롭게 유지되는 점을 눈여겨보자. 이류(移流)는 무엇도 뭉개지 않는다.

베르누이는 같은 답을 sine으로 다시 썼다#

다니엘 베르누이의 반론은 음악에서 나왔다. 현 하나는 기본음과 배음을 동시에 낸다. 그러니 해도 정상파의 중첩이어야 한다는 것이다.

u(x,t) = \sum_{n=1}^{\infty} b_n \sin\!\frac{n\pi x}{L}\cos\!\frac{n\pi c t}{L}, \qquad b_n = \frac{2}{L}\int_0^L f(x)\,\sin\!\frac{n\pi x}{L}\,dx

$b_n$ 은 $n$ 번째 모드의 진폭, $n\pi c/L$ 은 그 모드의 각진동수다. 모드 주파수가 $1:2:3$ 의 정수비라는 사실은 피타고라스 이래 알려져 있었다. 당시 파리에서는 라모의 화성론을 두고 논쟁이 한창이었고, 달랑베르 본인도 뉴턴 역학으로 화성학을 옹호하던 쪽이었다. 음계의 물리적 근거가 이 방정식에서 나온 셈이다.

오일러가 끼어든 지점은 "함수란 무엇인가"였다#

오일러는 달랑베르 편에 섰지만 이유가 달랐다. 현을 손가락으로 뜯으면 초기 모양은 꺾인 삼각형이다. 그 지점에서 $u_{xx}$ 가 존재하지 않는다. 달랑베르는 그런 곡선을 애초에 해로 인정하지 않으려 했다. 하나의 해석적 식으로 쓸 수 있는 곡선만 함수라고 보았기 때문이다. 오일러는 손으로 그린 임의의 곡선도 초기조건이 될 수 있다고 맞섰다.

베르누이는 한 발 더 나갔다. 어떤 곡선이든 sine의 무한합으로 쓸 수 있다는 주장이었다. 당시에는 근거 없는 낙관으로 보였다. 이 주장은 1822년 푸리에에 와서야 정리된다. 달랑베르가 원기둥 주위 포텐셜 유동에서 저항이 0으로 나오는 결과를 받아들여야 했던 것처럼, 여기서도 그의 직관은 절반만 맞았다.

Python으로 두 표현을 같은 시각에 세워 봤다#

세 가지 초기 변위를 준비했다. 점프가 있는 구간 함수, 꺾임만 있는 삼각형, 매끄러운 종 모양이다. 같은 시각 $tc/L = 0.15$ 에서 달랑베르 해와 $N$ 항 sine 급수를 비교한다.

import numpy as np
 
L, c, T = 1.0, 1.0, 0.15
X = np.linspace(0.0, L, 2001)
 
def hat_profile(x, a=0.35, b=0.55):
    """점프: [a,b] 구간만 1, 나머지 0 — 해머로 때린 구간"""
    return np.where((x >= a) & (x <= b), 1.0, 0.0)
 
def kink_profile(x, a=0.45):
    """꺾임: 연속이지만 x=a에서 기울기가 튄다 — 현을 뜯은 모양"""
    return np.where(x < a, x / a, (L - x) / (L - a))
 
def bell_profile(x, a=0.45, s=0.055):
    """매끄러움: 무한 번 미분 가능한 종 모양"""
    return np.exp(-((x - a) ** 2) / (2 * s ** 2))
 
def odd_extend(f0, xq):
    """고정단이 요구하는 홀함수·주기 2L 연장 위에서 보간"""
    xs = np.mod(xq, 2 * L)
    sign = np.where(xs > L, -1.0, 1.0)
    xs = np.where(xs > L, 2 * L - xs, xs)
    return sign * np.interp(xs, X, f0)
 
def dalembert_wave(f0, t):
    """달랑베르 해 — 절반씩 좌우로 옮긴 두 진행파의 합"""
    return 0.5 * (odd_extend(f0, X - c * t) + odd_extend(f0, X + c * t))
 
def modal_wave(f0, t, n_modes):
    """베르누이 해 — sine 모드 n_modes개의 중첩"""
    n = np.arange(1, n_modes + 1)[:, None]
    k = n * np.pi / L
    b = 2.0 / L * np.trapezoid(f0[None, :] * np.sin(k * X[None, :]), X, axis=1)
    return (b[:, None] * np.sin(k * X[None, :]) * np.cos(k * c * t)).sum(axis=0)
 
def overshoot_pct(u, exact):
    """참해 진폭 대비 최대 초과량 (%)"""
    return 100.0 * (u.max() - exact.max()) / (exact.max() - exact.min())
 
for name, f0 in (("jump  ", hat_profile(X)),
                 ("kink  ", kink_profile(X)),
                 ("smooth", bell_profile(X))):
    exact = dalembert_wave(f0, T)
    print(f"[{name}]  t*c/L = {T}")
    for n_modes in (8, 32, 128, 512):
        u = modal_wave(f0, T, n_modes)
        print(f"   N={n_modes:4d}   max|modal - dAlembert| = {np.abs(u - exact).max():.5f}"
              f"   overshoot = {overshoot_pct(u, exact):+6.2f} %")

[jump  ]  t*c/L = 0.15
   N=   8   max|modal - dAlembert| = 0.30480   overshoot = +21.01 %
   N=  32   max|modal - dAlembert| = 0.26454   overshoot = +12.07 %
   N= 128   max|modal - dAlembert| = 0.23847   overshoot =  +9.86 %
   N= 512   max|modal - dAlembert| = 0.18738   overshoot =  +8.94 %
[kink  ]  t*c/L = 0.15
   N=   8   max|modal - dAlembert| = 0.02073   overshoot =  -0.23 %
   N=  32   max|modal - dAlembert| = 0.00639   overshoot =  -0.19 %
   N= 128   max|modal - dAlembert| = 0.00157   overshoot =  -0.03 %
   N= 512   max|modal - dAlembert| = 0.00038   overshoot =  -0.01 %
[smooth]  t*c/L = 0.15
   N=   8   max|modal - dAlembert| = 0.04435   overshoot =  -6.25 %
   N=  32   max|modal - dAlembert| = 0.00000   overshoot =  -0.00 %
   N= 128   max|modal - dAlembert| = 0.00000   overshoot =  +0.00 %
   N= 512   max|modal - dAlembert| = 0.00000   overshoot =  -0.00 %

세 줄이 서로 다른 이야기를 한다. 매끄러운 종 모양은 32항에서 이미 배정밀도 바닥에 닿는다. 꺾인 삼각형은 항을 4배로 늘릴 때마다 오차가 4분의 1이 된다. 1차 수렴이다. 점프가 있는 경우만 오버슛이 사라지지 않고 9% 근처에 주저앉는다.

9%는 항을 늘려도 줄지 않는다#

이 값이 깁스 현상(Gibbs phenomenon)이다. 1848년 윌브러햄이 먼저 발견했고, 1899년 깁스가 다시 확인해 이름이 붙었다. 이론값은 점프 크기의 8.95%다. 위 실행에서 512항이 내놓은 8.94%가 그 값이다.

핵심은 오버슛이 좁아지기만 하고 낮아지지는 않는다는 점이다. 항을 늘리면 진동 구간의 폭은 $1/N$ 로 줄어든다. 그래서 적분값이나 $L^2$ 노름으로 재면 수렴한다. 하지만 최대값으로 재면 수렴하지 않는다. 이 두 노름의 차이가 실무에서는 음수 밀도로 나타난다.

modes N24speed0.25

Drag modes N to the right on jump: the red ringing narrows but its height parks near 9% (now 0.00%), and the right-hand curve flattens out. Switch to kink and then smooth with N untouched — the same slider that bought nothing now buys everything, and the error curve turns into a cliff (max error 0.0000).

modes N 슬라이더를 오른쪽 끝까지 밀어 보자. jump에서는 붉은 진동이 얇아질 뿐 키가 그대로다. 오른쪽 수렴 곡선도 바닥에서 평평해진다. 같은 슬라이더를 smooth에서 움직이면 곡선이 절벽처럼 떨어진다. 슬라이더는 그대로인데 초기조건만 바꾼 것이다.

계수의 감쇠율이 이 차이를 전부 설명한다. 점프가 있으면 $b_n \sim n^{-1}$ , 꺾임만 있으면 $b_n \sim n^{-2}$ , 매끄러우면 어떤 거듭제곱보다도 빠르게 줄어든다. 잘라 버린 꼬리가 오차이므로, 꼬리가 두꺼울수록 잘린 자리가 크게 남는다.

1755년, 같은 사람이 비선형을 만났을 때#

오일러는 1755년 유체 유동을 편미분방정식으로 처음 적었다. 오일러 방정식이다. 파동방정식과 달리 이쪽은 특성선의 기울기가 해에 의존한다. 초기조건이 아무리 매끄러워도 특성선이 교차하면 유한 시간 안에 불연속이 생긴다. 초음속 유동이 상류를 모르는 이유와 같은 구조다.

그래서 압축성 해석에서는 매끄러운 초기조건을 골라도 소용이 없다. 충격파가 스스로 점프를 만들어 내고, 그 순간부터 고차 스킴은 1747년의 문제로 되돌아간다. von Neumann이 1950년에 충격파를 일부러 흐려 놓은 것도 같은 진동을 다루기 위해서였다. TVD 제한자와 WENO 가중치는 충격파 근처에서만 차수를 떨어뜨린다. 9%를 지우는 유일한 방법이 국소적으로 매끄러움을 되찾는 것이기 때문이다.

초기조건의 매끄러움부터 확인하는 습관#

새 스킴에서 진동이 보이면 스킴을 의심하기 전에 초기조건과 경계자료를 먼저 본다. 초기 필드를 셀 단위 상수로 깔았는지, 계면을 점프로 넣었는지, 입구 프로파일이 시간에 대해 $C^0$ 인지 확인한다. 하나라도 걸리면 그 진동은 버그가 아니라 표현 방식의 대가다.

검증 케이스를 고를 때도 같은 기준이 쓰인다. 매끄러운 해로는 설계 차수가 그대로 나온다. 점프가 들어간 순간 최대 노름 수렴은 사라지고 $L^1$ 만 남는다. 1747년의 세 사람은 이 구분을 언어로 갖지 못했다. 우리는 그것을 노름이라고 부른다.