Why naively chopped finite difference matrix works for different ODE boundary conditions

Question

We know finite difference method (FDM) can replace $y''(x)$ as $\frac{1}{h^2}[y(x+h)+y(x-h)-2y(x)]$ or so. One naive way to write down the matrix of the differential operator is like the following, which somehow simply cuts off at the two edges. Similarly, we can write the matrix for $y'=\frac{1}{2h}[y(x+h)-y(x-h)]$ as

$$\begin{bmatrix} 2t_0 &-t_0 & & & \\ -t_0 &2t_0 &-t_0 & & \\ &\ddots &\ddots &\ddots & \\ & & -t_0 &2t_0 &-t_0\\ & & & -t_0 &2t_0 \end{bmatrix}$$

But what boundary condition (BC) correspond to this?

I first thought this means zero-value Dirichlet BC (like $y(a)=y(b)=0$). However, it turns out not to be the case. It automatically handles either zero or nonzero boundary values as shown in the solution below.

Based on such naive FDM matrics, I solve the eigenvalue ($\lambda$) problem of this linear ODE $y''+\frac{2}{x}y'+[2\lambda-x^2-\frac{l(l+1)}{x^2}]y=0$ with nonnegative integer $l$.
It can be analytically solved and we know at the boundary, $y(0)\neq0$ for $l=0$, $y(0)=0$ for $l>0$, and $y(\infty)=0$. Surely we use some large enough interval to mimic the infinity. The following code solves the two cases well. But again, why? Is this chop scheme generally correct for fixed (arbitrary) valued Dirichlet boundaries or anything else?

---

Following is my simple Mathematica code.

l = 0; a = 1; n = 1001; h = 16/(n - 1);
M1 = -2.0 IdentityMatrix[n] + DiagonalMatrix[Table[1, {n - 1}], 1] + 
  DiagonalMatrix[Table[1, {n - 1}], -1];(*//MatrixForm*)
M2 = DiagonalMatrix[Table[-1/(i + 1), {i, 1, n - 1}], -1] + 
  DiagonalMatrix[Table[1/i, {i, 1, n - 1}], 1];
M3 = DiagonalMatrix[Table[(i)^2, {i, 1, n}]];(*//MatrixForm*)
M4 = l (l + 1) DiagonalMatrix[Table[(i)^-2, {i, 1, n}]];
M = 1/h^2 M1 + h^-2 M2 - h^2 M3 - h^-2 M4;
qq = Eigenvectors[M];
ListPlot[qq[[n - a + 1]], PlotRange -> All]

$a$ means the $a$th smallest eigenvalue and certainly one can change the value of $l$ and $a$.

---

The exact solution looks like $$ y_{0,l}(x)\propto x^l e^{-\frac{1}{2}x^2}\\ y_{1,l}(x)\propto x^l e^{-\frac{1}{2}x^2}(\frac{2l+3}{2}-x^2)\\ y_{2,l}(x)\propto x^l e^{-\frac{1}{2}x^2}(\cdots)\\ \cdots $$

Answer 1

If you convert the first line of your matrix operator back into a difference equation you get $$2y(x)-y(x+h)=h^2b,$$ where $b$ is whatever is on the right hand side.

This is consistent with a non-centred finite difference approximation to $$ y - h y' = h^2 b.$$

This is an example of a Robin, or third type boundary condition, although a slightly unusual one, since it contains your mesh length inside it. Note that in the limit of $h\rightarrow 0$, this does tend towards a zero Dirichlet condition, so you would expect your scheme to tend towards that solution as your mesh converges, providing the gradient of your solution at the boundary stays well behaved.

Answer 2

When first discretizing equations, it's useful to draw a diagram showing where the data lives relative to physical boundaries. Here are some schematics for visual aid:

Applying finite difference to a differential equation takes a different form for cell-centered (CC) and node (N) data.

--------------------------------------------------------------------------------------------------------------------------------

Node data Consider N data.

Consider the Laplace stencil operating on data $u$ with $f$ on the right-hand side (RHS),

\begin{equation} \left(\frac{\partial^2 u}{\partial x^2}\right)_{boundary} = \frac{u_g - 2 u_b + u_i}{\Delta x^2} = f_b \end{equation}

Where $u_g,u_b,u_i$ are the ghost, boundary and first interior node points. This Laplacian operator may be written in matrix form as

$ A = \frac{1}{\Delta x^2} \left[\begin{array}{ccccccccc} 0 & 0 & 0 & & & & & & 0 \\ 1 & -2 & 1 & & & & & & \\ 0 & 1 & -2 & 1 & & & & & \\ & & & \ddots & \ddots & \ddots & & & \\ & & & & & 1 & -2 & 1 & 0 \\ & & & & & & 1 & -2 & 1 \\ 0 & & & & & & 0 & 0 & 0 \\ \end{array} \right] $

Note that this stencil reaches to the ghost points.

Dirichlet BCs

Consider Dirichlet BCs, $u_b$ is known, and the solve is only from $u_i$ onward. Therefore, the only equation we must consider is at location $i$:

\begin{equation} \left(\frac{\partial^2 u}{\partial x^2}\right)_{i} = \frac{u_b - 2 u_i + u_{i+1}}{\Delta x^2} = f_i \end{equation}

If we were to write an equation for the boundary point, we could write \begin{equation} \left(\frac{\partial^2 u}{\partial x^2}\right)_{b} = \frac{u_g - 2 u_b + u_{i}}{\Delta x^2} = f_b \rightarrow u_g = 2 u_b - u_{i} - \Delta x^2 f_b \end{equation} If we insist $f_b=0$ then we have a simplified version: \begin{equation} u_g = 2 u_b - u_{i} \end{equation} So the equation changes to \begin{equation} \left(\frac{\partial^2 u}{\partial x^2}\right)_{b} = 0 = 0 \end{equation} And we may remove it from our system. Looking back the first interior point, let's move boundary value to the RHS: \begin{equation} \left(\frac{\partial^2 u}{\partial x^2}\right)_{i} = \frac{u_b - 2 u_i + u_{i+1}}{\Delta x^2} = f_i \rightarrow \frac{- 2 u_i + u_{i+1}}{\Delta x^2} = f_i - \frac{u_b}{\Delta x^2} \end{equation}

This means that the ghost point should not enter the computations at all. They are a means to an end to apply the desired BCs. Correspondingly, the matrix $A$ is: $ A = \frac{1}{\Delta x^2} \left[\begin{array}{ccccccccc} 0 & 0 & 0 & & & & & & 0 \\ 1 & -2 & 1 & & & & & & \\ 0 & 1 & -2 & 1 & & & & & \\ & & & \ddots & \ddots & \ddots & & & \\ & & & & & 1 & -2 & 1 & 0 \\ & & & & & & 1 & -2 & 1 \\ 0 & & & & & & 0 & 0 & 0 \\ \end{array} \right] \\ \rightarrow A = \frac{1}{\Delta x^2} \left[\begin{array}{ccccccccc} 0 & 0 & 0 & & & & & & 0 \\ 0 & 0 & 0 & & & & & & \\ 0 & \textbf{0} & -2 & 1 & & & & & \\ 0 & 0 & 1 & -2 & & & & & \\ & & & \ddots & \ddots & \ddots & & & \\ & & & & & -2 & 1 & 0 & 0 \\ & & & & & 1 & -2 & \textbf{0} & 0 \\ 0 & & & & & & 0 & 0 & 0 \\ 0 & & & & & & 0 & 0 & 0 \\ \end{array} \right] $

This is nice because the equation has effectively become smaller. NOTE: the first two 0's in the above equation refer to the ghost and boundary equations.

Result

As you can see, this "truncation" of the first and last column in $A$ will effectively apply Dirichlet BCs to our system, so long $u_b$ is in fact zero, so this only works for a special case.

--------------------------------------------------------------------------------------------------------------------------------

Cell-Centered data

Consider the Laplace stencil operating on data $u$ with $f$ on the right-hand side (RHS),

\begin{equation} \left(\frac{\partial^2 u}{\partial x^2}\right)_{1} = \frac{u_g - 2 u_1 + u_{2}}{\Delta x^2} = f_1 \end{equation} Where $u_g,u_1,u_2$ are the ghost, first interior and second interior cells. This Laplacian operator may be written in matrix form as

$ A = \frac{1}{\Delta x^2} \left[\begin{array}{ccccccccc} 1 & 1 & & & & & & & 0 \\ 1 & -2 & 1 & & & & & & \\ & 1 & -2 & 1 & & & & & \\ & & & \ddots & \ddots & \ddots & & & \\ & & & & & 1 & -2 & 1 & \\ & & & & & & 1 & -2 & 1 \\ 0 & & & & & & & 1 & 1 \\ \end{array} \right] $

Let's substitute $u_g$ and adjust the stencil accordingly.

Dirichlet

Consider Dirichlet BCs, $u_g$ may be computed observing that the boundary value is the average of the neighboring two cell center values. \begin{equation} u_b = \frac{u_g + u_1}{2} \rightarrow u_g = 2u_b - u_1 \end{equation} So our Laplacian stencil becomes: \begin{equation} \left(\frac{\partial^2 u}{\partial x^2}\right)_{1} = \frac{(2 u_b - u_1) - 2 u_1 + u_{2}}{\Delta x^2} = f_1 \end{equation} This boundary point must be moved to the RHS to maintain a consistent matrix-vector multiplication. Therefore our equation changes: \begin{equation} \frac{u_g - 2 u_1 + u_{2}}{\Delta x^2} = f_1 \rightarrow \frac{- 3 u_1 + u_{2}}{\Delta x^2} = f_1 - \frac{2 u_b}{\Delta x^2} \end{equation} Correspondingly, the matrix $A$ changes: $ A = \frac{1}{\Delta x^2} \left[\begin{array}{ccccccccc} 1 & 1 & & & & & & & 0 \\ 1 & -2 & 1 & & & & & & \\ & 1 & -2 & 1 & & & & & \\ & & & \ddots & \ddots & \ddots & & & \\ & & & & & 1 & -2 & 1 & \\ & & & & & & 1 & -2 & 1 \\ 0 & & & & & & & 1 & 1 \\ \end{array} \right] \\ \rightarrow A = \frac{1}{\Delta x^2} \left[\begin{array}{ccccccccc} 0 & 0 & 0 & & & & & & 0 \\ 0 & -3 & 1 & & & & & & \\ 0 & 1 & -2 & 1 & & & & & \\ & & & \ddots & \ddots & \ddots & & & \\ & & & & & 1 & -2 & 1 & 0 \\ & & & & & & 1 & -3 & 0 \\ 0 & & & & & & 0 & 0 & 0 \\ \end{array} \right] $

Note that the first and last equations are identities ($0=0$) which are reserved for the ghost points.

Result

So, I stand corrected with my comment (on the OP's question) about assuming cell centered data. This "truncation" only works in a special case for Node data.