I am trying to solve the Poisson equation in a rectangular domain using a finite difference scheme with a rectangular mesh.

I have happily generated the matrix system of equations Ax = b which is required to be solved, but when I try to impose Neumann boundary conditions (using ghost cells) of zero gradient I am running in to problems, I believe because the matrix is singular.

I know that there are several similar threads on this topic: Discrete Poisson Equation with Pure Neumann Boundary Conditions , Writing the Poisson equation finite-difference matrix with Neumann boundary conditions , Solving a linear equation system with pure Neumann condition but I cannot find how to fix my problem. Some people suggest that I should impose a Dirichlet condition at one point on the boundary, but when I try this I am getting a spike at this point.

I was using a Gauss-Seidel procedure to solve the system, and at each iteration I would re-set the solution to have an average of 0. However, as this is very slow I moved to using the Bicongugate Gradient Stabilised Method which is much quicker. However, I cannot see a way in which to do something similar in this case, so any help or suggestions would be welcome.

One reason I think this may be the case is that the test problems I am using may be unphysical. For the case where I had b being a sinusoidal forcing with an average of zero over the domain (-2 pi^2 cos(pi x)cos(pi y) from 0-1 in x and y), setting one corner to a Dirichlet condition appeared to work fine. However, for a point source in the centre of the domain, setting a Dirichlet condition at one corner caused a spike there. In my head/reasoning physically a point source in the centre of the domain with Neumann conditions is just going to keep 'pushing' and the solution will float off to infinity with an ever growing arbitrary constant and cannot be anchored. So there is some kind of condition that says I need zero net forcing by b on my closed domain must be zero i.e. the integral of the source term across the domain must be zero. Is this correct?

If it helps, the application of this is solving a Poisson equation for pressure in a fluid dynamics pressure correction scheme. Here I cannot guarantee that b will have an average of zero, though it should hopefully be close.

I am happy to provide more information as necessary.

  • 1
    $\begingroup$ Welcome to SciComp.SE! You missed a question: scicomp.stackexchange.com/questions/19807/…. If this doesn't answer your question (which I think it should), please edit it to describe why. $\endgroup$ Commented Nov 9, 2016 at 10:01
  • $\begingroup$ My apologies, I always try to check for previous cases as it's irritating to have multiple questions, but sometimes miss things. I will have a read and get back to you. $\endgroup$
    – FluidFox
    Commented Nov 9, 2016 at 10:19
  • $\begingroup$ No worries, there are quite a lot of questions on this topic so it's easy to miss one. $\endgroup$ Commented Nov 9, 2016 at 10:27
  • $\begingroup$ I have updated with a more specific question The 'Dirichlet at one cell edge' solution works in some cases and not others, so I think I am maybe missing a mathematical constraint about the source term when I set my problem up. $\endgroup$
    – FluidFox
    Commented Nov 9, 2016 at 10:53
  • $\begingroup$ Yes, the compatibility condition must be satisfied, otherwise the problem you are trying to solve doesn't have a solution, and hence any attempt to calculate one inevitably fails. $\endgroup$ Commented Nov 9, 2016 at 10:58

1 Answer 1


One of the standard ways to handle singular problems like Poisson with Neumann boundary conditions which lead to a sinuglar problem $Ax=b$ is to make the obtained solution orthogonal to the kernel (which is for simple domains only the constant vector) if you use something like a Krylov subspace method.

For instance, you can solve your system with conjugated gradients (CG) and then make the solution orthogonal to the kernel, i.e., subtract the average from your solution. But before solving the linear system make sure that your righthand side vector is also orthogonal to the kernel (again, subtract the average of its values), which is the compatibility condition for singular systems.

The approach with 1-point Dirichlet condition is known to be bad because you spoil the underlying boundary value problem (which is singular, but nice) and one should usually avoid this way.

  • $\begingroup$ Hi, thanks for your reply. Is there somewhere I can read more about this? Yousef Saad's Iterative Methods perhaps? $\endgroup$
    – FluidFox
    Commented Dec 8, 2016 at 11:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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