Weyl's law

In this post we'll look at eigenfunctions and eigenvalues of the Laplace operator $\Delta$ on a domain $\Omega$ in $\mathbb{R}^d$. A function $\phi$ on $\Omega$ and a number $\lambda$ are an eigenfunction/eigenvalue pair if

$$-\Delta\phi = \lambda^2\phi$$

along with the Dirichlet boundary condition $\phi|_{\partial\Omega} = 0$. The operator $-\Delta$ is symmetric and positive-definite, so the eigenvalues are real and positive. I've chosen a slightly different way of writing things in terms of $\lambda^2$ because this makes the units of the eigenvalues an inverse length.

The Weyl asymptotic law describes how the eigenvalues grow as a function of the domain size and shape. Weyl proved in 1911 that, if $N(\lambda)$ is the number of eigenvalues of the Dirichlet Laplacian less than $\lambda$, that

$$N(\lambda) = (2\pi)^{-d}\omega_d\cdot\text{vol}(\Omega)\cdot\lambda^{d} + \mathscr{O}(\lambda^{d})$$

as $\lambda \to \infty$, where $\omega_d$ is the volume of the unit ball in $\mathbb{R}^d$. As a sanity check, note that $\lambda$ has units of length${}^{-1}$, so the formula above is dimensionless. As another sanity check, you can look at the analytical expression for the eigenvalues on a box or a sphere. The proof given in volume 1 of Courant and Hilbert is pretty easy to follow. Weyl conjectured that the second term could be expressed in terms of the area of the boundary:

$$N(\lambda) = (2\pi)^{-d}\omega_d\cdot\text{vol}(\Omega)\cdot\lambda^d - \frac{1}{4}(2\pi)^{1 - d}\omega_{d - 1}\cdot\text{area}(\partial\Omega)\cdot\lambda^{d - 1} + \mathscr{o}\left(\lambda^{d - 1}\right)$$

but this wasn't proved in his lifetime. Here we'll come up with a simple domain and show how you might verify this law numerically.

Making a mesh¶

First, we'll generate a mesh using the Python API for gmsh. The calls to add a plane surface and a physical plane surface are easy to forget but essential.

To make sure everything worked right, we'll visualize the mesh after loading it in.

Using SLEPc¶

To compute the eigenvalues and eigenfunctions of the Laplace operator, we'll use the Scalable Library for Eigenvalue Problem Computations (SLEPc). Both SLEPc and Firedrake are built on top of PETSc, so creating the eigenvalue problem is just a matter of assembling the Firedrake form objects representing our linear operators, extracting the underlying PETSc matrix objects, and passing them to SLEPc.

To solve the right problem and to help SLEPc get the right answer we'll pass it several options. First, we're solving a generalized Hermitian eigenproblem. Since the eigenproblem is Hermitian, all the eignevalues are real, which is a very convenient simplifying assumption.

For this problem we're going to use a spectral transformation. Rather than find the eigenvalues of a matrix $A$ directly, we'll instead find the eigenvalues of a matrix $f(A)$ where $f$ is invertible and holomorphic on a domain containing the spectrum of $A$. We can then compute the eigenvalues of $A$ as the function $f^{-1}$ aplied to the eigenvalues of $f(A)$. The advantage of spectral transformations is that, with a good choice of $f$, the eigenvalues of $f(A)$ can be easier to compute than those of $A$ itself. Since $A$ is positive-definite and we're looking for the smallest eigenvalues, a good choice is

$$f(z) = 1/(z - \sigma),$$

i.e. shifting and inverting. This spectral transformation is equivalent to finding the eigendecomposition of $(A - \sigma M)^{-1}$. Computing the inverse of a matrix is generally a bad idea, but under the hood it's enough to be able to solve linear systems.

Anything in SLEPc having to do with spectral transformations is prefixed with st. In our case, we're using the shift-and-invert transformation (sinvert). To solve these linear systems, we'll a Krylov subspace method (ksp_type) with some preconditioner (pc_type). Since $A$ is symmetric and positive-definite, we can use the conjugate gradient method (cg).

To check that everything worked right, we can see how many eigenvalues converged:

Just for fun, we can plot one of the eigenfunctions. The zero contours of eigenfunctions are a fascinating subject -- the Courant nodal domain theorem tells us that the $n$-th eigenfunction can have no more than $n$ nodal domains.

The following plot shows exact eigenvalue counting function and the order-1 and order-2 approximations from Weyl's law.

The accuracy difference is even more stark if we look at the relative error in the eigenvalue counting function.

The order-1 approximation is pretty good, but the order-2 approximation is startlingly accurate. Of course we've only looked at the first few hundred eigenvalues on a mesh with several thousand vertices. Once the corresponding wavelengths get close to the diameter of a triangle of our mesh, I'd expect the approximation to break down. The mesh is too coarse at that point to resolve the highly oscillatory eigenfunctions.

Conclusions¶

The Weyl asymptotic law has some interesting physical implications. The first-order version of the law tells us that you can hear the area of a drumhead by fitting the sequence of harmonic frequencies to the right power. The second-order version of the law tells us that you can, in the same way, hear the perimeter of the drumhead by fitting the remainder of the first-order approximation.

Victor Ivrii gave a proof in 1980 of the Weyl law up to second order, under some special conditions that are thought to hold for a wide class of domains. While proving the law up to first order is relatively elementary, Ivrii's proof used microlocal analysis, which is well and truly above my pay grade.