FTCS가 떨릴 때 업윈드는 무뎌진다 — 1차 이류 스킴, 인공 점성, 그리고 CFL의 등가교환

N=200, c=0.8이면 한 바퀴는 250스텝이다. 코드의 400스텝은 1.6바퀴라 시작 펄스와 자리가 다르다.

1960년대 초, von Neumann 옆의 사람들이 1차원 선형 이류에 중심차분을 넣었다. 사각 펄스가 조용히 흘러야 하는데, 뒤 가장자리에서 진동이 커지더니 몇 스텝 만에 격자 밖으로 나갔다.

원인은 차수가 아니다. 하류 값을 같은 무게로 끌어오면 인과가 깨진다. 한 칸 업윈드는 그 폭주를 잠재우는 대신 모서리를 죽인다. 그 경계가 CFL이다.

CFL 눈금을 0.05에서 1.5까지 움직이면, 같은 사각 펄스에서 FTCS는 계속 떨고 업윈드는 1 근처에서만 점선과 겹친다.

CFL number c = u·dt / dx: 0.80

centered (FTCS)1st-order upwindLax–Wendroffanalytic shifttry CFL > 1 to break the safe zone

눈금을 0.5에 두면 업윈드가 가장 넓게 퍼진다. 1.05만 넘겨도 업윈드와 Lax–Wendroff가 같이 날아간다. 스텐실이 특성선보다 짧아진 자리다.

가장 자연스러운 차분이 가장 먼저 무너진다#

1차원 선형 이류는 한 줄이다.

\partial_t q + u\,\partial_x q = 0

$q$ 는 흐르는 양, $u>0$ 는 일정한 속도다. 해는 초기 모양이 속도 $u$ 로 평행 이동하는 것이다.

시간에 forward Euler, 공간에 중심차분(2차)을 쓰면 갱신은 이렇다.

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

$c \equiv u\,\Delta t / \Delta x$ 는 Courant 수다. 한 스텝에 파가 건너는 칸 수다. 공간 2차, 시간 1차라 LTE(국소 절단 오차)는 $\mathcal{O}(\Delta t,\Delta x^2)$ 이다. 그런데 $c$ 를 아무리 줄여도 이 스킴은 폭주한다.

$q^n_i = g^n e^{i k x_i}$ 를 넣으면 증폭계수는

g = 1 - i\,c\,\sin(k\Delta x)

이고 $|g|^2 = 1 + c^2 \sin^2(k\Delta x) > 1$ 이다. 모든 모드가 매 스텝 자란다. $|g|$ 가 1을 넘는 모드가 하나면 전체가 죽는다.

코드의 $c=0.8$ 에서 $\sin=1$ 인 모드의 크기는 $\sqrt{1+c^2}=\sqrt{1.64}\approx 1.28$ 이다. 한 스텝에 1.28배, 400스텝은 그 거듭제곱이다. 아래 눈대중의 $1.78^{400}$ 은 이 모드의 인자가 아니다. 식은 1.28을 준다.

정보는 한쪽에서만 흘러온다#

$u>0$ 이면 시각 $t_{n+1}$ 의 $q_i$ 는 상류, 즉 $x<x_i$ 만 봐야 한다. 중심차분은 하류 값 $q_{i+1}^n$ 을 같은 무게로 가져온다. 흐름이 모르는 쪽을 미리 읽는 것이다.

업윈드는 그 방향을 고정한다.

q_i^{n+1} = q_i^{n} - c\,\big(q_i^{n} - q_{i-1}^{n}\big), \quad u>0

증폭계수는

g = 1 - c + c\,e^{-i k \Delta x}

이고 $|g|^2 = 1 - c(1-c)\cdot 4\sin^2(k\Delta x/2)$ 이다. $0 \le c \le 1$ 이면 $|g|\le 1$ 이다.

감쇠 항의 계수 $c(1-c)$ 는 $c=0.5$ 에서 0.25로 최대다. $c=0.8$ 이면 0.16이다. 최단 파장에서 $|g|^2 = 1 - 4\cdot 0.16 = 0.36$ , 즉 $|g|=0.6$ 이다. $c=0.5$ 에서는 그 모드의 $|g|$ 가 0까지 내려간다. $c=1$ 이면 빼는 항이 사라져 모든 모드가 $|g|=1$ 이다. 무뎌짐의 정점은 0.5이지, 안정 구간의 끝이 아니다.

업윈드 = 중심차분 + 보이지 않는 점성#

1차 업윈드 기울기를 쪼개면 항등식이다.

\frac{q_i - q_{i-1}}{\Delta x} = \frac{q_{i+1} - q_{i-1}}{2\Delta x} - \frac{\Delta x}{2}\,\frac{q_{i+1} - 2q_i + q_{i-1}}{\Delta x^2}

앞 항은 중심차분이다. 뒤 항은 2차 미분의 이산형이다. 그래서 업윈드가 실제로 푸는 식은

\partial_t q + u\,\partial_x q = \frac{u\,\Delta x}{2}\,\partial_x^2 q

이다. 이류를 풀려다 점성 $\nu_{\text{num}} = u\Delta x/2$ 를 끼워 넣은 수정 방정식이다. 인공 점성(artificial viscosity), 또는 numerical diffusion이다.

이 점성이 폭주를 죽인다. 대가로 사각 펄스의 모서리가 두께 $\sqrt{4\nu_{\text{num}} t}$ 로 번진다. $\Delta x$ 를 줄이면 점성도 줄지만, 같은 시각까지 스텝이 늘어 번짐이 그대로 남는다.

뒤에 나오는 계수 $u\Delta x(1-c)/2$ 가 이 식과 다른 이유는 $c$ 다. 위 식의 $u\Delta x/2$ 는 $c\to 0$ 쪽이다. $c=0.8$ 이면 $(1-c)=0.2$ 라, 실제로 남는 계수는 그 5분의 1, 즉 $u\Delta x/10$ 이다. $c=1$ 에서 0이 되는 것은 $(1-c)$ 가 붙은 쪽뿐이다.

CFL — 격자가 따라잡지 못하는 순간#

시각 $t_{n+1}$ 의 $q_i$ 는 시각 $t_n$ 의 위치 $x_i - u\Delta t$ 에 있던 값이다. 그 점이 $x_{i-1}$ 과 $x_i$ 사이에 있으면 두 격자로 보간된다. $u\Delta t > \Delta x$ 이면 발은 $x_{i-2}$ 너머다. 스텐실은 거기까지 안 본다.

Courant–Friedrichs–Lewy 조건은 그 기하다.

\Delta t \le \frac{\Delta x}{|u|}

수치 의존 영역(스텐실)이 물리 의존 영역(특성선)을 포함해야 한다. CFL은 explicit 스킴의 필요 조건이지 충분 조건이 아니다. FTCS는 CFL 안에서도 터진다. 폭발은 스텐실 길이 문제가 아니라 인과가 깨진 문제다.

NumPy 50줄로 따라가는 세 스킴#

같은 격자, 같은 사각 펄스, 같은 CFL이다.

import numpy as np
 
N = 200
dx = 1.0 / N
u = 1.0
cfl = 0.8
dt = cfl * dx / u
n_steps = 400
 
x = (np.arange(N) + 0.5) * dx
q0 = ((x > 0.2) & (x < 0.4)).astype(float)  # 사각 펄스
 
def ftcs_advect(q, c):
    return q - 0.5 * c * (np.roll(q, -1) - np.roll(q, 1))
 
def upwind_advect(q, c):
    return q - c * (q - np.roll(q, 1))  # u > 0
 
def lax_wendroff_advect(q, c):
    flux_diff = 0.5 * c * (np.roll(q, -1) - np.roll(q, 1))
    visc = 0.5 * c * c * (np.roll(q, -1) - 2 * q + np.roll(q, 1))
    return q - flux_diff + visc
 
q_ftcs, q_up, q_lw = q0.copy(), q0.copy(), q0.copy()
for _ in range(n_steps):
    q_ftcs = ftcs_advect(q_ftcs, cfl)
    q_up = upwind_advect(q_up, cfl)
    q_lw = lax_wendroff_advect(q_lw, cfl)
 
# 분석해: 한 바퀴 돌아 제자리
q_exact = q0
print("FTCS  L1 err:", np.abs(q_ftcs - q_exact).sum() * dx)
print("Upwind L1 err:", np.abs(q_up - q_exact).sum() * dx)
print("LW    L1 err:", np.abs(q_lw - q_exact).sum() * dx)

cfl=0.8, 400스텝에서 본문이 적어 둔 눈대중은 이렇다.

FTCS: 진폭이 약 $1.78^{400}$ 배로 커져 NaN 또는 inf가 된다.
Upwind: 사각 펄스가 가우시안으로 번져 L1 오차가 0.04 근처에 머문다.
Lax–Wendroff: 진폭은 보존되지만 모서리 뒤로 잔물결(Gibbs)이 남는다.

c=0.5가 업윈드를 가장 퍼뜨린다#

$c(1-c)$ 가 0.5에서 최대라는 사실이 화면의 번짐과 맞는다. $c=1$ 에서는 $(1-c)$ 가 0이라 1차 업윈드가 그 사각을 그대로 한 칸씩 민다. $c=0.8$ 은 그 사이의 0.16이다. 안정하다고 덜 무딘 것은 아니다. 0.8은 0.5보다 점성이 작을 뿐이다.

1차로는 부족하지만 2차는 또 다른 골치다#

업윈드 점성은 충격 근처의 가짜 진동을 눌러 준다. 매끈한 구간에서도 같은 항이 켜져 있어, 오차는 2차 미분 크기에 묶인다. 그래서 실무 코드는 둘 중 하나로 간다.

더 높은 차수(Lax–Wendroff, MUSCL, WENO)를 쓰되, 충격 근처에서만 차수를 내리는 플럭스 리미터 또는 슬로프 리미터를 붙인다. 계단에서 Lax–Wendroff가 떠는 자리를 리미터로 누르는 계산이 그 선택이다.
셀 면에서 Riemann 해결자(Godunov, HLLC, Roe)로 인과를 직접 푼다.

둘 다 "점성을 어디서는 받고 어디서는 거절할지"를 정하는 장치다. 1차 업윈드는 모든 칸에서 받는다. 단순하고, 안전하고, 가장 많이 퍼진다.

400스텝을 한 바퀴라고 적으면 실패한다#

$N=200$ , $u=1$ , $c=0.8$ 이면 $\Delta t = c/N = 0.004$ 다. 도메인 길이 1을 한 바퀴 도는 스텝 수는 $1/\Delta t = 250$ , 곧 $N/c$ 다. 400스텝의 이동 거리는 $400\times 0.004 = 1.6$ 이다. 펄스 폭 0.2짜리가 0.6만큼 밀리면 $[0.2,0.4]$ 와 $[0.8,1.0]$ 은 겹치지 않는다. 날카로운 두 사각의 L1은 0.4다.

본문의 업윈드 L1 0.04는 그 비교가 될 수 없다. 0.04를 한 바퀴의 번짐이라고 읽으려면 스텝이 250이어야 한다. 1차원 업윈드가 $c=1$ 까지 사는 이 상한이, 옮겨지는 것이 계면이 되면 왜 훨씬 작은 CFL에서 멈추는지가 다음 질문이다. 그 계면 이류에서 0.05가 남은 자리를 보면, 스텐실 길이와 계면 구속이 어디서 갈라지는지 보인다.