# Why ODE's naive finite difference matrix works well for different 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. The 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 for $y'$.

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, 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 is analytically solved in textbooks 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. It should work well for $l>0$, but why also for $l=0$ case???

Two things I noticed but not sure if relevant or not:
1. If we really somehow impose a nonzero left boundary value $y(0)$, the value you use doesn't matter for any eigenvalue problem. But we intuitively think it at least should be nonzero.
2. For the strange $l=0$ case, we have $y'(0)=0$ in the analytic solution.

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]


M2 is for this $\frac{2}{x}y'=\frac{2}{ih}\frac{y(x_{i+1})-y(x_{i-1})}{2h}=\frac{y(x_{i+1})-y(x_{i-1})}{i\,h^{2}}$. $a$ means the $a$th smallest eigenvalue and certainly one can change the value of $l$ and $a$.

• Seems to be more of a mathematics question than one about how Mathematica works. – Michael E2 Jan 11 '18 at 13:55
• @xzczd No, it's correct. M2 is for the $\frac{2}{x}y'$ term in the ODE. Then no need for the $\frac{1}{2}$ for the other three matrices since overall factor doesn't matter. In any case, it reproduces the analytic results in some textbook perfectly. – xiaohuamao Jan 12 '18 at 6:06
• @xzczd M2 is for this $\frac{2}{x}y'=\frac{2}{ih}\frac{y(x_{i+1})-y(x_{i-1})}{2h}=\frac{y(x_{i+1})-y(x_{i-1})}{i\,h^{2}}$. $2\lambda$ or $\lambda$ doesn't matter. – xiaohuamao Jan 12 '18 at 6:30
• Sorry, I made a simple mistake… So your questoin is, you think the code should not reproduce the correct result when $l=0$ if the matrix represents zero b.c.? – xzczd Jan 12 '18 at 6:30
• @xzczd Exactly! – xiaohuamao Jan 12 '18 at 6:31

Not a complete answer, but the following partly explains the 2 observations you made: The equation has a removable singularity at $x=0$, which leads to different implicit b.c.s depending on whether $l$ is $0$ or not:

Clear@l;
eq = y''[x] + 2/x y'[x] + (2 λ - x^2 - (l (l + 1))/x^2) y[x] == 0;
eq /. x -> 0
(* Indeterminate == 0 *)
neweq = x^2 # & /@ eq // Simplify;
neweq /. x -> 0
(* (l + l^2) y[0] == 0 *)

l = 0;
eq /. x -> 0
(* False *)
neweq2 = x # & /@ eq // Simplify;
neweq2 /. x -> 0
(* 0 == 2 Derivative[1][y][0] *)
Clear@l


Sadly I can't figure out why the matrix built in your way respects the implicit b.c. and isn't influenced by the wrong b.c..

BTW the following is my approach for calculating eigenvector with FDM, notice pdetoae is used for discretization:

domain = {bL, bR} = {eps, 16};
lhs = λ y[x] /. First@Solve[neweq2, \[Lambda]]
bcL = y'[bL] == 0;
bcR = y[bR] == 0;

points = 300;
grid = Array[# &, points, domain];
(*Definition of pdetoae is not included in this post,please find it in the link above.*)

ptoafunc = pdetoae[y[x], grid, 2];
ae = ptoafunc[lhs];
aebcL = Solve[ptoafunc@bcL, y[bL]][[1]]
aebcR = Solve[ptoafunc@bcR, y[bR]][[1]]
delete = #[[2 ;; -2]] &;
aewithbc = delete@ae /. aebcL /. aebcR;

{b, mat} = CoefficientArrays @@ ({aewithbc, y /@ grid // delete} /. eps -> 0);
ListLinePlot[Eigenvectors[N[mat], -1], PlotRange -> All]


I think this approach is a bit easier to understand.

• This sounds to make sense. My question is: for what reason implicit B.C. prevails over the wrong B.C. imposed by the matrix form? Why is this contradiction just insignificant? BTW, to analyse behavior at $x\rightarrow 0$, is it correct to say $y'=0$ for something like $y'+\frac{l}{x}y+f(x)y=\lambda y$ when $l=0$ and $f(0)$ is finite? – xiaohuamao Jan 12 '18 at 8:28
• Sadly I can't think out an explanation for the first question at the moment. As to the last question, the answer is no. Just repeat our analysis above: y'[x] + l/x y[x] + f[x] y[x] == \[Lambda] y[x] /. l -> 0/.x->0 gives f[0] y[0] + Derivative[1][y][0] == \[Lambda] y[0] so there's no implicit b.c. in this case. A simple counter example: test = y'[x] + l/x y[x] + f[x] y[x] == \[Lambda] y[x] /. l -> 0 /. f -> (0 &);D[DSolve[test, y[x], x], x] /. x -> 0 – xzczd Jan 12 '18 at 9:53
• Thanks a lot for the extra solution. Have you used this to solve some eigenvalue problem of several coupled ode or pde? (Just would like to learn some code example if you've ever shown somewhere.) Will it be very different for pde? $\qquad$ What I actually want to study is a few coupled 1D ode and a few coupled 2D pde with periodic b.c. in one direction. – xiaohuamao Jan 12 '18 at 10:58
• This is the first time I used pdetoae for eigenvalue problem (you see this type of problem is relatively infrequent in this site, and to be honest I'm not quite familiar with it), but I think the solving process won't be that different from solving ODE/PDE. (The only extra step seems to be transforming the equation(s) to $Lu=\lambda u$ form and take the left part. ) (Update: Oh I forgot I used to write this answer: mathematica.stackexchange.com/a/128275/1871 This method is more efficient than the one above, but less flexible and a bit harder to understand. ) – xzczd Jan 12 '18 at 12:42
• Thanks for the link. I found that in your code, replacing the left b.c. $y'(0)=0$ by $y(0)=\textbf{any value including 0}$ also works well if only we set eps nonzero but vanishingly small to avoid singularity, which reassures the 1st point I noted in the post. This crazily reminds me of the $\epsilon$-$\delta$ language of limit, which is actually the 1st part of my tentative understanding of why the wrong $y(0)=0$ in the naive FDM works. – xiaohuamao Jan 13 '18 at 9:03

Your matrix is not properly constructed. Your "boundary condition" is $f[x]/h=0$, where $h$ is a grid step.

• You caught some point and can take a look at these answers. But my actual question is why it still gives correct solutions. I cannot see the direct reason for this in those answers from that link. Any comment or suggestion is welcome. – xiaohuamao Mar 19 '18 at 3:16
• @xiaohuamiao I’m focusing your attention on the fact that there is no differential equation with any boundary conditions that lead to the matrix you use. It’s not "naive", it’s wrong. – Vsevolod A. Mar 19 '18 at 17:07