I've recently posted this question on Computational Science Stack Exchange in hopes of getting some insights on how to deal with ill-conditioned 2D Poisson partial differential equation (PDE) which arises when the boundaries imposed on the edges of the computational region are non-zero finite Neumann boundary conditions. The specifics of the problem revolve around the successive over-relaxation method (SOR) which is used as a proffered numerical solver. Unfortunately, there I received mostly generic answers and comments, or comments which mostly apply to a different type of solvers. It is for this reason I am reaching out to the Mathematics Stack Exchange community for some additional technical details.
Majority of the implementation details (code+discretization scheme) you can find in the link to my original question I provided above. I will here summarize briefly the setup and the problem statement:
The equation I am trying to solve (in 2D) is Poisson PDE $$\nabla(\epsilon\nabla\varphi)=\rho,$$ which is generally supposed to be solved in arbitrary 2D geometry. The discretization and code implementation is based on the finite differences. The PDE should be able to be solved for any combination of the Neumann's and the Dirichlet's boundary conditions, including a case where only Neumann's or only Dirichlet's boundaries exist. Out of the two options, the situation where I only have Neumann boundaries (which is possible in a lot of physical systems) is a problematic one, as it leads to a situation where there is no unique solution for the PDE, being that $\varphi+const$ solves it, and $const$ can be anything. Typically, the Neumann boundaries are defined as: $$\frac{\partial\varphi}{\partial x}(x=X_{0,1}, y)\to G_{0,1}$$ and $$\frac{\partial\varphi}{\partial y}(x, y=Y_{0,1})=H_{0,1},$$where $X_0, X_1, Y_0, Y_1$ are initial and final point along the numerically relevant x-axis and y-axis, respectively, and $G_0, G_1, H_0, H_1$ are the boundary conditions along x and y axis at those places. In general, $G_{0,1} = G_{0,1}(y)$ and $G_{0,1} = G_{0,1}(x)$, but for the sake of argument here, let's keep them constant.
Due to the way how I discretized the equation (in the link to the original post), it turns out that if $G_0=-G_1$ and $H_0=-H_1$ the SOR solver is not able to converge, and it starts stalling until it reaches the maximal allowed number of iterations. Under normal operation, it should converge for significantly lower number of iterations than the maximally allowed, but for this particular scenario this is not the case.
I am lead to believe that the reason for the convergence stall is the fact that the Neumann boundary conditions, which I store in the right-hand side (RHS) of the equation together with the source term, $\rho$, invalidate the requirement that the mean of the RHS should be zero.
For instance, if $G_0=G_1$ and $H_0=H_1$, stall does not arise (at least I have not seen it in the examples I've tested, it might as well be still present and hidden).
It is this problem of the convergence stall I have and I currently do not know how to solve, so I am looking here a possible answer to resolve the issue.
Up to now, I have seen suggestions to replace one of the Neumann boundaries with a Dirichlet one, which simply is not possible for some problems which demand all-Neumann conditions. Another thing was Tikhonov regularization in the form of small positive constant multiplying the $\varphi$ on the left-hand side (LHS) of the Poisson equation as: $$\nabla(\epsilon\nabla\varphi)-\varepsilon\varphi=\rho, \varepsilon>0.$$This seems like an interesting and implementable solution, but I am not sure if there is any special form it takes for SOR. In fact, straight-forwardly adding a small positive value (1e-6 – 1e-3) to the main diagonal term yielded no change in observed behavior.
Thus, any technical insights and fixes you can propose for ill-condition Poisson PDE with pure Neumann boundaries would be rather helpful.
Many thanks in advance!
Best Answer
The issue, as it seems, revolves around the compatibility condition which the right-hand side of the Poisson equation and the Neumann boundary conditions must satisfy. This is also addressed in this post. If this condition is not satisfied, then the problem is ill-posed, and it simply cannot be solved. Thus one needs to redefine the problem in some way to make it either well-posed or ill-conditioned.
As often suggested, adding a single point in which Dirichlet's boundary condition is specified is a way of redefining a system to a well-posed problem.
Interestingly, one may face the same problem with pure periodic (torus) boundary conditions, as the boundary value problem in that case should also satisfy the compatibility condition in order to be solvable.
So to summarize, the problem was not in the solver algorithm, it is in the problem definition.