수치분산

편집 역사 토론
수치해석 전산유체역학 마지막 수정: 2026-08-15 04:37:05

1. 개요[편집]

수치분산
Numerical dispersion
정의이산화가 만들어 낸 가짜 분산 — 위상속도가 파수에 의존하게 되는 현상
진단 도구수치 분산관계 $\omega_h(k)$
이류 중앙차분$\omega_h = \dfrac{c}{\Delta x}\sin(k\Delta x)$
실무 지표파장당 격자점 수 PPW $=\lambda/\Delta x$
흔한 지침2차 도식 10~20 PPW · 4차 8~10 · DRP/스펙트럴 5 내외
증상후행 잔물결, 역행 파속, 전파거리에 비례해 커지는 위상오차

진폭은 멀쩡한데 파가 엉뚱한 데 가 있다. 이건 물리가 아니라 격자가 한 짓이다.

수치분산(numerical dispersion)은 원래 방정식에는 분산이 없는데도 이산화 때문에 위상속도가 파수 kk에 의존하게 되어, 파형이 전파하면서 스스로 흩어지는 현상이다. 물리적으로 비분산인 이류·파동 방정식을 풀 때조차, 격자 위의 도식은 짧은 파를 긴 파보다 느리게(또는 빠르게) 보내 파속(wave packet)을 뭉개고 늘어뜨린다.

진폭이 아니라 위상이 틀린다는 점이 이 오차를 특별히 악질로 만든다. 잔차 그래프는 예쁘게 떨어지고, 에너지 보존도 잘 되고, 그림도 매끈한데, 파면이 있어야 할 자리에 없다. 전산음향학·FDTD 전자기·탄성파 계산처럼 파를 수십~수백 파장 거리로 실어 나르는 문제에서 계산 비용의 상한을 결정하는 것은 안정성도 정확도 차수도 아니고 대개 이 수치분산이다.

이 문서는 분산관계 ωh(k)\omega_h(k)의 관점에 집중한다. 이산 도식이 실제로 푸는 등가 PDE를 테일러 전개로 뽑아내는 절차는 수정 방정식 문서가, 소산과 분산의 대비 및 TVD·한계자 이야기는 수치 소산과 분산 문서가 각각 다룬다.

2. 수치 분산관계 유도[편집]

연속계의 분산 관계는 평면파 u=ei(kx−ωt)u = e^{i(kx-\omega t)}를 방정식에 대입해 ω\omega와 kk의 관계를 읽는 것이다. 이산 도식에도 똑같이 하되, 격자점 xj=jΔxx_j = j\Delta x에서만 대입한다.

2.1. 이류 방정식의 중앙차분[편집]

ut+c ux=0u_t + c\,u_x = 0을 공간만 2차 중앙차분으로 이산화하면

dujdt+c uj+1−uj−12Δx=0\frac{du_j}{dt} + c\,\frac{u_{j+1}-u_{j-1}}{2\Delta x} = 0

이다. uj=ei(kjΔx−ωt)u_j = e^{i(kj\Delta x - \omega t)}를 넣으면 e±ikΔxe^{\pm ik\Delta x}가 튀어나오고, 정리하면

ωh(k)=cΔx sin⁡(kΔx)\omega_h(k) = \frac{c}{\Delta x}\,\sin(k\Delta x)

가 된다. 정확한 관계 ω=ck\omega = ck와 비교하면, kΔx≪1k\Delta x \ll 1일 때 sin⁡(kΔx)≈kΔx\sin(k\Delta x)\approx k\Delta x라 잘 맞지만 파수가 커질수록 어긋난다. 수치 위상속도는

cp(k)c=sin⁡(kΔx)kΔx  ≈  1−(kΔx)26\frac{c_p(k)}{c} = \frac{\sin(k\Delta x)}{k\Delta x} \;\approx\; 1 - \frac{(k\Delta x)^2}{6}

이다. 모든 파수에서 1보다 작다 — 즉 모든 모드가 실제보다 느리게 간다. 그리고 kΔx=πk\Delta x = \pi(파장당 격자점 2개, 나이퀴스트 한계)에서는 ωh=0\omega_h = 0이라, 격자가 표현할 수 있는 가장 짧은 파는 아예 움직이지 않는다.

2.2. 완전 이산 — 파동 방정식과 마법의 시간 스텝[편집]

시간까지 이산화하면 Δt\Delta t가 관계에 들어온다. 1차원 파동 방정식 utt=c2uxxu_{tt}=c^2u_{xx}를 시간·공간 모두 2차 중앙차분(leapfrog)으로 놓고 같은 대입을 하면, 양쪽에서 −4sin⁡2(⋅/2)-4\sin^2(\cdot/2) 꼴이 떨어져 나와

sin⁡2(ωhΔt/2)(cΔt)2=sin⁡2(kΔx/2)Δx2⟺sin⁡ ⁣(ωhΔt2)=ν sin⁡ ⁣(kΔx2)\frac{\sin^2(\omega_h \Delta t/2)}{(c\Delta t)^2} = \frac{\sin^2(k\Delta x/2)}{\Delta x^2} \quad\Longleftrightarrow\quad \sin\!\left(\frac{\omega_h \Delta t}{2}\right) = \nu\,\sin\!\left(\frac{k\Delta x}{2}\right)

가 된다. 여기서 ν=cΔt/Δx\nu = c\Delta t/\Delta x는 CFL 조건의 쿠랑수다. 이 한 줄에서 세 가지가 한꺼번에 읽힌다.

  • ν≤1\nu \le 1이어야 우변의 절댓값이 1을 넘지 않아 ωh\omega_h가 실수로 남는다. 넘으면 ωh\omega_h가 복소수가 되어 지수 증폭 — 안정성 한계가 분산관계에서 그대로 떨어진다.
  • ν<1\nu < 1에서 전개하면 cp/c≈1+(kΔx)224(ν2−1)c_p/c \approx 1 + \frac{(k\Delta x)^2}{24}(\nu^2-1), 즉 역시 느려진다.
  • ν=1\nu = 1이면 ωhΔt/2=kΔx/2\omega_h\Delta t/2 = k\Delta x/2가 되어 모든 파수에서 ωh=ck\omega_h = ck가 정확히 성립한다. 시간 오차와 공간 오차가 정확히 상쇄되는 이 지점이 FDTD에서 말하는 “마법의 시간 스텝”(magic time step)이다.

물론 공짜는 아니다. ν=1\nu=1은 안정 한계 바로 그 지점이라 여유가 0이고, 무엇보다 2·3차원에서는 격자 대각선 방향의 유효 Δx\Delta x가 축 방향과 달라서 모든 전파 방향에 대해 동시에 마법을 걸 수 없다. 그래서 다차원 FDTD에는 위상속도가 전파 각도에 따라 달라지는 수치 이방성(numerical anisotropy)이 남고, 이건 격자를 줄여야만 줄어든다.

3. 뒤처짐, 후행 잔물결, 역행 위상[편집]

분산관계에서 두 가지 속도를 읽어야 한다. 개별 파의 마루가 움직이는 위상속도 cp=ωh/kc_p = \omega_h/k와, 파속(에너지)이 움직이는 군속도 cg=dωh/dkc_g = d\omega_h/dk다. 군속도 쪽이 훨씬 사납다. 이류 중앙차분에서

cg(k)=dωhdk=c cos⁡(kΔx)c_g(k) = \frac{d\omega_h}{dk} = c\,\cos(k\Delta x)

인데, kΔx>π/2k\Delta x > \pi/2(파장당 격자점 4개 미만)이면 **cg<0c_g < 0**이다. 즉 격자가 겨우 표현하는 짧은 파들은 물리와 반대 방향으로 흘러간다. 불연속이나 점 소스처럼 고파수가 잔뜩 들어간 초기조건을 주면, 상류 쪽으로 거슬러 올라가는 유령 파속이 생기고 경계에서 반사되어 계산 영역을 돌아다닌다.1

이 관점이 잘 알려진 현상 하나를 깔끔하게 설명한다. 급경사 부근의 잔물결(wiggle)이 앞에 생기느냐 뒤에 생기느냐는 도식의 위상속도 부호가 결정한다. 락스-벤드로프는 수치 위상속도가 정확값보다 느려서 진동이 불연속 뒤에 끌려오고, 빔-워밍(2차 상류)은 빨라서 진동이 불연속 앞으로 달려 나간다. 코드 결과 그림만 보고도 어느 계열의 도식인지 짐작할 수 있다는 뜻이다.

4. 파장당 격자점(PPW)[편집]

실무에서 쓰는 단위는 kΔxk\Delta x가 아니라 파장당 격자점 수다.

NPPW=λΔx=2πkΔxN_{\mathrm{PPW}} = \frac{\lambda}{\Delta x} = \frac{2\pi}{k\Delta x}

이걸로 위 오차식을 다시 쓰면 숫자가 손에 잡힌다. 2차 중앙차분 이류에서 위상속도 상대오차는

ε=1−cpc≈(kΔx)26=2π23N2≈6.6N2\varepsilon = 1 - \frac{c_p}{c} \approx \frac{(k\Delta x)^2}{6} = \frac{2\pi^2}{3N^2} \approx \frac{6.6}{N^2}

이다. N=10N=10이면 6.6%, N=20N=20이면 1.6%. “2차 도식은 파장당 10~20점”이라는 국룰의 출처가 바로 이 숫자다. 4차 중앙차분이면 ε≈(kΔx)4/30\varepsilon \approx (k\Delta x)^4/30이라 N=10N=10에서 0.5%로 한 자릿수 이상 좋아지고, 그래서 고차 도식은 같은 정확도를 훨씬 성긴 격자로 낸다 — 파를 멀리 보내는 문제에서 고차 도식이 “비싸 보이는데 실제로는 싼” 이유다.

주의할 점은 차수가 같아도 계수가 다르다는 것이다. 같은 2차라도 1차 미분 중앙차분은 ε≈(kΔx)2/6\varepsilon \approx (k\Delta x)^2/6인데, 2차 미분 중앙차분(파동 방정식형)은 ε≈(kΔx)2/24\varepsilon \approx (k\Delta x)^2/24로 네 배 작다. “몇 차 도식이냐”만 물어보고 PPW를 정하면 4배씩 틀린다.

5. 오차의 누적 — 왜 멀리 갈수록 나빠지나[편집]

수치분산이 무서운 진짜 이유는 오차가 전파 거리에 비례해 쌓이기 때문이다. 파수 kk인 모드가 거리 LL을 갔을 때 누적 위상오차는

Φ=k (c−cp) t=2π Lλ ε\Phi = k\,(c - c_p)\,t = 2\pi\,\frac{L}{\lambda}\,\varepsilon

라디안이다. 격자당 오차가 아니라 지나온 파장 수 L/λL/\lambda에 곱해진다. 2차 도식에 N=10N=10이면 ε=0.066\varepsilon=0.066이라 한 파장만 가도 24도가 틀어지고, 열 파장을 가면 위상이 통째로 뒤집힌다. 도메인이 100파장짜리인 항공기 소음이나 고주파 초음파 계산에서 10 PPW를 쓰면, 격자 수렴 시험은 통과하는데 결과는 쓰레기가 나온다.

이 누적을 조건 하나로 정리하면 Φ≲const\Phi \lesssim \text{const}, 즉

Lλ⋅(kΔx)p≲const\frac{L}{\lambda}\cdot(k\Delta x)^{p} \lesssim \text{const}

(pp차 도식). 도메인을 파장 단위로 키우면 kΔxk\Delta x를 고정하는 것만으로는 부족하고 더 줄여야 한다. 주파수 영역에서 이 사실이 나타나는 형태가 헬름홀츠 문제의 오염 효과(pollution effect)로, khkh를 상수로 잡아도 kk를 키우면 오차가 k(kh)2pk(kh)^{2p} 항 때문에 커진다. 정량적 논의와 바부슈카-자우터의 “2차원 이상에서는 오염을 피할 수 없다”는 결과는 헬름홀츠 방정식 문서에 있다. 시간영역의 누적 위상오차와 주파수영역의 오염은 같은 현상을 두 방향에서 본 것이다.

6. 소산과의 구분[편집]

폰 노이만 안정성 해석의 증폭계수 G(k)G(k)로 보면 둘의 분업이 명확해진다.

성분읽는 곳증상
수치소산∣G(k)∣<1\lvert G(k)\rvert < 1진폭 감소, 충격파·와류가 뭉개짐
수치분산arg⁡G(k)≠−ωΔt\arg G(k) \ne -\omega\Delta t진폭은 살고 위치·파형이 틀어짐

상류차분처럼 소산이 있는 도식은 고파수를 눌러 죽여서 분산의 증상(잔물결)을 가려 준다. 그래서 “1차 상류는 분산 문제가 없다”고 착각하기 쉬운데, 사실은 분산도 있고 소산이 그걸 덮은 것뿐이다. 반대로 스펙트럴 방법이나 무소산 중앙차분처럼 진폭을 정확히 보존하는 도식일수록 잘못된 위상이 그대로 살아남아 진동이 화려해진다. 정확한 도식이 더 지저분한 그림을 준다는 이 역설이 초보자를 제일 많이 혼란시킨다.2

7. 어떻게 줄이나[편집]

  • 고차 도식. 가장 정직한 답. ε∼(kΔx)p\varepsilon \sim (k\Delta x)^p에서 pp를 올린다. 다만 스텐실이 넓어져 경계 처리와 병렬 유령셀 폭이 같이 커진다.
  • 컴팩트(파데) 도식. 미분값을 음함수 관계로 묶어 좁은 스텐실로 넓은 스텐실급 해상도를 얻는다. 렐레(Lele, 1992)의 “spectral-like resolution” 계열이 표준 참고문헌이고, 대가는 매 스텝 삼중대각 계를 푸는 것.
  • DRP 도식. 탐과 웹(Tam & Webb, 1993)이 전산음향학용으로 제안한 분산관계 보존(Dispersion-Relation-Preserving) 유한차분. 발상의 전환이 재미있다 — 테일러 차수를 최대로 만드는 대신, 관심 파수 대역 [0,kmax⁡Δx][0, k_{\max}\Delta x]에서 ωh(k)\omega_h(k)와 ckck의 오차 적분을 최소화하도록 스텐실 계수를 최적화한다. 형식 차수는 오히려 내려가지만(7점 스텐실이 6차가 아니라 4차), 쓸 수 있는 파수 대역은 훨씬 넓어져 PPW 5~7에서도 실용 정확도가 나온다. “차수를 팔아 해상도를 산다”는 이 트레이드오프는 이후 최적화 계수 도식 전체의 원형이 됐다.3
  • 스펙트럴/고차 요소. 주기 문제의 푸리에 스펙트럴 방법은 표현 가능한 모든 파수에서 분산관계가 정확하다. 복잡 형상에서는 스펙트럴 요소·불연속 갈러킨이 같은 역할을 한다.
  • 질량행렬 혼합. 유한요소법에는 이 계열만의 우아한 수법이 있다. 1차 요소를 일관 질량행렬로 쓰면 cp/c≈1+(kh)2/24c_p/c \approx 1 + (kh)^2/24로 빨라지고, 집중 질량(lumped)으로 쓰면 1−(kh)2/241 - (kh)^2/24로 느려진다. 부호가 반대이므로 둘을 반반 섞으면 2차 오차항이 상쇄되어 분산 정확도가 4차로 올라간다. 지구물리 탄성파 계산에서 오래 쓰여 온 기법이고, 요즘은 최적 혼합 계수를 스펙트럴 요소 차수별로 구한 결과들이 나와 있다.
  • 격자를 줄인다. 언제나 통하지만, 3차원에서 Δx\Delta x를 절반으로 줄이면 비용이 16배(공간 8배 × 시간 스텝 2배)다. 위의 다섯 가지가 존재하는 이유.

8. 관련 문서[편집]

9. Footnotes[편집]

  1. 트레페선(Trefethen)이 1982년 SIAM Review에 쓴 “Group velocity in finite difference schemes”가 이 관점을 정리한 고전이다. 경계조건이 왜 어떤 도식에서만 불안정을 일으키는지도 결국 “군속도가 어느 방향을 향하느냐”로 설명된다. 위상속도만 보고 있으면 절대 안 보이는 것들이 군속도로 바꿔 보면 한꺼번에 보인다. ↩

  2. “격자를 늘렸더니 진동이 심해졌어요”라는 문의의 절반은 여기서 온다. 격자를 늘려 표현 가능한 파수 대역이 넓어졌는데, 그 대역의 고파수 성분이 죄다 틀린 속도로 굴러다니는 것이다. 나머지 절반은 그냥 초기조건이 격자에 안 맞는 경우다. ↩

  3. 이 발상은 처음 들으면 이단처럼 들린다. 수치해석 수업 내내 “차수가 높을수록 좋다”고 배웠는데 차수를 일부러 낮추라니. 하지만 테일러 차수는 kΔx→0k\Delta x \to 0 극한의 성질이고, 우리가 실제로 계산하는 파는 kΔx≈1k\Delta x \approx 1 근처에 있다. 안 쓰는 극한에서 최적인 것과 쓰는 대역에서 최적인 것 중 무엇을 고를지의 문제일 뿐이다. ↩