Poisson equation finite-difference with pure Neumann boundary conditions

Question

I'm trying to solve a 1D Poisson equation with pure Neumann boundary conditions. I've found many discussions of this problem, e.g.

1) Poisson equation with Neumann boundary conditions

2) Writing the Poisson equation finite-difference matrix with Neumann boundary conditions

3) Discrete Poisson Equation with Pure Neumann Boundary Conditions

4) Finite differences and Neumann boundary conditions

And many more.

Some of the answers seem unsatisfactory though. For example, the answer in 2) contains a document that explains how the matrix $A$ changes when applying Neumann boundary conditions, but does not explain how to solve the singular system. The comment left in 2) by @Evgeni Sergeev is a reference to a problem with mixed boundary conditions, and not pure Neumann boundary conditions.

That said, there are some useful things I've found. I do like the second comment given in the 1) by @Sumedh Joshi. It suggests to subtract the mean from the RHS. Other references I've read suggest using Dirichlet BCs at one location in the computational domain which should result in a unique solution, however, I prefer removing the mean since this seams more elegant.

The problem setup is:

$\frac{\partial^2 u}{\partial x^2} = f, \qquad f = -(2\pi)^2 cos(2\pi x), \qquad 0 \le x\le 1, \qquad \left(\frac{\partial u}{\partial x}\right)_{0} = 0, \qquad \left(\frac{\partial u}{\partial x}\right)_{1} = 0$

I have ghost points in my computational domain, so my matrix $A$ has some extra zeros (I remove them later). Using central differences everywhere in the computational domain, matrix $A$ takes the following form:

$A = \frac{1}{\Delta x^2} \left[ \begin{array}{ccccccccc} 0 & 0 & 0 & & & & & & \\ 1 & -2 & 1 & & & & & & \\ 0 & 1 & -2 & 1 & & & & & \\ & & & \ddots & \ddots & \ddots & & & \\ & & & & & 1 & -2 & 1 & 0 \\ & & & & & & 1 & -2 & 1 \\ 0 & & & & & & 0 & 0 & 0 \\ \end{array} \right] $

Once applying the ghost points and adjusting the matrix $A$, I end up with

$A = \frac{1}{\Delta x^2} \left[ \begin{array}{ccccccccc} 0 & 0 & 0 & & & & & & \\ 0 & -2 & 2 & & & & & & \\ 0 & 1 & -2 & 1 & & & & & \\ & & & \ddots & \ddots & \ddots & & & \\ & & & & & 1 & -2 & 1 & 0 \\ & & & & & & 2 & -2 & 0 \\ 0 & & & & & & 0 & 0 & 0 \\ \end{array} \right] $

Correspondingly, the RHS has changed to

$\frac{- 2u_b + 2u_i}{\Delta x^2} = f_b - \frac{-2\theta \Delta x}{\Delta x^2}$

Where $u_b,u_i,f_b,\theta$ are the $u$ at the boundary, $u$ at the first interior point, $f$ at the boundary and the slope of the solution at the boundary (zero in this case) respectively. The interior of this matrix is the same as the one discussed in the document given in the accepted answer in 2).

Here I show the setup and results of a small MATLAB script where I print out $A, f, mean(f)$ and the (attempted) solution $u$ using MATLAB's backslash operator.

  A =

   -100.0000  100.0000         0         0         0         0         0         0         0         0         0
    100.0000 -200.0000  100.0000         0         0         0         0         0         0         0         0
           0  100.0000 -200.0000  100.0000         0         0         0         0         0         0         0
           0         0  100.0000 -200.0000  100.0000         0         0         0         0         0         0
           0         0         0  100.0000 -200.0000  100.0000         0         0         0         0         0
           0         0         0         0  100.0000 -200.0000  100.0000         0         0         0         0
           0         0         0         0         0  100.0000 -200.0000  100.0000         0         0         0
           0         0         0         0         0         0  100.0000 -200.0000  100.0000         0         0
           0         0         0         0         0         0         0  100.0000 -200.0000  100.0000         0
           0         0         0         0         0         0         0         0  100.0000 -200.0000  100.0000
           0         0         0         0         0         0         0         0         0  100.0000 -100.0000


  F =

    -19.7392
    -31.9387
    -12.1995
     12.1995
     31.9387
     39.4784
     31.9387
     12.1995
    -12.1995
    -31.9387
    -19.7392


  meanF =

   -9.6892e-016

  Warning: Matrix is singular to working precision.
  > In main at 121

  u =

     NaN
     NaN
     NaN
     NaN
     NaN
     NaN
     NaN
     NaN
     NaN
    -Inf
    -Inf

My question is why is this not working? I thought that removing the mean of the RHS would work as suggested in @Sumedh Joshi's comment. Any help is greatly appreciated.

UPDATE:

I found out that this same exact problem setup works perfectly fine for 12 unknowns (and 15 and much larger, e.g. 100), but does not work for 10 unknowns (as posted). So I suppose that the problem is setup correctly. This still begs the question though, what is going on here? It seems that this may be more of a question regarding numerical analysis.

Here are the matlab results for 12 unknowns:

 A =
 0         0         0         0         0         0         0         0         0         0         0         0         0         0         0
 0   -1.0    1.0         0         0         0         0         0         0         0         0         0         0         0         0
 0    1.0   -2.0    1.0         0         0         0         0         0         0         0         0         0         0         0
 0         0    1.0   -2.0    1.0         0         0         0         0         0         0         0         0         0         0
 0         0         0    1.0   -2.0    1.0         0         0         0         0         0         0         0         0         0
 0         0         0         0    1.0   -2.0    1.0         0         0         0         0         0         0         0         0
 0         0         0         0         0    1.0   -2.0    1.0         0         0         0         0         0         0         0
 0         0         0         0         0         0    1.0   -2.0    1.0         0         0         0         0         0         0
 0         0         0         0         0         0         0    1.0   -2.0    1.0         0         0         0         0         0
 0         0         0         0         0         0         0         0    1.0   -2.0    1.0         0         0         0         0
 0         0         0         0         0         0         0         0         0    1.0   -2.0    1.0         0         0         0
 0         0         0         0         0         0         0         0         0         0    1.0   -2.0    1.0         0         0
 0         0         0         0         0         0         0         0         0         0         0    1.0   -2.0    1.0         0
 0         0         0         0         0         0         0         0         0         0         0         0    1.0   -1.0         0
 0         0         0         0         0         0         0         0         0         0         0         0         0         0         0

 f =

        0
 -19.7392
 -34.1893
 -19.7392
  -0.0000
  19.7392
  34.1893
  39.4784
  34.1893
  19.7392
   0.0000
 -19.7392
 -34.1893
 -19.7392
        0

 mean(f)
 ans =
 -9.4739e-016

 A(2:end-1,2:end-1)\f(2:end-1)
 ans =
  0.9167
  0.7796
  0.4051
 -0.1065
 -0.6181
 -0.9926
 -1.1297
 -0.9926
 -0.6181
 -0.1065
  0.4051
  0.7796
  0.9167

Answer 1

I know this question is old but I found an error that I wanted to clear up in @DrHansGruber 's answer. It took me a few hours to spot the issue and I wanted to post the solution to save anyone else from working through this.

First, as has been stated by Hans and Peter, there are infinitely many solutions to the 1D Poisson equation with Neumann BCs at both end points. They differ by a constant of integration.

In order to determine the constant of integration, the value of the function at a given point or the average value over an interval is required. It's important to note that, even after specifying the function's value at a point (or its average over an interval), the resulting solution is not the same as the solution obtained from the 1D Poisson equation with Dirchlet on one side and Neumann on the other.

If you look at Gruber's Matlab code you can see that when he solves the problem with Neumann boundaries on both sides, he enforces the discrete compatibility condition: $$ \sum f = 0 $$ by subtracting the mean of $f$ from the RHS vector before feeding it to the conjugate gradient method.

When he specifies the value of $u$ at the left boundary by enforcing a Dirichlet boundary there, however, he does not enforce the discrete compatibility condition. The discrete compatibility condition "encodes" the second Neumann boundary condition into the solution. Without this compatibility condition we are solving the following problem: $$ \frac{d^2u}{dx^2}=f $$ $$ u(0)=a $$ $$ \frac{du}{dx}(1)=0 $$

If, in Gruber's code, you make the following changes, you will obtain identical solutions using Matlab's backslash command with Neumann BCs and a given point specified, and using CG with only Neumann BCs. I'm not familiar with CG method but my guess is that somewhere in the iterations a given constant is being selected for you.

clear;

n = 64;
lx = 1.0;
% BCA = 0; 
% BCA = -0.1533;
BCA = 0.0883;   % CHANGED TO MATCH VALUE OF NEUMANN soltuion

% problem setup
dx=lx/(n-1);
x=0:dx:1;

% set up A's to compute Laplacian on internal points
nx=n-2;
e=ones(nx-2,1);
% left Dirichelt, right Neumann
AD = - spdiags([1 -2 1; e*[1 -2 1]; 1 -1 1],-1:1,nx,nx)/dx^2;
% both sides Neumann
AN = - spdiags([1 -1 1; e*[1 -2 1]; 1 -1 1],-1:1,nx,nx)/dx^2;

% set up RHS for internal points
b=zeros([nx 1]);
b(1:nx)=4*sin(33/pi.*x(2:nx+1));
%b(1:nx)=-2*pi*cos(2*pi.*x(2:nx+1));

% incorporate left Dirichelt BC conditon into RHS
% e.g. p.2 of https://web.stanford.edu/class/cs205b/lectures/lecture16.pdf
bD=b-mean(b);  %CHANGED SO THAT RHS ENCODES COMPATIBILITY CONDITION!
bD(1)=-AD(2,1)*BCA; % fix xD(1) to random value BCA by adjusting rhs

% ensure Discrete Compatibility for Neumann case, i.e. mean(bN)=0
bN=b-mean(b);

% solve both with CG
xD=AD\bD; %pcg(AD,bD,1e-6,2*n);
xN=pcg(AN,bN,1e-6,2*n);

% solutions are obtained on internal points only and boundary conditions
% are incorporated into the RHS. For plotting values need to be appended.
xD=[ BCA ; xD; xD(end)];
xN=[xN(1); xN; xN(end)];

% compute first derivatives, frequently of interest for potential problems
% and also to see the Neumann boundary conditions
dxDdx=gradient(xD,dx);
dxNdx=gradient(xN,dx);

% plot solutions
f=figure(1); set(f,'Position', [200, 300, 720, 310]);
subplot(1,2,1);
plot(x,xD,'r-'); hold on
plot(x,xN,'b--'); hold off
title('Solutions'); grid on;

subplot(1,2,2);
plot(x,dxDdx,'r-'); hold on
plot(x,dxNdx,'b--'); hold off
title('First Derivatives of Solutions'); grid on;
% legend('left Dirichelt','both Neumann','location','SouthEast');
legend('left Dirichelt','both Neumann','location','NorthEast');

I hope this helps someone else. I don't have enough reputation or I would have left a comment.

Answer 2

First a short answer - it did not work because your matrix is singular. The matrix is singular, because you did not do all steps necessary for the case of zero Neumann boundary condition.

A longer answer - I guess your expectation it should work (I might be wrong) were due to the fact that the mean of right hand side is in your case (numerically almost) zero. There exists a step in this kind of problems where a mean of right hand side is subtracted from it, but this is not enough, and in your case it is not necessary because it is equal 0.

In fact, the zero mean of $f$ (or zero sum of entries in your $f$) means in your case that there exists a solution to your problem, in fact, there are infinitely many solutions. You still have to do some additional work to find at least one solution from those infinitely many ones.

Let me first note how to see easily that your matrix is singular. I mean the matrix in your output from Matlab, i.e. without the first and last row and column with zero entries only. If you sum all rows in your matrix $A$ you get the row with zeros only. This can happen only for singular matrix. In fact it means that you can delete e.g. the last row in your matrix, because it is redundant, it can be expressed from other rows.

Luckily, if you do the same with the right hand side $f$ then you get also the zero (it means also the mean of $f$ is zero). If this would not be case, your system has no solution (a brief explanation - it would be like you have the zero row in the matrix with non zero right hand side, such equation can not be fulfilled).

So now what to do to find at least one representative solution. In general you must add some additional equation that makes the solution unique. It means you have to add one row in the matrix and corresponding one right hand side.

One option you mentioned in your question (and it is the most typical one) is to say e.g. that the last unknown (for which in fact you erased the row in $A$) you prescribe an arbitrary value. It is like you replace the zero Neumann boundary condition with Dirichlet one. So it is like to add the row $(0 \, \ldots \, 0 \, 1)$ to replace the erased one.

Another one is that you require that the mean of unknowns has a prescribed value, e.g. zero. It is like you have to add a row with ones only, i.e. $(1 \, \ldots 1)$ to replace the erased row. That is a reason that such option is not very popular, because it changes the structure of your matrix.

Answer 3

As stated by Peter the matrix of the all Neumann boundary Poisson problem is indeed singular and therefore the Matlab Backslash operator, which attempts to use direct solvers, cannot deal with it. Conjugate Gradient on the other hand can, you just need to, as you have, ensure the mean of f equals 0 before handing it to e.g. PCG.

I would like to stress that arbitrarily replacing a Neumann boundary with a fixed (arbitrary) valued Dirichelt boundary is not the correct way of dealing with the problem.

To demonstrate this, here is some code for anyone to verify in Matlab or Octave:

clear;

n = 64;
lx = 1.0;
BCA = 0; %-0.1533;

% problem setup
dx=lx/(n-1);
x=0:dx:1;

% set up A's to compute Laplacian on internal points
nx=n-2;
e=ones(nx-2,1);
% left Dirichelt, right Neumann
AD = - spdiags([1 -2 1; e*[1 -2 1]; 1 -1 1],-1:1,nx,nx)/dx^2;
% both sides Neumann
AN = - spdiags([1 -1 1; e*[1 -2 1]; 1 -1 1],-1:1,nx,nx)/dx^2;

% set up RHS for internal points
b=zeros([nx 1]);
% b(1:nx)=4*sin(33/pi.*x(2:nx+1));
b(1:nx)=-2*pi*cos(2*pi.*x(2:nx+1));

% incorporate left Dirichelt BC conditon into RHS
% e.g. p.2 of https://web.stanford.edu/class/cs205b/lectures/lecture16.pdf
bD=b;
bD(1)=-AD(2,1)*BCA; % fix xD(1) to random value BCA by adjusting rhs

% ensure Discrete Compatibility for Neumann case, i.e. mean(bN)=0
bN=b-mean(b);

% solve both with CG
xD=AD\bD; %pcg(AD,bD,1e-6,2*n);
xN=pcg(AN,bN,1e-6,2*n);

% solutions are obtained on internal points only and boundary conditions
% are incorporated into the RHS. For plotting values need to be appended.
xD=[ BCA ; xD; xD(end)];
xN=[xN(1); xN; xN(end)];

% compute first derivatives, frequently of interest for potential problems
% and also to see the Neumann boundary conditions
dxDdx=gradient(xD,dx);
dxNdx=gradient(xN,dx);

% plot solutions
f=figure(1); set(f,'Position', [200, 300, 720, 310]);
subplot(1,2,1);
plot(x,xD,'r-'); hold on
plot(x,xN,'b--'); hold off
title('Solutions'); grid on;

subplot(1,2,2);
plot(x,dxDdx,'r-'); hold on
plot(x,dxNdx,'b--'); hold off
title('First Derivatives of Solutions'); grid on;
% legend('left Dirichelt','both Neumann','location','SouthEast');
legend('left Dirichelt','both Neumann','location','NorthEast');

This will produce the following output:

One might be inclined to think that fixing the left boundary and thus being able to use the Backslash operator has done the job. The agreement between the two is good enough, however for a different rhs the results can differ significantly. Changing the rhs to e.g.:

b(1:nx)=4*sin(33/pi.*x(2:nx+1));

Produces the following output:

Here both the solutions and their first derivatives are very different, especially near where the boundary condition of the problem has been changed.

There are many different ways of dealing with singular matrices, but using an arbitrary Dirichlet boundary condition is not one. For the problem at hand it is easiest to switch to an iterative solver like PCG and ensure Discrete Compatibility.