Computing Geometric Flows

6D.9 Finite Difference Methods and 6D.10 Monte Carlo, Fourier and Finite Element Methods built the numerical methods of Book 6D, finite differences, finite elements and Monte Carlo, on the heat equation and the equations of finance, and ended with two forward glances: simulations of Ricci flow neckpinches discretise a parabolic equation in one space variable, and the stiffness matrix of piecewise linear elements on a triangulated surface is the cotangent Laplacian of discrete geometry. This chapter, the last of Course 6, follows those glances to the end. It describes how the geometric flows of this guide are computed in practice, which features of the analysis survive discretisation and which must be built in by hand, and what computation has contributed to the theory: conjectures suggested by pictures, singularities seen before they were proved, and, through discrete versions of the Ricci flow, tools used in computer graphics, medical imaging and network science.

Computation has a particular place in the story of the Poincaré conjecture. Perelman's proof is not computational, and nobody has simulated a three-dimensional Ricci flow with surgery on a general manifold. But the pictures that guided the field, the dumbbell neckpinch of mean curvature flow (6A.10 Curve Shortening and the First Geometric Flows), the shape of a Ricci flow neck near pinching (11B.6 Rotationally Symmetric Flows and the Standard Solution), the degenerate neckpinch modelled by the Bryant soliton (11B.5 Singularities), were computed long before or alongside their proofs, and the algorithmic consequences of geometrisation now include a solution of the homeomorphism problem for three-manifolds. The chapter ends there.

By the end of this chapter you will be able to:

• discretise curve shortening by a polygon flow, show that regular polygons shrink exactly like circles, and explain the role of tangential redistribution;
• describe the three main approaches to mean curvature flow (parametric finite elements, level sets, phase fields) and their trade-offs, including behaviour through singularities and topology change;
• derive formally why Merriman–Bence–Osher threshold dynamics moves interfaces by curvature;
• discretise the rotationally symmetric Ricci flow, explain the time-step restriction near a neckpinch, and describe what the simulations of Garfinkle–Isenberg and Rubinstein–Sinclair found;
• define circle-packing metrics and the Chow–Luo combinatorial Ricci flow, prove discrete Gauss–Bonnet, and describe discrete uniformisation and its uses;
• define Ollivier's Ricci curvature of a graph and compute it in simple cases;
• describe the numerical computation of Ricci-flat metrics by Ricci–DeTurck flow, and the algorithmic consequences of geometrisation.

11.1Curves: polygons that shrink like circles

The curve shortening flow $\partial_t\gamma = \kappa N$ (6A.10 Curve Shortening and the First Geometric Flows) is a quasilinear heat equation for the position vector: with arclength $s$, $\partial_t\gamma = \partial_s^2\gamma$. The obvious discretisation replaces the curve by a closed polygon with vertices $x_1, \dots, x_N$ and the second derivative by a difference quotient. With $\ell_j = |x_{j+1} - x_j|$, a standard choice is

the difference of unit tangent vectors divided by the length of the "dual edge": a discrete curvature vector. For a regular polygon inscribed in a circle of radius $R$, this vector points to the centre with length exactly $\frac1R$ times a factor tending to $1$ as $N \to \infty$, and with the cruder choice $\frac{x_{j+1} - 2x_j + x_{j-1}}{\ell^2}$ the factor is exactly $1$: a regular polygon shrinks self-similarly with $R^2 = R_0^2 - 2t$, exactly like the circle (Exercise 11.1). So the polygon flow reproduces the model solution.

Two numerical facts distinguish geometric flows from the heat equation. First, the equation $\partial_t\gamma = \kappa N$ does not specify how points move along the curve, since tangential motion is a reparametrisation (the gauge freedom of 6A.4 Fourier Methods and Scale Space in its simplest form). A discretisation chooses a gauge implicitly, and a bad choice lets vertices bunch up where the curve shrinks fastest and spread out elsewhere, which destroys accuracy and eventually makes the polygon cross itself. Practical schemes add tangential redistribution, a tangential velocity that keeps vertices equally spaced or spaced according to curvature; this is the numerical counterpart of DeTurck's trick, a choice of gauge that makes the discrete problem well behaved without changing the geometry. Second, the explicit time step is limited by the smallest edge, $\Delta t \lesssim \ell_{\min}^2$ (the parabolic stability condition of 6D.9 Finite Difference Methods), and edges shrink as the curve does; implicit schemes, or semi-implicit ones in which the curvature is taken at the new time and the geometry at the old, remove the restriction. Gerhard Dziuk (1994) proved convergence, with error estimates, for a semi-discrete finite element scheme of this kind.

11.2Surfaces: parametric, level set and phase field

For mean curvature flow of surfaces there are three approaches, each matching one of the mathematical frameworks of the guide.

Parametric finite elements. Represent the surface by a triangulation and move its vertices by a discrete mean curvature vector. Since the mean curvature vector is the Laplace–Beltrami operator applied to the position, $\vec H = \Delta_\Sigma x$ (9A.8 Submanifolds and Minimal Surfaces), the discretisation is the finite element Laplacian of the triangulated surface, the cotangent Laplacian of 6D.10 Monte Carlo, Fourier and Finite Element Methods:

with $\alpha_{ij}$, $\beta_{ij}$ the angles opposite the edge $ij$ and $A_i$ an area associated with vertex $i$. Dziuk (1990) introduced the method; it is accurate and efficient while the surface is smooth, but it cannot pass through singularities: when a neck pinches, the triangulation degenerates, and the topology change must be performed by hand, the numerical analogue of surgery.

Level sets. Stanley Osher and James Sethian (1988) represented a moving hypersurface as the zero level set of a function $u(x, t)$ on all of space, and moved every level set by mean curvature at once:

the degenerate parabolic equation whose viscosity solutions (6C.9 Hamilton–Jacobi Equations and Viscosity Solutions) were constructed by Evans and Spruck and by Chen, Giga and Goto (1991). The level-set method passes through singularities and topology changes automatically (a dumbbell's zero level set simply separates into two components), and its solutions are unique, but it pays for this: it computes in one dimension more than the surface, and it can fatten, the zero level set acquiring interior when the flow is genuinely non-unique.

Phase fields. The Allen–Cahn equation of 6A.9 Reaction–Diffusion: Fronts, Patterns and Pinching, $u_t = \Delta u - \varepsilon^{-2}W'(u)$, has interfaces of width $\varepsilon$ that move by mean curvature as $\varepsilon \to 0$. Solving it with $\varepsilon$ a few grid spacings computes mean curvature flow, again through singularities, with the convergence theorems of Chen, de Mottoni–Schatzman, Evans–Soner–Souganidis and Ilmanen as justification. Phase-field models are the standard tool in materials science for solidification, grain growth and crack propagation, where interfaces with extra physics (anisotropy, latent heat, elasticity) are modelled by extra terms.

Threshold dynamics. The simplest scheme of all, due to Barry Merriman, James Bence and Stanley Osher (1992), alternates two steps: given a set $\Omega_k$, solve the heat equation for time $h$ with initial data the indicator $1_{\Omega_k}$, and let $\Omega_{k+1}$ be the set where the solution exceeds $\frac12$. Diffusion rounds off corners and threshold sharpens the result. A short computation (Exercise 11.2) shows that each step moves a curved boundary inward by about $h$ times its mean curvature, so the scheme approximates mean curvature flow with time step $h$; Lawrence Evans (1993) and Guy Barles and Christine Georgelin (1995) proved convergence to the viscosity solution. It is the heat equation of 6A.3 The Heat Equation on ℝⁿ and the maximum principle of 6A.5 Maximum Principles, composed into a geometric flow, and it can be computed with fast Fourier transforms (6D.9 Finite Difference Methods).

What computation suggested. Computed flows have repeatedly shown the way for theory. David Chopp (1994) computed numerous self-shrinkers, including the torus whose existence Sigurd Angenent had proved (1992), and families that suggested which shrinkers are stable. Angenent, Tom Ilmanen and Chopp (1995) computed an example of a smooth initial surface in $\mathbb{R}^3$ whose mean curvature flow is non-unique after a singular time, the level-set flow fattening; rigorous examples of fattening from smooth surfaces in $\mathbb{R}^3$ were constructed by Ilmanen and White in 2024. Image processing used the same flows in the other direction: Vicent Caselles, Ron Kimmel and Guillermo Sapiro's geodesic active contours (1997) and Tony Chan and Luminita Vese's level-set segmentation (2001) segment medical and natural images by curvature-driven level-set flows (6A.4 Fourier Methods and Scale Space).

11.3The rotationally symmetric Ricci flow

The Ricci flow on $S^{n+1}$ for metrics $g = ds^2 + \psi(s)^2g_{S^n}$, invariant under rotations about an axis, reduces to the system (11B.6 Rotationally Symmetric Flows and the Standard Solution)

where $ds = \varphi\,dx$ in a fixed coordinate $x$, so that $\partial_s$ and $\partial_t$ do not commute. It is a one-dimensional reaction–diffusion system (6A.9 Reaction–Diffusion: Fronts, Patterns and Pinching), and finite differences in $x$ solve it (6D.9 Finite Difference Methods), with three difficulties that are the difficulties of the analysis seen through a grid.

• The poles. Where the metric closes off, $\psi \to 0$ and the term $\frac{1 - \psi_s^2}{\psi}$ is a ratio of two vanishing quantities, finite only because smoothness forces $\psi_s \to 1$ (11B.6 Rotationally Symmetric Flows and the Standard Solution). Schemes must impose the pole conditions exactly and evaluate the ratio by a Taylor expansion, or they produce spurious cone points.
• The gauge. The coordinate $x$ is fixed while arclength changes, and grid points in the neck are stretched or compressed by the factor $\varphi$. As in the parametric curve flow, a good scheme chooses a gauge: a DeTurck-type vector field, or a reparametrisation that keeps grid points evenly spaced in arclength, or a grid adapted to the neck.
• The time step. The diffusion coefficient in $x$ is $\varphi^{-2}$ and the reaction term has rate $(n - 1)\psi^{-2}$, so an explicit scheme needs $\Delta t \lesssim \min(\varphi^2\Delta x^2, \psi^2)$ (6D.9 Finite Difference Methods). As a neck pinches, $\psi_{\min}^2 \approx 2(n - 1)(T - t)$, so the admissible time step shrinks linearly as the singular time approaches; with adaptive steps, reaching time $T - \varepsilon$ takes a number of steps of order $\log\frac1\varepsilon$ (Exercise 11.3), which is why simulations can follow a neckpinch very close to the singular time.

What the simulations showed. Hyam Rubinstein and Robert Sinclair (2005) computed rotationally symmetric Ricci flows of spheres with necks and visualised them as surfaces of revolution, showing neckpinches forming and confirming the shape predicted by Angenent and Knopf (11B.6 Rotationally Symmetric Flows and the Standard Solution): the neck becomes a round cylinder, with curvature blowing up at the rate $\frac{1}{T - t}$. David Garfinkle and James Isenberg (2003) studied one-parameter families of initial metrics, from dumbbells with thick necks (which become round and shrink to points) to dumbbells with thin necks (which pinch), and found at the threshold between the two a critical solution: a degenerate neckpinch at a pole, in which a small cap pinches off, with the rescaled geometry near the pole looking like the Bryant soliton, the steady soliton of 11B.1 Ricci Solitons. In 2008 they gave more detailed numerical evidence for this Bryant-soliton picture. Sigurd Angenent, James Isenberg and Dan Knopf later proved (2015) that degenerate neckpinches exist and have exactly this structure: type II singularities modelled on the Bryant soliton, as the computations suggested. The degenerate neckpinch is the Ricci flow counterpart of the critical phenomena that Matthew Choptuik had found numerically in 1993 in gravitational collapse (with self-similar solutions at the threshold of black-hole formation), and Garfinkle and Isenberg's work was explicitly motivated by that analogy.

11.4Discrete Ricci flow on surfaces

On a triangulated surface, there is a way to define curvature and a Ricci flow without any limit $h \to 0$, by making the metric itself discrete. William Thurston proposed in 1985 to approximate conformal maps by circle packings: place a circle of radius $r_i$ at each vertex $i$ of a triangulation, tangent to the circles at neighbouring vertices, so that each triangle $ijk$ becomes a Euclidean triangle with side lengths $r_i + r_j$, $r_j + r_k$, $r_k + r_i$. The discrete Gauss curvature at vertex $i$ is the angle deficit

with $\theta_i^{jk}$ the angle at $i$ in triangle $ijk$: zero when the triangles around $i$ fit flat, positive when they close up like the tip of a cone. These curvatures satisfy a discrete Gauss–Bonnet theorem, $\sum_iK_i = 2\pi\chi$, for the Euler characteristic $\chi = V - E + F$ (Exercise 11.4, 8A.9 The Curvature of Surfaces). Burton Rodin and Dennis Sullivan (1987) proved Thurston's conjecture that hexagonal circle packings of a domain, as the circles shrink, converge to the Riemann map (5A.4 Harmonic Functions and Conformal Mapping).

Bennett Chow and Feng Luo (2003) defined the combinatorial Ricci flow

the analogue of the two-dimensional Ricci flow $\partial_tg = -2Kg$, which shrinks the metric where curvature is positive and expands it where it is negative, with a normalised version that preserves the total area. They proved that it is the gradient flow of a convex function (in the variables $\log r_i$), that it exists for all time, and that it converges exponentially fast to a circle-packing metric of constant curvature exactly when the conditions of Thurston's circle-packing theorem hold; this reproves Thurston's theorem by a flow, as Hamilton's flow on surfaces proves uniformisation (11A.12 Ricci Flow on Surfaces). Feng Luo (2004) introduced a discrete Yamabe flow that changes edge lengths by vertex scaling, $\ell_{ij} = e^{u_i + u_j}\ell^0_{ij}$; Alexander Bobenko, Ulrich Pinkall and Boris Springborn (2015) found the convex variational principle behind it, through the volumes of ideal hyperbolic polyhedra; and Xianfeng David Gu, Feng Luo, Jian Sun and Tianqi Wu (2018) proved a discrete uniformisation theorem: every polyhedral surface is discretely conformal to one of constant curvature, after allowing changes of triangulation, the discrete counterpart of surgery-free uniformisation.

11.5Ricci curvature of graphs and networks

Discrete curvature has also been defined directly on graphs, without any triangulation, by the transport characterisation of 6C.11 Optimal Transport and Ricci Curvature. Yann Ollivier (2009) defined the coarse Ricci curvature of an edge $xy$ of a graph as

where $m_x$ is the uniform probability measure on the neighbours of $x$ and $W_1$ the transport distance with cost the graph distance. Positive curvature means that the neighbourhoods of $x$ and $y$ are closer than $x$ and $y$ are: a random step from each brings them together, the discrete form of the von Renesse–Sturm coupling of 6D.5 Diffusions on Manifolds and Curvature. Complete graphs have positive curvature, cycles have zero curvature, and trees have negative curvature (Exercise 11.5). Ollivier proved that positive coarse curvature implies concentration inequalities and spectral gaps, as $\operatorname{Ric} \geq K$ does on manifolds.

Network scientists have used these curvatures and discrete Ricci flows on graphs, which shrink positively curved edges and stretch negatively curved ones, to find communities in social and biological networks (edges between communities are bottlenecks, with negative curvature, and the flow stretches them until the communities separate; Chien-Chun Ni, Yu-Yao Lin, Feng Luo and Jie Gao, 2019), to measure the robustness of financial and internet networks, and, in machine learning, to diagnose over-squashing in graph neural networks, where negatively curved bottleneck edges limit how information propagates (Jake Topping and collaborators, 2022). These are discrete analogues of the guide's central theme: curvature controls how diffusion spreads, and a flow driven by curvature exposes the structure of a space.

11.6Computing Ricci-flat metrics

The Ricci flow is also a numerical method for finding special metrics. Calabi–Yau manifolds, which carry Ricci-flat Kähler metrics by Yau's theorem (6C.10 Fully Nonlinear Equations and Monge–Ampère), are needed explicitly in string theory, where the geometry of the extra dimensions determines the physics of the four-dimensional world, but Yau's proof is not constructive and no Ricci-flat metric on a compact Calabi–Yau manifold of this kind is known in closed form. Matthew Headrick and Toby Wiseman (2005) computed the Ricci-flat metric on the K3 surface numerically by discretising it on a lattice and running a Ricci–DeTurck-type flow (and its Kähler variant) to convergence, and Headrick and Wiseman (2006) used Ricci flow to study black-hole metrics in Euclidean quantum gravity, where unstable directions of the flow correspond to instabilities of the black holes. Simon Donaldson (2009) gave an alternative algorithm through balanced embeddings into projective space, and from about 2020 machine-learned approximations of Calabi–Yau metrics have been developed. The DeTurck gauge (6B.10 Harmonic Coordinates and the Regularity of Metrics, 11A.4 Short-Time Existence and Uniqueness) is what makes these flows strictly parabolic, and therefore computable with the schemes of 6D.9 Finite Difference Methods.

11.7From flows to algorithms

Geometrisation turned three-dimensional topology into a subject with algorithms. Before Perelman, Hyam Rubinstein (1992) and Abigail Thompson (1994) had found an algorithm to recognise the three-sphere, using normal surfaces rather than geometry; Jeffrey Weeks's program SnapPea (from the 1980s, now SnapPy) computes hyperbolic structures on three-manifolds from triangulations by solving Thurston's gluing equations, and it showed experimentally that "most" three-manifolds are hyperbolic, as the geometrisation conjecture predicts (10A.9 Hyperbolic Three-Manifolds). After Perelman, the homeomorphism problem for closed three-manifolds is decidable: given two triangulated three-manifolds, an algorithm decides whether they are homeomorphic. Greg Kuperberg (2019) gave a careful proof of this as a corollary of geometrisation, assembling the geometric decomposition, the recognition of each geometric piece and Mostow rigidity (10A.9 Hyperbolic Three-Manifolds) into an algorithm. The Ricci flow itself is not the algorithm; but the theorem it proves is what makes the algorithm exist. The flow of this guide, a heat equation for the shape of space, ends in a statement that a computer can use.

11.8History

Thurston proposed circle packing in 1985, and Rodin and Sullivan proved convergence in 1987. Osher and Sethian introduced level-set methods in 1988; Dziuk computed mean curvature flow by finite elements in 1990; Evans and Spruck and Chen, Giga and Goto constructed level-set solutions in 1991; Merriman, Bence and Osher proposed threshold dynamics in 1992, and Evans (1993) and Barles and Georgelin (1995) proved its convergence. Brakke released the Surface Evolver in 1992, and Weaire and Phelan found their foam in 1993. Choptuik discovered critical collapse in 1993. Chopp computed self-shrinkers in 1994, and Angenent, Ilmanen and Chopp their non-uniqueness example in 1995. Rubinstein (1992) and Thompson (1994) found three-sphere recognition. Chow and Luo's combinatorial Ricci flow appeared in 2003, Garfinkle and Isenberg's simulations from 2003, Luo's discrete Yamabe flow in 2004, Rubinstein and Sinclair's visualisations in 2005, Headrick and Wiseman's numerical K3 metric in 2005, Ollivier's coarse Ricci curvature in 2009, Bobenko, Pinkall and Springborn's variational principle in 2015, Angenent, Isenberg and Knopf's degenerate neckpinches in 2015, Gu, Luo, Sun and Wu's discrete uniformisation in 2018, and Kuperberg's account of the homeomorphism algorithm in 2019.

11.9Exercises

Solution

(a) Writing $x_j$ as the complex number $Re^{2\pi ij/N}$, $x_{j+1} + x_{j-1} = x_j(e^{2\pi i/N} + e^{-2\pi i/N}) = 2\cos\frac{2\pi}{N}x_j$, so the second difference is $(2\cos\frac{2\pi}N - 2)x_j = -4\sin^2\frac\pi Nx_j$. (b) $\frac{dx_j}{dt} = -\frac{4\sin^2(\pi/N)}{4R^2\sin^2(\pi/N)}x_j = -\frac{x_j}{R^2}$. So every vertex moves radially with $\frac{d}{dt}x_j = -\frac{x_j}{R^2}$; the polygon stays regular, and $\frac{dR}{dt} = -\frac1R$. The circle under curve shortening has curvature vector $-\frac{x}{R^2}$ too.

Solution

The value is $\frac12$ when $y_0$ is the median of $-\frac\kappa2X^2 - Y$. For small $\kappa\sqrt h$, the distribution of $-\frac\kappa2X^2 - Y$ is that of $-Y$ (symmetric, median $0$) shifted by a small amount; to first order the median shifts by the mean of the perturbation, $-\frac\kappa2\mathbb{E}X^2 = -\frac\kappa2\cdot2h = -\kappa h$. So the $\frac12$-level passes through $(0, -\kappa h)$: the boundary has moved a distance $\kappa h$ into the set. Over time $h$ this is normal velocity $\kappa$, curve shortening (and in higher dimensions mean curvature flow, with $\kappa$ replaced by the sum of the principal curvatures).

Solution

(a) The step at time $t$ is $\Delta t = 2(n - 1)c(T - t)$, so the number of steps is approximately $\int_0^{T-\varepsilon}\frac{dt}{2(n - 1)c(T - t)} = \frac{1}{2(n - 1)c}\log\frac{T}{\varepsilon}$: each step removes a fixed fraction of the remaining time. (b) A fixed step $\Delta t$ violates stability once $2(n - 1)(T - t) < \frac{\Delta t}{c}$, so it cannot get closer to $T$ than a distance of order $\Delta t$; reaching $T - \varepsilon$ requires $\Delta t \lesssim \varepsilon$, hence of order $\frac T\varepsilon$ steps, exponentially more than the adaptive scheme.

Solution

(a) $\sum_iK_i = \sum_i2\pi - \sum_i\sum\theta = 2\pi V - (\text{sum of all angles of all triangles}) = 2\pi V - \pi F$. (b) Each triangle has three edges and each edge lies in two triangles, so $3F = 2E$ and $\pi F = 2\pi(E - F)$; then $2\pi V - \pi F = 2\pi(V - E + F)$. (c) Tetrahedron: three equilateral angles $\frac\pi3$ at each vertex, deficit $\pi$, total $4\pi = 2\pi\cdot2$. Cube: three right angles at each vertex, deficit $\frac\pi2$, eight vertices, total $4\pi$. (The diagonals add no curvature: the triangles in each face fit flat.)

Solution

(a) $m_x$ and $m_y$ both put mass $\frac{1}{n-1}$ on each of the $n - 2$ common neighbours; $m_x$ also puts $\frac{1}{n-1}$ on $y$ and $m_y$ puts $\frac{1}{n-1}$ on $x$. Moving that mass from $y$ to $x$ costs $\frac{1}{n-1}$, and no plan can do better (the two measures differ by $\frac{1}{n-1}(\delta_y - \delta_x)$ and these points are at distance $1$), so $W_1 = \frac{1}{n-1}$ and $\kappa = \frac{n - 2}{n - 1}$. (b) For $y = x + 1$: $m_x = \frac12(\delta_{x-1} + \delta_{x+1})$, $m_y = \frac12(\delta_x + \delta_{x+2})$. Moving $x - 1 \to x$ and $x + 1 \to x + 2$ costs $\frac12 + \frac12 = 1$; by the $1$-Lipschitz test function $f$ taking the values $-1, 0, 1, 2$ at $x - 1, x, x + 1, x + 2$ (extendable to a $1$-Lipschitz function on the cycle when $N \geq 6$), $W_1 \geq \int f\,dm_y - \int f\,dm_x = 1$. So $\kappa = 0$. (c) $m_x$ puts $\frac13$ on $y$ and on two other neighbours $a$, $b$; $m_y$ puts $\frac13$ on $x$ and two others $c$, $d$. The plan $a \to x$, $y \to c$, $b \to d$ costs $\frac13(1 + 1 + 3) = \frac53$. Conversely, the function equal to $-1$ at $a$, $b$, $0$ at $x$, $1$ at $y$ and $2$ at $c$, $d$ is $1$-Lipschitz, and $\int f\,dm_y - \int f\,dm_x = \frac43 - (-\frac13) = \frac53$, so $W_1 = \frac53$ and $\kappa = 1 - \frac53 = -\frac23 < 0$. The tree's branching pulls neighbourhoods apart, the discrete form of negative curvature.