# Solving the Poisson equation with Neumann Boundary Conditions - Finite Difference, BiCGSTAB

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.

• 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. Nov 9, 2016 at 10:01
• 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. Nov 9, 2016 at 10:19
• No worries, there are quite a lot of questions on this topic so it's easy to miss one. Nov 9, 2016 at 10:27
• 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. Nov 9, 2016 at 10:53
• 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. Nov 9, 2016 at 10:58

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.