# preconditioner for Laplace “without” boundary values

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:

import matplotlib.pyplot as plt
import numpy as np
from dolfin import (
EigenMatrix,
FacetNormal,
FunctionSpace,
TestFunction,
TrialFunction,
UnitIntervalMesh,
UnitSquareMesh,
assemble,
dot,
ds,
dx,
)

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)

out = np.linalg.eigvals(A.toarray())
plt.plot(out.real, out.imag, "o")
plt.show()

• Have you tried multi-grid? – Alone Programmer Mar 29 at 2:41
• 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. – Wolfgang Bangerth Mar 29 at 3:32
• @AloneProgrammer The matrix isn't even symmetric so that's going to be tough. Unless you have a hint? – Nico Schlömer Mar 29 at 7:43
• @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. – Nico Schlömer Mar 29 at 7:44
• @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. – Alone Programmer Mar 29 at 14:24

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.