Solution of an underdetermined system stemming from a PDE with Neumann BC

Consider the Poisson’s equation in 1D with homogeneous b.c.’s $mathrm{d} phi/mathrm{d} x=0$ with the seven point Laplacian (1 -54 783 -1460 783 -54 1 / 576 on a uniform grid). The resulting system is underdetermined because the solution is known only up to an additive constant. It only exists if the compatibility condition is fulfilled (sum of the RHS is 0).

When being solved with LAPACK (dgbtrf and dgbtrs) on the uniform grid an accurate solution is recovered even when the system is left underdetermined.

However, I am having problems with the non-uniform grid (the second derivative in difference form is derived by applying a discrete gradient to field of a discrete gradients in intermediate grid points).

When the grid is perturbed, the very same procedure fails. It either returns NaN or the solution is so large in magnitude that there is a large numerical error in each point. E.g. this system:

$frac{1}{576overline{Delta x}^2} begin{pmatrix}-624.5 &&675 && -50.5
1350 &&-2700 && 1350
-50.5 && 675 && -624.5end{pmatrix} begin{pmatrix}x_1 x_2 x_3end{pmatrix} = begin{pmatrix}0.216 -0.864 0.216end{pmatrix}, overline{Delta x} = 1/3$

When $n=3$ and grid is symmetric, like here, it is enough to add another equation that fixes the solution in one point, e.g., $x_1 = 0$. For that I just multiply $A_{1,1}$ by 2. I get the solution

$frac{1}{576overline{Delta x}^2} begin{pmatrix}-1249 &&675 && -50.5
1350 &&-2700 && 1350
-50.5 && 675 && -624.5end{pmatrix} begin{pmatrix}0 0.02048 0end{pmatrix} = begin{pmatrix}0.216 -0.864 0.216end{pmatrix}, overline{Delta x} = 1/3$

that is also a solution of the original system.

However, when the grid is asymmetric or just larger, the same procedure yields solution that is accurate, except for the value of $x_1$, where the residue is typically 1e-4 or a similar number a few orders more or less.

The RHS is adjusted to zero mean by subtracting the weighted average

$b_0 = b – frac{sum_1^n b(i)Delta x_i}{sum_1^n Delta x_i}$.

What is the proper procedure to solve such a system? Why do I get the non-zero residuum after adding $x_1=0$ to the first row in a non-trivial case?

It is hard to show a complete example due to rounding, but something like

$frac{1}{576overline{Delta x}^2} begin{pmatrix}-535.58897243 && 582.46753247 && -46.87856004
774.43609023 && -1671.42857143 && 896.99248120
-57.92258424 && 895.23809524 && -837.31551100end{pmatrix} begin{pmatrix}x_1 x_2 x_3end{pmatrix} = begin{pmatrix}0.198 -0.852 0.588end{pmatrix}, overline{Delta x} = 1/3$

giving, after the modification,
$x^T = begin{pmatrix}1.0558979027522729E-004 &&2.0058136568986706E-002 &&-2.3505247184514636E-002end{pmatrix}^T$
with residuum
$begin{pmatrix}8.8363636363628095E-004&& 0 &&0end{pmatrix}^T$.

Maybe the non-zero $x_1$ is a que?