Solving the Schrödinger Equation with Numerov's Algorithm

The Schrödinger equation describes the energy and time-evolution of a particle or system of particles, and is one of the fundamental building blocks of modern physics. In it’s general form, the (time-independent) Schrödinger equation for a one-dimensional harmonic oscillator reads thus:¹

$$\begin{array}{}\text{(1)}& \frac{-{\hslash }^2}{2m}\frac{{\partial }^2}{\partial z^2}\psi (z)+\frac{m{z}^2}{2}\psi (z)=E\psi (z)\end{array}$$

There are relatively few situations in which the Schrödinger equation can be solved analytically, and numerical methods and approximations are one way around that analytical limitation. To demonstrate how this is possible and how a numerical solution works, what better way than to solve a system which can be solved analytically and comparing the results.

In solving the Schrödinger equation, we will start with one of the simplest interesting quantum mechanical systems, the quantum mechanical harmonic oscillator.² This is described in $\text{(1)}$ above.

We will make use of the Numerov algorithm which is particularly suited to solving second order differential equations of the form $y\prime \prime (x)+k(x)y(x)=0$. You can read more about it elsewhere, including its derivation etc., but its form is quite simple and easy to code:

$$\begin{array}{}\text{(2)}& (1+\frac{1}{12}{h}^2{k}_{n+1})y_{n+1}=2(1-\frac{5}{12}{h}^2{k}_n)y_n-(1+\frac{1}{12}{h}^2{k}_{n-1})y_{n-1}+O(h^6)\end{array}$$

As you can see, it provides $6^{\text{th}}$ order accuracy which is pretty impressive for such a simple algorithm. In the above equation, $h$ is the step size between each iteration, and $n$ is the index of iteration; $y$ and $k$ relate to those in the formula in the paragraph above. Thus we need to manipulate $\text{(1)}$ into a (dimensionless) form which the Numerov algorithm can solve: using a substitution $E=\epsilon \hslash \omega$ and $z=x\sqrt{\frac{\hslash }{m\omega }}$ we can rearrange $\text{(1)}$ into the form:

$$\begin{array}{}\text{(3)}& \psi \prime \prime (x)+(2\epsilon –x^2)\psi (x)=0\end{array}$$

Now the Schrödinger equation is in the correct form, we can simply plug it into the Numerov algorithm and see what comes out.

 Finding the Eigenvalues Numerically

To determine the eigenvalues, the program scans a range of energies, denoted by the Greek letter $\epsilon$ in the above equations, and tests for when the tail of the graph changes from $+\mathrm{\infty }$ to $-\mathrm{\infty }$ or vice versa. When that happens, the tail must have crossed zero, and therefore it must have stepped over a solution.³ The programme then goes backwards and so on with increased resolution, honing in until it finds all of the solutions we want.

Given the substitution above, we should expect the eigenvalues to be $\epsilon =n+\frac{1}{2}$ where $n$ is an integer from zero (representing the ground state) upwards.⁴ Hit F12 to pull up the web console before you run the simulation to see what results you actually get.

You can find the code for this in JavaScript or a really bare-bones version in Python.