Solving one dimensional Schrodinger equation with finite difference method

Question

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)

Answer 1

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.

Answer 2

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.