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

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!

• Are you sure that you do not mean the Stokes eigenvalue problem $-\Delta u + \operatorname{grad} p = \lambda u$, $\operatorname{div} u = 0$ with $\int p(x) \, \mathrm{d} x =0$? – Henrik Schumacher Oct 15 '19 at 15:22
• See also here for an ansatz that utilizes a stream function $\psi$ with $u = \operatorname{curl} \psi$ to recast the eigenvalue problem into an eigenvalue problem for the bi-Laplacian. However, if I am not mistaken, this works only in simply-connected domains. – Henrik Schumacher Oct 15 '19 at 15:45

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.

• Henrik, thank you for the detailed answer! I've checked the article you mentioned - it is rather useful. I'm new to numerical methods and now I see that the problem is not so simple as I expected. – all Oct 17 '19 at 4:53