3
$\begingroup$

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'$.
enter image description here
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$.

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ Seems to be more of a mathematics question than one about how Mathematica works. $\endgroup$
    – Michael E2
    Commented Jan 11, 2018 at 13:55
  • $\begingroup$ @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. $\endgroup$
    – xiaohuamao
    Commented Jan 12, 2018 at 6:06
  • $\begingroup$ @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. $\endgroup$
    – xiaohuamao
    Commented Jan 12, 2018 at 6:30
  • $\begingroup$ 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.? $\endgroup$
    – xzczd ♦
    Commented Jan 12, 2018 at 6:30
  • $\begingroup$ @xzczd Exactly! $\endgroup$
    – xiaohuamao
    Commented Jan 12, 2018 at 6:31

2 Answers 2

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
7
  • $\begingroup$ 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? $\endgroup$
    – xiaohuamao
    Commented Jan 12, 2018 at 8:28
  • $\begingroup$ 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 $\endgroup$
    – xzczd ♦
    Commented Jan 12, 2018 at 9:53
  • $\begingroup$ 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. $\endgroup$
    – xiaohuamao
    Commented Jan 12, 2018 at 10:58
  • $\begingroup$ 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. ) $\endgroup$
    – xzczd ♦
    Commented Jan 12, 2018 at 12:42
  • $\begingroup$ 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. $\endgroup$
    – xiaohuamao
    Commented Jan 13, 2018 at 9:03
0
$\begingroup$

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

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – xiaohuamao
    Commented Mar 19, 2018 at 3:16
  • $\begingroup$ @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. $\endgroup$
    – Vsevolod A.
    Commented Mar 19, 2018 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.