On numerically solving the Schrödinger equation

Question

I just read a paper 'A pocket calculator determination of energy eigenvalues' by J Killingbeck (1979).

Link: http://iopscience.iop.org/0305-4470/10/6/001

I have some questions about section 2.

1. Why $R<0$ implies the $E$ estimate was too high and $R>1$ implies too low? And is it possible that $0<R<1$ remains?

2. Why "if the starting point is chosen at $r=h(l+1)$, then $R(hl)$ can be assigned any finite value without disturbing the calculation at larger $r$ values" ?

3. How to guess a trial energy $E$ (especially the ground state energy)?

(Actually I'm interested in a later paper by Killingbeck (1982) 'Finite-difference methods for eigenvalues'(Link: http://iopscience.iop.org/0022-3700/15/6/009), in which the author applies the method in 'A pocket calculator determination of energy eigenvalues' to some potentials with highly singularity at origin)

---

mod note: copying the text of whoplisp's post here since it isn't actually an answer to the question

They convert the Schroedinger equation into a finite difference equation (instead of a variational approach):

$-\frac{1}{2}\Delta\psi +V\psi=E\psi$

$V$ is a central potential

Ansatz with solid harmonic $\xi_l$ of degree $l$ (e.g. $x$ or $xy$)

$\psi=\xi_l\phi(r)$

leads to

$2r(V-E)\phi=rD^2\phi+(2l+2)D\phi$

Notation: $D=\frac{d}{dr}$

finite difference approximation

$h^2D^2\phi \rightarrow \phi(r+h)+\phi(r-h)-2\phi(r)$

$2hD\phi \rightarrow \phi(r+h)-\phi(r-h)$

crucial step: introduce ratio variable $R(r)$:

$\phi(r+h)=R(r)\phi(r)$

convert Schroedinger equation in recursive equation to calculate $R(r)$:

$2r[1+h^2(V-E)]=R(r)[r+h(l+1)]+[r-h(l+1)](R(r-h))^{-1}$

the corresponding 1D Schroedinger equation is:

$2[1+h^2(V-E)]=R(x)+(R(x-h))^{-1}$

if $r=h(l+1)$ is chosen then $R(hl)$ can be assigned any finite value without disturbing the calculation at large $r$ values

to find the ground states with angular momentum $l$ start with an estimate for $E$, some non-zero $R(hl)$ and some small $h$. $R(r)$ as $r$ increases will either become negative or pass through unity from below.

first thing happens when guessed $E$ was too high

second thing indicates $E$ was too low (because $\psi$ increases with $r$)

you just try various values and 'sandwhich' the solution

example for $V=\lambda r-1/r$:

$\lambda=1$, $h=.1$, state=2p: 1.9759

Answer 1

I can't entirely make sense of some the claims in that letter, but this looks like a version of solving for the bound states of potential which is fairly standard. It's in Griffiths Quantum Mechanics I think. I'll write it down here for you and maybe you can piece together the rest. Anyway it goes like this:

We're trying to find the values of $E$ for which the radial Schrodinger equation $-\frac{1}{r^2}\frac{\partial}{\partial r}(r^2\frac{\partial\psi}{\partial r} )+ \frac{\ell(\ell+1)}{r^2}\psi + V(r)\psi = E\psi$ has normalizable solutions. On the other hand that's a second order ODE for any $E$ so it has to have two solution - they're just not going to be normalizable. Specifically they are going to blow up quite badly as $r\rightarrow\infty$. So we are just going to guess a value of $E$, numerically solve that equation by finite differences and see if the solutions blow up at infinity or not. When we found a value of $E$ where they don't blow up we know we found a bound state.

As for the particular issues that come up when we are looking for the ground state.

1. $R < 0$ means your wavefunction has a node (it changes sign) and $R > 1$ means your wavefunction is increasing. It is the case that the ground state wavefunction has no nodes so you must be looking in the vicinity of one of the excited states. As for $R > 1$, while asymptotically your wavefunction has to decrease, I don't see why its not possible that it could increase at some point. So this does not seem to be valid termination rule in general (though given a particular potential you might be able to say more).

2. See that $R(h(l+1))$ does not depend on $R(hl)$ since the prefactor is zero. After that $R(x)$ depends only on $R(x-h)$ so $R(hl)$ never enters.

3. This method is very computationally cheap so why not just try a lot of energies? You can always bound the energy from above by the standard variational methods. Also just dimensional analysis could give you the order of magnitude.

Answer 2

Bebop answered all the original questions already. However it might be useful to see an implementation. The following code calculates the energies of the 1s state with lambda=0 and lambda=1 with varying values of $h$ (see table 1 in the paper).

I start with a value of 0.5 for $r=h(l+1)$ and then use the recursive equation to calculate R(r) for increasing r in the function calc-domain. The function returns, when R is either smaller then 0 or bigger than 1.

The function bracket-E calls calc-domain with different guesses for E and tries to home in onto the value where R stays within $[0,1]$ for the longest distance.

In table 1 the paper gives -0.4815 as the value for E for lambda=0 and h=0.4. This code returns -0.48145604506134987. Note that the calculation was done using rational bignums and should be correct to 42 digits. I rounded the result into a double for convenience.

;see 0305-4470_10_6_001.pdf, runs in Steel Bank Common Lisp http://www.sbcl.org/

(defun Rr (r h l E Rr-h lam)
  "Caculate R at position r using a value of R at r-h and a guess E."
  (let* ((V (- (* lam r) (/ r)))
     (lhs (* 2 r (+ 1 (* h h (- V E))))))
    (/ (- lhs (/ (- r (* h (1+ l)))
         Rr-h))
       (+ r (* h (1+ l))))))

(defun calc-domain (E h l lam)
  "Calculate ratio function over the domain from r0 to big r,
if Rhl becomes negative, the guess E was too high"
  (let* ((r0 (* h (1+ l)))
     (Rhl .5s0))
    (loop while (< (- (* Rhl Rhl) Rhl) 0) 
       do
     (setf Rhl (Rr r0 h l E Rhl lam))
     (incf r0 h))
    Rhl))

(defun bracket-E (left right h l lam)
  (let* ((rleft (calc-domain left h l lam))
     (rright (calc-domain right h l lam))
     (center (/ (+ left right) 2))
     (rcenter (calc-domain center h l lam)))
    (when (< (abs (- left right)) 1e-42)
      (return-from bracket-E (* 1d0 center)))
    (cond 
      ((< rcenter 0 rleft)
       (bracket-E left center h l lam))
      ((< rright 0 rcenter)
       (bracket-E center right h l lam))
      (t (break "error: There is no zero between ~a, try increasing the interval" 
        (list left right))))))

;; reproduce values for the 1s state of table 1:
(loop for h in '(4/10 2/10 1/10 5/100 1/1000) collect
     (bracket-E -2 2 h 0 0))

;; lam=0
;; (-0.48145604506134987d0 -0.49509748816490173d0 -0.49875694513320923d0
;;     -0.4996879608376049d0 -0.49949414802320075d0)

;; lam=1
;; (0.6055558981315698d0 0.5841694474220276d0 0.5794421434402466d0
;;            0.5782957077026367d0 0.5777209093467737d0)