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

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$$

• Would you mind adding the exact solution to the eigenproblem on $[0,\infty]$ to your question? – GoHokies Jan 16 '18 at 7:54
• @GoHokies Added. – xiaohuamao Jan 16 '18 at 8:15
• Can you please specify whether you are referring to cell centered data, or cell cornered (node) data. This is important. It seems that you are considering cell centered data. – Charles Jan 25 '18 at 5:40
• @Charlie Sorry, what does that mean? Why is it cell centered in my case? Thank you. – xiaohuamao Jan 25 '18 at 5:42
• @origimbo has the right idea with writing down the equation for the boundary and seeing exactly what equation holds on the boundary. First, I think, you need to answer the question “where are the discrete unknowns relative to the boundary?” (Cell centered/node centered). – Charles Jan 25 '18 at 5:43

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.

• But as I increase the mesh density, there's completely no sign of approaching a zero Dirichlet condition for the ODE I showed. It's a very stable nonzero boundary value as expected from analytical results. Anything wrong? – xiaohuamao Jan 11 '18 at 14:12
• And your claim of Robin type seems to contradict page 3 of this document. Their claim of Dirichlet is not in the limit of vanishing $h$. – xiaohuamao Jan 11 '18 at 15:04
• That depends on whether your first data point is on the boundary, or one grid point inside. In the latter case, then the boundary value at the boundary ends up on the right hand side, as per their eqn (2), scaled by $1/h^2$, which is akin to replacing the condition at the boundary with $y=y_L$ etc. – origimbo Jan 11 '18 at 15:39

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.