# Solving a linear equation system with pure Neumann condition

I am trying to solve a linear equation system $\textbf{A}\textbf{x}=\textbf{b}$, e.g. a Poisson equation discretized in strong form, using biCGstab method.

Since there are only natural Neumann boundary conditions, at least one Dirichlet condition must be given. I simply add condition $\sum_{j} x=0$ to row $i$ in $\textbf{A}$. The result converges, but the error in the gradient field $\nabla x$ at $i$ is very large.

What is the reason for this, and how should I imposition this condition?

• Look at Bochev and Lehoucq, On the Finite Element Solution of the Pure Neumann Problem, SIAM Review 2005, epubs.siam.org/doi/abs/10.1137/S0036144503426074. Jun 1 '15 at 14:05
• Depending on your preconditioner, you may get away with not specifying such a condition. Krylov solvers will, in some sense, pick a pressure for you. Jun 1 '15 at 14:31

Pure Neumann problem is unique up to a constant. My two favourite solution strategies:

1. Modifying the equation $-\Delta u = f$ to $-\Delta u + \varepsilon u = f$ for some small $\varepsilon>0$. If you perform reduced integration this corresponds roughly to adding $\varepsilon$ to the diagonal of your system matrix.

2. Imposing $\int_\Omega u \,\mathrm{d}x=0$ using Lagrange multiplier. Then your system would be $$\begin{bmatrix} A & 1 \\ 1^T & 0 \end{bmatrix} \begin{bmatrix} x \\ \lambda \end{bmatrix} = \begin{bmatrix} f \\ 0 \end{bmatrix}.$$

• Just a note: BiCGstab won't care, only CG requires symmetry and positive definiteness. (On the other hand, CG will work for the unmodified system, as long as the (discrete) system is consistent -- at least in exact arithmetic, and up to a constant.) Jan 28 '16 at 19:33
• Thanks for the clarification. I made a faulty assumption that this property is shared by CG and BiCGstab although I don't know the latter method well.
– knl
Jan 28 '16 at 20:12
• @knl I tried to follow your approach in my last question, however it's not apparently working Oct 3 '21 at 10:52

Including even one Dirichelt condition changes the problem you are trying to solve and will not give you the correct solution! You must fulfil the Discrete Compatibility Criteria, see e.g. first and second pages of: http://eprints.ma.man.ac.uk/894/2/covered/MIMS_ep2007_156_Sample_Chapter.pdf

It basically states that the summation of each element of b must be equal to zero for all-Neumann problems, otherwise your solution will drift as there is nothing holding it at the boundaries.

You simply need to add a line stating that b = b - mean(b) in your code before solving, if you are solving in double precision. In single precision it might also be necessary to ensure, at every other iteration, that the residual in biCGstab also meets the Discrete Compatibility Criteria.

I attempted to solve an all-Neumann problem using CG by imposing a Dirichlet condition at a single point, as one is usually told to do anecdotally by e.g. some Prof., and then compared it to using a Discrete Cosine Transform approach. It does not yield the same result, using the Discrete Compatibility Criteria instead does.

• Could you clarify what you mean by: "In single precision it might also be necessary to ensure, at every other iteration, that the residual in biCGstab also meets the Discrete Compatibility Criteria."? May 16 '21 at 15:06

I remember that the matrices from the Poisson equations with pure Neumann boundary conditions is singular. Is it right for your case ? If so, you cannot use the BiCGSTAB to solve the resulting linear systems.

• Yes, it is singular, and that is why I have to add one Dirichlet bc. Because I only care about the gradient field. Maybe I can use QR solvers, but they are too slow. Jun 2 '15 at 4:59
• Well, if the matrix is singular, maybe I suggest you use some regularization methods, such as Tikhonov regularization. Jun 2 '15 at 7:20
• The resulting matrix will still be singular, even with the additional condition added. Jun 3 '15 at 7:44
• @shuhalo but the new system converged fast, does that happen to a singular matrix? Jun 3 '15 at 13:48
• I case someone still stumbles upon this answer, it is incorrect and Conjugate Gradients works for pure Neumann problems, see e.g. cs.sandia.gov/~pbboche/papers_pdf/2005SIREV.pdf Jan 21 '17 at 23:54