An Optimization Journey Toward Minimal Surfaces
Mean Curvature Flow (MCF) is one of the most beautiful geometric flows in differential geometry. It deforms a surface in the direction of its mean curvature vector, continuously decreasing area until the surface converges to a minimal surface — or develops singularities along the way. In this post, we’ll explore the mathematics, implement a concrete numerical example from scratch in Python, and visualize the entire evolution in 3D.
1. The Mathematics
Given a surface $\mathcal{M}_t$ parametrized by $\mathbf{X}(\cdot, t)$, the mean curvature flow is governed by:
$$\frac{\partial \mathbf{X}}{\partial t} = H \mathbf{n}$$
where $H$ is the mean curvature and $\mathbf{n}$ is the unit outward normal. The mean curvature itself is defined as the divergence of the normal:
$$H = \nabla \cdot \mathbf{n} = \kappa_1 + \kappa_2$$
where $\kappa_1, \kappa_2$ are the principal curvatures. For a graph surface $z = u(x,y,t)$, the flow reduces to a quasilinear PDE:
$$\frac{\partial u}{\partial t} = \left(1 + |\nabla u|^2\right)^{-1} \left[ \left(1 + u_y^2\right) u_{xx} - 2 u_x u_y u_{xy} + \left(1 + u_x^2\right) u_{yy} \right]$$
In the small-gradient regime $|\nabla u| \ll 1$, this simplifies to the linear heat equation:
$$\frac{\partial u}{\partial t} \approx \Delta u = u_{xx} + u_{yy}$$
The area functional along the flow is:
$$A(t) = \iint \sqrt{1 + |\nabla u|^2} , dx , dy$$
and it satisfies the monotone decrease:
$$\frac{dA}{dt} = -\int H^2 , d\mu \leq 0$$
2. The Concrete Example
We take an initial surface that is a bumpy perturbation of the flat plane:
$$u_0(x, y) = A_0 \sin!\left(\frac{2\pi x}{L}\right) \sin!\left(\frac{2\pi y}{L}\right)$$
on the domain $[-L/2, L/2]^2$ with periodic boundary conditions. Under MCF, this surface will smooth out — the mean curvature drives the bumps flat — and converge to the minimal surface $u \equiv 0$.
We solve this with both the full nonlinear MCF and the linearized heat equation and compare them.
3. Source Code
1 | # ============================================================ |
4. Code Walkthrough
4-1. Grid & Initial Condition
We discretize $[-L, L]^2$ with $N=64$ points per axis with periodic boundary conditions throughout. The initial surface is:
$$u_0(x,y) = 0.5,\sin!\left(\frac{\pi x}{L}\right)\sin!\left(\frac{\pi y}{L}\right)$$
This is a single-mode double-sine bump with amplitude $A_0=0.5$.
4-2. Spectral Laplacian (Speed Optimization)
For the linearized version we use the FFT-based Laplacian, which is $O(N^2 \log N)$ vs $O(N^2)$ per grid point for finite differences:
$$\widehat{\Delta u}(\mathbf{k}) = -|\mathbf{k}|^2 \hat{u}(\mathbf{k})$$
This gives spectral accuracy (exponential convergence) for smooth periodic functions.
4-3. Nonlinear MCF Right-Hand Side
The full nonlinear operator:
$$F(u) = \frac{(1+u_y^2)u_{xx} - 2u_xu_yu_{xy} + (1+u_x^2)u_{yy}}{1 + u_x^2 + u_y^2}$$
is discretized using second-order periodic finite differences with np.roll for all stencil shifts. The mixed derivative $u_{xy}$ uses the standard 4-point cross stencil.
4-4. Area & $\int H^2,d\mu$ Diagnostics
The area functional:
and the Willmore-type integral $\int H^2,d\mu$ are tracked at every save_every=50 steps. By the first variation formula, $dA/dt = -\int H^2,d\mu$, so this acts as a consistency check.
4-5. Time Integration & Stability
We use explicit Euler with $\Delta t = 2\times10^{-4}$. The CFL-like stability condition for MCF on a graph is roughly:
$$\Delta t \lesssim \frac{\Delta x^2}{2}$$
With $\Delta x = 2L/N \approx 0.0625$ and $\Delta x^2/2 \approx 0.002$, our step $2\times10^{-4}$ is safely inside the stability region.
4-6. Linear Decay Theory
For the linearized equation $\partial_t u = \Delta u$ with initial mode $\sin(\pi x/L)\sin(\pi y/L)$, the exact solution is:
$$u(x,y,t) = A_0 e^{-\lambda t}\sin!\left(\frac{\pi x}{L}\right)\sin!\left(\frac{\pi y}{L}\right), \quad \lambda = 2!\left(\frac{\pi}{L}\right)^2$$
We overlay this theoretical curve on the $\max|u|$ plot to validate the numerics.
5. Figure-by-Figure Analysis
Figure 1 — 3D Surface Evolution
Six 3D snapshots from $t=0$ to $t=T_{\rm end}=0.30$ show the surface progressively flattening. The initial double-sine bump smooths monotonically — there is no singularity formation for this convex initial data. The plasma colormap encodes height, making the decay of amplitude visually striking.
Figure 2 — Diagnostics
Left — Area $A(t)$: Area decreases monotonically toward the flat-plane area $(2L)^2 = 16$, confirming the fundamental inequality $dA/dt \leq 0$.
Center — $\int H^2,d\mu$ (log scale): The $L^2$ norm of mean curvature decays exponentially. This is consistent with the analytic eigenmode decay:
$$\int H^2,d\mu \sim e^{-2\lambda t}$$
Right — Amplitude $\max|u|$ (log scale): The nonlinear MCF (solid) and the linear theory (dashed white) overlap almost perfectly for small amplitude, confirming the linearization is accurate for $A_0 = 0.5$ and that the nonlinear corrections are higher order in $A_0$.
Figure 3 — Mean Curvature Heatmap
At $t=0$, the mean curvature field $H(x,y)$ mirrors the surface shape — it is positive at the peak and negative in the valleys (for our sign convention). By $t=T_{\rm end}$, $H$ is essentially zero everywhere: the surface has reached the minimal surface $u \equiv 0$.
Figure 4 — Cross-Section $y=0$
The 1D slice $u(x, 0, t)$ confirms pointwise convergence to zero. The profile transitions from a full sine arch to a flat line, with the decay rate governed by the leading eigenvalue $\lambda = 2\pi^2/L^2$.
6. Execution Results
Steps: 1499 | Final max|u|: 1.1767e-01 Initial area: 18.25907 | Final area: 16.13568
[Fig 1 saved]
[Fig 2 saved]
[Fig 3 saved]
[Fig 4 saved]
7. Singularity Analysis Notes
For more complex initial data — for example a dumbbell surface (two spheres connected by a thin neck) — the neck pinches in finite time, forming a Type-I singularity modeled by the shrinking cylinder $\mathbb{S}^1 \times \mathbb{R}$:
$$H \sim \frac{1}{\sqrt{2(T-t)}}$$
Our example is chosen to be singularity-free (convex, small amplitude), making it the ideal setting for demonstrating the convergence to a minimal surface. The full singularity analysis — including blowup profile extraction and rescaling — is the subject of the Huisken–Sinestrari and Brendle programs in geometric analysis.
The key theoretical result (Huisken 1984) is:
If $\mathcal{M}_0$ is a compact, convex hypersurface, then MCF contracts it to a round point in finite time $T$, and the rescaled surface $\tilde{\mathcal{M}}_s = (2(T-t))^{-1/2}\mathcal{M}_t$ converges smoothly to a round sphere.
For graph surfaces converging to a flat plane (our case), the analogous statement is that $u(\cdot, t) \to 0$ in $C^\infty$ as $t\to\infty$, which our numerics confirm beautifully.