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}

for some 1-form A_h, and then B_h := dA_h satisfies

\displaystyle \|B_h - B\|_{L^2} \leq \|A - A_h\| \lesssim h(\|A\|_{W^{1, 2}} + \|B\|_{W^{1, 2}})

by the Bramble-Hilbert lemma. By Ladyzhenskaya’s inequality ad the fact that A is in Coulomb gauge,

\displaystyle h(\|A\|_{W^{1, 2}} + \|B\|_{W^{1, 2}}) \lesssim h\|B\|_{W^{1, 2}}

Finally, by elliptic regularity,

\displaystyle h\|B\|_{W^{1, 2}} \lesssim h\|J\|_{L^2}

so we conclude

\displaystyle \|B_h - B\|_{L^2} \lesssim h\|J\|_{L^2}

and so B is a good approximation to J in L^2, which is in addition divergence conforming.


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s