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:1

\begin{equation} \label{eq:sch} \frac{-\hbar^2}{2m} \frac{\partial^2}{\partial z^2}\psi(z) + \frac{mz^2}{2} \psi(z) = E\psi(z) \end{equation}

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.2 This is described in \eqref{eq:sch} 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{equation} \label{eq:numerov} \left(1+\frac{1}{12}h^2k_{n+1}\right)y_{n+1} = 2\left(1-\frac{5}{12}h^2k_n\right)y_n -\left(1+\frac{1}{12}h^2k_{n-1}\right)y_{n-1}+O(h^6) \end{equation}

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 \eqref{eq:sch} into a (dimensionless) form which the Numerov algorithm can solve: using a substitution \(E=\varepsilon \hbar \omega\) and \(z=x\sqrt{\frac{\hbar}{m\omega}}\) we can rearrange \eqref{eq:sch} into the form:

\begin{equation} \label{eq:solve} \psi\prime\prime(x) + (2\varepsilon – x^2)\psi(x) = 0 \end{equation}

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 \(\varepsilon\) in the above equations, and tests for when the tail of the graph changes from \(+\infty\) to \(-\infty\) or vice versa. When that happens, the tail must have crossed zero, and therefore it must have stepped over a solution.3 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 \(\varepsilon = n + \frac 12\) where \(n\) is an integer from zero (representing the ground state) upwards.4 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.

  1. where \(m\) is the mass of the particle, \(x\) is the position, \(\hbar\) is the Plank constant, \(V(x)\) is the potential the particle is in, \(E\) is the particle’s energy, and \(\psi(x)\) is the wavefunction. ↩︎

  2. A harmonic oscillator is simply an object which is moving in a constant field of some kind. For example a gravitational or electric field. A good example is a pendulum, or a ball bouncing on a piece of elastic. Harmonic oscillators are vitally important in physics and physical chemistry, because they can be used to model the complex vibrations and oscillations of molecules, atoms, and sub-atomic particles to a reasonable degree of accuracy, and are really rather simple to solve. ↩︎

  3. This is because wavefunctions have to be normalizable, and so a wavefunction which goes to infinity as x increases is not a physically relevant one. ↩︎

  4. Don’t confuse this \(n\) with the index \(n\) from the definition of the Numerov algorithm! ↩︎