Solving one dimensional Schrodinger equation with finite difference method

11,898

Solution 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.

Solution 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.

Share:
11,898

Related videos on Youtube

NGY
Author by

NGY

Updated on April 10, 2020

Comments

  • NGY
    NGY over 3 years

    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)

    • BebopButUnsteady
      BebopButUnsteady over 12 years
      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?
    • Ted Bunn
      Ted Bunn over 12 years
      I bet that @BebopButUnsteady is right. Often when numerically solving systems with singularities, you have to artifically "soften" the singularity.
    • Marek
      Marek over 12 years
      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.
    • NGY
      NGY over 12 years
      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?
    • Marek
      Marek over 12 years
      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.
    • NGY
      NGY over 12 years
      What do you mean by "bounded from below"? The value $-1/2$ is a classical result(see, for example, iopscience.iop.org/0305-4470/10/6/001), in fact, it is analogous to the ground state of hydrogen atom.
    • Marek
      Marek over 12 years
      @NGY: "bounded from below" means precisely that: quantity $b$ is b.f.b. if there is quantity $a$ such that $b \geq a$. You are right about the hydrogen. | Ah, right, thanks. I thought that it is bounded from below because of the centrifugal term $L^2 \over r^2$ but that only occurs for $L > 0$. Spherically symmetric states are indeed described by your equation.
    • NGY
      NGY over 12 years
      Can anyone calculate it and compare the result with mine?
    • Ruslan
      Ruslan almost 9 years
      This equation won't work because it has kinetic energy operator of flat space instead of that of radial kinetic energy operator, which will also contain first derivative. In one dimension the Coulomb potential has unbounded from below spectrum.
  • Keenan Pepper
    Keenan Pepper over 12 years
    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.
  • Marek
    Marek over 12 years
    +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).
  • BebopButUnsteady
    BebopButUnsteady over 12 years
    @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?
  • NGY
    NGY over 12 years
    @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.
  • Qmechanic
    Qmechanic over 12 years
    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...
  • NGY
    NGY over 12 years
    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.