Explicitly including boundary points in a set of finite-difference equations

Question

Consider the 1D poisson equation

$$ \frac{d^2 u}{dx^2} = -\rho $$

with Dirichlet boundary conditions $u(0) = u(l) = g$. Using a finite difference scheme, with a 5-point grid $u_1,u_2,u_3,u_4,u_5$ (excluding boundary points $u_0$ and $u_l$), we get the set of linear equations

$$ \left( \begin{array}{ccc} 2 & -1 & & & \\ -1 & 2 & -1 & & \\ & -1 & 2 & -1 & \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \\ \end{array} \right)\left( \begin{array}{c} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5\end{array} \right) = \left( \begin{array}{c} \rho_1+g \\ \rho_2 \\ \rho_3 \\ \rho_4 \\ \rho_5+g\end{array} \right) $$

My question is: What would the matrix look like if I made $u_3$ a boundary point too?

Would it look like

$$ \left( \begin{array}{ccc} 2 & -1 & & & \\ -1 & 2 & 0 & & \\ & 0 & 1 & 0 & \\ & & 0 & 2 & -1 \\ & & & -1 & 2 \\ \end{array} \right)\left( \begin{array}{c} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5\end{array} \right) = \left( \begin{array}{c} \rho_1+g \\ \rho_2+g \\ g \\ \rho_4+g \\ \rho_5+g\end{array} \right) $$

I ask because I have a 3d system with small but irregular internal boundary regions, and it would be currently more convenient for my purposes to leave them in the matrix (even if it means extra computational cost).

Answer 1

If you include the boundary condition directly in the matrix, you will only get the g value at the points where the boundary is prescribed. If we use 5 nodes with the following BCS: $$u_1=g_1$$ and $$u_5=g_5$$

Then the matrix resulting from the finite difference will be:

$$ \left( \begin{array}{ccc} 1 & 0 & & & \\ -1 & 2 & -1 & & \\ & -1 & 2 & -1 & \\ & & -1 & 2 & -1 \\ & & & 0 & 1 \\ \end{array} \right)\left( \begin{array}{c} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5\end{array} \right) = \left( \begin{array}{c} g_1 \\ \rho_2 \\ \rho_3 \\ \rho_4 \\ g_5\end{array} \right) $$

That is, the value of $$u_1$$ and $$u_5$$ follow the Dirichlet BCs and within the domain, the Poisson equation applies.

If you have Neumann boundary conditions or Robin Boundary condition, than line 1 and 5 will have to change to respect the application of these boundary conditions.

If the third point is a boundary point and following my example, then the resulting matrix should be :

$$ \left( \begin{array}{ccc} 1 & 0 & & & \\ -1 & 2 & -1 & & \\ & 0 & 1 & 0 & \\ & & -1 & 2 & -1 \\ & & & 0 & 1 \\ \end{array} \right)\left( \begin{array}{c} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5\end{array} \right) = \left( \begin{array}{c} g_1 \\ \rho_2 \\ g_3 \\ \rho_4 \\ g_5\end{array} \right) $$ But this is impossible in our case, since this would imply prescribing 3 conditions on a problem which requires 2. In 2D, however, you will get lines with strictly diagonal term that appear within the matrix due to Dirichlet BCs.

Answer 2

Another option to the answer @BlaisB is to incorporate the constraints via Lagrange multipliers, which basically leads to a linear system which is augmented by two variables,

$$ \left( \begin{array}{ccc} 2 & -1 & & & & 1 &\\ -1 & 2 & -1 & & \\ & -1 & 2 & -1 & \\ & & -1 & 2 & -1 \\ & & & -1 & 2 & & 1 \\ 1 & & & & & \\ & & & & 1 &\\ \end{array} \right)\left( \begin{array}{c} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5 \\ \lambda_1 \\ \lambda_2 \\ \end{array} \right) = \left( \begin{array}{c} \rho_1 \\ \rho_2 \\ \rho_3 \\ \rho_4 \\ \rho_5 \\ g_1\\ g_2 \end{array} \right) $$

The found Lagrange multipliers $\lambda_i$ itself are not important and can be discarded after solution. Their only purpose is to keep up the boundary condition.

The advantage of this method is that you don't need to bother with adjustments of the original system. Moreover, it's straightforward to incorporate also more sophisticated boundary conditions, e.g. of Neumann-type (--try that with the explicit method!)