I'm looking at solving systems with the FEM discretization $$ -\int_\Omega (\Delta u) v = \int_\Omega \nabla u \cdot \nabla v - \int_\Gamma (n\cdot\nabla u) v. $$ without applying Dirichlet- or Neumann-type boundary conditions. The resulting matrix is generally not self-adjoint (except in for 1D meshes).

The kernel consists of all linear functions over the domain, so there are always $n+1$ eigenvalues 0 (where $n$ is the dimensionality of the mesh).

The right-hand side is such that the system is consistent and I would like to find a solution of the system. The idea is to use GMRES, starting off with an all-zero initial guess.

Is there a good preconditioner for this known in literature? I played around with Dirichlet- and Neumann-Laplace without much success.

To get a feel, here's how to create the matrix and plot the spectrum in FEniCS/SciPy:

enter image description here

import matplotlib.pyplot as plt
import numpy as np
from dolfin import (

def _assemble_eigen(form):
    L = EigenMatrix()
    assemble(form, tensor=L)
    return L

k = 20
# mesh = UnitIntervalMesh(n)
mesh = UnitSquareMesh(k, k)

degree = 1
V = FunctionSpace(mesh, "CG", degree)

u = TrialFunction(V)
v = TestFunction(V)
mesh = V.mesh()
n = FacetNormal(mesh)

A = _assemble_eigen(dot(grad(u), grad(v)) * dx - dot(n, grad(u)) * v * ds).sparray()

out = np.linalg.eigvals(A.toarray())
plt.plot(out.real, out.imag, "o")
  • $\begingroup$ Have you tried multi-grid? $\endgroup$ – Alone Programmer Mar 29 at 2:41
  • $\begingroup$ The question is as ill-posed as your problem. If you say you want to solve a linear system with this matrix, then this will not work: You have a singular matrix, so there will be infinitely many solutions. If you have a specific solution in mind, you have to add this information either to your system, or at least to the question here. $\endgroup$ – Wolfgang Bangerth Mar 29 at 3:32
  • $\begingroup$ @AloneProgrammer The matrix isn't even symmetric so that's going to be tough. Unless you have a hint? $\endgroup$ – Nico Schlömer Mar 29 at 7:43
  • $\begingroup$ @WolfgangBangerth I clarified the question: I'm just looking for a solution, and I usually do that by starting off with Krylov and using an all-0 initial guess. $\endgroup$ – Nico Schlömer Mar 29 at 7:44
  • 2
    $\begingroup$ @NicoSchlömer OK, got it. I suggest you to use PETSC with a multi-grid preconditioner. The fact that your matrix is not symmetric is not a big problem. At least giving it a shot worth to try. By the way, if you realized multi-grid is not a good choice here there are plenty of other choices available in PETSC such as SOR. $\endgroup$ – Alone Programmer Mar 29 at 14:24

Here is at least an idea, whether it works is a different question.

Let's say you sort unknowns so that you have the ones in the interior of the domain first, and then all those at the boundary. Then the matrix that corresponds to your problem decomposes in the following way: $$ A = \begin{pmatrix} A^{\circ\circ} & B^{\circ\partial} \\ C^{\partial\circ} & D^{\partial\partial} \end{pmatrix} $$ where $\circ$ indicates shape functions in the interior and $\partial$ on the boundary of the domain. $A^{\circ\circ}$ is simply the Dirichlet-boundary condition matrix for the Laplace operator and we know how to invert it efficiently. One could then think about using a Schur complement approach where you first solve for the boundary unknowns and then the interior unknowns, or maybe the other way around.

If you follow the arguments of the Silvester-Wathen preconditioner for the Stokes system, you will also be able to construct a good preconditioner for the whole matrix based on the Schur complement approach; it will probably look something like this: $$ P^{-1} = \begin{pmatrix} A^{\circ\circ} & B^{\circ\partial} \\ 0 & S^{\partial\partial} \end{pmatrix}^{-1} $$ where $S^{\partial\partial}=D^{\partial\partial}-D^{\partial\circ}[A^{\circ\circ}]^{-1}B^{\circ\partial}$ is the Schur complement. It would be worthwhile thinking about whether one can approximate $S^{\partial\partial}$ by a simpler matrix, in the same way as for the Stokes equations, one approximates the Schur complement by the mass matrix.


If you are not enforcing any boundary conditions, you should remove the nullspace of the problem (to make sure that there is a unique solution, MatSetNullSpace can be used to achieve this in PETSc) and then use multigrid on the problem (for example, PCMG of PETSc).

What surprises me is that there are n+1 zero eigenvalues. I would expect only one. Are you solving the Laplace's equation on surfaces (like a 2D surface embedded in 3D space) and the boundary is the boundary of that surface (like 1D boundary embedded in 3D space)? That is the only scenario I can imagine that the nullspace is that populous.

  • $\begingroup$ The kernel consists of all linear functions over the domain, that's how you can see that dimensionality. (I've added this to the original question.) Q: Why should you remove the null space before applying multigrid? And what makes you think multigrid is a good choice here? The system has complex-valued eigenvalues after all. $\endgroup$ – Nico Schlömer Mar 30 at 7:00
  • 1
    $\begingroup$ Heuristically, by removing the nullspace you are preventing false divergence which results from the accumulation of coefficients corresponding to the basis vectors of the kernel. More mathematically, but still in a handwavy way, what you need for multigrid convergence is ellipticity on the kernel of the operator -which you have- (I couldn't find a direct citation for this, but check Boffi, Brezzi and Fortin's Mixed FEM and Applications). Removing the nullspace allows for effective restriction and prolongation. Lastly, I don't know whether multigrid is the best choice. But it is a must-try. $\endgroup$ – Abdullah Ali Sivas Mar 30 at 15:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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