Numerical methods for the Schrodinger equation¶

by Justin Cai and Sam Bateman¶

Solving the time-independent Schrodinger equation¶

Recall the that time-independent Schrodinger equation (TISE) is

$$-\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2} + V \psi = E\psi.$$

When solving the time-independet Schrodinger equation (TISE), we are faced with solving a boundary value problem (BVP), or finding a solution to a differential equation with the solution being specified at two or more points. This is contrasted with an initial value problem (IVP), where we want to find a solution with the solution only specified at one point. For the TISE, if we restrict our particle to the finite domain $[a, b]$, then we'll have that $\psi(a)=0=\psi(b)$ as our boundary conditions.

Shooting method¶

One method to solving boundary value problems is called the shooting method. The shooting method turns the problem into an initial value problem by starting at one of the boundaries, guessing an initial free parameter, solving up to the boundary (that is, "shooting" to the boundary), and checking if our solution matches the boundary condition. If it does not, we update our guess on our free parameter and then repeat this process until our solution is sufficiently close to the boundary condition.

Let's see this in a simple scenario with an infinite square well. Recall that the potential of an infinite square is

$$V(x) = \begin{cases}0 & 0 \leq x \leq L \\ \infty & \text{otherwise} \end{cases},$$

so we get the boundary conditions $\psi(0)=0=\psi(L)$ discussed earlier. We'll have to make a guess for the initial slope, so turning this into an IVP, we have $\psi(0)=0, \psi'(0)=s$, where $s$ is a free parameter. For arbitrary potentials, we also don't know $E$, the total energy of the system, but for now, we'll plug in the analytic solution for the energy of an infinite square well,

$$E_n = \frac{n^2 \pi^2 \hbar^2}{2mL^2}.$$

For convenience, we'll also set $m=1$ and $h=1$ (so $\hbar = \frac{1}{2\pi}$). If we solve this IVP with the correct energy, then we see that the initial slope doesn't really matter:

As we can see, assuming we have the right energy, the initial energy only rescales the wave function. And it turns out that if you satisify the TISE and the boundary conditions, you only need to scale your solution by

$$C = \left(\int_{-\infty}^{\infty} \psi(x)^2 dx\right)^{-1}$$

to satisfy normalization.

However, what if we do not know $E$? The choice of $E$ does matter, so an incorrect choice of $E$ might not satisfy the boundary conditions. In the example below, we solve the IVP but with the energy given perturbed a small amount from the actual energy.

Instead of guessing the initial slope (which in general, any value will do), we can guess the energy. If $\psi^{(E)}$ is a solution with $E$ as the energy, then we can rephrase the shooting problem as finding the roots of the function

$$ g(E) = \psi^{(E)}(L). $$

In the example above, we have $g(E_1-0.05) = \psi^{(E_1-0.05)}(L)=0.0027$ and $g(E_1-0.05)=\psi^{(E_1+0.05)}(L)=-0.0015$.

A simple root finding method for a continuous function (which we are assuming $g$ is) comes from the intermediate value theorem: if $f(a)$ and $f(b)$ have opposite signs, then there exists some $x$ such that $f(x)=0$. The bisection method takes advantage of this. Given a interval and a (continuous) function taking opposite signs at the end of the interval, the method calculates the function at the midpoint of the interval, checks the sign, and then updates the interval with this new midpoint so that the function still has opposite signs with the new interval, and repeating this until the interval is sufficiently small.

From the example above, since $g(E_1-0.05)$ and $g(E_1+0.05)$ have opposite signs, we know that there exists a root within that interval (and of course, the root is at $E_1$). Let's see the shooting method in action now. Using the bisection method, we're going to need an inital energy range, so for simplicity, we'll just use the $[E_1-0.05, E_1+0.05]$ interval from above.

As we can see, the calculated energy is very close to the original energy, and there is almost no difference in the solutions when using either energy.

Matching method¶

The shooting method, however, will not be sufficient for other potentials. Let's see why with the finite square well. We'll be using the potential

$$V(x)=\begin{cases} -V_0 & -L \leq x \leq L \\ 0 & -L' \leq x < L \text{ or } L < x \leq L'\\ \infty & \text{otherwise}\end{cases}$$

for the finite square well (only difference from class is the infinite height barrier at a distance $L'$ from the origin, keeping the particle in a finite domain). Let's set $V_0=1$, $L=1$, $L'=5$, and then see what a solution looks like with $E=-0.5$.

In the classically forbiden region (that is, where $V(x)>E$), the solution blows up due to numerical error. As the solutions in these regions are exponetial functions, the "opposite" exponential takes over (that is, instead of the decaying solution $e^{-\alpha x}$, we get the increasing solution $e^{\alpha x}$). However, note that the solution is correct in the classically forbidden region on the left, which is where we start solving from. What does the solution look like if we start solving from the right?

The solution from the right has the correct functional form, so we could "splice" the two solutions together to make a new solution. If we do this for this particular energy by taking the solution from the left all the way up to the classically beginning of the forbidden region on the right, and then taking the solution from the right, we get a discontinous solution, which does not satisfy the Schrodinger equations. However, if we change the energy to be close to an allowed energy...

the solutions "match", and so the spliced solution is continous. This is the essence of the matching method. We shoot from the left and right up to some matching point, which should happen before any divergent behavior for both solutions, and then aim to match the function and the derivative at that point. As a root finding problem, if $\overrightarrow{\psi^{(E)}}$ is a solution from the left and $\overleftarrow{\psi^{(E)}}$ is a solution from the right, both with energy $E$, and $x_m$ is the match point, then we want to find the root of

$$g(E) = \overrightarrow{\psi^{(E)}}\,'(x_m)-\frac{\overrightarrow{\psi^{(E)}}(x_m)}{\overleftarrow{\psi^{(E)}}(x_m)} \overleftarrow{\psi^{(E)}}\,'(x_m)\quad,$$

which will find the energy when the slope difference between the matched solutions is 0. For the matching method, we'll use the secant method, which is similar to Newton's root finding method, but we don't need to supply the derivative of $g$ (and is instead approximated). Below is a solution that we've matched.

Direct matrix method¶

The matching works well for general potentials, but we still need to search for the energy values. However, there is a method that can solve for the approximate energies directly, called the direct matrix method. Recall that we can write the TISE as

$$\hat{H}\psi(x) = E \psi(x)$$

where $\hat{H} = -\frac{\hbar}{2m}\frac{d^2}{dx^2} + V(x)$. So finding the energies reduces to finding the eigenvalues of $\hat{H}$. However, we'll need to discretize $\hat{H}$ in order to use numerical methods for solving eigenvalues. By discretizing our domain into $N$ points as $x_j$ where $j ={0, 1, \dots, N-1}$, we can use the central difference formula to approxmiate the TISE as

$$-\frac{\hbar^2}{2m}\frac{\psi_{j+1} - 2\psi_{j} + \psi_{j-1}}{\Delta x^2} + V_j\psi_j = E\psi_j + O(\Delta x^2) \quad \text{for $j=1, 2, \dots, N-2$} $$

where $x_{j+1}-x_{j} = \Delta x$, $\psi_j = \psi(x_j)$, and $V_j = V(x_j)$. Rearranging with $k = \frac{\hbar^2}{2m\Delta x^2}$, $\psi_0=0$, we have the following system of equations:

\begin{align} -k\psi_2 + (2k+V_1)\psi_1 &= E\psi_1 \\ -k\psi_3 + (2k+V_2)\psi_2 -k\psi_1 &= E\psi_2\\ -k\psi_4 + (2k+V_3)\psi_3 -k\psi_2 &= E\psi_3\\ \vdots\\ (2k+V_{N-2})\psi_{N-2}-k\psi_{N-3} &= E\psi_{N-2} \end{align}

We can write this in matrix form,

\begin{align} \begin{bmatrix} 2k+V_1 & -k & 0 & \dots & 0 \\ -k & 2k+V_2 & -k & \ddots & \vdots \\ 0 & -k & 2k+V_3 & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -k \\ 0 & \dots & 0 & -k & 2k+V_{N-2} \end{bmatrix} \begin{bmatrix}\psi_1 \\ \psi_2 \\ \psi_3 \\ \vdots \\ \psi_{N-2} \end{bmatrix} = E \begin{bmatrix}\psi_1 \\ \psi_2 \\ \psi_3 \\ \vdots \\ \psi_{N-2} \end{bmatrix} \end{align}

or more compactly, $H\psi = E\psi + O(\Delta x^2)$. Hence, we can approximate the eigenvalues of $\hat{H}$ by finding the eigenvalues of the discretized Hamiltonian $H$. In practice, the eigenvalues accumlate error for the higher energy states, but for small enough energy states and the matching method, we can now solve for the $n$th exited state of arbitrary potentials.

One system we can use to check our calculations is the quantum harmonic oscillator, whose analytic solutions are well known. Its potential is given by $V(x) = \frac{1}{2}m\omega^2x^2$, and its energy levels are $E_n = \hbar\omega (n+\frac{1}{2})$. Using the direct matrix method and letting $\omega=1$, $m=1$ and $h=1$, the plot below shows the difference between the analytic and diect matrix method. We can see that increasing the subdivision (and thus decreasing $\Delta x$) increases the accuracy of the estimate.

Using the energies from the direct matrix method, below are some of the excited states for the quantum harmonic oscillator.