파동방정식

편집 역사 토론
수치해석 물리 시뮬레이션 마지막 수정: 2026-09-02 04:13:27

1. 개요[편집]

라플라스 방정식은 퍼지고, 확산 방정식은 잊어버리고, 파동방정식은 기억한다.

파동방정식(wave equation)은 교란이 유한한 속도 cc 로 형태를 유지한 채 전파하는 현상을 기술하는 2계 선형 쌍곡형 편미분방정식 utt=c2∇2uu_{tt} = c^2 \nabla^2 u 이다. 현의 진동, 공기 중의 소리, 지반의 탄성파, 진공 속 전자기파가 전부 이 한 줄로 환원된다 — 물리량 uu 가 무엇인지와 cc 가 무엇으로 만들어지는지만 다를 뿐이다.

수치해석 입장에서 이 방정식이 특별한 이유는 소산이 없다는 것이다. 확산 방정식은 오차를 스스로 뭉개 없애 주지만, 파동방정식은 격자가 저지른 위상 오차를 파가 지나가는 내내 고스란히 누적해서 들고 간다. 그래서 파동 문제의 계산 비용 상한을 결정하는 것은 대개 안정성이 아니라 수치분산이다.

시간 조화(harmonic) 가정을 넣어 주파수 영역으로 옮긴 형태는 헬름홀츠 방정식이고, 그쪽이 주파수 하나만 필요한 정상상태 문제에서는 훨씬 싸다. 이 문서는 시간영역 쪽을 다룬다.

2. 어디서 나오는가[편집]

같은 방정식이 완전히 다른 세 출발점에서 나온다.

팽팽한 현. 장력 TT, 선밀도 μ\mu 인 현의 미소 요소에 뉴턴 2법칙을 쓰고 기울기가 작다고 가정(sin⁡θ≈tan⁡θ=ux\sin\theta \approx \tan\theta = u_x)하면 곧바로 μutt=Tuxx\mu u_{tt} = T u_{xx}, 즉 c=T/μc = \sqrt{T/\mu} 다. 기타 줄을 조이면(=TT 증가) 음이 올라가고 감아 놓은 굵은 줄이 낮은 음을 내는(=μ\mu 증가) 이유가 이 한 줄에 다 있다.

음향. 정지 매질에서 오일러 방정식을 작은 요동에 대해 선형화하고 등엔트로피 상태식을 쓰면 압력 요동에 대해 ptt′=c02∇2p′p'_{tt} = c_0^2 \nabla^2 p', c0=γp0/ρ0=Ks/ρ0c_0 = \sqrt{\gamma p_0/\rho_0} = \sqrt{K_s/\rho_0} 를 얻는다. 비압축 가정을 넣는 순간 Ks→∞K_s \to \infty, 즉 c0→∞c_0 \to \infty 가 되어 파동이 사라진다 — 비압축 CFD가 소리를 못 내는 것은 모델의 한계가 아니라 정의다. 유동이 있을 때의 지저분함은 전산음향학 쪽 이야기.

전자기. 소스 없는 맥스웰 방정식에서 ∇×\nabla\times 를 한 번 더 취하고 ∇⋅E=0\nabla\cdot\mathbf{E}=0 을 쓰면 각 성분이 Ett=(1/με)∇2E\mathbf{E}_{tt} = (1/\mu\varepsilon)\nabla^2\mathbf{E} 를 만족한다. 맥스웰이 여기서 튀어나온 1/μ0ε01/\sqrt{\mu_0\varepsilon_0} 이 이미 측정돼 있던 광속과 일치하는 것을 보고 “빛은 전자기파다”라고 선언한 것이 1865년이다.1

탄성체. 나비에-코시 운동방정식은 스칼라 하나로 안 떨어지고 헬름홀츠 분해를 거쳐 탄성파의 두 종류, 즉 종파(P파, cp=(λ+2μ)/ρc_p=\sqrt{(\lambda+2\mu)/\rho})와 횡파(S파, cs=μ/ρc_s=\sqrt{\mu/\rho})로 갈라진다. 지진계에 P파가 먼저, S파가 나중에 찍히는 시간차가 진원 거리를 주고, 그 도달시간 데이터를 뒤집어 지구 내부 속도 구조를 복원하는 것이 지진 토모그래피다. 광섬유 안의 유도파도 결국 경계가 있는 매질에서 푼 같은 방정식이다.

3. 달랑베르 해와 특성선[편집]

1차원에서는 연산자가 그냥 인수분해된다.

∂tt−c2∂xx=(∂t−c ∂x)(∂t+c ∂x)\partial_{tt} - c^2\partial_{xx} = \left(\partial_t - c\,\partial_x\right)\left(\partial_t + c\,\partial_x\right)

따라서 일반해는 모양을 유지한 채 오른쪽으로 가는 파와 왼쪽으로 가는 파의 합이다. 초기조건 u(x,0)=f(x)u(x,0)=f(x), ut(x,0)=g(x)u_t(x,0)=g(x) 에 대한 달랑베르 공식(1747)은

u(x,t)=f(x−ct)+f(x+ct)2+12c∫x−ctx+ctg(s) dsu(x,t) = \frac{f(x-ct)+f(x+ct)}{2} + \frac{1}{2c}\int_{x-ct}^{x+ct} g(s)\,ds

이다. 여기서 두 가지가 바로 읽힌다.

  • 의존 영역. (x,t)(x,t) 의 값은 초기 시각의 구간 [x−ct, x+ct][x-ct,\,x+ct] 에만 의존한다. 그 바깥에서 무슨 일이 일어나든 상관없다. 이 유한 전파 속도가 CFL 조건의 물리적 원본이다 — 수치 도식의 의존 영역이 물리의 의존 영역을 덮지 못하면 그 도식은 정답을 알 방법 자체가 없다.
  • 특성선. x∓ct=constx \mp ct = \text{const} 를 따라 정보가 흐른다. 이 선을 좌표축으로 삼아 방정식을 푸는 관점이 특성곡선법이고, 1계 시스템 형태에서는 리만 불변량으로 재등장한다.

차원이 올라가면 이야기가 미묘해진다. 3차원에서는 점 폭발이 정확히 두께가 유지되는 껍질로 퍼져 나가 지나간 뒤엔 아무것도 남기지 않지만(강한 하위헌스 원리, 홀수 차원 3 이상에서만 성립), 2차원에서는 파면이 지나간 뒤에도 꼬리(wake)가 영원히 남아 천천히 잦아든다. 우리가 대화를 할 수 있는 것은 공간이 3차원이기 때문이라는 진술은 농담이 아니다.2

4. 에너지 보존[편집]

유계 영역 Ω\Omega 에서 에너지를 운동 항과 변형 항의 합으로 정의하면

E(t)=12∫Ω(ut2+c2∣∇u∣2)dΩE(t) = \frac{1}{2}\int_\Omega \left( u_t^2 + c^2 |\nabla u|^2 \right) d\Omega

이고, 방정식에 utu_t 를 곱해 부분적분하면

dEdt=c2∮∂Ωut ∂u∂n dS\frac{dE}{dt} = c^2 \oint_{\partial\Omega} u_t \,\frac{\partial u}{\partial n}\,dS

가 된다. 즉 에너지 변화는 전부 경계에서만 일어난다. 고정단(u=0⇒ut=0u=0 \Rightarrow u_t=0)이든 자유단(∂u/∂n=0\partial u/\partial n = 0)이든 우변이 0이라 에너지는 정확히 보존된다. 이것이 파동 솔버를 짤 때의 첫 번째 점검 항목이다 — 감쇠를 안 넣었는데 에너지가 줄면 도식이 소산을 넣고 있는 것이고, 늘면 곧 터진다.

5. 경계조건과 반사[편집]

한쪽 끝에 도달한 펄스가 어떻게 되돌아오는지는 경계조건이 결정한다. 거울상 파원(image source)을 놓고 생각하면 답이 즉시 나온다.

끝단조건거울상반사파계수
고정단디리클레 u=0u=0반대 부호위상 반전R=−1R=-1
자유단노이만 ux=0u_x=0같은 부호그대로R=+1R=+1
임피던스 정합ut+c ux=0u_t + c\,u_x=0없음없음R=0R=0

고정단에서 위상이 뒤집히는 것은 “벽이 현을 붙잡고 있으려면 반대 방향으로 밀어야 한다”는 뉴턴 3법칙의 파동판이다. 두 매질이 임피던스 Z=ρcZ=\rho c 로 만나면 일반적으로 R=(Z2−Z1)/(Z2+Z1)R=(Z_2-Z_1)/(Z_2+Z_1) 이고, 위 세 줄은 Z2→∞Z_2 \to \infty, Z2→0Z_2 \to 0, Z2=Z1Z_2 = Z_1 인 극한이다. 양 끝이 고정된 유한 구간에서 반사가 반복되면 특정 파장만 살아남아 정상파(=모드 해석의 고유모드)가 되고, 구면 경계에서 같은 일을 하면 구면조화 함수가 나온다.

6. 수치해법 — 중앙차분과 CFL[편집]

가장 흔한 도식은 시간·공간 모두 2차 중앙차분을 쓴 립프로그(leapfrog)다. 1차원에서

ujn+1=2ujn−ujn−1+ν2(uj+1n−2ujn+uj−1n),ν=c ΔtΔxu_j^{n+1} = 2u_j^{n} - u_j^{n-1} + \nu^2\left(u_{j+1}^{n} - 2u_j^{n} + u_{j-1}^{n}\right), \qquad \nu = \frac{c\,\Delta t}{\Delta x}

3-레벨 명시적 도식이라 값이 싸고, 이산 에너지를 거의 보존한다는 좋은 성질이 있다.3 폰 노이만 안정성 해석을 넣으면 증폭인자가 크기 1인 복소수 쌍이 될 조건으로

sin⁡2 ⁣(ωΔt2)=ν2sin⁡2 ⁣(kΔx2)  ≤  1⟹ν≤1\sin^2\!\left(\frac{\omega \Delta t}{2}\right) = \nu^2 \sin^2\!\left(\frac{k\Delta x}{2}\right) \;\le\; 1 \quad\Longrightarrow\quad \nu \le 1

이 나온다. dd 차원 등방 격자에서는 ν≤1/d\nu \le 1/\sqrt{d} 로 조여진다. CFL 조건이 부등식으로 정확히 이 형태다.

재미있는 것은 1차원에서 ν=1\nu = 1 로 정확히 맞추면 위 관계식이 ωΔt=kΔx\omega\Delta t = k\Delta x 가 되어 이산해가 정확해와 완전히 일치한다는 점이다(“매직 타임스텝”). 도식이 특성선을 정확히 한 칸씩 밟기 때문인데, 물론 2차원 이상이나 비균일 매질에서는 모든 방향에 대해 동시에 ν=1\nu=1 을 맞출 수 없어 이 마법은 깨진다.

ν<1\nu < 1 이면 이산 위상속도가 cc 보다 느려지고, 그 정도가 kΔxk\Delta x 에 의존한다. 원래 비분산인 방정식이 격자 위에서 분산성이 되는 것이다. 위상속도뿐 아니라 군속도도 틀어지므로 파속(wave packet) 전체가 뒤처지면서 뒤로 잔물결을 흘린다. 실무 국룰인 “파장당 격자점 10~20개”(PPW)는 이 위상 오차를 허용 범위로 누르기 위한 것이며, 오차가 전파 거리에 비례해 쌓이므로 멀리 보낼수록 격자를 더 촘촘히 해야 한다.4 자세한 정량은 수치분산과 수정 방정식 문서에.

고차 정확도로 이 비용을 낮추는 길은 여럿이다. 공간에 4차·6차 유한차분법 스텐실을 쓰거나, 아예 스펙트럴 방법으로 공간 미분을 푸리에 변환에 맡기거나, 복잡 형상이 걸리면 스펙트럴 요소법 계열의 유한요소법을 쓴다. 지진파 계산 커뮤니티가 스펙트럴 요소법으로 대거 넘어간 이유가 정확히 이 분산 예산 문제다.

7. 1계 시스템 — 속도-응력 형식[편집]

2계 방정식을 그대로 푸는 대신, 변수를 늘려 1계 쌍곡형 시스템으로 바꾸는 것이 현대 파동 솔버의 표준이다. 1차원 음향/탄성의 경우 입자속도 vv 와 응력 σ\sigma(또는 압력)에 대해

ρ ∂v∂t=∂σ∂x,∂σ∂t=K ∂v∂x\rho\,\frac{\partial v}{\partial t} = \frac{\partial \sigma}{\partial x}, \qquad \frac{\partial \sigma}{\partial t} = K\,\frac{\partial v}{\partial x}

이고, 두 식을 결합하면 c=K/ρc=\sqrt{K/\rho} 인 파동방정식이 복원된다. 이 형식의 이점은 실무적이다.

  • 엇갈림 격자와 궁합이 좋다. vv 와 σ\sigma 를 반 칸·반 스텝 어긋나게 배치하면(엇갈림 격자) 2차 정확도를 가장 조밀한 스텐실로 얻는다. 전자기의 Yee 격자(FDTD)가 정확히 같은 구조이고, 탄성파 쪽 표준인 비리외(Virieux, 1984·1986) 도식도 마찬가지다.
  • 경계조건이 자연스럽다. 자유표면은 σ=0\sigma=0, 강체면은 v=0v=0 으로 물리량에 직접 걸린다. 2계 형식에서 uxu_x 조건을 유령점으로 구현하는 것보다 깔끔하다.
  • 불연속 물성을 다룰 수 있다. 특성 변수 v±σ/(ρc)v \pm \sigma/(\rho c) 로 리만 문제를 풀면 임피던스가 급변하는 계면에서 반사·투과가 자동으로 나온다. 유한체적·불연속 갤러킨 계열이 이 길을 쓴다.
  • 감쇠 모델을 붙이기 쉽다. 점탄성 완화 메커니즘(SLS)이나 뒤에 나올 PML의 보조 변수가 전부 1계 항으로 들어간다.

8. 무반사 경계 — ABC와 PML[편집]

무한 매질을 유한 격자에 담으면 격자 끝에서 파가 되돌아온다. 계산 영역을 넓히는 것은 언제나 지는 싸움이라, 인공 경계에서 파를 흡수해야 한다.

단방향 파동 방정식(ABC). 위 표의 셋째 줄 ut+cux=0u_t + c u_x = 0 은 오른쪽으로 나가는 파만 통과시키는 1차 흡수 경계조건이다. 엔퀴스트-마이다(1977)가 이를 체계화했고, 전자기 쪽 구현으로 유명한 것이 뮤어(Mur, 1981) 조건이다. 구현은 몇 줄로 끝나지만 치명적 약점이 있다 — 비스듬히 입사하는 파에 대해 반사가 급격히 커진다. 입사각 θ\theta 에서 1차 조건의 반사계수는 대략 ∣R∣≈tan⁡2(θ/2)|R| \approx \tan^2(\theta/2) 수준이라, 정면 입사는 완벽하지만 60∘60^\circ 만 돼도 10% 가까이 되돌아온다.

완전정합층(PML). 베랑제(Bérenger, 1994)의 해법은 경계조건을 고치는 대신 흡수하는 가짜 매질 층을 덧대는 것이다. 핵심 트릭은 좌표를 복소 평면으로 늘이는 것(x→x+iω∫σx dx′x \to x + \frac{i}{\omega}\int \sigma_x\,dx')인데, 그러면 층 안에서 파가 지수적으로 감쇠하면서도 입사각과 주파수에 무관하게 계면 반사가 이론상 0 이 된다. 연속계에서 완벽하다는 이 성질이 이산 격자에서는 유한한 값으로 열화되므로 실무에서는 흡수 계수 σx\sigma_x 를 층 두께에 걸쳐 다항식으로 서서히 키우고, 저주파·스치는 입사에서의 문제는 CFS-PML로 보완한다. 자세한 내용은 완전정합층 문서에 있다.

주의할 점 하나. PML은 무조건 안전한 부품이 아니다. 이방성 탄성 매질이나 역방향 군속도를 갖는 매질에서는 고전 PML이 불안정해져 층 안에서 스스로 발산하는 사례가 잘 알려져 있다. 흡수층을 켜 놨는데 에너지가 늘어나면 도식이 아니라 PML을 먼저 의심하는 게 맞다.

9. 관련 문서[편집]

10. Footnotes[편집]

  1. 정확히는 맥스웰이 1862년 논문에서 이미 “빛은 같은 매질의 횡파동으로 보인다”고 적었고, 1865년 A Dynamical Theory of the Electromagnetic Field 에서 정식화했다. 당시 그가 쓴 광속 값은 피조와 베버-콜라우슈의 측정치였는데, 두 완전히 다른 실험에서 같은 숫자가 나오는 것을 보고 물리학자가 무슨 기분이었을지는 상상에 맡긴다. ↩

  2. 강한 하위헌스 원리는 홀수 차원 n≥3n\ge 3 에서만 성립한다. 2차원(예: 수면파나 얇은 판)에서 잔향이 남는 것은 매질이 지저분해서가 아니라 차원 탓이다. 짝수 차원 우주였다면 모든 소리가 잔향으로 뭉개져 언어가 성립하기 어려웠을 거라는 이야기를 교과서들이 꼭 한 번씩 하고 지나간다. ↩

  3. 립프로그는 ν≤1\nu\le 1 에서 어떤 이산 에너지 범함수를 정확히 보존한다. 다만 그 범함수는 연속 에너지와 O(Δt2)O(\Delta t^2) 만큼 다른 물건이라, 진짜 에너지는 아주 작은 폭으로 영원히 진동한다. 심플렉틱 적분기에서 에너지가 드리프트 없이 떨리기만 하는 것과 같은 이유다 — 실제로 립프로그는 파동 해밀토니안에 대한 슈퇴르머-베를레와 같은 물건이다. ↩

  4. “그냥 Δt\Delta t 를 줄이면 되지 않나”라는 유혹은 파동 문제에서 특히 비싸다. ν\nu 를 줄이면 안정성은 좋아지지만 위상 오차는 오히려 커지고(1차원에서 ν=1\nu=1 이 최적이라는 것을 기억하자), 스텝 수는 늘어난다. 확산 문제의 직관을 그대로 들고 오면 계산비만 태우고 정확도는 나빠지는, 보기 드물게 배신당하는 경우다. ↩