I've recently started looking into writing a finite-differences-based solver for the Poisson equation of the form $$\nabla\left(\varepsilon\nabla\varphi\right)=\rho$$ in 2D for arbitrary geometries (some additional details can be found here). However, I've soon discovered that upon imposing pure Neumann boundaries the system no longer can be solved in a straight forward manner, as the problem becomes ill-posed. First I thought the failure to solve the system arises perhaps due to some kind of a limitation of the solver I use (successive relaxation method), but from answer here I learned that the compatibility conditions are not satisfied in my problem (in particular the discrete compatibility condition), on top of the fact that the all Neumann boundaries make it an ill-posed problem (Edit: not sure about the nomenclature here, but the lack of satisfied compatibility condition could be classified as an ill-conditioned problem).

Now, how I understood one can impose the compatibility condition is just to modify the right-hand side (RHS) of the Poisson equation (including both the original $\rho$ and the Neumann boundary conditions that are stored in the RHS vector) by subtracting its mean as: $$\rho_{h}=\rho-mean(\rho),$$ where mean is calculated over all points in the vector $\rho$. This should be trivial enough, right? Well, my implementation of this in C++ did not yield any results and the convergence issues continued. In fact, the solvers I tried out with this (SOR, biCGstat, CG, CGsquare) all failed to converge as their respective residuals either diverged or just get stuck at a particular value and looped around it until the iteration reached maximal allowed iteration count.

I have recently encountered this post, where a user complains about a similar issue. He, luckily, managed to resolve it, but did not go too much into the details (he only explained that something was incorrect with the counter for the mean value calculation/subtraction).

Thus I wanted to ask if somebody can explain the practical implementation details of the discrete compatibility condition. What needs to be taken into account when calculating the average of the RHS vector, and what needs to be taken into account when subtracting the said average from the original RHS?


2 Answers 2


The problem $$ \nabla (\varepsilon \nabla u) = \rho $$ with pure Neumann conditions admits a family of solutions that differ by a constant value (i.e. $u + c$ is also a solution, for any constant $c$).

For most numerical methods, after discretizing, you will obtain a linear system $$ L \boldsymbol{u} = \boldsymbol{\rho}, $$ where $L$ is a symmetric matrix with a one-dimensional nullspace spanned by the constants, i.e. $$ L \boldsymbol{1} = \boldsymbol{0}, $$ where $\boldsymbol{1}$ is the vector with 1 in every entry.

In order for the discrete problem to be well posed, you need $\boldsymbol{\rho}$ to be in the range of $L$. Note that for any $\boldsymbol{y} \in \operatorname{Range}(L)$ (i.e. with $\boldsymbol{y} = L \boldsymbol{x}$, we have $$ \boldsymbol{1}^T \boldsymbol{y} = 0, $$ since, by symmetry of $L$, $$ \boldsymbol{1}^T \boldsymbol{y} = \boldsymbol{1}^T L \boldsymbol{x} = \boldsymbol{1}^T L^T \boldsymbol{x} = \boldsymbol{0}^T \boldsymbol{x} = 0. $$

Therefore, the discrete compatibility condition is that the vector $\boldsymbol{1}$ should be orthogonal to your right-hand side.

Assume you are given $\tilde{\boldsymbol{\rho}}$, which may not satisfy this condition. You can set $$ \boldsymbol{\rho} = \tilde{\boldsymbol{\rho}} - \left( \frac{\boldsymbol{1}^T \tilde{\boldsymbol{\rho}}}{n} \right) \boldsymbol{1}, $$ where $n$ is the dimension of the discrete space (i.e. $L \in \mathbb{R}^{n \times n}$). This is what is meant by subtracting off the mean of $\tilde{\boldsymbol{\rho}}$. It is easy to check that this modified vector $\boldsymbol{\rho}$ satisfies the discrete compatibility condition. Then, the linear system $$ L \boldsymbol{u} = \boldsymbol{\rho} $$ will have a solution.

The conjugate gradient method applied to this problem will converge to a solution without performing any special modifications (see this answer).

Note that if you are using preconditioned conjugate gradient, you may need to ensure that the output of the preconditioner remains orthogonal to $\boldsymbol{1}$. This can be achieved by performing an orthogonalization step every preconditioner application (i.e, if $M$ is the preconditioner, first compute $\tilde{\boldsymbol{z}} = M^{-1} \boldsymbol{r}$, and then set $\boldsymbol{z} = \tilde{\boldsymbol{z}} - \left( \frac{\boldsymbol{1}^T \tilde{\boldsymbol{z}}}{n} \right) \boldsymbol{1}$). If you are not using a preconditioner, this is not necessary.

  • $\begingroup$ Your answer only covers the homogeneous case of Neumann BCs. Anyway, thank you for your hint. See my updated answer. $\endgroup$
    – ConvexHull
    Commented Nov 15, 2022 at 19:53
  • $\begingroup$ Will P., so my problem comes exactly here - I implement the same condition you gave in the form of $\rho_h=\rho-mean(\rho)$, where the mean value is a sum over all points in $\rho$ and scaled by the total number of points in $\rho$ matrix. However, the solver then fails to converge. But if I pick such NBCs so that their mean is zero, and apply them to e.g. Laplace PDF, everything works perfectly, I do not even need to subtract any average from the right hand side vector. From there I am convinced I am doing something wrong when imposing this condition, but I have no idea what it is. $\endgroup$
    – Akhaim
    Commented Nov 15, 2022 at 20:46
  • $\begingroup$ @Akhaim, can you check that the matrix that results from your discretization is symmetric positive-semidefinite, and has a one-dimensional nullspace spanned by the constants? If those conditions hold, then the described procedure should work. Perhaps something nonstandard in your discretization is causing an issue with the matrix. $\endgroup$
    – Will P.
    Commented Nov 16, 2022 at 19:43
  • $\begingroup$ @WillP. that could be a good test and I think I can try it out fairly straightforwardly. My discretization scheme is nothing fancy - central differences, applied in a rectangular box with equal grid spacing in $x$ and $y$ directions. For the sake of simplicity, I even set $\varepsilon$ to a constant value, so that it does not contribute to the problem. $\endgroup$
    – Akhaim
    Commented Nov 18, 2022 at 11:39

It is actually not that complicated, try to think less analytically.

Assume you have already given the ill-posed & ill-conditioned system $A u=f$ of the original Poisson Neumann problem ($\varphi\equiv u,~\rho\equiv f, \epsilon =1).$ All you have to do is to introduce one single constraint

\begin{align} \frac{1}{b-a}\int_a^b u \,\mathrm{d}x&=0, \\ \int_a^b u \,\mathrm{d}x&=0, \end{align}

which can be discretisized with the most simple quadrature rule (rectangle rule)

\begin{align} \sum_{i=1}^N h u_i &=0, \\ u_1 +u_2 + \dots + u_N &= 0 \equiv \lambda. \end{align}

Now add this single line to your system

$$\begin{bmatrix} A \\ 1^T \end{bmatrix} \begin{bmatrix} u \end{bmatrix} = \begin{bmatrix} f \\ 0\end{bmatrix}.$$

At this point your system is still irregular because you have more lines than columns. In order to fix this add a "clever" zero using $\lambda$ resulting in

$$\begin{bmatrix} A & 1 \\ 1^T & 0 \end{bmatrix} \begin{bmatrix} u \\ \lambda \end{bmatrix} = \begin{bmatrix} f \\ 0 \end{bmatrix}$$

which is invertible and solvable by most standard solvers.

Edit: Credits to @Will P. which however is not general at all.

The Poisson equation with Neumann BCs is given by

\begin{align} \Delta u(x) &= f(x),\qquad x\in\Omega\\ \frac{\partial u(x)}{\partial\mathbf n} &= g(x),\qquad x\in\partial\Omega. \end{align}

To actually be well posed the divergence theorem has to be fullfilled

\begin{align} \int_\Omega f(x)\,dV = \int_\Omega \operatorname{div}\nabla u(x) \,dV = \oint_{\partial\Omega}\nabla u(x)\cdot \mathbf n(x)\,dS = \oint_{\partial\Omega} g(x)\,dS, \end{align}

which ultimately results in following condition

\begin{align} \int_\Omega f(x)\,dV = \oint_{\partial\Omega} g(x)\,dS. \end{align}

Note that the answer by @Will P. only covers the case $\oint_{\partial\Omega} g(x) = 0$.

  • $\begingroup$ thanks for the answer. Just to check one thing - the constraint you suggest, it is to ensure the uniqueness of the solution, right? As far as I understand, if the discrete compatibility condition is not satisfied, Poisson-Neumann problem cannot be solved at all, so the uniqueness of the solution is not really a problem at that stage. My understanding is that only when compatibility condition is satisfied one should be concerned with the uniqueness of the solution, and even then depending on the initial condition some of the iterative solvers should still converge to a solution. $\endgroup$
    – Akhaim
    Commented Nov 13, 2022 at 21:35
  • $\begingroup$ I am not quite sure about your comment. You may confuse the topics ill-posed and ill-conditioned. $\endgroup$
    – ConvexHull
    Commented Nov 13, 2022 at 22:01
  • $\begingroup$ I am not sure you understood my question properly. Please refer to this link scicomp.stackexchange.com/questions/34493/… in order to understand what is meant by the discrete compatibility condition and how it differs from uniqueness constraint you use in your post. $\endgroup$
    – Akhaim
    Commented Nov 13, 2022 at 23:31
  • $\begingroup$ I am not sure what the link is supposed to tell me. The extra condition for $f$ is not needed (from a numerical point of view) and even wrong for other arbitrary source terms. I removed the term compatibilty condition for $u$ in my answer to avoid confusion. $\endgroup$
    – ConvexHull
    Commented Nov 14, 2022 at 18:34
  • $\begingroup$ Also note that pure Neumann BCs ist not the same as pure periodic BCs. $\endgroup$
    – ConvexHull
    Commented Nov 14, 2022 at 20:29

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.