# 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}{576\overline{\Delta x}^2} \begin{pmatrix}-624.5 &&675 && -50.5 \\ 1350 &&-2700 && 1350 \\ -50.5 && 675 && -624.5\end{pmatrix} \begin{pmatrix}x_1 \\ x_2 \\ x_3\end{pmatrix} = \begin{pmatrix}0.216 \\ -0.864 \\ 0.216\end{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}{576\overline{\Delta x}^2} \begin{pmatrix}-1249 &&675 && -50.5 \\ 1350 &&-2700 && 1350 \\ -50.5 && 675 && -624.5\end{pmatrix} \begin{pmatrix}0 \\ 0.02048 \\ 0\end{pmatrix} = \begin{pmatrix}0.216 \\ -0.864 \\ 0.216\end{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}{576\overline{\Delta x}^2} \begin{pmatrix}-535.58897243 && 582.46753247 && -46.87856004 \\ 774.43609023 && -1671.42857143 && 896.99248120 \\ -57.92258424 && 895.23809524 && -837.31551100\end{pmatrix} \begin{pmatrix}x_1 \\ x_2 \\ x_3\end{pmatrix} = \begin{pmatrix}0.198 \\ -0.852 \\ 0.588\end{pmatrix}, \overline{\Delta x} = 1/3$$

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

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

• Trying to understand this statement "with LAPACK (dgbtrf and dgbtrs) on the uniform grid an accurate solution is recovered anyway" - does it mean that the compatibility condition is not satisfied, so it is an ill-posed problem in this case but the solver does not show it? Apr 30 at 0:28
• @MaximUmansky No, the compatibility condition is satisfied, but the system remains underdetermined (any solution plus a constant is also a solution) but LAPACK returns an accurate solution. Apr 30 at 8:12