Numerical Analysis via Discrete Exterior Calculus

0. Why geometry matters for computation

Numerical analysis is the art of replacing continuous objects with finite ones and controlling the error. Derivatives become differences, integrals become sums, differential equations become linear systems. The standard toolkit for making these replacements (finite differences, finite elements, spectral methods) was developed over three centuries, and each tool comes with its own convergence theorems, its own stability conditions, and its own failure modes.

This article presents a different perspective. The classical tools are all, in a precise sense, incomplete versions of a single geometric construction: discrete exterior calculus (DEC). When a classical method works, it is because it accidentally preserves some fragment of the geometric structure of the underlying continuous problem. When it fails, the failure can be diagnosed as a geometric structure violation. DEC preserves all of this structure by construction.

To make this claim tangible, we begin with a failure.

0.1. The Poisson equation from physics

Consider a thin metal plate lying flat in the plane . The plate occupies a bounded region , and its edges are held at fixed temperatures. Heat flows from hot to cold. After waiting long enough, the temperature distribution reaches a steady state where nothing changes. What equation governs this steady state?

Two physical facts are sufficient. The first is Fourier's law: heat flux is proportional to the negative temperature gradient. If is the temperature at position , the heat flux vector is

where is the thermal conductivity (a material property we can set to 1 without loss of generality). Heat flows from high temperature to low temperature, in the direction of steepest descent.

The second is conservation of energy: in steady state, the net heat flux out through the boundary of any sub-region equals the total heat generated inside it. For a sub-region with boundary curve , this is

where is the outward unit normal to , is the arc length element along the boundary, is the area element, and is the heat source density (energy generated per unit area per unit time). The symbol denotes a closed line integral: an integral taken along a curve that forms a closed loop (here, the boundary of the sub-region). It is computed exactly like an ordinary line integral, but the circle emphasises that the path returns to its starting point.

The bridge between this integral statement and a pointwise equation is the divergence theorem. Since our thin plate lives in , we state the theorem in this setting first, then generalise.

Applying the divergence theorem to [eq:conservation] with , and noting that the equation holds for every sub-region , we conclude that the integrand itself must vanish:

The operator is called the Laplacian. Writing for the unknown function (since this equation appears far beyond heat conduction), we obtain the Poisson equation:

This derivation generalises immediately to . The same argument (a flux law plus conservation plus the divergence theorem) produces the Poisson equation in any dimension, provided we use the correct form of the divergence theorem.

The key point is that the boundary is always one dimension lower than the region , and the theorem relates an integral over the interior to an integral over this lower-dimensional boundary.

The Poisson equation in is then

In the Laplacian is ; in it adds . When (no sources), the equation reduces to Laplace's equation , and the solution is called a harmonic function.

0.2. The five-point stencil

To solve [eq:poisson-intro] on a computer, we need to approximate the Laplacian with a finite number of operations. The oldest and most natural approach is finite differences. Place a uniform rectangular grid over the domain with spacing in both directions. At each interior grid point , approximate the second derivatives using the Taylor expansion:

and similarly for . Adding the two gives the five-point stencil for the discrete Laplacian:

This approximation is accurate, meaning the error between and at each grid point is bounded by a constant times , provided is sufficiently smooth.

On a square domain with interior grid points per side, the discrete Poisson equation is a sparse linear system of equations in unknowns, solvable in operations via fast methods.

For rectangular domains, this works beautifully. The trouble begins when the domain is not rectangular.

0.3. The failure on an L-shaped domain

Consider the L-shaped domain formed by removing the upper-right quadrant from the unit square: . This domain has a re-entrant corner at the origin, where the interior angle is .

The exact solution to the Laplace equation on this domain (with appropriate boundary conditions) has a singularity at the re-entrant corner. In polar coordinates centred at the origin, the leading term of the solution behaves as

The exponent means the gradient diverges as near the corner. The solution is continuous but not smooth.

Now apply the five-point stencil [eq:five-point] to this domain. Two things go wrong:

1. Staircase boundary. The rectangular grid cannot represent the re-entrant corner exactly. Grid points near the corner are either inside the domain or outside it, and the corner itself falls between grid lines. The discrete domain is a staircase approximation that introduces an geometric error at every boundary cell.

2. Regularity failure. The five-point stencil's error estimate requires (four continuous derivatives). But the exact solution has unbounded derivatives at the corner. The stencil does not know this: it treats the corner the same as any interior point, and the local error there is much larger than .

The result is a computed solution that is globally degraded. The convergence rate drops from to roughly , and the error contaminates the entire domain, not just the neighbourhood of the corner. Refining the grid helps, but much more slowly than the theory for smooth problems predicts.

This is not a minor inconvenience. Non-rectangular domains with corners and curved boundaries are the rule in applications, not the exception. The L-shaped domain is the simplest possible example of a universal problem: the rectangular grid imposes a geometry on the computation that may have nothing to do with the geometry of the problem.

The question this article answers is: what geometric structure did the rectangular grid destroy, and how do we build numerical methods that preserve it? The answer involves replacing the rectangular grid with a simplicial complex, replacing sampled function values with cochains, and replacing the Laplacian with the Hodge Laplacian. We begin with the simplicial complex.

1. Simplicial complexes

The rectangular grid is a special case of a more general idea: dividing a domain into simple geometric pieces. The trouble with rectangles is that they are rigid. They cannot conform to curved boundaries without staircase artifacts, and they impose a preferred set of coordinate directions that may have nothing to do with the problem. Triangles (and their higher-dimensional analogues) are flexible: any domain in can be triangulated, with the mesh conforming exactly to the boundary. This section develops the combinatorial machinery for working with triangulations.

1.1. Simplices

The starting point is the simplest possible geometric building block: the generalisation of a triangle to arbitrary dimensions. To define it, we first need the notion of points in "general position."

1.2. Orientation

To define the boundary operator, we need to distinguish between the two "directions" a simplex can be traversed. An edge can be traversed from to or from to . A triangle can be traversed clockwise or counter-clockwise.

1.3. Simplicial complexes

A single simplex is a building block. A simplicial complex is a collection of simplices glued together along shared faces, subject to two compatibility conditions that prevent overlaps and dangling pieces.

1.4. Chains and the boundary operator

The central algebraic fact about the boundary operator is:

1.5. Boundary matrices

Once we fix an ordering of the simplices in each dimension, the boundary operator is a matrix. The entry in row (a -simplex) and column (a -simplex) is if appears in with a positive sign, if it appears with a negative sign, and otherwise.

2. Cochains and the coboundary operator

Chains represent discrete pieces of a domain: vertices, edges, triangles, and their formal combinations. To do numerical analysis, we need to assign numerical values to these pieces. A temperature at each vertex. A flux across each edge. A charge density on each triangle. These assignments are cochains, and they are the discrete analogues of differential forms.

2.1. Cochains

A cochain formalises the idea of assigning measurements to the simplices of a mesh. The definition is dual to that of a chain: where chains are formal sums of simplices, cochains are linear functions that evaluate on chains to produce real numbers.

2.2. The coboundary operator

The boundary operator maps -chains to -chains. Its dual, the coboundary operator, maps -cochains to -cochains. It is the discrete analogue of the exterior derivative.

2.3. Finite differences recovered

The coboundary operator on a 1D mesh is the finite difference operator. This is the first instance of a classical numerical method emerging from the DEC framework.

2.4. Interpolation error of the coboundary

The coboundary of a sampled function approximates the integral of the derivative. How well?

2.5. Example: the 1D simple harmonic oscillator

We now apply the coboundary operator to a genuine eigenvalue problem. The one-dimensional quantum harmonic oscillator is governed by the time-independent Schrodinger equation:

with boundary conditions as . The exact eigenvalues are for , and the eigenfunctions are Hermite functions.

Discretisation. Take a uniform mesh on with vertices at , . Impose (the wavefunction vanishes at the boundary, which is a good approximation for large enough). The unknowns are the interior values for .

The second derivative is approximated by the coboundary composed with its transpose. With being the coboundary matrix (transpose of the boundary matrix from Section 1), the kinetic energy operator is

which is the familiar tridiagonal matrix with on the diagonal and on the super- and sub-diagonals. The potential is a diagonal matrix . The discrete Hamiltonian is

Solving the eigenvalue problem (a symmetric tridiagonal eigenproblem, cost ) gives the approximate eigenvalues.

Results. With and (), the first five eigenvalues are:

The eigenvalue error grows with (higher modes oscillate faster, requiring finer meshes to resolve) and decreases as under mesh refinement (consistent with the second-order accuracy of the second-difference stencil ).

use cartan::simplicial::SimplicialComplex1D;
use cartan::dec::CoboundaryOperator;
use nalgebra::{DMatrix, DVector, SymmetricEigen};

// Build 1D mesh on [-8, 8] with 201 vertices
let n = 200;
let l = 8.0_f64;
let h = 2.0 * l / n as f64;
let mesh = SimplicialComplex1D::uniform(-l, l, n);

// Coboundary d0: N x (N+1) matrix (transpose of boundary)
let d0 = mesh.coboundary_matrix(0);

// Kinetic operator: T = -(1/2h^2) d0^T d0, restricted to interior
let d0_int = d0.columns(1, n - 1); // drop boundary columns
let kinetic = &d0_int.transpose() * &d0_int * (-0.5 / (h * h));

// Potential: diagonal V_ii = x_i^2 / 2
let potential = DMatrix::from_diagonal(&DVector::from_fn(n - 1, |i, _| {
  let x = -l + (i + 1) as f64 * h;
  x * x / 2.0
}));

// Hamiltonian and eigensolve
let hamiltonian = kinetic + potential;
let eigen = SymmetricEigen::new(hamiltonian);
let mut eigenvalues: Vec<f64> = eigen.eigenvalues.iter().copied().collect();
eigenvalues.sort_by(|a, b| a.partial_cmp(b).unwrap());

for (i, &e) in eigenvalues.iter().take(5).enumerate() {
  let exact = i as f64 + 0.5;
  println!("E_{i} = {e:.6}  (exact {exact:.4}, error {:.2e})", (e - exact).abs());
}

import numpy as np
from scipy.linalg import eigh_tridiagonal

N = 200
L = 8.0
h = 2 * L / N
x = np.linspace(-L + h, L - h, N - 1)  # interior vertices

# Tridiagonal Hamiltonian: T + V
diag = 1.0 / h**2 + x**2 / 2         # main diagonal
off  = -0.5 / h**2 * np.ones(N - 2)   # off-diagonals

eigenvalues, _ = eigh_tridiagonal(diag, off)
exact = np.arange(5) + 0.5

for i in range(5):
    err = abs(eigenvalues[i] - exact[i])
    print(f"E_{i} = {eigenvalues[i]:.6f}  (exact {exact[i]:.4f}, error {err:.2e})")

3. Differential forms, the wedge product, and the exterior derivative

The Whitney map (Section 4) turns cochains into differential forms, and it uses two operations we have not yet defined: the wedge product and the exterior derivative . These are the fundamental operations of the calculus of differential forms. This section constructs both from scratch, assuming only multivariable calculus and linear algebra.

3.1. Differential forms

A function assigns a number to each point. A differential form generalises this: it assigns a number not just to a point, but to a point together with a collection of tangent vectors at that point. The number of tangent vectors it requires determines the degree of the form.

In coordinates on , the coordinate differentials are the basic 1-forms. Each is the linear map that reads off the -th component of a tangent vector:

3.2. The wedge product

The wedge product is the operation that builds higher-degree forms from lower-degree ones. It is the exterior algebra analogue of multiplication.

The full permutation formula [eq:wedge] is useful for proofs, but in practice the wedge product is computed using the following rules.

3.3. The exterior derivative

The exterior derivative is the universal differentiation operator for differential forms. It generalises the gradient, curl, and divergence of vector calculus into a single operation.

These three properties determine uniquely. For a general -form (where the sum is over increasing multi-indices ), the exterior derivative is

The nilpotency gives us a cochain complex of smooth forms, the de Rham complex:

where denotes the space of smooth -forms on . This is the continuous analogue of the discrete cochain complex from Section 2.

3.4. Integration of forms and Stokes' theorem

Differential forms are the natural objects to integrate. A -form is integrated over a -dimensional oriented submanifold, producing a number. The general Stokes' theorem is the master identity that connects integration and differentiation of forms.

The discrete coboundary from Section 2 was defined by . Stokes' theorem shows this is the discrete version of : the coboundary is the discrete exterior derivative, and it satisfies the discrete Stokes' theorem exactly, with no approximation.

3.5. A note on integration notation

4. Whitney forms and discrete-to-continuous

With differential forms, the wedge product, and the exterior derivative in hand, we can now construct the maps between the continuous and discrete worlds. Cochains assign numbers to simplices. Differential forms assign numbers to tangent vectors at every point. The de Rham map turns a continuous form into a cochain by integrating over simplices. The Whitney map goes the other direction: it turns a cochain into a continuous differential form by interpolation. Together, they form a bridge between the continuous and discrete worlds, and the error of this bridge controls the accuracy of every DEC computation.

4.1. Barycentric coordinates

The Whitney map requires a coordinate system that is intrinsic to each simplex. Barycentric coordinates provide exactly this: they express any point inside a simplex as a weighted average of its vertices, with the weights summing to one.

4.2. The Whitney map

The Whitney map constructs a continuous differential form from a cochain by interpolating with barycentric coordinates. The construction uses the wedge product and exterior derivative from Section 3: the resulting form is piecewise polynomial, smooth on the interior of each simplex, and continuous across shared faces.

This formula is dense, so let us unpack it for each degree.

4.3. The de Rham map

The Whitney map goes from discrete to continuous. The de Rham map goes the other direction: it converts a continuous differential form into a cochain by integrating over each simplex.

4.4. The round trip: R, W, and the de Rham theorem

The Whitney and de Rham maps are approximate inverses of each other. The following results make this precise.

4.5. Interpolation error

The other direction of the round trip, , is not the identity on smooth forms. It is an approximation, and the error is controlled by the mesh size.

5. The discrete Hodge star

Everything we have built so far, simplicial complexes, cochains, the coboundary operator, Whitney forms, is purely topological. None of it depends on distances, angles, or volumes. The coboundary matrix is the same whether the simplices are equilateral or wildly deformed. The identity is combinatorial, not geometric.

But to solve PDEs, we must measure things. The Poisson equation [eq:poisson-intro] involves the Laplacian, which depends on the metric (distances and angles). The energy of a heat distribution involves the norm, which requires integration against a volume form. The Hodge star is where the metric enters the DEC framework.

5.1. The dual mesh

Every simplicial complex has a natural dual decomposition of the same domain, constructed from the circumcentres of its simplices. Where the primal mesh decomposes the domain into simplices (triangles, tetrahedra), the dual mesh decomposes it into cells that are roughly perpendicular to the primal elements. The ratio of dual cell volume to primal simplex volume is the Hodge star.

5.2. Well-centredness

The dual mesh construction assumes that circumcentres lie inside their simplices. When this condition fails (for obtuse triangles, for example), the dual cell volumes can become negative, and the Hodge star breaks down. The following condition guarantees that this does not happen.

5.3. The discrete Hodge star

With the dual mesh in place, the Hodge star has a remarkably simple form: it is a diagonal matrix whose entries are ratios of dual cell volumes to primal simplex volumes.

5.4. The cochain inner product

The Hodge star turns the space of cochains into an inner product space, which is essential for defining norms, energies, and adjoint operators. The inner product of two -cochains weights each simplex's contribution by the corresponding Hodge star entry.

5.5. Convergence of the Hodge star

5.6. Example: the 2D simple harmonic oscillator

We extend the 1D SHO computation from Section 2 to two dimensions. The 2D quantum harmonic oscillator has the Schrodinger equation

with exact eigenvalues and degeneracies: (non-degenerate), (doubly degenerate), (triply degenerate), and so on.

Discretisation. Triangulate the square with a Delaunay mesh. The Hodge Laplacian on 0-forms (which we will define formally in Section 6, but preview here) is

The Hamiltonian is where is the diagonal matrix with . The eigenvalue problem is the generalised system .

Results. With and a Delaunay mesh with approximately 2500 interior vertices, the first six eigenvalues are:

The degeneracies are resolved correctly: the two eigenvalues near agree to , and the three near agree to . On a log-log plot of eigenvalue error versus , the slope is approximately 2, confirming the convergence rate predicted by the theory.

use cartan::simplicial::SimplicialComplex2D;
use cartan::dec::{CoboundaryOperator, HodgeStar};
use cartan::mesh::delaunay_triangulate_square;
use nalgebra_sparse::CsrMatrix;

// Delaunay triangulation of [-6, 6]^2 with ~2500 interior vertices
let mesh = delaunay_triangulate_square(-6.0, 6.0, 50);

// DEC operators
let d0 = mesh.coboundary_matrix(0);   // edges x vertices
let star0 = mesh.hodge_star(0);       // diagonal, vertices
let star1 = mesh.hodge_star(1);       // diagonal, edges

// Hodge Laplacian: L = star1_inv * d0^T * star1 * d0
let star1_d0 = &star1 * &d0;
let laplacian = &d0.transpose() * &star1_d0;

// Potential: V_ii = (x_i^2 + y_i^2) / 2
let potential = mesh.vertex_diagonal(|v| (v.x * v.x + v.y * v.y) / 2.0);

// Generalised eigenproblem: (-0.5 * L + star0 * V) psi = E * star0 * psi
let (h_int, star0_int) = mesh.restrict_to_interior(
  &(-0.5 * &laplacian + &(&star0 * &potential)),
  &star0,
);

let eigenvalues = sparse_generalised_eigsh(&h_int, &star0_int, 6);
for (i, &e) in eigenvalues.iter().enumerate() {
  println!("E_{i} = {e:.4}");
}

import numpy as np
from scipy.spatial import Delaunay
from scipy.sparse import lil_matrix, diags, csr_matrix
from scipy.sparse.linalg import eigsh

# --- Grid and Delaunay triangulation ---
n_side = 40
L = 6.0
x = np.linspace(-L, L, n_side)
xx, yy = np.meshgrid(x, x)
points = np.column_stack([xx.ravel(), yy.ravel()])
tri = Delaunay(points)
n_verts = len(points)

# --- Build DEC operators from triangle geometry ---
# Collect unique oriented edges
edge_map = {}
for simplex in tri.simplices:
    for i in range(3):
        e = tuple(sorted((simplex[i], simplex[(i + 1) % 3])))
        if e not in edge_map:
            edge_map[e] = len(edge_map)
n_edges = len(edge_map)

# d0: coboundary operator (n_edges x n_verts)
d0 = lil_matrix((n_edges, n_verts))
for (v0, v1), idx in edge_map.items():
    d0[idx, v0] = -1
    d0[idx, v1] = 1
d0 = csr_matrix(d0)

# Hodge stars from circumcentric duals
star0 = np.zeros(n_verts)  # Voronoi dual areas
star1 = np.zeros(n_edges)  # dual-edge / primal-edge length ratios
for simplex in tri.simplices:
    i, j, k = simplex
    pi, pj, pk = points[i], points[j], points[k]
    area = 0.5 * abs(np.cross(pj - pi, pk - pi))
    # Circumcentre
    D = 2 * (pi[0] * (pj[1] - pk[1]) + pj[0] * (pk[1] - pi[1])
             + pk[0] * (pi[1] - pj[1]))
    if abs(D) < 1e-14:
        cc = (pi + pj + pk) / 3
    else:
        ux = ((pi @ pi) * (pj[1] - pk[1]) + (pj @ pj) * (pk[1] - pi[1])
              + (pk @ pk) * (pi[1] - pj[1])) / D
        uy = ((pi @ pi) * (pk[0] - pj[0]) + (pj @ pj) * (pi[0] - pk[0])
              + (pk @ pk) * (pj[0] - pi[0])) / D
        cc = np.array([ux, uy])
    for v in [i, j, k]:
        star0[v] += area / 3
    for (v0, v1) in [(i,j), (j,k), (i,k)]:
        ek = tuple(sorted((v0, v1)))
        eidx = edge_map[ek]
        mid = 0.5 * (points[v0] + points[v1])
        star1[eidx] += np.linalg.norm(cc - mid) / np.linalg.norm(points[v1] - points[v0])

# --- Assemble Hamiltonian ---
S1 = diags(star1)
Lap = d0.T @ S1 @ d0                       # stiffness matrix
V = 0.5 * (points[:, 0]**2 + points[:, 1]**2)
S0 = diags(star0)
H = -0.5 * Lap + S0 @ diags(V)

# --- Restrict to interior vertices ---
edge_count = {}
for simplex in tri.simplices:
    for i in range(3):
        e = tuple(sorted((simplex[i], simplex[(i + 1) % 3])))
        edge_count[e] = edge_count.get(e, 0) + 1
bnd = sorted({v for (v0, v1), c in edge_count.items() if c == 1 for v in (v0, v1)})
int_mask = np.ones(n_verts, dtype=bool)
int_mask[bnd] = False
int_idx = np.where(int_mask)[0]
H_int = csr_matrix(H.toarray()[np.ix_(int_idx, int_idx)])
S0_int = diags(star0[int_idx])

# --- Solve generalised eigenproblem ---
eigenvalues, _ = eigsh(H_int, k=6, M=S0_int, sigma=0, which='LM')
eigenvalues.sort()
exact = sorted(i + j + 1.0 for i in range(4) for j in range(4))[:6]
for i in range(6):
    err = abs(eigenvalues[i] - exact[i])
    print(f"E_{i} = {eigenvalues[i]:.4f}  (exact {exact[i]:.1f}, error {err:.2e})")

6. The discrete Laplacian and Poisson's equation

We now assemble the operators from the preceding sections into a discrete Laplacian that solves the problem posed in Section 0: a geometrically faithful discretisation of the Poisson equation [eq:poisson-intro] that works on arbitrary domains, including the L-shaped domain where the five-point stencil failed.

6.1. The Hodge Laplacian

In the smooth setting, the Laplacian on functions (0-forms) is : apply the exterior derivative, then the Hodge star, then the exterior derivative again, then the Hodge star once more. The general Hodge Laplacian on -forms also includes a term from the codifferential , but on 0-forms this term vanishes, so we focus on the 0-form case.

Each factor has a clear geometric role. The coboundary computes differences along edges (the discrete gradient). The Hodge star converts these edge differences into dual-edge fluxes, weighting each edge by the ratio of dual-edge length to primal-edge length. The transpose sums the incoming fluxes at each vertex (the discrete divergence). The inverse Hodge star normalises by the Voronoi cell area at each vertex. The composition is the discrete divergence of the discrete gradient: the discrete Laplacian.

6.2. The discrete Poisson equation

The discrete Poisson equation is: given a 0-cochain (the source), find a 0-cochain such that

subject to boundary conditions. For Dirichlet boundary conditions ( prescribed on boundary vertices), we partition the vertices into interior and boundary sets, restrict [eq:discrete-poisson] to interior vertices, and move the known boundary values to the right-hand side. The resulting system is sparse, symmetric, and positive definite (on the interior, by the well-centredness of the mesh), and can be solved by a sparse Cholesky factorisation or a conjugate gradient iteration.

6.3. The L-shaped domain revisited

Return to the L-shaped domain from Section 0. The exact solution to the Laplace equation has a singularity at the re-entrant corner: .

The five-point stencil on a uniform grid converges at in the energy norm because it cannot adapt to the singularity. The DEC approach, using a triangulation that can be refined near the corner, does significantly better. On a graded mesh with element sizes proportional to the distance from the corner (so that the smallest elements are concentrated where the gradient is largest), the DEC solution achieves the optimal rate in and in , matching the smooth-domain rates. Even the uniform triangular mesh, thanks to the Aubin-Nitsche duality argument, converges at in , already a substantial improvement over the stencil.

The key is that the simplicial complex conforms exactly to the re-entrant corner (no staircase approximation), and the graded mesh resolves the singularity without wasting degrees of freedom in the smooth interior. The Hodge Laplacian on this mesh inherits all the structural properties (symmetry, positive definiteness, exactness of the coboundary) that guarantee convergence, regardless of the domain geometry.

use cartan::simplicial::SimplicialComplex2D;
use cartan::dec::{CoboundaryOperator, HodgeStar};
use cartan::mesh::delaunay_triangulate_polygon;
use nalgebra_sparse::CsrMatrix;

// L-shaped domain vertices (counter-clockwise)
let boundary = vec![
  (-1.0, -1.0), (1.0, -1.0), (1.0, 0.0),
  (0.0, 0.0), (0.0, 1.0), (-1.0, 1.0),
];

// Graded Delaunay triangulation: refine near the re-entrant corner (0,0)
let mesh = delaunay_triangulate_polygon(&boundary, |x, y| {
  let r = (x * x + y * y).sqrt().max(0.01);
  0.02 + 0.15 * r  // element size grows with distance from corner
});

// DEC operators
let d0 = mesh.coboundary_matrix(0);
let star0 = mesh.hodge_star(0);
let star1 = mesh.hodge_star(1);

// Hodge Laplacian: L = star0_inv * d0^T * star1 * d0
let laplacian = mesh.hodge_laplacian_0();

// RHS: f = 0 (Laplace equation)
// Boundary: u(r, theta) = r^(2/3) sin(2*theta/3)
let (l_int, rhs) = mesh.apply_dirichlet(&laplacian, |x, y| {
  let r = (x * x + y * y).sqrt();
  let theta = y.atan2(x);
  let theta = if theta < 0.0 { theta + 2.0 * std::f64::consts::PI } else { theta };
  r.powf(2.0 / 3.0) * (2.0 * theta / 3.0).sin()
});

// Solve: sparse Cholesky
let u_h = sparse_cholesky_solve(&l_int, &rhs);

// Compute L2 error against exact solution
let l2_err = mesh.l2_error(&u_h, |x, y| {
  let r = (x * x + y * y).sqrt();
  let theta = y.atan2(x);
  let theta = if theta < 0.0 { theta + 2.0 * std::f64::consts::PI } else { theta };
  r.powf(2.0 / 3.0) * (2.0 * theta / 3.0).sin()
});
println!("L2 error: {l2_err:.4e}");

import numpy as np
from scipy.spatial import Delaunay
from scipy.sparse import lil_matrix, diags, csr_matrix
from scipy.sparse.linalg import spsolve

def build_dec_2d(points, simplices):
    """Build d0, star0, star1 from a 2D triangulation."""
    n_v = len(points)
    edge_map = {}
    for s in simplices:
        for i in range(3):
            e = tuple(sorted((s[i], s[(i+1)%3])))
            if e not in edge_map:
                edge_map[e] = len(edge_map)
    n_e = len(edge_map)
    d0 = lil_matrix((n_e, n_v))
    for (v0, v1), idx in edge_map.items():
        d0[idx, v0] = -1; d0[idx, v1] = 1
    d0 = csr_matrix(d0)
    star0 = np.zeros(n_v); star1 = np.zeros(n_e)
    for s in simplices:
        i, j, k = s; pi, pj, pk = points[i], points[j], points[k]
        area = 0.5 * abs(np.cross(pj - pi, pk - pi))
        D = 2*(pi[0]*(pj[1]-pk[1]) + pj[0]*(pk[1]-pi[1]) + pk[0]*(pi[1]-pj[1]))
        cc = (pi+pj+pk)/3 if abs(D)<1e-14 else np.array([
            ((pi@pi)*(pj[1]-pk[1])+(pj@pj)*(pk[1]-pi[1])+(pk@pk)*(pi[1]-pj[1]))/D,
            ((pi@pi)*(pk[0]-pj[0])+(pj@pj)*(pi[0]-pk[0])+(pk@pk)*(pj[0]-pi[0]))/D])
        for v in [i,j,k]: star0[v] += area/3
        for (v0,v1) in [(i,j),(j,k),(i,k)]:
            ek = tuple(sorted((v0,v1))); eidx = edge_map[ek]
            mid = 0.5*(points[v0]+points[v1])
            star1[eidx] += np.linalg.norm(cc-mid)/np.linalg.norm(points[v1]-points[v0])
    return d0, star0, star1

def l_shaped_exact(x, y):
    r = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y, x)
    theta = np.where(theta < 0, theta + 2*np.pi, theta)
    return r**(2/3) * np.sin(2*theta/3)

def l_shaped_points(n, graded=False):
    """Points in [-1,1]^2 \\ (0,1]^2."""
    if graded:
        t = np.linspace(0, 1, n)
        v = np.unique(np.concatenate([-1 + t, t]))
    else:
        v = np.linspace(-1, 1, n)
    xx, yy = np.meshgrid(v, v)
    pts = np.column_stack([xx.ravel(), yy.ravel()])
    return pts[~((pts[:,0] > 1e-10) & (pts[:,1] > 1e-10))]

def solve_l(n, graded=False):
    pts = l_shaped_points(n, graded)
    tri = Delaunay(pts)
    d0, star0, star1 = build_dec_2d(pts, tri.simplices)
    n_v = len(pts)
    L = d0.T @ diags(star1) @ d0
    # Boundary vertices
    ec = {}
    for s in tri.simplices:
        for i in range(3):
            e = tuple(sorted((s[i], s[(i+1)%3]))); ec[e] = ec.get(e,0)+1
    bnd = sorted({v for (v0,v1),c in ec.items() if c==1 for v in (v0,v1)})
    mask = np.ones(n_v, bool); mask[bnd] = False; ii = np.where(mask)[0]
    u_ex = l_shaped_exact(pts[:,0], pts[:,1])
    Ld = L.toarray()
    rhs = -Ld[np.ix_(ii, bnd)] @ u_ex[bnd]
    u_int = spsolve(csr_matrix(Ld[np.ix_(ii,ii)]), rhs)
    u_h = u_ex.copy(); u_h[ii] = u_int
    return 2.0/n, np.sqrt(np.sum(star0 * (u_h - u_ex)**2))

for label, graded in [("Uniform", False), ("Graded", True)]:
    print(f"\n{label} mesh:")
    for n in [15, 20, 30, 40, 60]:
        h, err = solve_l(n, graded)
        print(f"  h = {h:.4f}, L2 error = {err:.4e}")

7. Parabolic problems: the heat equation

The Poisson equation is elliptic: it describes a steady state with no time dependence. The heat equation introduces time. Rather than asking "what is the temperature after all transients have died?", we ask "how does the temperature evolve from a given initial distribution?" The DEC discretisation of the heat equation is a direct consequence of the Hodge Laplacian.

7.1. The continuous heat equation

The heat equation on a domain is

with initial condition and boundary conditions (Dirichlet, Neumann, or mixed) on . The solution describes how an initial temperature distribution diffuses over time. Heat flows from hot to cold at a rate proportional to the Laplacian, which measures the local average deviation of from its surroundings.

7.2. The semidiscrete system

Replacing the spatial Laplacian with the Hodge Laplacian [eq:hodge-laplacian] and leaving time continuous gives the semidiscrete heat equation: a system of ODEs for the cochain :

The mass matrix on the left accounts for the fact that different vertices control different amounts of the domain (their Voronoi cells). On a uniform mesh in 1D, and the system reduces to the familiar method of lines for the heat equation.

7.3. Time integration

To integrate [eq:semidiscrete-heat] in time, we discretise the time derivative. Let be the time step and .

Forward Euler (explicit). The simplest choice is

which gives . This is cheap (one matrix-vector multiply per step) but conditionally stable: the time step must satisfy for a mesh-dependent constant , or the solution blows up.

Backward Euler (implicit). Evaluating the right-hand side at time level gives

This requires solving a sparse linear system at each step, but is unconditionally stable: any time step produces a bounded solution. The matrix on the left is symmetric positive definite (by well-centredness), so a single sparse Cholesky factorisation, reused at every time step, makes the cost per step comparable to the explicit method.

Crank-Nicolson (trapezoidal). Averaging the forward and backward Euler gives

This is unconditionally stable and second-order accurate in time (error ), compared to first-order () for forward and backward Euler.

7.4. Example: heat diffusion on a disk

Consider the heat equation on the unit disk with homogeneous Dirichlet boundary conditions ( on ) and initial condition . The exact solution is a sum of Bessel functions, and the temperature decays exponentially to zero as the heat escapes through the boundary.

On a Delaunay triangulation of the disk with approximately 1000 interior vertices, backward Euler with produces the following decay of the maximum temperature:

The solution preserves positivity (the discrete maximum principle) and monotonically decreases, consistent with the continuous problem. The error is dominated by the spatial discretisation () at early times and by the temporal discretisation () at late times.

use cartan::simplicial::SimplicialComplex2D;
use cartan::dec::{CoboundaryOperator, HodgeStar};
use cartan::mesh::delaunay_triangulate_disk;

// Delaunay triangulation of the unit disk, ~1000 interior vertices
let mesh = delaunay_triangulate_disk(1.0, 0.05);
let d0 = mesh.coboundary_matrix(0);
let star0 = mesh.hodge_star(0);
let star1 = mesh.hodge_star(1);

// Stiffness matrix: A = d0^T * star1 * d0
let stiffness = &d0.transpose() * &(&star1 * &d0);

// Backward Euler: (star0 + dt * A) u^{n+1} = star0 * u^n
let dt = 0.001;
let lhs = &star0 + &(&stiffness * dt);
let chol = sparse_cholesky(&lhs);  // factor once, reuse

// Initial condition: u0 = 1 - x^2 - y^2, restricted to interior
let mut u = mesh.interior_cochain(|v| 1.0 - v.x * v.x - v.y * v.y);

for step in 0..500 {
  let rhs = &star0 * &u;
  u = chol.solve(&rhs);  // reuse factorisation
  if step % 50 == 49 {
      let t = (step + 1) as f64 * dt;
      println!("t = {t:.2}, max u = {:.4}", u.max());
  }
}

import numpy as np
from scipy.spatial import Delaunay
from scipy.sparse import lil_matrix, diags, csr_matrix
from scipy.sparse.linalg import splu

# --- Mesh: concentric rings in the unit disk ---
pts = [[0.0, 0.0]]
n_rad, n_ang = 20, 24
for i in range(1, n_rad + 1):
    r = i / n_rad
    nt = max(6, int(n_ang * r))
    for j in range(nt):
        theta = 2 * np.pi * j / nt
        pts.append([r * np.cos(theta), r * np.sin(theta)])
points = np.array(pts)
tri = Delaunay(points)
n_v = len(points)

# --- DEC assembly (same helper as previous cells) ---
edge_map = {}
for s in tri.simplices:
    for i in range(3):
        e = tuple(sorted((s[i], s[(i+1)%3])))
        if e not in edge_map: edge_map[e] = len(edge_map)
n_e = len(edge_map)
d0 = lil_matrix((n_e, n_v))
for (v0, v1), idx in edge_map.items():
    d0[idx, v0] = -1; d0[idx, v1] = 1
d0 = csr_matrix(d0)
star0 = np.zeros(n_v); star1 = np.zeros(n_e)
for s in tri.simplices:
    i, j, k = s; pi, pj, pk = points[i], points[j], points[k]
    area = 0.5 * abs(np.cross(pj - pi, pk - pi))
    D = 2*(pi[0]*(pj[1]-pk[1]) + pj[0]*(pk[1]-pi[1]) + pk[0]*(pi[1]-pj[1]))
    cc = (pi+pj+pk)/3 if abs(D)<1e-14 else np.array([
        ((pi@pi)*(pj[1]-pk[1])+(pj@pj)*(pk[1]-pi[1])+(pk@pk)*(pi[1]-pj[1]))/D,
        ((pi@pi)*(pk[0]-pj[0])+(pj@pj)*(pi[0]-pk[0])+(pk@pk)*(pj[0]-pi[0]))/D])
    for v in [i,j,k]: star0[v] += area/3
    for (v0,v1) in [(i,j),(j,k),(i,k)]:
        ek = tuple(sorted((v0,v1))); eidx = edge_map[ek]
        mid = 0.5*(points[v0]+points[v1])
        star1[eidx] += np.linalg.norm(cc-mid)/np.linalg.norm(points[v1]-points[v0])

# --- Stiffness and boundary restriction ---
A = d0.T @ diags(star1) @ d0
ec = {}
for s in tri.simplices:
    for i in range(3):
        e = tuple(sorted((s[i], s[(i+1)%3]))); ec[e] = ec.get(e,0)+1
bnd = sorted({v for (v0,v1),c in ec.items() if c==1 for v in (v0,v1)})
mask = np.ones(n_v, bool); mask[bnd] = False
int_idx = np.where(mask)[0]
S0_int = diags(star0[int_idx])
A_int = csr_matrix(A.toarray()[np.ix_(int_idx, int_idx)])

# --- Backward Euler time stepping ---
dt = 0.001
lhs = S0_int + dt * A_int
lu = splu(lhs.tocsc())
u = 1.0 - points[int_idx, 0]**2 - points[int_idx, 1]**2

for step in range(500):
    u = lu.solve(S0_int @ u)
    if (step + 1) % 50 == 0:
        print(f"t = {(step+1)*dt:.3f}, max u = {u.max():.4f}")

8. Maxwell's equations and computational electromagnetics

The Poisson and heat equations involve only 0-forms (scalar fields). Maxwell's equations demand the full de Rham complex: the electric field is a 1-form, the magnetic field is a 2-form, charge density is a 3-form, and the equations relating them are the exterior derivative and the Hodge star. DEC discretises Maxwell's equations with the same operators we have already built, and the identity automatically enforces the constraint equations (Gauss's laws) at the discrete level.

8.1. Maxwell's equations in differential forms

On a 3-dimensional domain, Maxwell's equations in the language of differential forms are:

where is the electromagnetic 2-form (the Faraday tensor), is its Hodge dual, and is the current 3-form. In components on (the static case, separating electric and magnetic fields), these become the familiar four equations:

The equation (no magnetic monopoles) and (Faraday) are topological: they involve only the exterior derivative and hold on any manifold, with any metric. The equations involving are metric-dependent: the Hodge star encodes the constitutive relations of the medium (permittivity, permeability).

8.2. DEC discretisation

The DEC discretisation assigns each field to the appropriate cochain degree:

Faraday's law becomes , which is exact (the coboundary is the transpose of the boundary, no approximation). Gauss's law becomes , which is automatically satisfied by if for some 1-cochain (the vector potential). The constitutive relations (, ) are encoded in the Hodge star, which now carries the material parameters.

8.3. Example: electrostatics of a biological cell

Consider the electrostatic potential inside a spherical biological cell of radius , modelled as a dielectric sphere with permittivity surrounded by a medium with permittivity . A point charge at position inside the cell creates a potential that satisfies the Poisson equation , with continuity conditions at the membrane.

In DEC, the variable permittivity enters through the Hodge star: is multiplied by evaluated at the midpoint of each edge (or, more precisely, by the average of over the dual cell of the edge). The discrete Poisson equation becomes

where is the material Hodge star and is the discrete charge distribution (a 0-cochain that is nonzero only at the vertex closest to ).

On a tetrahedral mesh of the cell with approximately 5000 interior vertices and a refined layer near the membrane (where the permittivity jumps), the DEC solution captures the image charge effects that arise from the dielectric mismatch. The potential inside the cell is enhanced relative to the free-space Coulomb potential, because the lower permittivity of the membrane partially reflects the electric field lines back into the cell.

9. Elliptic eigenvalue problems: atomic orbitals

The SHO examples in Sections 2 and 5 solved the Schrodinger equation in 1D and 2D with a harmonic potential. This section extends to 3D with a Coulomb potential, computing the hydrogen atom's energy levels and orbital shapes. The Hodge Laplacian discretises the kinetic energy; the potential is a diagonal 0-cochain; and the eigenvalue problem is a generalised sparse eigenproblem whose solutions are the atomic orbitals.

9.1. The hydrogen atom

The time-independent Schrodinger equation for the hydrogen atom (in atomic units, ) is

on , with as . The exact eigenvalues are for , with degeneracy . The eigenfunctions (orbitals) are products of radial functions and spherical harmonics: the familiar , , , , etc.

9.2. DEC discretisation in 3D

We work on a tetrahedral mesh of a large ball with chosen large enough that the wavefunctions have decayed to negligible values at the boundary (typically Bohr radii for the first few orbitals). The mesh is refined near the origin, where the Coulomb singularity demands higher resolution.

The discrete Hamiltonian is

where is the 3D Hodge Laplacian on 0-forms [eq:hodge-laplacian] (now built from a tetrahedral mesh, with the vertex-to-edge coboundary, encoding edge-to-dual-face volume ratios, and encoding Voronoi cell volumes). The eigenvalue problem is

a generalised sparse eigenproblem. We seek the smallest (most negative) eigenvalues.

9.3. Results: energy levels and degeneracies

On a graded tetrahedral mesh with approximately 15,000 interior vertices (element size near the origin, near ), the first eigenvalues are:

The degeneracies are resolved correctly: the 4 eigenvalues in the shell (one and three states) agree to , and the 9 eigenvalues in the shell agree to . The error grows with because higher orbitals extend further and are less well-resolved by the finite mesh.

use cartan::simplicial::SimplicialComplex3D;
use cartan::dec::{CoboundaryOperator, HodgeStar};
use cartan::mesh::delaunay_triangulate_ball;

// Graded tetrahedral mesh of B_20, refined near origin
let mesh = delaunay_triangulate_ball(20.0, |r| 0.1 + 0.1 * r);

let d0 = mesh.coboundary_matrix(0);
let star0 = mesh.hodge_star(0);
let star1 = mesh.hodge_star(1);

// Hodge Laplacian on 0-forms
let laplacian = mesh.hodge_laplacian_0();

// Coulomb potential: V_vv = -1/r
let potential = mesh.vertex_diagonal(|v| {
  let r = (v.x * v.x + v.y * v.y + v.z * v.z).sqrt().max(1e-10);
  -1.0 / r
});

// Hamiltonian: H = -0.5 * L + star0 * V
let hamiltonian = -0.5 * &laplacian + &(&star0 * &potential);

// Restrict to interior, solve generalised eigenproblem
let (h_int, s_int) = mesh.restrict_to_interior(&hamiltonian, &star0);
let eigenvalues = sparse_generalised_eigsh(&h_int, &s_int, 14);

for (i, &e) in eigenvalues.iter().enumerate() {
  println!("E_{i} = {e:.6}");
}

import numpy as np
from scipy.spatial import Delaunay
from scipy.sparse import lil_matrix, diags, csr_matrix
from scipy.sparse.linalg import eigsh

# --- Graded tetrahedral mesh of a ball of radius R ---
R = 15.0
pts = [[0, 0, 0]]
# Concentric shells with grading: finer near origin
for i in range(1, 16):
    r = R * (i / 15)**1.5          # grade toward origin
    n_phi = max(4, int(8 * (i / 15)))
    for ip in range(n_phi):
        phi = np.pi * (ip + 0.5) / n_phi
        n_theta = max(6, int(12 * np.sin(phi) * (i / 15)))
        for it in range(n_theta):
            theta = 2 * np.pi * it / n_theta
            pts.append([r*np.sin(phi)*np.cos(theta),
                        r*np.sin(phi)*np.sin(theta), r*np.cos(phi)])
points = np.array(pts)
tri = Delaunay(points)
n_v = len(points)

# --- DEC assembly (3D: tetrahedra) ---
edge_map = {}
for s in tri.simplices:
    for i in range(4):
        for j in range(i+1, 4):
            e = tuple(sorted((s[i], s[j])))
            if e not in edge_map: edge_map[e] = len(edge_map)
n_e = len(edge_map)
d0 = lil_matrix((n_e, n_v))
for (v0, v1), idx in edge_map.items():
    d0[idx, v0] = -1; d0[idx, v1] = 1
d0 = csr_matrix(d0)

# Hodge stars: lumped (barycentric) for simplicity
star0 = np.zeros(n_v)
star1 = np.zeros(n_e)
for s in tri.simplices:
    verts = points[s]
    # Tetrahedron volume
    mat = np.column_stack([verts[1]-verts[0], verts[2]-verts[0], verts[3]-verts[0]])
    vol = abs(np.linalg.det(mat)) / 6
    for v in s:
        star0[v] += vol / 4
    for i in range(4):
        for j in range(i+1, 4):
            ek = tuple(sorted((s[i], s[j])))
            eidx = edge_map[ek]
            primal_len = np.linalg.norm(points[s[j]] - points[s[i]])
            # Barycentric dual: vol / (6 * edge_length)
            star1[eidx] += vol / (6 * primal_len)

# --- Hamiltonian ---
Lap = d0.T @ diags(star1) @ d0
V = np.array([-1.0 / max(np.linalg.norm(p), 1e-8) for p in points])
S0 = diags(star0)
H = -0.5 * Lap + S0 @ diags(V)

# --- Restrict to interior ---
ec = {}
for s in tri.simplices:
    for i in range(4):
        for j in range(i+1, 4):
            for k in range(j+1, 4):
                f = tuple(sorted((s[i], s[j], s[k])))
                ec[f] = ec.get(f, 0) + 1
bnd = sorted({v for face, c in ec.items() if c == 1 for v in face})
mask = np.ones(n_v, bool); mask[bnd] = False
int_idx = np.where(mask)[0]
H_int = csr_matrix(H.toarray()[np.ix_(int_idx, int_idx)])
S0_int = diags(star0[int_idx])

# --- Solve generalised eigenproblem ---
eigenvalues, _ = eigsh(H_int, k=10, M=S0_int, sigma=-0.6, which='LM')
eigenvalues.sort()
exact = [-0.5 / n**2 for n in range(1, 4)]  # -0.5, -0.125, -0.0556
print("Hydrogen eigenvalues (Hartree):")
for i, e in enumerate(eigenvalues[:10]):
    print(f"  E_{i} = {e:.6f}")
print(f"\nExact: E_1s = {exact[0]:.4f}, E_2 = {exact[1]:.4f}, E_3 = {exact[2]:.4f}")

The visualisation below shows the exact hydrogen wavefunctions on the -plane. Select quantum numbers to explore the orbital geometries: the number of radial nodes grows with , the angular structure is governed by the spherical harmonic , and the sign changes (blue to red) reflect the nodal surfaces that the DEC eigensolver must resolve.

10. Nonlinear problems: the Navier-Stokes equations

Every problem so far has been linear: the Poisson equation, the heat equation, Maxwell's equations, the Schrodinger equation. The Navier-Stokes equations for viscous incompressible flow are nonlinear, because the fluid velocity appears in both the transport and the transported quantity. DEC handles the nonlinearity by discretising the vorticity-stream function formulation, where the topological structure (the de Rham complex) enforces the incompressibility constraint exactly.

10.1. The vorticity-stream function formulation

The incompressible Navier-Stokes equations in 2D are

where is the velocity, is the pressure, and is the kinematic viscosity. The divergence-free constraint is a topological condition: it says the velocity 1-form is coclosed.

In the vorticity-stream function formulation, we define the vorticity (a scalar in 2D, the curl of the velocity) and the stream function (satisfying ). The incompressibility constraint is automatically satisfied, and the Navier-Stokes equations reduce to a single scalar equation:

This is an advection-diffusion equation for the vorticity, with the velocity field recovered from the vorticity by solving (a Poisson equation) and computing .

10.2. DEC discretisation

In the DEC framework, the vorticity is a 0-cochain (a scalar at each vertex), and the stream function is also a 0-cochain. The velocity is a 1-cochain (circulations along edges), recovered from the stream function by the coboundary: , where the Hodge star rotates the gradient 1-form by to give the perpendicular velocity.

The discrete vorticity equation is

where is the discrete advection operator, assembled from the edge circulations of the velocity. At each time step:

1. Solve for the stream function (a Poisson solve).
2. Compute edge velocities .
3. Assemble the advection operator from edge velocities.
4. Advance by one time step using [eq:dec-vorticity].

The figure below shows the result of this DEC solver on the standard lid-driven cavity benchmark at . The top wall moves with unit velocity to the right, driving a primary recirculation vortex that fills the cavity. The vorticity concentrates at the top corners (where the moving lid meets the stationary walls) and diffuses into the interior, exactly as the continuous solution predicts.

10.3. Energy estimates

The cochain inner product provides the natural energy functional.

The figure below verifies the energy dissipation theorem numerically. Starting from a Taylor-Green-like initial vorticity (sinusoidal in both directions, zero on the boundary), the DEC solver evolves the vorticity under pure diffusion (no advection forcing). The normalised enstrophy is plotted for two viscosities. In both cases the enstrophy is strictly monotone decreasing, as the theorem guarantees. Higher viscosity dissipates faster, exactly as the continuous estimate predicts.

10.4. Time integration and solver comparison

The DEC framework separates the spatial discretisation (which is geometric and structure-preserving) from the time integration (which is a choice). The spatial operator is the same regardless of how we advance in time, but the choice of time integrator determines the stability, accuracy, and energy behaviour of the full scheme.

Consider the DEC heat equation on the unit disk with on the boundary. Three classical time integrators applied to this semi-discrete system produce strikingly different behaviour:

• Forward Euler (explicit, first-order): . Unconditionally unstable for the heat equation unless , where is the largest eigenvalue of . The solution blows up within a few time steps at the time step used here.

• Backward Euler (implicit, first-order): . Unconditionally stable (the implicit solve damps all modes), but introduces temporal diffusion that smears the solution.

• Classical RK4 (explicit, fourth-order): applies four explicit stages per step. Much more accurate than either Euler variant at the same , but still subject to a CFL stability constraint (satisfied here). The RK4 error is orders of magnitude below backward Euler.

The same spatial discretisation also allows a direct comparison with classical finite differences. On the unit square with the exact Poisson solution , both the 5-point stencil and the DEC cotangent Laplacian converge at in , but the DEC method achieves lower absolute error because the cotangent weights adapt to the Delaunay triangulation's geometry.

11. What was preserved

We began with a failure: the five-point stencil on an L-shaped domain, converging at instead of because the rectangular grid destroyed the geometry of the domain. We asked what geometric structure was lost and how to preserve it.

The answer, developed over the preceding ten sections, is a separation principle: topology is exact, geometry is approximate.

The topological layer consists of the simplicial complex (Section 1), the cochain complex (Section 2), the coboundary operator (Section 2), and the identity . These structures are combinatorial. They do not depend on the positions of the vertices, the shapes of the simplices, or the metric of the ambient space. They are preserved exactly in the DEC discretisation, with no approximation at any stage. The consequences are:

• Stokes' theorem holds exactly at the discrete level: .
• Gauss's laws are automatically satisfied: prevents magnetic monopoles from appearing and charge from being created or destroyed.
• Incompressibility of fluid flow is exact in the stream function formulation.
• Conservation laws that follow from (charge conservation, mass conservation, vorticity conservation) hold to machine precision.

The geometric layer consists of the Hodge star (Section 5), which encodes distances, angles, volumes, and material properties. This is where the approximation lives. The Hodge star converges at , and the Whitney interpolation error also contributes , giving overall in and in for the Poisson equation.

The classical numerical methods of the standard curriculum are special cases:

Each classical method preserves a fragment of the geometric structure by accident, because it was designed for a specific problem on a specific mesh. DEC preserves all of it by construction, because it discretises the entire de Rham complex, not just the particular operator appearing in the PDE at hand.

The price of this generality is small: a well-centred simplicial mesh (which any mesh generator can produce) and a diagonal Hodge star (which costs nothing to invert). The reward is a unified framework in which convergence theorems, energy estimates, and conservation laws hold for every problem simultaneously, from the Poisson equation on an L-shaped domain to the Navier-Stokes equations on a biological cell.