Motion Involving Mean Curvature

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

3장에서는 인터페이스가 외부에서 주어진 속도장 $\vec{V}(\vec{x},t)$를 따라 움직이는 상황(convection)을 다뤘습니다. 이 챕터는 완전히 다른 종류의 운동, 즉 속도장 자체가 $\phi$로부터 계산되는 자기생성(self-generated) 속도장을 다룹니다. 그 대표적인 예가 곡률에 비례해서 움직이는 mean curvature motion입니다. 이 챕터가 중요한 이유는 단순히 새로운 응용 예제를 하나 추가하는 데 있지 않습니다. mean curvature motion을 지배하는 방정식은 3장의 level set equation과 수학적 성질(쌍곡형 vs 포물형)이 근본적으로 다르고, 그래서 3장에서 애써 구축한 upwind 기반 수치 스킴을 그대로 쓸 수 없습니다. 즉 이 챕터는 “같은 $\phi$를 진화시키는 문제라도 방정식의 성질에 따라 완전히 다른 수치해법이 필요하다”는 것을 보여주는 좋은 대조 사례입니다.

Introduction: 곡률을 따라 움직이는 인터페이스

mean curvature motion은 인터페이스가 법선 방향으로, 그 지점의 곡률 $\kappa$에 비례하는 속도로 움직이는 운동입니다.

여기서 $b>0$은 상수입니다. $b>0$일 때 인터페이스는 오목한(concave) 방향, 즉 곡률을 줄이는 방향으로 움직입니다. 그 결과 2차원에서 원은 계속 작아지다가 한 점으로 수축해서 사라집니다. 반대로 $b<0$이면 인터페이스가 볼록한 방향으로 움직여서 원이 점점 커지는데, 책은 이 경우를 아예 다루지 않습니다. 이유가 흥미로운데, $b<0$인 상황에서는 인터페이스 위의 아주 작은 섭동(예를 들어 반올림 오차로 생긴 미세한 요철)조차 점점 자라나서 결국 $O(1)$ 크기의, 즉 무시할 수 없는 크기의 왜곡으로 성장해버립니다. 작은 오차가 시간이 지날수록 억제되기는커녕 오히려 증폭되는 이런 상황을 ill-posed하다고 부르며, 이는 수치적으로 다룰 수 없는 문제입니다. 반대로 $b>0$인 경우는 작은 섭동이 곡률에 의해 스스로 매끄러워지면서 사라지는 방향이므로, 안정적으로 계산할 수 있습니다.

책이 예시로 든 그림 두 개가 이 운동의 성격을 잘 보여줍니다. 감겨 있는 나선(spiral) 모양의 인터페이스를 이 흐름에 놓으면, 곡률이 큰 나선의 양 끝부분이 곡률이 작은 몸통 부분보다 훨씬 빠르게 움직여서, 시간이 지날수록 나선이 풀리며 점점 짧아지다가 결국 사라집니다. 별 모양의 인터페이스에 적용한 다른 예시에서는, 곡률이 큰 별의 뾰족한 끝 부분은 안쪽으로 빠르게 들어오는 반면 곡률이 작은(오목한) 골 부분은 바깥쪽으로 밀려 나가서, 결과적으로 뾰족한 부분과 오목한 부분이 서로 상쇄되며 인터페이스가 점점 원에 가까운 매끄러운 모양으로 수렴합니다. 이 두 예시 모두 mean curvature motion이 본질적으로 인터페이스를 매끄럽게 만드는(smoothing) 확산적 운동이라는 사실을 보여줍니다.

운동 방정식 유도

이 속도장은 법선 방향 성분만 가지고 접선(tangential) 방향 성분이 정확히 0이라는 점이 중요합니다. 사실 이는 mean curvature motion만의 특별한 성질이 아니라, level set 프레임워크 전반에 적용되는 일반적인 사실입니다. 접선 벡터 $\vec{T}$는 정의상 $\vec{T}\cdot\nabla\phi=0$을 만족하므로, 속도장의 접선 성분은 level set equation에 대입했을 때 자동으로 사라집니다. 2차원에서 $\vec{V}=V_n\vec{N}+V_t\vec{T}$로 속도를 법선 성분과 접선 성분으로 나누어 쓰면, level set equation

은 $\vec{T}\cdot\nabla\phi=0$이므로

로 줄어듭니다. 여기서 핵심은 $\vec{N}\cdot\nabla\phi$가 사실 아주 단순한 형태로 정리된다는 점입니다.

이므로, 최종적으로

를 얻습니다. 이것이 normal velocity 형태의 level set equation입니다. 3장의 $\phi_t+\vec{V}\cdot\nabla\phi=0$과 비교하면, 벡터 내적 $\vec{V}\cdot\nabla\phi$가 스칼라 곱 $V_n\lvert\nabla\phi\rvert$로 바뀌었을 뿐인데, 이 작은 표기상의 차이가 실은 이후 수치해석 전략을 완전히 갈라놓는 지점입니다.

여기에 mean curvature motion의 정의인 $V_n=-b\kappa$를 대입하면

를 얻습니다. 책은 이 우변 $b\kappa\lvert\nabla\phi\rvert$가 포물형(parabolic) 항이라서 3장에서 다룬 upwind 방식으로는 이산화할 수 없다고 명확히 짚습니다. 특히 $\phi$가 signed distance function이라면($\lvert\nabla\phi\rvert=1$, $\kappa=\Delta\phi$), 이 식은 정확히 열방정식(heat equation)

이 됩니다. 여기서 $\phi$는 온도, $b$는 열전도율에 대응하는 역할을 합니다. 열방정식은 포물형 방정식의 가장 기본적인 형태이므로, mean curvature motion은 본질적으로 “인터페이스라는 도형에 열방정식을 적용한 것”이라고 이해할 수 있습니다. 다만 이 두 식이 서로 바꿔 쓸 수 있는 것은 $\phi$가 signed distance function일 때뿐입니다. forward Euler로 한 스텝을 전진시키고 나면 새로운 $\phi$ 값은 더 이상 signed distance function이 아니게 되므로, 그다음 스텝에서 두 식은 더 이상 같은 값을 주지 않습니다. 이 문제를 해결하려면 매 스텝(혹은 주기적으로) $\phi$를 signed distance function으로 되돌리는 재초기화(reinitialization)가 필요한데, 이는 7장에서 다룬다고 예고합니다. 중요한 점은, 이 재초기화가 인터페이스의 실제 위치 자체를 바꾸는 것이 아니라 그 위치를 표현하는 embedding 함수 $\phi$의 형태만 다시 다듬는 작업이라는 것입니다.

수치 이산화

열방정식 같은 포물형 방정식은 어느 한 방향에서만 정보가 흘러오는 쌍곡형 방정식과 달리, 모든 공간 방향의 정보에 동시에 의존합니다(domain of dependence가 전 방향입니다). 따라서 특정 방향으로 치우친 upwind 근사는 애초에 맞지 않고, 중심차분(central differencing)을 써야 합니다. 열방정식의 $\Delta\phi$ 항은 1장의 2차 정확도 공식으로, mean curvature motion의 곡률 $\kappa$와 $\nabla\phi$ 항도 마찬가지로 1장에서 다룬 2차 정확도 중심차분으로 근사합니다. 이 공간 이산화는 2차 정확도에 그치지만, 책은 포물형 방정식 특유의 소산적(dissipative) 성질 때문에 이 정도 정확도로도 실용적으로 충분하다고 설명합니다.

forward Euler의 가혹한 시간 스텝 제약

문제는 시간 이산화입니다. $\Delta\phi$를 중심차분으로 근사하고 forward Euler로 시간 전진시키면, 안정성을 위해

이라는 조건이 필요합니다. 여기서 $\Delta t$가 $O((\Delta x)^2)$ 스케일로 제약된다는 점이 3장의 쌍곡형 방정식(CFL 조건에서 $\Delta t\sim O(\Delta x)$)과 결정적으로 다릅니다. 격자를 두 배 촘촘하게 하면 쌍곡형 문제에서는 필요한 스텝 수가 두 배로 늘지만, 포물형 문제에서는 네 배로 늘어난다는 뜻입니다. 흥미롭게도 이 $O((\Delta x)^2)$ 시간 스텝 제약을 그대로 받아들이면, 공간 2차·시간 1차인 forward Euler 조합 전체가 결과적으로 $O((\Delta x)^2)$ 정확도를 갖게 됩니다. 시간 스텝 자체가 이미 $(\Delta x)^2$ 스케일이므로, 1차 정확도인 시간 오차 $O(\Delta t)$가 곧 $O((\Delta x)^2)$이 되기 때문입니다.

암시적(implicit) 방법으로 제약을 완화하기

이 가혹한 스텝 제약을 피하는 표준적인 방법은 안정성 영역이 더 넓은 시간 적분기, 즉 implicit 방법을 쓰는 것입니다. 1차 정확도의 backward Euler를 열방정식에 적용하면

이 되는데, 이 식은 $\Delta t$ 크기에 대한 안정성 제약이 전혀 없습니다. 즉 $\Delta t$를 순수하게 정확도만 고려해서 고를 수 있고, 보통 $\Delta t=O(\Delta x)$ 정도로 선택합니다. 다만 이렇게 하면 시간 방향 정확도가 1차에 그쳐서 전체 정확도가 $O(\Delta x)$로 떨어집니다. 이를 개선한 것이 사다리꼴 공식(trapezoidal rule)입니다.

이 식은 시간에 대해 $O((\Delta t)^2)$ 정확도를 가지므로, $\Delta t=O(\Delta x)$로 두더라도 전체 정확도가 $O((\Delta x)^2)$로 유지됩니다. 이렇게 중심차분(공간)과 사다리꼴 공식(시간)을 결합한 방식을 흔히 Crank-Nicolson 스킴이라고 부릅니다.

물론 공짜는 없습니다. backward Euler나 Crank-Nicolson 모두 매 시간 스텝마다 $\phi^{n+1}$을 구하기 위한 선형 연립방정식을 풀어야 합니다. 다행히 열방정식의 경우 $\Delta\phi^{n+1}$이 $\phi^{n+1}$에 대해 선형이므로 이 선형계를 푸는 것 자체는 어렵지 않습니다. 하지만 (4.5)식, 즉 $\kappa\lvert\nabla\phi\rvert$ 형태를 implicit으로 이산화하려 하면 상황이 훨씬 복잡해집니다. $\kappa^{n+1}\lvert\nabla\phi^{n+1}\rvert$는 $\phi^{n+1}$에 대해 비선형이기 때문입니다.

implicit 방법을 쓸 때는 두 형태가 더 이상 호환되지 않는다

여기서 책이 강하게 경고하는 대목이 있습니다. explicit(forward Euler) 방법에서는 $\phi$가 signed distance function을 유지하는 한 $\Delta\phi$와 $\kappa\lvert\nabla\phi\rvert$를 서로 바꿔 써도 괜찮았습니다. 그런데 implicit 시간 이산화에서는 이 호환성이 깨집니다. $\phi^n$이 정확히 signed distance function이었다고 해도, 선형계를 풀어서 얻은 $\phi^{n+1}$은 일반적으로 signed distance function이 아닙니다. 그 결과 $\Delta\phi^{n+1}$은 더 이상 $\kappa^{n+1}\lvert\nabla\phi^{n+1}\rvert$의 좋은 근사가 아니게 됩니다. 즉 $\Delta\phi^n$이 $\kappa^n\lvert\nabla\phi^n\rvert$과 정확히 같았다는 사실이, 다음 스텝에서도 $\Delta\phi^{n+1}$이 $\kappa^{n+1}\lvert\nabla\phi^{n+1}\rvert$과 같으리라는 것을 보장하지 않습니다. 책은 이것이 “$\phi$가 signed distance function일 때 얻어지는 여러 단순화(법선을 $\nabla\phi$로, 곡률을 $\Delta\phi$로 대체하는 것 등)가 매우 유용하지만, $\phi$가 signed distance function이 아닌 경우에는 그 단순화 사이에 중요한 차이가 생긴다”는 책 전체를 관통하는 주의사항의 구체적인 사례라고 볼 수 있습니다.

Convection-Diffusion 방정식

실제 응용에서는 외부 속도장에 의한 이류(convection)와 곡률에 의한 확산(diffusion)이 함께 작용하는 경우가 흔합니다. 표준적인 convection-diffusion 방정식

의 level set 버전은

이며, 이 둘도 $\phi$가 signed distance function을 유지하는 한 서로 바꿔 쓸 수 있습니다. 수치적으로는 지금까지 다룬 두 가지 기법을 그대로 결합합니다. 이류를 나타내는 $\vec{V}\cdot\nabla\phi$ 항은 3장의 upwind 방법으로, 확산을 나타내는 $b\Delta\phi$(또는 $b\kappa\lvert\nabla\phi\rvert$) 항은 이 챕터의 중심차분으로 이산화합니다. 시간 이산화로 TVD Runge-Kutta를 쓴다면, 안정성을 위한 시간 스텝 제약은 쌍곡형 조건과 포물형 조건을 단순히 더한 형태가 됩니다.

이는 자연스러운 결과입니다. 이류 항이 요구하는 스텝 제약과 확산 항이 요구하는 스텝 제약 중 어느 하나라도 어기면 전체 방법이 불안정해지므로, 두 제약을 동시에 만족시켜야 하고 그 결과 조건이 합쳐지는 형태로 나타납니다.

인공 점성(artificial viscosity)과의 연결

이 챕터는 흥미로운 이론적 연결 고리 하나를 짚으며 마무리됩니다. 만약 $O(1)$ 크기의 $b$ 대신, 격자를 세밀화할수록 0으로 사라지는 $O(\Delta x)$ 크기의 작은 계수 $\epsilon$을 쓴다면

라는 식을 얻고, 이는 $\epsilon\to 0$일수록 3장의 순수 convection 방정식 $\phi_t+\vec{V}\cdot\nabla\phi=0$에 점점 가까워집니다. 우변에 이렇게 작은 확산 항을 인위적으로 추가하는 것을 artificial viscosity(인공 점성) 방법이라고 부릅니다. 이는 전산유체역학(CFD)에서 유래한 아이디어로, 불연속(충격파 등)이 있어 고전적인(classical) 해가 존재하지 않는 상황에서도, $\epsilon\to 0$ 극한으로 물리적으로 올바른 약해(weak solution)를 골라내는 vanishing viscosity solution 개념과 연결됩니다.

여기서 저자들이 짚는 통찰이 특히 흥미롭습니다. 3장에서 다룬 upwind 이산화 자체가, 사실 이 $\epsilon\Delta\phi$ 항과 똑같은 역할을 하는 절단 오차(truncation error)를 내재적으로 갖고 있다는 것입니다. 1차 정확도 upwind는 $O(\Delta x)$ 크기의 내재적 인공 점성을, $r$차 정확도의 고차 upwind 방법(HJ ENO/WENO 등)은 $O((\Delta x)^r)$ 크기의 인공 점성을 스스로 만들어냅니다. 즉 3장에서 “정보가 오는 방향을 존중한다”는 직관으로 도입했던 upwind differencing이, 사실은 CFD의 인공 점성 이론과 같은 수학적 근거(적절한 크기의 소산을 주입해서 올바른 약해를 선택한다)를 공유하고 있다는 것입니다.

다만 인터페이스 문제에는 CFD와 다른 특수성이 있습니다. Sethian은 곡선이 코너(꺾인 지점)로 흘러 들어가야 한다는 entropy condition을 제안하고, 이것이 스스로 교차하는(self-intersecting) 곡선에 대해서도 올바른 약해를 만들어낸다는 수치적 증거를 제시했습니다. 이 entropy condition이 시사하는 바는, 인터페이스의 진화에서는 균일한 $\epsilon\Delta\phi$보다 곡률로 가중된 $\epsilon\kappa\lvert\nabla\phi\rvert$가 더 적절한 형태의 vanishing viscosity라는 것입니다. Osher와 Sethian은 이를 엄밀하게 정리해서,

가 (4.13)의 $\epsilon\Delta\phi$ 버전보다 level set 방법에 더 자연스러운 선택이라고 밝혔습니다. 물론 이 두 형태도 $\phi$가 signed distance function일 때만 서로 같습니다.

Comment

이 챕터를 3장과 나란히 놓고 보면 하나의 대조가 뚜렷하게 드러납니다. 3장의 방정식(순수 convection)은 쌍곡형이라 upwind가 필요했고, 이 챕터의 방정식(mean curvature motion)은 포물형이라 중심차분이 필요합니다. 같은 $\phi$, 같은 level set 프레임워크를 쓰고 있지만 방정식의 수학적 분류가 다르면 요구되는 수치기법이 정반대가 될 수 있다는 것을, 아주 명확한 대비를 통해 보여주는 챕터라고 생각합니다. 그리고 convection-diffusion처럼 두 성질이 섞인 방정식에서는 두 기법을 항 별로 나누어 적용하고, 안정성 조건도 단순히 더해서 만족시키면 된다는 사실은 실용적으로 매우 유용한 지침입니다.

개인적으로 가장 인상 깊었던 부분은 “implicit 방법을 쓸 때는 $\Delta\phi$와 $\kappa\lvert\nabla\phi\rvert$가 더 이상 호환되지 않는다”는 경고입니다. 앞선 챕터들에서 signed distance function이 주는 단순화(법선 $\to\nabla\phi$, 곡률 $\to\Delta\phi$)는 거의 공짜 혜택처럼 소개되었는데, 이 챕터는 그 혜택에 숨은 전제 조건(φ가 계속 signed distance function이어야 한다는 것)이 명시적으로 깨지는 첫 번째 구체적인 상황을 보여줍니다. explicit 방법에서는 이 전제가 매 스텝 살짝만 깨지고 재초기화로 회복 가능하지만, implicit 방법은 선형계를 푸는 과정 자체가 $\phi$를 signed distance function에서 상당히 멀어지게 만들 수 있어서 문제가 더 근본적입니다. 이는 왜 실무에서 mean curvature 항을 implicit으로 처리하려는 시도가 조심스러운지를 이해하는 데 중요한 단서입니다.

다만 이 챕터가 다루지 않는 부분도 있습니다. implicit 방법에서 비선형 $\kappa^{n+1}\lvert\nabla\phi^{n+1}\rvert$ 항을 실제로 어떻게 푸는지(뉴턴 반복법 등 구체적인 알고리즘)는 설명되지 않고, $\phi$를 signed distance function으로 되돌리는 재초기화 절차도 7장으로 미뤄져 있습니다. 즉 이 챕터는 “포물형 방정식이라는 다른 종류의 문제가 있고, 이를 위해서는 중심차분과 (선택적으로) implicit 시간 적분이 필요하다”는 개념적 틀과 그 함정(SDF 호환성 붕괴)을 제시하는 데 집중하고, 실제 구현 디테일은 다른 장으로 넘긴 것으로 보입니다. 지금 단계에서 기억해야 할 핵심은 다음과 같습니다: (1) mean curvature motion은 $V_n=-b\kappa$인 자기생성 속도장이며 본질적으로 인터페이스를 매끄럽게 만드는 확산적 운동이다, (2) 이 운동을 지배하는 방정식 $\phi_t=b\kappa\lvert\nabla\phi\rvert$는 포물형이라 upwind가 아니라 중심차분이 필요하다, (3) forward Euler의 시간 스텝 제약은 $O((\Delta x)^2)$로 쌍곡형($O(\Delta x)$)보다 훨씬 가혹하며, 이를 완화하려면 backward Euler나 Crank-Nicolson 같은 implicit 방법을 쓰되 비선형 항과 SDF 호환성 붕괴라는 대가를 치러야 한다, (4) 이류와 확산이 섞인 문제는 upwind와 중심차분을 항별로 나누어 적용하고 두 안정성 조건을 더해서 만족시킨다, (5) upwind differencing의 내재적 절단 오차는 CFD의 인공 점성과 같은 역할을 하며, 인터페이스 문제에서는 균일한 $\epsilon\Delta\phi$보다 곡률로 가중된 $\epsilon\kappa\lvert\nabla\phi\rvert$가 더 자연스러운 vanishing viscosity 형태다.