3
$\begingroup$

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!

Improve this question
$\endgroup$
2
  • $\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

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\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$ Commented Oct 17, 2019 at 4:53

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.