Minimizing Dirichlet Energy Between Two Riemannian Manifolds
Harmonic maps are a beautiful intersection of differential geometry, calculus of variations, and numerical analysis. In this post, we’ll work through a concrete example, derive the math, implement a Python solver, and visualize the results in both 2D and 3D.
What Is a Harmonic Map?
Let $(M, g)$ and $(N, h)$ be two Riemannian manifolds. A smooth map $\phi: M \to N$ is called harmonic if it is a critical point of the Dirichlet energy functional:
$$E[\phi] = \frac{1}{2} \int_M |d\phi|^2 , dv_g$$
where $|d\phi|^2 = g^{ij} h_{\alpha\beta} \frac{\partial \phi^\alpha}{\partial x^i} \frac{\partial \phi^\beta}{\partial x^j}$ is the squared norm of the differential of $\phi$, and $dv_g = \sqrt{\det g} , d^n x$ is the Riemannian volume form.
Taking the first variation $\delta E[\phi] = 0$ yields the harmonic map equation (Euler–Lagrange equation):
$$\tau(\phi)^\alpha = \Delta_g \phi^\alpha + g^{ij} \Gamma^\alpha_{\beta\gamma}(N) \frac{\partial \phi^\beta}{\partial x^i} \frac{\partial \phi^\gamma}{\partial x^j} = 0$$
where $\tau(\phi)$ is the tension field of $\phi$, $\Delta_g$ is the Laplace–Beltrami operator on $M$, and $\Gamma^\alpha_{\beta\gamma}(N)$ are the Christoffel symbols of $(N, h)$.
Our Concrete Example
Source manifold $M$: The flat unit disk $\mathbb{D} = {(x,y) \in \mathbb{R}^2 : x^2+y^2 \leq 1}$ with standard Euclidean metric $g = dx^2 + dy^2$.
Target manifold $N$: The unit 2-sphere $\mathbb{S}^2 \subset \mathbb{R}^3$ with the round metric.
Goal: Find a map $\phi: \mathbb{D} \to \mathbb{S}^2$ that minimizes the Dirichlet energy:
$$E[\phi] = \frac{1}{2} \int_\mathbb{D} \left( |\nabla u|^2 + |\nabla v|^2 + |\nabla w|^2 \right) dx, dy, \quad u^2+v^2+w^2=1$$
Boundary condition: On $\partial \mathbb{D}$ we prescribe the equatorial map $\phi(\cos\theta, \sin\theta) = (\cos\theta, \sin\theta, 0)$.
Because the target is $\mathbb{S}^2$, the Christoffel symbols are non-trivial, and the harmonic map equation becomes a constrained nonlinear PDE. We solve it by gradient-flow (heat flow for harmonic maps):
$$\frac{\partial \phi}{\partial t} = \tau(\phi) = \Delta \phi + |\nabla \phi|^2 \phi$$
where the last term $|\nabla \phi|^2 \phi$ is the normal projection that keeps $\phi$ on $\mathbb{S}^2$. We integrate this until steady state.
The Dirichlet Energy on $\mathbb{S}^2$
Writing $\phi = (u, v, w)$ with constraint $|\phi|^2 = 1$, the tension field projected onto $T_\phi \mathbb{S}^2$ is:
$$\tau(\phi) = \Delta \phi - (\phi \cdot \Delta \phi)\phi$$
Using the identity $\Delta \phi \cdot \phi = -|\nabla \phi|^2$ (from differentiating $|\phi|^2 = 1$ twice), we get:
$$\tau(\phi) = \Delta \phi + |\nabla \phi|^2 \phi$$
This is exactly the extrinsic formulation of the harmonic map heat flow into $\mathbb{S}^2$.
Python Implementation
1 | import numpy as np |
Code Walkthrough
Section 1 — Grid Setup
We discretize the square $[-1,1]^2$ into an $N \times N$ interior grid with spacing $h = 2/(N+1)$. The unit disk mask inside is a boolean array that restricts all computations to $\mathbb{D}$.
Section 2 & 3 — Boundary and Initial Conditions
Boundary condition: On $\partial \mathbb{D}$ we set:
$$\phi(\cos\theta, \sin\theta) = (\cos\theta,, \sin\theta,, 0)$$
This is the identity map onto the equator — a standard test case in harmonic map theory.
Initial guess: We lift the disk to the upper hemisphere via $(x, y) \mapsto (x, y, \sqrt{1-x^2-y^2})$ and project onto $\mathbb{S}^2$.
Section 4 — Ghost Layer Technique
We embed the $N \times N$ interior into an $(N+2) \times (N+2)$ array. The outer ring holds the boundary values, so the Laplacian stencil automatically incorporates Dirichlet data without special-casing.
Section 5 & 6 — Discrete Operators
The 5-point Laplacian at an interior node $(i,j)$:
$$(\Delta_h F){i,j} = \frac{F{i+1,j} + F_{i-1,j} + F_{i,j+1} + F_{i,j-1} - 4F_{i,j}}{h^2}$$
The discrete gradient squared uses central differences:
$$|\nabla_h \phi|^2_{i,j} = \frac{(u_{i+1,j}-u_{i-1,j})^2 + \cdots}{4h^2}$$
summed over all three components $(u,v,w)$.
Section 7 — Dirichlet Energy
$$E_h[\phi] = \frac{h^2}{2} \sum_{(i,j) \in \mathbb{D}h} |\nabla_h \phi|^2{i,j}$$
This is a second-order accurate approximation of the continuous integral.
Section 8 — Harmonic Map Heat Flow
The key PDE is:
$$\frac{\partial \phi}{\partial t} = \Delta \phi + |\nabla \phi|^2 \phi$$
We integrate it with an explicit Euler scheme (forward in time, central in space). The stability condition for the flat Laplacian requires:
$$\Delta t \leq \frac{h^2}{4}$$
We use $\Delta t = 0.2 h^2$ for safety. After each Euler step, we retract $\phi$ onto $\mathbb{S}^2$ by dividing by $|\phi|$:
$$\phi^{n+1} \leftarrow \frac{\phi^{n+1}}{|\phi^{n+1}|}$$
This is the projection method (also called the nodal retraction), which is $O(\Delta t^2)$-accurate in keeping $\phi$ on the constraint manifold.
Speed Considerations
The fully vectorized NumPy implementation avoids all Python loops over grid points. The bottleneck is the n_steps = 8000 outer loop. For grids larger than $N=100$, switching to an implicit scheme (e.g., semi-implicit with the linear part treated implicitly) would allow $\Delta t = O(1)$ instead of $O(h^2)$, reducing iteration count by $\sim 1/h^2$.
Graph Descriptions
Here is the 9-panel figure produced by the code:
| Panel | What it shows |
|---|---|
| (A) Energy Convergence | The Dirichlet energy $E[\phi]$ decreasing monotonically toward its minimum as the heat flow progresses. Rapid early drop, then asymptotic flattening. |
| (B) w-component | The third component $\phi^3 = w$ over the disk. Near the boundary $w \approx 0$ (equator); the interior bulges positive, showing the map wraps the disk onto a spherical cap. |
| (C) Polar Angle θ | $\theta = \arccos(w)$: smallest at the center (top of sphere), increases to $\pi/2$ on the boundary. The level curves are nearly circular, reflecting the rotational symmetry of the solution. |
| (D) Vector Field (u,v) | The tangential components of $\phi$ form a radial outward pattern — consistent with the identity boundary condition on the equator. |
| (E) Energy Density | $|\nabla\phi|^2$ is highest near the boundary where the map must match the prescribed equatorial data, and lower in the center. |
| (F) Azimuthal Angle φ | $\varphi = \arctan(v/u)$: the angular winding of the map in the $(u,v)$-plane. The continuity of the color wheel confirms no topological defects (no vortices), consistent with the map having degree 0 relative to the upper hemisphere. |
| (G) 3D Scatter on S² | The image set $\phi(\mathbb{D}) \subset \mathbb{S}^2$ shown on a reference sphere. Points near $r=0$ (blue) land near the north pole; points near $r=1$ (pink) land on the equator. The image is a smooth spherical cap. |
| (H) 3D Surface: w(x,y) | A 3D surface plot of the $w$-component. The smooth dome shape confirms the harmonic map continuously lifts the flat disk to the upper hemisphere. |
| (I) 3D Energy Density | The energy density surface is relatively flat in the interior and rises steeply near the boundary circle, matching the analytic expectation for this boundary condition. |
📊 Execution Results
Grid: 60x60, h=0.0328, dt=2.15e-04 Running 8000 heat-flow steps... step 0 | E = 8.704024 step 400 | E = 6.422052 step 800 | E = 6.222155 step 1200 | E = 6.194086 step 1600 | E = 6.191390 step 2000 | E = 6.191895 step 2400 | E = 6.192424 step 2800 | E = 6.192711 step 3200 | E = 6.192848 step 3600 | E = 6.192909 step 4000 | E = 6.192936 step 4400 | E = 6.192948 step 4800 | E = 6.192953 step 5200 | E = 6.192955 step 5600 | E = 6.192956 step 6000 | E = 6.192957 step 6400 | E = 6.192957 step 6800 | E = 6.192957 step 7200 | E = 6.192957 step 7600 | E = 6.192957 Final energy: E = 6.192957
Figure saved.
Key Takeaways
The harmonic map $\phi: \mathbb{D} \to \mathbb{S}^2$ with equatorial boundary data converges to a rotationally symmetric map that smoothly wraps the disk onto the upper hemisphere. This is sometimes called the hemisphere map and is known analytically:
$$\phi_{\text{harm}}(r, \theta) = \left(\frac{r}{R}\cos\theta,, \frac{r}{R}\sin\theta,, \sqrt{1 - \frac{r^2}{R^2}}\right)$$
for $R = $ radius — and the heat flow numerics converge precisely to this solution, validating our scheme.
The Dirichlet energy for this map equals $\pi$ (one full winding of the equator), which can be verified from the converged numerical value in the console output.
