Finding eigenvalues of the Laplacian on solenoidal (divergence-free) vector fields

Question

In Mathematica it is easy to find eigenvalues of the Laplacian in simple cases. For example, on $\Omega\in \mathbb{R}^2$:

{vals, funs} = NDEigensystem[{-Laplacian[u[x, y], {x, y}], DirichletCondition[u[x, y] == 0, True]}, u[x, y], {x, y} \[Element] \[CapitalOmega], 6]

But is it possible to restrict the abovementioned problem to solenoidal (divergence-free) vector fields? That is, how could I numerically find $(u, \lambda)$, such that:

$-\Delta u=\lambda u\\ \nabla\cdot u=0$

With Dirichlet boundary conditions on $\partial \Omega$ and $u=u(x,y)$, $(x,y)\in \Omega\subset\mathbb{R}^2$. Thank you very much!

Answer 1

This is not really an answer, but may help to find a solution.

Iterative approach

What might work is solving the Stokes eigenvalue problem

$$- P \, A \, P \, v = \lambda \, P \, M \, P\, v$$

with the naive projector $P = I - B^T\, (B B^T)^{-1} \, B$ onto $\operatorname{ker}(B)$. Here $A$ is the stiffness matrix, $M$ is the mass matrix, and $B$ is the finite-element discretization of $\operatorname{div}$, all with respect a suitable (stable!) finite element discretization. Then $u = P \, v$ should be what you are looking for.

The matrix $B^T\, (B B^T)^{-1} \, B$ is dense, so I would not recommend to assemble it; instead, only its action should be implemented in a matrix-free way (by exploiting a sparse $LU$-factorization of $B B^T$). While the matrices $M$, $A$, $B$ might be obtainable from Mathematica, Mathematica's Arnoldi method does (as far as I know) not support matrix-free methods. Also one certainly wants to use a good preconditioner (e.g., geometric multigrid), which is also not available out of the box. (Notice that we cannot use ILU-preconditioner due to the absence of any concrete matrix.)

Alternate approach

Alternatively, one could study the generalized eigenvalue problem $$ \begin{pmatrix} A &B^T\\ B &0 \end{pmatrix} \begin{pmatrix} u \\ p \end{pmatrix} = \lambda \, \begin{pmatrix} M &0\\ 0 & \varepsilon \, I \end{pmatrix} \begin{pmatrix} u \\ p \end{pmatrix} $$ with some small $\varepsilon>0$.

That could be set up by

\[Epsilon] = 10^-12;
AA = ArrayFlatten[{{A, B\[Transpose]}, {B, 0.}}];
MM = ArrayFlatten[{{M, 0. B\[Transpose]}, {0. B, \[Epsilon] IdentityMatrix[Length[B],SparseArray, WorkingPrecision->MachinePrecision]}}];

and solved with

{\[CapitalLambda], U} = Eigensystem[{AA, MM}, -6, Method -> "Arnoldi"];

A positive $\varepsilon$ is necessary, otherwise, the Arnoldi solver will complain. But don't ask me how stable or how accurate this might be in practice.