# On numerically solving the Schrödinger equation

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

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)$

$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

-
I would like to point out that both of those articles are behind paywalls. People without institutional subscriptions are unlikely to be able to access them. – Richard Terrett Jul 22 '11 at 15:01
This might amuse you: nibot-lab.livejournal.com/86063.html – nibot Jul 22 '11 at 17:04
What /Who is Schordinger? – Georg Jul 22 '11 at 17:40
Is it possible to extract the relevant information from the paper(s) for those of us who do not have access to it ? – Frédéric Grosshans Jul 22 '11 at 17:56

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.

-
2. What about the one-dimensional case? How to start to calculate $R(r)$? – NGY Jul 23 '11 at 1:02
3. Variational method can give an upper bound of $E$, then how to give a lower bound? – NGY Jul 23 '11 at 1:04

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)

-
Sorry, how to run the code? I just installed SBCL on my Ubuntu, but dont't know how to run your program. – NGY Jul 23 '11 at 12:40
The easiest way would be to run sbcl and paste the code in there. This defines the two functions and runs the loop with different values for h. After a few seconds it should return a list starting with -0.481 and ending with -0.4994. You can also search for a solution with a call to (bracket-E -2 2 1/10 0 0). Note that no one actually runs Lisp like this. Normally one uses the SLIME environment within Emacs as I described here stackoverflow.com/questions/6546834/…. – whoplisp Jul 23 '11 at 13:01
I opened a channel chat.stackexchange.com/rooms/902/lisp. You can direct your questions there. – whoplisp Jul 23 '11 at 15:09
Thanks. And I have a question about the algorithm(not the code), which I comment on Behop's answer(but it seems that he didn't see the comment). That is, how to guess a value of $E$? It really puzzled me. – NGY Jul 23 '11 at 17:08
You can use any value for $E$. If the function calc-domain returns a negative value, the guess was too high. If the calc-domain returns a value bigger than 1, the guess was too low. By construction it never returns a value in [0,1]. bracket-E is called with -2 and 2, but you could also call it with -100 and 100. Just make sure the real value is inbetween. – whoplisp Jul 23 '11 at 17:40