Solving one dimensional Schrodinger equation with finite difference method

Consider the one-dimensional Schrodinger equation

$$-\frac{1}{2}D^2 \psi(x)+V(x)\psi(x)=E\psi(x)$$

where $D^2=\dfrac{d^2}{dx^2},V(x)=-\dfrac{1}{|x|}$.

I want to calculate the ground state energy(which is known to be $-0.5$) with finite difference method.

Suppose $\psi(-a)=\psi(a)=0$ for sufficiently large $a$.

Split $[-a,a]$ into $N$ sub-intervals with equal length $h$: $$-a=x_1<x_2<\cdots<x_{N+1}=a$$ $x_{i+1}-x_i=h.$

Now use

$$D^2\psi(x_n)\to\frac{\psi(x_{n-1})-2\psi(x_{n})+\psi(x_{n+1})}{2h^2}$$

and

$$-\frac{1}{2}D^2 \psi(x_n)+V(x_n)\psi(x_n)=E\psi(x_n)$$

we get

$$-y_{n-1}+2(1+V_n h^2)y_n-y_{n+1}=2h^2 Ey_n$$

where $y_n=\psi(x_n)$.

Recall that $y_1=y_{N+1}=0$, we get

$$H \begin{pmatrix} y_2 \\ \vdots \\ y_N \end{pmatrix} = 2h^2E \begin{pmatrix} y_2 \\ \vdots \\ y_N \end{pmatrix}$$

where

$$H=\begin{pmatrix} 2(1+V_2h^2) & -1 & & \\ -1 & 2(1+V_3h^2) & \ddots & \\ & \ddots & \ddots & -1 \\ & & -1 & 2(1+V_Nh^2) \end{pmatrix}$$

is a $N-1$ tridiagonal matrix.

If the least eigenvalue of $H$ is $\lambda$, then I expect $\dfrac{\lambda}{2h^2}$ to approximate the ground state energy $E$.

So I write a MATLAB program:

function m = onedimen1(a,M)

N = 2*M+1;

h = 2*a/N;
t = -a;

for i = 1 : N-2
t = t+h;
H(i,i) = 2*(1+vp(t)*h*h);
H(i+1,i) = -1;
H(i,i+1) = -1;
end

H(N-1,N-1) = 2*(1+vp(t+h)*h*h);

m1=eig(H)/(2*h*h);
m = m1(1:5);

function y = vp(x)
y = - 1/(abs(x)) ;
end
end


but when I try

onedimen1(10,100)
onedimen1(10,1000)


the least eigenvalue of $H$ divided by $2h^2$ appears to be too small(and tends to negative infinity when $N$ increases), but the second largest divided by $2h^2$ is always near $-0.5$.

Can anyone explain why?

(P.S. If I change $V(x)$ to $V(x)=\dfrac{1}{|x|}$ then the least eigenvalue of $H$ is stable)

-
If you use a non-singular potential like $V(x) = -\exp(-x^2)$ does it work in that case? The eigenvector corresponding to your diverging eigevalue: does it look like it diverges as well? Maybe it corresponds to a non-renormalizable state? –  BebopButUnsteady Jul 25 '11 at 14:24
I bet that @BebopButUnsteady is right. Often when numerically solving systems with singularities, you have to artifically "soften" the singularity. –  Ted Bunn Jul 25 '11 at 15:13
This is obviously connected to the fact that $-{1 \over |x|}$ is not bounded from below. Can you support your claim that the least eigenvalue is $-1/2$? If you use $1 \over |x|$ then there is only continuous spectrum bounded from below. Naturally, this must be stable. –  Marek Jul 25 '11 at 15:43
You mean "non-normalizable state"? Is it acceptable in physics? And Ted, what do you mean by "soften the singularity"? Could you give an example? –  NGY Jul 25 '11 at 15:44
Another thing to check: a bound eigenstate (i.e. corresponding to the discrete part of the spectrum) must be from your Hilbert space. It is quite possible that there is some strange state corresponding $-\infty$ energy that can't be normalized. But when considering discretized system, this state will reappear in the finite-dimensional Hilbert space (consisting of piece-wise constant functions with compact support [-a, a]) since there can be no divergence in this finite approximation; it only shows up in the $N \to \infty$ limit. –  Marek Jul 25 '11 at 15:46

I'm making this an answer because its too long to be a comment but its just an expansion of the things already stated in comments:

Non-normalizable states: The Schroedinger equation has an infinity of solutions but almost all of them do not have a finite norm ($\int|\psi(x)|^2dx$ is not finite). These are not phyiscally acceptable, since there would not be a probabilistic interpretation, amongst other issues. It is possible that when you discretize your system you are getting solutions that would have infinite norm in the limit $h\rightarrow 0$ and so are not physically acceptable. I strongly suspect this is the issue. To see if this is case look at the wavefunction associated with your bad eigenvalue ( this is the eigenvector associated with the eigenvalue). If in the limit $h\rightarrow0$, your solution goes to increases faster than $\frac{1}{\sqrt{x}}$ as $x\rightarrow0$ than it is non-normalizable and unphysical.

'Softening the singularity': the issue is that you potential goes to infinity at some point and that general numerical simulations don't like it when things go to infinity. Especially having your potential go to negative infinity since the wavefunction wants to gather around the singularity. The way around this would be to replace your potential $\frac{1}{|x|}$ with say $\frac{1}{\sqrt{x^2+a^2}}$ where $a$ is some number you pick, and solve for the states. This new potential is the same as the old when $x\gg a$ but it doesn't have any singularity anymore, so it shouldn't be a problem. To get back the real answer solve for the states with smaller and smaller values of $a$ until you see that if a limit is reached. It does not necessarily reach a limit and in this case you will still get the state that diverges.

-
Except that if you solve this with any softening of the singularity, such as $1/\sqrt{x^2+a^2}$. you'll get a no-node ground state, and the energy of that state goes to $-\infty$ as $a$ goes to zero. NGY wants the one-node state to be the ground state, but for any softening of the singularity there is a no-node ground state. Only for $a=0$ does that no-node state not exist. –  Keenan Pepper Jul 25 '11 at 18:33
+1, just a sidenote: the states corresponding to the continuous part of the spectrum are also not normalizable but they are physically acceptable as idealizations (plane-waves, Dirac deltas). –  Marek Jul 25 '11 at 18:35
@Keenan Pepper: I agree with you. I was just trying to explain the idea - it does not solve his problem here on its own. I have made that more explicit. I think in this case you can just apply the boundary condition $\psi \rightarrow 0$ as $x\rightarrow0$. Do you agree? –  BebopButUnsteady Jul 25 '11 at 19:05
@BebopButUnsteady I just tried to compute in the interval $[0,a]$, using the boundary condition $\psi(0)=\psi(a)=0$ and get the right answer $-0.5$. But I wonder whether it is reasonable to suppose $\psi(0)=0$ if we don't know the really wave function. –  NGY Jul 26 '11 at 0:57

Mathematically, the problems with the potential $1/|x|$ are due to the fact, that the resulting Hamiltonian is not continuous on $L^2(R)$. You can expect anything with such operators - e.g., for this particular example, there is a continuum of normalizable "states" for every negative energy, which are not mutually orthogonal w.r.t. the standard norm in $L^2(R)$. One is better to avoid such anomailies at all.

One possible way to regularize the problem is to consider the odd solutions only (which is equivalent to a 3D hydrogen atom problem).

Concerning the issues of numerical stability, I'd try to change the variable for the spatial coordibnate to make the potential less singular.

-
The singularity itself seems not to explain all of OP's problem. Note that the number of sub intervals $N = 2*M+1$ is chosen in the MATLAB code to be odd, so for fixed $N$, the MATLAB program doesn't "know" that the potential is singular... –  Qmechanic Jul 28 '11 at 9:37
How to change the variable for the spatial coordinate properly in this case? I have tried $r\to t:\exp(Kt)=1+Kr$, but the numerical solution seems to be more unstable. –  NGY Jul 30 '11 at 16:20