Finite elements for the highbrow mathematician

In recent years, study of the finite element method for numerically solving PDE has taken an increasingly homological turn. However, at least in my circles, “pure” mathematicians have largely not taken notice of this trend. I suspect that there’s an opportunity for a mutually beneficially partnership here — and so I shall here try to translate some of the ideas of the FEM to the language spoken by pure mathematicians.

The Dirichlet Laplacian, and the basic idea. Let $Pu = f$ be a partial differential equation, which for simplicity I will take to be a linear, elliptic equation of second order with homogeneous Dirichlet data on a compact polytope $\Omega$ . The canonical example is that P is the Laplacian, and then this system reads

$\displaystyle \begin{cases} -\Delta u = f \\ u|_{\partial \Omega} = 0 \end{cases}$

The key thing to observe is that since P is differential, it is local: we can determine $Pu(x)$ from the germ of u at x. Motivated by this observation, we divide and conquer by fixing a family $\mathcal T_h$ of partitions of $\Omega$ into simplices, where each simplex in $\mathcal T_h$ has edge lengths comparable to h. We then fix the space $V_h$ of PL shape functions, which are continuous functions on $\Omega$ whose restriction to each simplex is a linear polynomial. This is a finite-dimensional, hence closed, subspace of the Sobolev space $W^{1, 2}_0(\Omega)$ . Of course, other Sobolev spaces admit better choices of shape functions, which in some cases can be easily read off of the Sobolev embedding theorem.

We now take the orthogonal projection of the equation $Pu = f$ to $V_h$ to obtain a new equation $P_h u_h = f_h$ . If we can solve this equation, then by the orthogonality of the projection, the Bramble-Hilbert lemma for Sobolev spaces, and elliptic regularity, we will have

$\displaystyle \|u - u_h\|_{W^{1, 2}} \lesssim h \|u\|_{\dot W^{2, 2}} \lesssim h \|f\|_{L^2}$

thus for $h \ll 1$ we have a good solution in $W^{1, 2}$ . By modifying the choice of Sobolev spaces and shape functions we can obtain bounds at other regularities.

There are three obstructions to the ability to invert a matrix $P_h$ in a practical amount of time, however:

Unisolvence. The matrix must be square. This is no issue for the Laplacian, since $P_h: V_h \to V_h$ .
Conditioning. The size of the smallest eigenvalue must be bounded from below (so the matrix is not too close to a singular matrix). The smallest eigenvalue $\lambda_{1, h}$ of $P_h$ is controlled from below by the Dirichlet spectral gap of $\Omega$ , which is bounded from below as long as $\Omega$ is not too fractalline or degenerate:
$\displaystyle |\lambda_{1, h}| \gtrsim |\lambda_1| \gtrsim 1$ .
Sparsity. Most of the entries of $P_h$ should be zero, and the nonzero entries should be easy to pick out, so that we don’t have to iterate over the whole matrix when we invert it. This is ensured by locality: if $(i, j)$ is a pair of indices corresponding to a shape function on a simplex K, then the only nonzero entries in row i or column j are for shape functions on simplices which are adjacent to K.

Elliptic Maxwell, and dealing with systems. Unfortunately the Laplacian is positive-definite, not a system, and otherwise generally just about the nicest PDE there is. We want to solve other PDE, too, and they will probably be of greater interest.

So let us consider the Maxwell equations. They are hyperbolic, however, and this adds a lot of complications outside the scope of this blog post. So we assume time-independence, and in that case the Maxwell equations decouple into two elliptic systems: the electric and magnetic systems. The electric system can be rewritten as the Laplacian, so we focus on the magnetic system

$\displaystyle \begin{cases} \nabla \cdot B = 0 \\ \nabla \times B = J \\ B|_{\partial \Omega} \cdot \vec n_\Omega = 0 \end{cases}$

where B is the magnetic field, $\vec n_\Omega$ is the normal field to $\partial \Omega$ , and the electric current J is given. The right-hand side of the equation $\nabla \cdot B = 0$ is not data, so that zero is not secretly an epsilon: it asserts that B has no monopoles. Therefore it is natural to place B in a space of divergence-free vector fields. Doing so will impose that $\nabla \cdot B = 0$ as opposed to simply $\|\nabla \cdot B\|_{L^2} \ll 1$ or a similar “approximate monopole inequality”. I will refer to such a B as divergence conforming; a similar condition appears when one tries to solve the Stokes or Navier-Stokes equations for a liquid.

This can be done directly, but we’re going to have to go homological eventually anyways, so let’s rewrite the magnetic Maxwell equations as

$\displaystyle \begin{cases} dB = 0 \\ d * B = * J \\ \iota_{\partial \Omega}^* B = 0 \end{cases}$

where B is now viewed as a 2-form, J is now viewed as a 1-form, and $\iota: \partial \Omega \to \Omega$ is the inclusion map. This interpretation is quite natural, because it views B as something that we want to integrate over a surface $\Sigma$ , so that $\int_\Sigma B$ is now the magnetic flux.

Let $\Omega^\ell$ be the space of $\ell$ -forms, $K^\ell$ the kernel of $d: \Omega^\ell \to \Omega^{\ell + 1}$ , and observe that

$\displaystyle K^2 = d\Omega^1 \oplus H^1(\Omega)$

where the de Rham cohomology $H^1(\Omega)$ is a finite-dimensional space that does not depend on h, so we can solve it separately. So we may assume that $B = dA$ for some 1-form A. Without loss of generality, we may assume that A is in Coulomb gauge, in which case

$\displaystyle \begin{cases} d * dA = J \\ d * A = 0 \\ A|_{\partial \Omega} = 0 \end{cases}$

Here we don’t actually care about the fact that $d * A = 0$ , so if we pick up some numerical error there, it doesn’t matter. After all, we actually care about $dA$ , not A itself, and $dA$ is left invariant by gauge transformations of A.

The Whitney complex. It remains to find suitable shape $\ell$ -forms for all $\ell \in \{0, \dots, 3\}$ . We cannot choose them arbitrarily, however, because we must preserve unisolvence of the maps $d_h: \Omega^\ell \to K^\ell/H^\ell(\Omega)$ . The natural norm to put on $\Omega^\ell$ turns out to be

$\displaystyle \|\varphi\| := \|\varphi\|_{L^2} + \|d\varphi\|_{L^2}$

so that d is continuous and is elliptic modulo its kernel and cokernel. One can then use this fact, plus the theory of Fortin operators, to show that the discrete d on the space of shape forms is well-conditioned.

A natural way to preserve unisolvence is to find a family of morphisms of chain complexes $\Pi_h: \Omega^\bullet \to V_h^\bullet$ where $V_h^\bullet$ is a family of chain complexes associated to the triangulations $\mathcal T_h$ . This was more or less accomplished by Nédélec in the 1980s, but not in this language, and after quite a lot of work. However, the modern viewpoint is that Nédélec forms are nothing more than Whitney forms, which were discovered by Whitney decades prior:

Definition. Let $\sigma$ be a d-simplex, and let C be a simplicial k-chain in $\sigma$ . The Whitney form dual to C is the unique k-form $\varphi$ on $\sigma$ such that:

$\varphi_I$ is a linear function on $\sigma$ for every increasing multiindex I.

If $\iota: \tau \to \sigma$ is a k-face of $\sigma$ , then $(\iota^* \varphi)_I$ is a constant function on $\tau$ for every increasing multiindex I.

If $\iota: \tau \to \sigma$ is a k-face of $\sigma$ , then
$\displaystyle \int_\tau \varphi = \langle C, \tau\rangle$ .

The Whitney complex associated to a triangulation $\mathcal T_h$ is the chain complex of all $L^2$ k-forms $\varphi$ on $\Omega$ , such that:

For every d-simplex $\sigma \in \mathcal T_h$ , $\varphi|_\sigma$ is a Whitney k-form on $\sigma$ .

For every k-face $\iota: \tau \to \Omega$ in $\mathcal T_h$ , the pullback $\iota^* \sigma$ exists in the sense of the Sobolev trace theorem. In particular, there is no jump discontinuity.

The Whitney projection $\Pi_h$ associated to a triangulation $\mathcal T_h$ maps a k-form $\varphi$ to the unique k-form $\Pi_h \varphi$ in the Whitney complex associated to $\mathcal T_h$ such that for every k-face $\tau$ ,

$\displaystyle \int_\tau \Pi_h \varphi = \int_\tau \varphi$ .

See the paper of Dodziuk for more on the construction of Whitney forms. Anyways, we have a version of the Bramble-Hilbert lemma, as proven by Bramble, Falk, and Whinther:

Theorem. For any k-form $\varphi$ ,

$\displaystyle \|(1 - \Pi_h)\varphi\| \lesssim h (\|\varphi\|_{W^{1, 2}} + \|d\varphi\|_{W^{1, 2}})$ .

Putting it all together. To solve the elliptic Maxwell system, we can now proceed as when we solved the Laplacian.
Let $P_h$ be the projection of the matrix of differential operators

$\displaystyle P_h = \begin{bmatrix}d * d \\ d * \end{bmatrix}$

to the first space in the Whitney complex. Then $P_h$ is clearly unisolvent and sparse, and one can show that it is well-conditioned as well. We conclude that we can solve

$\displaystyle P_h A_h = \begin{bmatrix}\Pi_h J \\ 0\end{bmatrix}$