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 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 . The canonical example is that P is the Laplacian, and then this system reads
The key thing to observe is that since P is differential, it is local: we can determine from the germ of u at x. Motivated by this observation, we divide and conquer by fixing a family of partitions of into simplices, where each simplex in $\mathcal T_h$ has edge lengths comparable to h. We then fix the space of PL shape functions, which are continuous functions on whose restriction to each simplex is a linear polynomial. This is a finite-dimensional, hence closed, subspace of the Sobolev space . 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 to to obtain a new equation . 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
thus for we have a good solution in . 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 in a practical amount of time, however:
- Unisolvence. The matrix must be square. This is no issue for the Laplacian, since .
- 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 of is controlled from below by the Dirichlet spectral gap of , which is bounded from below as long as is not too fractalline or degenerate:
- Sparsity. Most of the entries of 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 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
where B is the magnetic field, is the normal field to , and the electric current J is given. The right-hand side of the equation 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 as opposed to simply 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
where B is now viewed as a 2-form, J is now viewed as a 1-form, and is the inclusion map. This interpretation is quite natural, because it views B as something that we want to integrate over a surface , so that is now the magnetic flux.
Let be the space of -forms, the kernel of , and observe that
where the de Rham cohomology is a finite-dimensional space that does not depend on h, so we can solve it separately. So we may assume that for some 1-form A. Without loss of generality, we may assume that A is in Coulomb gauge, in which case
Here we don’t actually care about the fact that , so if we pick up some numerical error there, it doesn’t matter. After all, we actually care about , not A itself, and is left invariant by gauge transformations of A.
The Whitney complex. It remains to find suitable shape -forms for all . We cannot choose them arbitrarily, however, because we must preserve unisolvence of the maps . The natural norm to put on turns out to be
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 where is a family of chain complexes associated to the triangulations . 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 be a d-simplex, and let C be a simplicial k-chain in . The Whitney form dual to C is the unique k-form on such that:
- is a linear function on for every increasing multiindex I.
- If is a k-face of , then is a constant function on for every increasing multiindex I.
- If is a k-face of , then
The Whitney complex associated to a triangulation is the chain complex of all k-forms on , such that:
- For every d-simplex , is a Whitney k-form on .
- For every k-face in , the pullback exists in the sense of the Sobolev trace theorem. In particular, there is no jump discontinuity.
The Whitney projection associated to a triangulation maps a k-form to the unique k-form in the Whitney complex associated to such that for every k-face ,
Theorem. For any k-form ,
Putting it all together. To solve the elliptic Maxwell system, we can now proceed as when we solved the Laplacian.
Let be the projection of the matrix of differential operators
to the first space in the Whitney complex. Then is clearly unisolvent and sparse, and one can show that it is well-conditioned as well. We conclude that we can solve
for some 1-form , and then satisfies
by the Bramble-Hilbert lemma. By Ladyzhenskaya’s inequality ad the fact that A is in Coulomb gauge,
Finally, by elliptic regularity,
so we conclude
and so B is a good approximation to J in , which is in addition divergence conforming.