The Problem
Suppose I have an equation of the form $\nabla^2 \phi(x) = f(x)$ on the interval $A \le x \le B$, where $f(x)$ is known and $\phi(x)$ is unknown. I have Neumann-type boundary conditions: $\frac{\partial \phi}{\partial x}\big|_{x=A} = g_A$ and $\frac{\partial \phi}{\partial x}\big|_{x=B} = g_B$, where $g_A$ and $g_B$ are known.
The compatibility constraint (derived via the divergence theorem) requires that $\int_A^B f(x) dx = g_B-g_A$. Let's assume for this problem that this is satisfied exactly such that a solution is possible.
Ultimately, I don't care about solving for $\phi$, but rather $\vec{E} = -\vec{\nabla} \phi$, but I want to use $\phi$ as the electrostatic potential because it enforces $\vec{\nabla} \times \vec{E} = 0$.
We know and expect that there is not a unique solution, but a family of solutions exist where $\phi(x) = \Phi(x) + C$, where C is any number. In other words, $\phi(x)$ is unique up to an additive constant. I don't care what this constant $C$ is, but it would be nice to be able to constrain it, and I think it is required for a numerical solution to be possible.
Finite Difference Approach
I could set up a 2nd-order finite difference approach over $N_x$ points as follows, where $f_i \equiv f(x_i)$ and $\phi_i \equiv \phi(x_i)$:
$x_i = a + \left(\frac{b-a}{N_x-1}\right) (i-1) = a + \delta x (i-1) \qquad \text{for $i=1, 2, …,N_x$}$,
where $\delta x \equiv \frac{b-a}{N_x-1}$.
The governing equations become $\frac{\phi_{i+1} – 2\phi_{i} + \phi_{i-1}}{\delta x^2} = f_i$.
The two boundary conditions become $\frac{\phi_2 – \phi_{0}}{\delta x} = g_A$ and $\frac{\phi_{N_x+1} – \phi_{N_x-1}}{\delta x} = g_B$.
This brings us to my issue. The PDE does not have a unique solution!
I can use ghost points ($x_0$ and $x_{N_x+1}$) and combine each boundary condition with the governing equation at each boundary. Doing so gives me $N_x$ equations and satisfies the boundary conditions. I can create a $N_x \times N_x$ matrix to solve this problem, but it is singular because $\phi(x)$ is not unique. I can add a uniqueness constraint to the matrix ($\phi_C = C$) to make a $(N_x+1) \times N_x$ matrix, but I'm afraid this will lead to a minimum-length solution.
Concrete Example
For $N_x=7$, I have the following system of equations:
$
\begin{align}
\phi_{2} – 2\phi_{1} + \phi_{0} \;=&\; \delta x^2 \; f_{1} \quad \text{and} \quad \phi_{2} – \phi_{0} \;=\; \delta x\; g_A \rightarrow 2 \phi_{2} – 2 \phi_{1} = \delta x g_A + \delta x^2 \; f_{1} \\
\phi_{3} – 2\phi_{2} + \phi_{1} \;=&\; \delta x^2 \; f_{2} \\
\phi_{4} – 2\phi_{3} + \phi_{2} \;=&\; \delta x^2 \; f_{3} \\
\phi_{5} – 2\phi_{4} + \phi_{3} \;=&\; \delta x^2 \; f_{4} \\
\phi_{6} – 2\phi_{5} + \phi_{4} \;=&\; \delta x^2 \; f_{5} \\
\phi_{7} – 2\phi_{6} + \phi_{5} \;=&\; \delta x^2 \; f_{6} \\
\phi_{8} – 2\phi_{7} + \phi_{6} \;=&\; \delta x^2 \; f_{7} = \quad \text{and} \quad \phi_{8} – \phi_{6} \;=\; \delta x\; g_B \rightarrow 2 \phi_{6} – 2 \phi_{7} = -\delta x \; g_B + \delta x^2 \; f_{7} \\
\end{align}
$
$
\begin{bmatrix}
-2 & 2 & 0 & 0 & 0 & 0 & 0 \\
1 & -2 & 1 & 0 & 0 & 0 & 0 \\
0 & 1 & -2 & 1 & 0 & 0 & 0 \\
0 & 0 & 1 & -2 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 & -2 & 1 & 0 \\
0 & 0 & 0 & 0 & 1 & -2 & 1 \\
0 & 0 & 0 & 0 & 0 & 2 & -2
\end{bmatrix}
\begin{bmatrix}
\phi_{1} \\
\phi_{2} \\
\phi_{3} \\
\phi_{4} \\
\phi_{5} \\
\phi_{6} \\
\phi_{7}
\end{bmatrix}
=
\begin{bmatrix}
\delta x^2 \; f_{1} + \delta x \; g_A \\
\delta x^2 \; f_{2} \\
\delta x^2 \; f_{3} \\
\delta x^2 \; f_{4} \\
\delta x^2 \; f_{5} \\
\delta x^2 \; f_{6} \\
\delta x^2 \; f_{7} – \delta x; g_B
\end{bmatrix}
$
The problem with the above matrix equation is that it is singular. This is to be expected because $\phi$ is non-unique and is only defined up to a scalar constant.
We could specify this constant by setting the value of $\phi$ at any given point. In the above example, we could try for a unique solution by adding the constraint $\phi_4 = C$ (where $C$ is a given finite number):
$
\begin{bmatrix}
-2 & 2 & 0 & 0 & 0 & 0 & 0 \\
1 & -2 & 1 & 0 & 0 & 0 & 0 \\
0 & 1 & -2 & 1 & 0 & 0 & 0 \\
0 & 0 & 1 & -2 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 & -2 & 1 & 0 \\
0 & 0 & 0 & 0 & 1 & -2 & 1 \\
0 & 0 & 0 & 0 & 0 & 2 & -2 \\
0 & 0 & 0 & 1 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
\phi_{1} \\
\phi_{2} \\
\phi_{3} \\
\phi_{4} \\
\phi_{5} \\
\phi_{6} \\
\phi_{7}
\end{bmatrix}
=
\begin{bmatrix}
\delta x^2 \; f_{1} + \delta x \; g_A \\
\delta x^2 \; f_{2} \\
\delta x^2 \; f_{3} \\
\delta x^2 \; f_{4} \\
\delta x^2 \; f_{5} \\
\delta x^2 \; f_{6} \\
\delta x^2 \; f_{7} – \delta x \; g_B \\
C
\end{bmatrix}
$
The problem with this approach is that it is no longer a square matrix and looks like it's an over-determined system. I'm concerned that matrix-solving algorithms for non-square matrices would treat this as a least-squares problem. This isn't an over-determined system, however: the rank of the matrix equals the number of unknowns, even though the matrix isn't square.
Alternatively, I could remove one of the equations to make the matrix square.
If I were to take out the uniqueness constraint ($\phi_4 = C$), then the matrix is square but singular. (MATLAB can solve this anyway, and it gets an answer with the least error, but I don't know how it does it, and I'm concerned about it not specifying $C$.)
If I were to remove one of the governing equations ($\frac{\phi_{i+1} – 2\phi_{i} + \phi_{i-1}}{\delta x^2} = f_i$), then the system does not account for the curvature at $x=x_i$, and the solution is inaccurate.
My Question
Is there a good algorithm for solving this type of problem, where I need more constraints than the number of unknowns, but the problem is not overdetermined?
In particular:
- Is there an algorithm that can solve the underdetermined case where $\phi_i = C$ is not included? (Is it absolutely necessary to include this constraint, or an equivalent one?)
- Is there a system of equation that is solved reliably, with minimal numerical error? I've tried rearranging the above equations and using different matrix solvers, but most of the matrices I've used give condition numbers that scale with $N_x^2$, which is concerning for the problem sizes I'm using for my research.)
Best Answer
A standard way how to make the solution unique is to look for the zero-mean solution $$ \int_A^B\phi(x)\,\mathrm{dx}=0. $$ Let $Mu=f$ be the linear system corresponding to your finite difference approximation (with uniform discretization, constant vectors form nullspace of $M$) and $e=[1,\ldots,1]^T$. Since $M$ is symmetric and semi-definite, any solution of $Mu=f$ also minimizes $\frac{1}{2}u^TMu-u^Tf$. Imposing the zero-mean condition gives the constrained minimization problem $$ \min_u \frac{1}{2}u^TMu-u^Tf\quad\text{such that}\quad e^Tu=0 $$ which by introducing a Lagrange multiplier $c$ leads to a system $$ \pmatrix{M&e\\e^T&0}\pmatrix{u\\c}=\pmatrix{f\\0}. $$