Hamilton-Jacobi Equations

Reference: Stanley Osher and Ronald Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Springer, 2003, Ch. 5

3장의 convection 방정식과 4장의 normal-velocity 형태 level set equation은 겉보기에 서로 다른 문제처럼 보이지만, 사실 둘 다 같은 수학적 틀 안에 들어갑니다. 이 챕터는 그 공통 틀인 Hamilton-Jacobi 방정식을 명시적으로 정의하고, 이 틀 아래에서 수치 스킴을 설계하는 일반적인 방법론(numerical Hamiltonian)을 소개합니다. 흥미롭게도 이 일반 프레임워크로 다시 살펴보면, 3장에서 다룬 단순한 upwind differencing이 사실은 이 챕터에서 소개하는 Godunov’s scheme의 특수한 경우였다는 사실이 드러납니다. 즉 이 챕터는 앞의 두 챕터를 사후적으로 통합해서 정리하는 역할을 합니다.

Introduction: Hamilton-Jacobi 방정식이란

일반적인 Hamilton-Jacobi 방정식은 다음과 같은 형태입니다.

3차원으로 풀어 쓰면 $\phi_t+H(\phi_x,\phi_y,\phi_z)=0$이 됩니다. $H$는 공간과 시간의 함수일 수도 있습니다. 지금까지 다룬 두 가지 핵심 방정식이 모두 이 형태에 들어맞습니다. 3장의 convection 방정식(외부 속도장에 의한 이류)은 $H(\nabla\phi)=\vec{V}\cdot\nabla\phi$인 경우이고, 4장에서 소개한 normal-velocity 형태의 level set equation $\phi_t+V_n\lvert\nabla\phi\rvert=0$은 $H(\nabla\phi)=V_n\lvert\nabla\phi\rvert$인 경우입니다. 여기서 $V_n$은 $\vec{x}$, $t$, 심지어 $\nabla\phi/\lvert\nabla\phi\rvert$(즉 법선 방향)에 의존할 수 있습니다.

다만 4장에서 다룬 mean curvature motion 자체(방정식 $\phi_t=b\kappa\lvert\nabla\phi\rvert$)는 Hamilton-Jacobi 방정식이 아닙니다. Hamilton-Jacobi 방정식은 $\phi$의 1차 도함수까지만 사용해야 하는데, 곡률 $\kappa$는 $\phi$의 2차 도함수로 계산되는 양이기 때문입니다. 이 구분이 중요한 이유는 앞선 챕터에서 반복해서 강조된 사실과 정확히 대응하기 때문입니다. Hamilton-Jacobi 방정식은 1차 도함수만 쓰므로 쌍곡형(hyperbolic)이고 upwind 계열 방법이 필요하며, mean curvature motion처럼 2차 도함수가 등장하는 방정식은 포물형(parabolic)이고 중심차분이 필요합니다. 즉 이 챕터가 다루는 대상은 3·4장에서 이미 구분해 두었던 “쌍곡형” 쪽 문제들 전체를 아우르는 일반화라고 이해하면 됩니다.

Conservation Law와의 연결

이 챕터의 수치 스킴들은 대부분 conservation law(보존법칙) 이론에서 빌려온 것입니다. 그래서 본격적인 수치해법에 들어가기 전에, 왜 이 둘이 그렇게 밀접하게 연결되어 있는지부터 살펴볼 필요가 있습니다.

1차원 스칼라 conservation law는 다음과 같은 형태입니다.

여기서 $u$는 보존되는 양이고 $F(u)$는 flux 함수입니다. 질량 보존을 나타내는 연속 방정식 $\rho_t+\partial_x(\rho u)=0$이 잘 알려진 예이고, 여기에 운동량·에너지 보존 방정식을 더하면 압축성 Navier-Stokes 방정식이, 점성을 무시하면 비압축성이 아닌 압축성 비점성 Euler 방정식이 나옵니다. Euler 방정식은 불연속(충격파, 접촉 불연속)을 포함할 수 있어서 고전적인 미분 가능한 해가 존재하지 않는 weak solution을 다뤄야 하고, 해가 유일하지 않을 수 있어 물리적으로 올바른 해를 골라내는 entropy condition(4장에서 다룬 vanishing viscosity solution)이 필요합니다. 이런 비선형적 성질을 가장 단순하게 압축해서 보여주는 예가 Burgers’ 방정식 $u_t+\partial_x(u^2/2)=0$이며, 매끄러운 초기 데이터에서도 불연속 충격파가 저절로 발생하고, vanishing viscosity를 강제하지 않으면 비물리적인 expansion shock이 나타난다는 점에서 훨씬 복잡한 Euler 방정식의 축소판 역할을 합니다.

1차원에서의 정확한 대응 관계

이제 1차원 Hamilton-Jacobi 방정식 $\phi_t+H(\phi_x)=0$을 생각해봅시다. 이 식 전체를 $x$에 대해 미분하면

을 얻습니다. 여기서 $u=\phi_x$로 치환하면

이 되는데, 이는 정확히 앞서 본 conservation law의 형태입니다! 즉 1차원에서는 Hamilton-Jacobi 방정식과 conservation law 사이에 정확한 대응 관계가 성립합니다. conservation law의 해 $u$는 Hamilton-Jacobi 방정식의 해 $\phi$를 미분한 것이고, 거꾸로 $\phi$는 $u$를 적분한 것입니다.

이 대응 관계는 단순한 수학적 유희가 아니라 실용적으로 매우 유용한 통찰을 줍니다. 불연속의 적분은 kink(1차 도함수가 불연속인 지점)가 된다는 사실을 이용하면, $u$가 불연속(충격파)을 가지더라도 그 적분인 $\phi$는 데이터가 매끄럽게 시작했더라도 kink를 가질 수는 있지만 일반적으로는 $\phi$ 자체가 불연속이 되지는 않는다는 것을 알 수 있습니다. $\phi$가 불연속이 되려면 대응하는 conservation law의 해가 델타 함수를 가져야 하는데, 이는 일반적인 상황이 아니기 때문입니다. 따라서 Hamilton-Jacobi 방정식 (5.2)의 해 $\phi$는 보통 연속입니다. 또한 conservation law가 해의 유일성을 보장하지 못하고 entropy condition이 필요했던 것처럼, Hamilton-Jacobi 방정식 역시 “물리적으로” 올바른 해를 고르기 위해 entropy condition이 필요합니다.

간단한 역사적 배경

Crandall과 Lions는 Hamilton-Jacobi 방정식에 대한 viscosity solution 개념과, 이를 수렴시키는 최초의 monotone 1차 정확도 수치 방법을 제안했습니다. 이후 Osher와 Sethian은 방금 살펴본 conservation law와의 대응 관계를 이용해서 더 높은 차수의 “artifact-free”(비물리적 잡음이 없는) 방법을 구축했습니다. 다만 이 1차원 대응 관계는 다차원으로 그대로 확장되지는 않습니다. 그럼에도 많은 Hamilton-Jacobi 방정식은 차원별로(dimension by dimension) 독립적으로 이산화하는 방식으로 다룰 수 있고, 이 아이디어가 Osher와 Shu의 일반적인 프레임워크로 정리되었습니다. 이 챕터의 나머지 부분은 이 Osher-Shu 프레임워크를 따라갑니다.

수치 이산화: Numerical Hamiltonian

Hamilton-Jacobi 방정식의 forward Euler 시간 이산화는 다음과 같이 씁니다.

여기서 $\hat{H}$는 $H(\phi_x,\phi_y,\phi_z)$의 수치적 근사이며 numerical Hamiltonian이라고 부릅니다. $\hat{H}$가 갖춰야 할 최소한의 요구조건은 consistency입니다. 즉 $\phi_x^-=\phi_x^+=\phi_x$(양쪽 근사가 정확한 값과 일치하는 극한)일 때 $\hat{H}(\phi_x,\phi_x;\phi_y,\phi_y;\phi_z,\phi_z)=H(\phi_x,\phi_y,\phi_z)$가 성립해야 합니다. $\phi_x^-$와 $\phi_x^+$ 자체는 1차 정확도의 one-sided differencing이나, 3장에서 다룬 고차 정확도 HJ ENO/HJ WENO로 근사할 수 있습니다. 이 챕터는 3차원으로의 확장이 직관적이라는 이유로 2차원 $H(\phi_x,\phi_y)$의 경우를 중심으로 설명을 진행합니다.

여기서 중요한 개념 하나가 monotone scheme입니다. $\phi^{n+1}$이 모든 $\phi^n$ 값들에 대해 단조증가하는(nondecreasing) 함수일 때 그 스킴을 monotone하다고 부릅니다. Crandall과 Lions는 이런 monotone scheme이 (비록 1차 정확도에 그치더라도) 올바른 해로 수렴한다는 것을 증명했습니다. 이 챕터에서 소개하는 여러 numerical Hamiltonian들이 monotone인지 아닌지가 계속 언급되는 이유가 여기에 있습니다. forward Euler 시간 이산화는 3장에서 다룬 TVD Runge-Kutta로 그대로 확장할 수 있습니다.

CFL 조건의 일반형

이 프레임워크의 CFL 조건은

이며, 여기서 $H_1,H_2,H_3$은 각각 $\phi_x,\phi_y,\phi_z$에 대한 $H$의 편미분입니다. 이 일반형이 앞선 챕터들의 특수한 CFL 조건들을 모두 포함한다는 사실이 이 프레임워크의 통합적 성격을 잘 보여줍니다. 예를 들어 3장의 convection 방정식 $H(\nabla\phi)=\vec{V}\cdot\nabla\phi$에서는 $H_1=u$, $H_2=v$, $H_3=w$이므로 이 조건은 정확히 3장에서 본 CFL 조건 (3.10)으로 환원됩니다. 4장의 level set equation $H(\nabla\phi)=V_n\lvert\nabla\phi\rvert$에서는 (단 $V_n$이 $\phi_x,\phi_y,\phi_z$에 의존하지 않는다는 가정 하에) $H_1=V_N\phi_x/\lvert\nabla\phi\rvert$, $H_2=V_N\phi_y/\lvert\nabla\phi\rvert$, $H_3=V_N\phi_z/\lvert\nabla\phi\rvert$로 조금 더 복잡해집니다.

Lax-Friedrichs 계열 스킴

$\hat{H}$을 설계하는 첫 번째 방법은 Lax-Friedrichs(LF) 스킴입니다.

이 식의 구조를 뜯어보면, 앞부분 $H(\cdots)$는 $\phi_x^-$와 $\phi_x^+$를 그냥 평균 내서(중심차분 스타일로) $H$에 대입한 항이고, 뒤에 붙은 두 항은 $\phi_x^+-\phi_x^-$, $\phi_y^+-\phi_y^-$(양방향 근사값의 차이, 즉 국소적으로 매끄럽지 않은 정도를 나타내는 양)에 비례하는 인공 점성(dissipation) 항입니다. 순수하게 평균만 쓰면(중심차분) 3장에서 이미 살펴봤듯 불안정해지므로, 이 인공 점성 항이 안정성을 위해 필수적입니다. 계수 $\alpha^x,\alpha^y$는

로, $H$의 편미분(특성 속도에 해당)의 최댓값으로 정합니다.

α를 어디까지 넓혀서 찾을 것인가

문제는 이 최댓값을 “어느 범위에서” 찾느냐입니다. 전통적인 LF 스킴은 Cartesian grid 전체에서 $\phi_x^-,\phi_x^+$의 최댓값과 최솟값을 모아 구간 $I^x$를 만들고, 그 구간 전체에서 $\lvert H_1\rvert$의 최댓값을 $\alpha^x$로 씁니다. 예를 들어 3장의 convection 방정식처럼 $H_1=u$가 $\phi_x,\phi_y$와 무관한 경우라면, $\alpha^x$는 그냥 격자 전체에서의 $\max\lvert u\rvert$가 됩니다.

여기서 실용적으로 중요한 트레이드오프가 등장합니다. $\alpha$를 크게 잡을수록 인공 점성이 늘어나 방법은 안정적이 되지만 해의 디테일이 뭉개집니다. 반대로 $\alpha$가 너무 작으면 진동이나 비물리적인 현상이 생깁니다. 문제는 도메인 전체에서 하나의 $\alpha$ 값을 쓰면, 속도가 큰 영역에서 필요한 만큼 $\alpha$를 키우는 순간 속도가 작은 다른 영역에서는 불필요하게 큰 점성이 적용되어 그 영역의 미세한 특징이 뭉개진다는 것입니다. 이 문제를 해결하기 위해 $\alpha$를 계산할 때 전역이 아니라 국소적인 이웃만 살펴보자는 아이디어가 등장합니다.

• Stencil Lax-Friedrichs (SLF): HJ WENO 근사에 실제로 쓰이는 stencil 이웃 격자점들(전형적으로 $x$ 방향은 $x_{i-3}$부터 $x_{i+3}$까지)만으로 $I^x,I^y$를 구성합니다.
• Local Lax-Friedrichs (LLF): 한 걸음 더 나아가, $\alpha^x$는 해당 격자점의 $\phi_x^-,\phi_x^+$ 값만으로 결정합니다. ($\alpha^y$는 여전히 LF 혹은 SLF 방식으로 정하며, 이 조합은 SLLF라고 구분해서 부릅니다.)
• Local Local Lax-Friedrichs (LLLF): $\alpha^x$와 $\alpha^y$ 모두 해당 격자점의 값만으로 결정합니다.

$H$가 $H(\phi_x,\phi_y)=H^x(\phi_x)+H^y(\phi_y)$처럼 분리 가능(separable)하다면, $\alpha^x$가 애초에 $\phi_y$와 무관하고 $\alpha^y$가 $\phi_x$와 무관하므로 LLLF와 LLF는 같아집니다. $H$가 분리 불가능할 때만 둘이 진짜로 다른 스킴이 됩니다. 저자들은 실무 경험상 LF와 SLF는 대체로 지나치게 소산적(dissipative)이고, LLLF는 대체로 소산이 부족해서 식 (5.11)의 중심 평균 근사가 만드는 문제를 완전히 억누르지 못하며, LLF가 실용적으로 가장 균형 잡힌 선택이라고 평가합니다. 참고로 LLF는 monotone scheme입니다.

Roe-Fix 스킴

LF 계열의 근본적인 어려움은 적절한 $\alpha$를 고르는 일이 상당히 까다롭다는 데 있습니다. 그래서 대안으로, 애초에 인공 점성을 손으로 조절할 필요가 없는 upwind 기반 방법을 쓰는 편이 낫다는 아이디어가 등장합니다. Shu와 Osher는 conservation law를 위해 Roe의 upwind 방법에 sonic point(특성 속도가 0을 지나가는 지점)에서만 LLF 보정을 더하는 방식을 제안했고, 이를 Hamilton-Jacobi 방정식에 옮긴 것이 Roe-Fix (RF) 스킴입니다.

평소에는 $\alpha^x=\alpha^y=0$으로 두어 인공 점성 항을 완전히 제거합니다. 대신 $H$의 편미분 $H_1,H_2$의 부호를 살펴서 정보가 흐르는 방향을 판단합니다. $H_1$이 (해당 구간 전체에서) 항상 양수면 정보가 왼쪽에서 오른쪽으로 흐르므로 $\phi_x^\star=\phi_x^-$를, 항상 음수면 $\phi_x^\star=\phi_x^+$를 씁니다. $H_2$에 대해서도 같은 논리가 적용됩니다. $H_1,H_2$ 모두 부호가 바뀌지 않는다면 완전히 upwind로만 처리하면 되고, 이는 3장에서 다룬 순수 upwind differencing과 본질적으로 같습니다.

문제는 $H_1$이나 $H_2$의 부호가 구간 안에서 바뀌는 경우, 즉 sonic point 근방입니다. 이는 특성 속도가 0이 되는 지점으로, 해가 유일하지 않을 수 있는 위험 구간이라 물리적으로 올바른 vanishing viscosity 해를 골라내기 위한 인공 점성이 다시 필요해집니다. 이때만 국소적으로 LLF 방식으로 전환해서 그 방향에만 감쇠를 추가합니다. 한 방향에서만 sonic point가 발견되면 그 방향에만 damping을 추가하고 다른 방향은 순수 upwind를 유지하는 식으로, 정말 필요한 곳에만 인공 점성을 쓰는 것이 이 스킴의 핵심 아이디어입니다.

이 방식은 계산 비용 면에서도 실용적인 이점이 있습니다. $\phi_x^-$와 $\phi_x^+$를 HJ WENO 같은 고차 스킴으로 계산하는 것은 비용이 큰데, RF 스킴에서는 upwind로 처리되는 방향에서는 둘 중 하나만 있으면 됩니다. 그래서 먼저 값싼 1차 정확도의 one-sided difference로 sonic point 여부를 먼저 판단하고, 실제로 필요한 값만 HJ WENO로 정밀하게 계산하는 전략을 씁니다. sonic point는 실무에서 드물게 발생하므로, 이 전략은 비용이 큰 HJ WENO 호출 횟수를 대략 절반으로 줄여줍니다.

Godunov’s 스킴

Godunov는 원래 조각별 상수(piecewise constant) 초기 데이터를 갖는 1차원 conservation law의 Riemann 문제를 정확히 풀어내는 방법을 제안했습니다. 이를 다차원 Hamilton-Jacobi 방정식으로 옮긴 형태는

로 쓸 수 있고, 이것이 가장 표준적인(canonical) monotone scheme입니다. 여기서 $\text{ext}^x H$는 $\phi_x^-<\phi_x^+$이면 $\phi_x\in I^x$ 구간에서 $H$의 최솟값을, $\phi_x^->\phi_x^+$이면 최댓값을, $\phi_x^-=\phi_x^+$이면 그 값을 그대로 $H$에 대입한 값을 뜻합니다. $\text{ext}^y H$도 마찬가지로 정의됩니다. 일반적으로 $\text{ext}^x\,\text{ext}^y H \neq \text{ext}^y\,\text{ext}^x H$이므로 연산 순서에 따라 서로 다른 버전의 Godunov’s method가 나올 수 있지만, $H$가 분리 가능한 경우를 포함한 많은 실용적인 경우에는 두 순서가 같은 결과를 줍니다.

Godunov’s method가 곧 upwind differencing이었다

이 챕터가 주는 가장 통합적인 통찰은 여기서 나옵니다. 3장의 convection 방정식 $H(\phi_x,\phi_y)=u\phi_x+v\phi_y$는 $x$와 $y$ 방향에 대해 분리 가능하므로 $\text{ext}^x\,\text{ext}^y H = \text{ext}^x(u\phi_x)+\text{ext}^y(v\phi_y)$로 각 방향을 독립적으로 처리할 수 있습니다. $\phi_x^-<\phi_x^+$인 경우 $u\phi_x$의 최솟값을 원한다면, $u>0$일 때는 $\phi_x^-$를, $u<0$일 때는 $\phi_x^+$를 쓰면 됩니다. $\phi_x^->\phi_x^+$인 경우(최댓값을 원하는 경우)도 똑같이 $u>0$이면 $\phi_x^-$, $u<0$이면 $\phi_x^+$가 답이 됩니다. 즉 어느 경우든 결론은 “$u>0$이면 $\phi_x^-$를 쓰고, $u<0$이면 $\phi_x^+$를 쓴다”로 요약되는데, 이는 3장에서 처음부터 직접 도입했던 표준 upwind differencing과 정확히 동일합니다. 다시 말해, convection 방정식에 대해서는 Godunov’s method가 곧 단순 upwind differencing입니다.

Comment

이 챕터를 다 읽고 나면, 3장과 4장에서 각각 독립적으로 등장했던 여러 스킴들이 사실 하나의 질문에 대한 서로 다른 답이었다는 것이 분명해집니다. “numerical Hamiltonian $\hat{H}$을, consistency를 지키면서도 안정적이 되도록 어떻게 근사할 것인가?” 이 질문에 대한 답이 upwind differencing(Godunov’s method의 특수한 경우)이었고, LF/SLF/LLF/LLLF는 그 답을 “인공 점성을 얼마나, 어디서 계산해서 넣을 것인가”의 관점에서 체계화한 것이며, Roe-Fix는 “웬만하면 upwind로 공짜로 처리하고, 정말 위험한 sonic point에서만 인공 점성을 쓰자”는 절충안이었습니다. 이런 관점에서 보면 이 챕터는 새로운 내용을 추가한다기보다, 앞선 두 챕터의 도구들이 왜 그런 형태를 취했는지에 대한 이론적 정당화를 사후적으로 제공하는 챕터에 가깝습니다.

1차원에서 Hamilton-Jacobi 방정식과 conservation law가 $u=\phi_x$라는 치환 하나로 정확히 대응한다는 사실도 개인적으로 상당히 우아하다고 생각합니다. 이 대응 관계 덕분에 “왜 level set 방법에서 $\phi$가 매끄러운 데이터에서 시작해도 kink를 가질 수 있는가”라는, 직접 증명하려면 까다로울 수 있는 질문에 “적분된 불연속은 kink가 된다”는 아주 직관적인 답을 즉시 얻을 수 있습니다. 다만 이 대응 관계가 다차원에서는 깨진다는 점, 그리고 그럼에도 불구하고 실용적으로는 차원별(dimension-by-dimension) 처리로 우회할 수 있다는 점은 이 챕터가 솔직하게 인정하는 한계입니다.

아쉬운 점이 있다면, 이 챕터는 각 스킴(LF, RF, Godunov)을 소개하고 실무적인 우열관계(LLF가 대체로 최선)를 언급하는 데 그치고, 왜 LLLF가 부족한 소산을 주는지, 왜 LLF가 정확히 그 중간 지점인지에 대한 정량적인(오차 추정이나 안정성 증명 같은) 근거는 제시하지 않는다는 것입니다. 저자들의 표현대로 이는 “실무 경험(practice)”에 근거한 권고이지 정리(theorem)가 아닙니다. 지금 단계에서 기억해야 할 핵심은 다음과 같습니다: (1) Hamilton-Jacobi 방정식 $\phi_t+H(\nabla\phi)=0$은 convection과 normal-velocity 형태의 level set equation을 포괄하지만 곡률처럼 2차 도함수에 의존하는 방정식은 배제한다, (2) 1차원에서는 $u=\phi_x$ 치환으로 conservation law와 정확히 대응하며 이로부터 $\phi$가 kink는 가질 수 있어도 불연속은 잘 갖지 않는다는 사실이 따라 나온다, (3) numerical Hamiltonian은 consistency를 만족해야 하고, monotone하면 (1차 정확도로) 수렴이 보장된다, (4) LF 계열은 인공 점성을 어디서(전역/stencil/국소) 계산하느냐로 갈리며 실무적으로 LLF가 가장 균형 잡혀 있다, (5) 단순 convection 방정식에서는 Godunov’s method가 곧 3장의 upwind differencing과 동일하다.