# Linear solvers: How to deal with a singular system? (Poisson equation with Neumann boundary conditions)

The question is in the context of iterative numerical solution of large PDE systems with Finite Differences or Finite Elements:

Stating the Poisson equation with Neumann boundary conditions will lead to a singular system because it is invariant when adding a constant function. An often proposed option is to just "pin" the system to a fixed value in one point, which leads to a well-defined system.

In my case, I'm rather interested in a solution where the constant part vanishes, so "pinning" a certain point does not help.

Furthermore, I'd like to use an iterative solver, and those seem to be significantly unhappy with a singular system. In the PETSc documentation (I may or may not use PETSc then) I found the option to set a null space of the system.

So how should I deal with that singluar system? I'm interested in some theoretical views as well as practical solutions using linear solver libraries.

• have you done a search on the site? there are quite a few questions (and answers) dealing with pure-Neumann conditions for Poisson systems, such as this one or this one. The latter q is concerned with biCGstab for a pure-Neumann problem, which seems very close to what you want. Jul 9 '17 at 12:45
• You can always post-process your solution to make the constant part vanish; simply subtract the mean.
– cfh
Jul 10 '17 at 8:49

Krylov solvers have no problem converging if the system is consistent, i.e., if the right-hand side is in the image of the system matrix $A$; see, e.g., Iterative Krylov Methods for Large Linear Systems by van der Vorst. PETSc's CG solver should work.

You can even use a multigrid preconditioner if you make sure that the coarse solver preserves the null space. PETSc's default ILU solver is not sufficient, better use Jacobi; see http://lists.mcs.anl.gov/pipermail/petsc-users/2012-February/012139.html.

The problem I imagine you are trying to solve is a diffusion equation with source term with homogeneous Neumann BC. To do so it must be well posed in order to obtain physical and good results.

The resultant system leads to a singular differential operator matrix $A$, once the BC have been applied. This is so, as you already mentioned, because the following problem (for example I use a linear equation because for the nonlinear case the same would apply ):

$$\left\{\begin{array}{ll}% -\partial_x^2u=f & x\in(0,1)\\ \partial_xu=0 & x=0,x=1 \end{array}\right. \tag{*}$$

is invariant under the transformation $u\to u+constant$.

And therefore some value for $u$ must be specify $\textbf{inside}$ the domain. Many people impose a Dirichlet BC instead one of the Neuman BC causing the solution to differ from the actual one, which now solves:

$$\left\{\begin{array}{ll}% -\partial_x^2u=f & x\in(0,1)\\ u=0 & x=0\\ \partial_xu=0 & x=1\\ \end{array}\right.$$

The problem $(*)$ is well posed if $\int_{0}^{1}{f\,dx}=0$, in fact the source function I propose:

$$f=cos(2\pi x)$$ fits well for the well-posedness of $(*)$, only some reference value for $u$ is required to manage the uniqueness, for example we require that $u(1/2)=1$ $\textbf{without any loss of generality}$. I put these words in bold due to your requirement that your constant must be set to $0$.

N = 100; % #nodes
uref = 1; %Ref value for u for centre node
deltax = 1/(N-1); % step
x = (0:deltax:1)';
f = cos(x*2*pi); %distributed source

b = zeros(N,1); % Source vector
A=zeros(N,N); % Stiffness matrix

for i = 2 : N-1
A(i,i-1:i+1) = -1/deltax^2*[1, -2, 1];
b(i) = f(i);
end

%Neumann BC
A(1,1:2) = 2/deltax^2*[1,-1];
b(1) = f(1);
A(N,N-1:N) = 2/deltax^2*[-1, 1];
b(N) = f(N);

%Uniqueness condition e.g. central node
idx = floor(N/2);
A(idx,:) = 0;
A(idx,idx) = 1;
b(idx) = uref;

% System solution
u = A\b;

% Verify the diff equation
ddu = zeros(N,1);
for i = 2 : N-1
ddu(i) = -(u(i+1)-2*u(i)+u(i-1))/(deltax^2);
end
ddu(1) = f(1);
ddu(N) = f(N);

%Solution residual
res = norm(ddu-f);

%Plots
plot(x,u,'r','linewidth',2); %Solution
figure
plot(x,ddu-f,'k','linewidth',1) % Error


The above code produces for $(*)$ a solution and the pointwise residual given by $e=\partial_x^2 u+f$ is given below for reference: You now are free to choose the arbitrary constant (forget the fact that $u(1/2) = 1$) to your problem, just add it.

Formally the uniqueness condition that you mention is simply imposed by the restriction: $$\int_{0}^{1}{u\,dx}=0$$ which constrains $u$ in a way that the transformation $u\to u+constant$ is not valid any more, and therefore the constant must be zero.

This condition would be imposed in an analogous manner when solving the system:

%Uniqueness condition
idx = N-1;
A(idx,:) = 1;
b(idx) = 0;


Or at the end, when $u_{calc}$ has been obtained, you do: $$u=u_{calc}-\overline{u}_{calc}=u_{calc}-\int_{0}^{1}{u_{calc}\,dx}$$

• (remark) It is not pointwise error, it is pointwise residual. Jul 10 '17 at 17:41
• I thought for a while about your idea - it will work badly in multidimensional case, won't it? You will make a solution constant on a line where actually it is not necessarily constant. Jul 10 '17 at 19:13
• No. You only set the value of one point... not of a line. @VorKir
– HBR
Jul 10 '17 at 21:49