Closest Point Method for Surface PDEs

• Closest Point Method is an embedding-based numerical framework that discretizes PDEs on smooth, open or closed manifolds by extending surface functions via the closest-point mapping.
• It converts intrinsic surface operators into standard Cartesian derivatives, enabling the use of finite-difference methods and interpolation on a narrow band around the manifold.
• CPM supports both explicit and implicit time-stepping schemes with rigorous convergence and stability, efficiently addressing a variety of elliptic, parabolic, and variational problems.

The Closest Point Method (CPM) is an embedding-based numerical framework for the discretization and solution of partial differential equations (PDEs) defined on smooth, possibly open or closed, manifolds immersed in Euclidean space. CPM transfers the complexity of intrinsic surface operators to the embedding space by exploiting the closest-point extension, allowing surface derivatives to be replaced by standard Cartesian derivatives operating on closest-point-extended data. This enables the application of classical finite-difference methods, interpolation, and time-stepping schemes on narrow uniform Cartesian bands surrounding the surface, thus facilitating high-order, efficient, and robust solvers for a broad class of PDEs including elliptic, parabolic, reaction–diffusion, and variational problems. CPM supports both explicit and implicit algorithms, is straightforward to implement, parallelizes efficiently, and accommodates nonlinear, variable-coefficient, and higher-order surface operators while maintaining rigorous convergence and stability guarantees (Glehn et al., 2013).

1. Mathematical Foundations and Closest-Point Extension

Let $\Gamma \subset \mathbb{R}^n$ be a smooth manifold (e.g., surface or curve). The method is constructed around the closest-point mapping

$$\operatorname{cp}: B(\Gamma) \to \Gamma, \quad \operatorname{cp}(x) = \arg\min_{y \in \Gamma} \|x - y\|$$

where $B(\Gamma)$ is a narrow tubular neighborhood ("band") around $\Gamma$. For any function $u$ defined on $\Gamma$, the closest-point extension operator $E$ is defined as

$$(Eu)(x) := u(\operatorname{cp}(x)),\quad x \in B(\Gamma).$$

By construction, $(E^2 = E)$ and $Eu$ is constant in the normal direction to $\Gamma$.

Critical analytic identities underlie the CPM:

• Gradient/Divergence/Laplacian principles: For $v = Eu$, and $x$ on the surface $y = \operatorname{cp}(x)$,
  • $\nabla (Eu)(y) = \nabla_\Gamma u(y)$,
  • $$\operatorname{div}(Eg)(y) = \operatorname{div}_\Gamma g(y)$$,
  • For Euclidean cp: $\Delta (Eu)(y) = \Delta_\Gamma u(y)$.

These equivalence properties justify replacing surface differential operators with Cartesian derivatives of extended functions for all intrinsic surface operators $L_\Gamma$, including nonlinear or higher-order cases (Glehn et al., 2013).

2. Embedding PDEs and Penalized Modified Equations

Given a time-dependent surface PDE $u_t(y, t) = L_\Gamma u(y,t)$ on $\Gamma$, CPM constructs an embedded formulation in the band:

• Extend to $v(x, t)= Eu(x, t)$. On $\Gamma$, $v_t = L_\Gamma u = [L (Eu)]|_\Gamma$. Extension via $E$ to $B(\Gamma)$ gives:

$v_t = E L v, \qquad v = E v.$

• Rather than enforce $v=Ev$ as a constraint, CPM introduces a penalized single PDE:

$$v_t = E L v - \sigma (v - Ev), \quad x \in B(\Gamma)$$

for penalty parameter $\sigma > 0$. The side-by-side or constraint formulation and this penalized system are equivalent as $v - Ev \to 0$ exponentially with time (Glehn et al., 2013).

3. Discretization: Method of Lines and Semi-Discrete System

To discretize, one introduces a uniform Cartesian grid of spacing $h$ covering the active band:

• Let $v \in \mathbb{R}^N$ denote the vector of unknowns $v_i \approx v(x_i)$ at grid nodes $x_i$.
• Cartesian finite-difference matrices approximate operators (e.g., $L$ for Laplacian with 5-point or 7-point stencil).
• Discrete extension/interpolation operator $E$ at each grid node implements $(Ev)(x_i) = \sum_j w_{ij} v_j$, where $w_{ij}$ are polynomial interpolation weights for $cp(x_i)$ over a local Cartesian stencil.

The semi-discrete method-of-lines (MOL) system is:

$v'(t) = E L v - \sigma (I - E) v.$

Time-stepping can involve explicit Euler:

$$v^{n+1} = (I - \Delta t\,\sigma(I - E)) v^n + \Delta t\, E L v^n$$

or implicit schemes (e.g., implicit Euler, IMEX, BDF2, Runge–Kutta), with appropriate treatment of stiff/nonlinear/penalty terms.

The solution on $\Gamma$ at any given time is obtained by restricting $v^n$ back to the surface grid, either at grid points $x_{i(j)}$ with $\operatorname{cp}(x_{i(j)}) = y_j$ or via further interpolation (Glehn et al., 2013).

4. Closest-Point Construction, Interpolation, and Implementation

For surfaces given as level sets ($\phi(x)=0$), the closest point map is:

$$\operatorname{cp}(x) = x - \phi(x)\, \frac{\nabla\phi(x)}{|\nabla\phi(x)|}.$$

For parameterized or triangulated $\Gamma$, local search or iterative projection (e.g., Newton's method) is applied to find $cp(x)$ onto the mesh (Glehn et al., 2013).

The interpolation stencil for $E$ typically uses polynomial interpolation of degree $p$ (commonly $p=3$ for second-order methods). The band width must cover all finite-difference and interpolation stencils.

• $E$ is a sparse matrix (with per-row sparsity determined by stencil size).
• No explicit "reinitialization" is required: the penalty term $\sigma (I-E)$ enforces normal constancy.

5. Consistency, Stability, and Convergence

5.1 Consistency

• Finite-difference truncation in $L$ is $O(h^2)$ for the 5-point Laplacian.
• Interpolation error in $E$ is $O(h^{p+1})$.
• Penalty term error is $O(h^{p+1-\alpha})$ if $\sigma = O(h^{-\alpha})$.
• To retain second order, one typically sets $p \geq 2$ and $\alpha \leq 2$, e.g., $\sigma = 2d/h^2$ for embedding dimension $d$.

5.2 Stability

• The penalty $-\sigma (v-Ev)$ drives $v-Ev$ to zero exponentially: $v_t = -\sigma (v-Ev)$, so $w_t = -\sigma w$, ensuring decay for $\sigma > 0$.
• Stability restrictions for Forward Euler: $$|1-\Delta t\,\sigma| \leq 1 \implies \Delta t \leq 2/\sigma$$.
• Diffusion (from $E L$) imposes $\Delta t \leq h^2/(2d)$; with $\sigma = 2d/h^2$, both restrictions coincide.

5.3 Convergence

• For Laplace–Beltrami and canonical problems (on the circle, sphere, etc.), second-order convergence in the maximum norm is observed with $p=3, \sigma=2d/h^2, \Delta t=h^2/(4d)$, e.g., for the heat equation.
• For higher-order PDEs (biharmonic), implicit time-stepping and higher-order interpolation ($p \geq 4$) are employed, still yielding second-order convergence (Glehn et al., 2013).

6. Model Problems and Representative Applications

Canonical equations and their embedded CPM forms include:

• Heat equation: $u_t = \Delta_\Gamma u$ becomes $v_t = E \Delta v - \sigma (v-Ev)$.
• Surface biharmonic: $$u_t = -\Delta_\Gamma^2 u \to v_t = -E E \Delta \Delta v - \sigma (v - E v)$$.
• Reaction–diffusion systems (e.g., Gray–Scott): Coupled CPM-embedded PDEs for each field, with IMEX time-stepping, directly applied to triangulated and genus-3 surfaces.
• Curvature-dependent diffusion: $$u_t = \operatorname{div}_\Gamma(a(\kappa)\nabla_\Gamma u)$$, with $\kappa$ (mean curvature) computed as the two-norm of the Cartesian Laplacian applied to each cp-component; standard variable-coefficient stencils used.

The approach uniformly applies to both closed and open geometries, as well as to non-linear or variable-coefficient surface differential operators (Glehn et al., 2013).

7. Comparison With Alternative CPM Formulations

CPM as a method-of-lines unifies several earlier strands:

• Ruuth–Merriman (2008): Explicit two-step alternating evolution and extension, restricted to explicit time-stepping and limited higher-order capabilities.
• Merriman–Ruuth implicit CPM (2009): Introduces diagonal stabilization and Laplacian-specific penalty, less general.
• RBF-based and mesh-free variants (Piret, Flyer–Wright): Avoid regular grids in favor of scattered-data interpolation to increase order, at the cost of dense matrix operations.

The method-of-lines CPM:

• Permits explicit, implicit, and high-order time-integration schemes in a unified PDE framework.
• Supports nonlinear, higher-order, and variable-coefficient surface operators.
• Only requires standard sparse finite-difference stencils and sparse interpolation.
• Achieves $O(N)$ per-time-step complexity for $N$ band nodes.
• Admits rigorous stability and convergence analysis paralleling classical MOL theory (Glehn et al., 2013).

Summary table comparing CPM variants:

Variant	Time-Stepping	Surface Operators	Implementation Features
Ruuth–Merriman (2008)	Explicit only	Laplace-Beltrami	Two-step extension
MOL-CPM (this method)	Explicit/implicit/HO	General ($L_\Gamma$)	Penalized single PDE, O(N)
RBF (Piret et al.)	Varied	High-order, all $L_\Gamma$	Mesh-free, dense ops

The CPM in the method-of-lines framework provides a robust, parallelizable, and theoretically sound strategy for solving a broad spectrum of PDEs on arbitrary smooth surfaces using only local Cartesian operations and standard ODE solvers (Glehn et al., 2013).