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!

  • $\begingroup$ 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$? $\endgroup$ Commented Oct 15, 2019 at 15:22
  • $\begingroup$ 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. $\endgroup$ Commented Oct 15, 2019 at 15:45

1 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.

  • $\begingroup$ 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. $\endgroup$
    – all
    Commented Oct 17, 2019 at 4:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.