# 1.5: Numerically Solving the Schrödinger Equation

$$\newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$$ $$\newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$$$$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\kernel}{\mathrm{null}\,}$$ $$\newcommand{\range}{\mathrm{range}\,}$$ $$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$ $$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}[1]{\| #1 \|}$$ $$\newcommand{\inner}[2]{\langle #1, #2 \rangle}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\kernel}{\mathrm{null}\,}$$ $$\newcommand{\range}{\mathrm{range}\,}$$ $$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$ $$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}[1]{\| #1 \|}$$ $$\newcommand{\inner}[2]{\langle #1, #2 \rangle}$$ $$\newcommand{\Span}{\mathrm{span}}$$

Often the bound potentials that we encounter are complex, and the time-independent Schrödinger equation will need to be evaluated numerically. There are two common numerical methods for solving for the eigenvalues and eigenfunctions of a potential. Both methods require truncating and discretizing a region of space that is normally spanned by an infinite dimensional Hilbert space. The Numerov method is a finite difference method that calculates the shape of the wavefunction by integrating step-by-step across along a grid. The DVR method makes use of a transformation between a finite discrete basis and the finite grid that spans the region of interest.

## The Numerov Method

A one-dimensional Schrödinger equation for a particle in a potential can be numerically solved on a grid that discretizes the position variable using a finite difference method. The TISE is

$[ T + V (x) ] \psi (x) = E \psi (x) \label{128}$

with

$T = - \dfrac {\hbar^{2}} {2 m} \dfrac {\partial^{2}} {\partial x^{2}},$

which we can write as

$\psi^{\prime \prime} (x) = - k^{2} (x) \psi (x) \label{129}$

where

$k^{2} (x) = \dfrac {2 m} {\hbar^{2}} [ E - V (x) ].$

If we discretize the variable $$x$$, choosing a grid spacing $$\delta x$$ over which $$V$$ varies slowly, we can use a three point finite difference to approximate the second derivative:

$f _ {i}^{\prime \prime} \approx \dfrac {1} {\delta x^{2}} \left( f \left( x _ {i + 1} \right) - 2 f \left( x _ {i} \right) + f \left( x _ {i - 1} \right) \right) \label{130}$

The discretized Schrödinger equation can then be written in the form

$\psi \left( x _ {i + 1} \right) - 2 \psi \left( x _ {i} \right) + \psi \left( x _ {i - 1} \right) = - k^{2} \left( x _ {i} \right) \psi \left( x _ {i} \right) \label{131}$

Using the equation for $$\psi \left( x _ {i + 1} \right)$$, one can iteratively solve for the eigenfunction. In practice, you discretize over a range of space such that the highest and lowest values lie in a region where the potential is very high or forbidden. Splitting the space into N points, chose the first two values $$\psi \left( x _ {1} \right) = 0$$ and $$\psi \left( x _ {2} \right)$$x to be a small positive or negative number, guess $$E$$, and propagate iteratively to $$\psi \left( x _ {N} \right)$$. A comparison of the wavefunctions obtained by propagating from $$x_1$$ to $$x_N$$ with that obtained propagating from $$x_N$$ to $$x_1$$ tells you how good your guess of $$E$$ was.

The Numerov Method improves on Equation \ref{131} by taking account for the fourth derivative of the wavefunction $$\Psi^{( 4 )}$$, leading to errors on the order $$O \left( \delta x^{6} \right)$$. Equation \ref{130} becomes

$f _ {i}^{\prime \prime} \approx \dfrac {1} {\delta x^{2}} \left( f \left( x _ {i + 1} \right) - 2 f \left( x _ {i} \right) + f \left( x _ {i - 1} \right) \right) - \dfrac {\delta x^{2}} {12} f _ {i}^{( 4 )} \label{132}$

By differentiating Equation \ref{129} we know

$\psi^{( 4 )} (x) = - \left( k^{2} (x) \psi (x) \right)^{\prime \prime}$

and the discretized Schrödinger equation becomes

\left.\begin{aligned} \psi^{\prime \prime} \left( x _ {i} \right) & = \dfrac {1} {\delta x^{2}} \left( \psi \left( x _ {i + 1} \right) - 2 \psi \left( x _ {i} \right) + \psi \left( x _ {i - 1} \right) \right) + \\ & \dfrac {1} {12} \left( k^{2} \left( x _ {i + 1} \right) \psi \left( x _ {i + 1} \right) - 2 k^{2} \left( x _ {i + 1} \right) \psi \left( x _ {i} \right) + k^{2} \left( x _ {i + 1} \right) \psi \left( x _ {i - 1} \right) \right) \end{aligned} \right. \label{33}

This equation leads to the iterative solution for the wavefunction

$\psi \left( x _ {i + 1} \right) = \dfrac{\psi \left( x _ {i} \right) \left( 2 + \dfrac {10 \delta x^{2}} {12} k^{2} \left( x _ {i} \right) \right) - \psi \left( x _ {i - 1} \right) \left( 1 - \dfrac {\delta x^{2}} {12} k^{2} \left( x _ {i - 1} \right) \right)}{1 - \dfrac {\delta x^{2}} {12} k^{2} \left( x _ {i + 1} \right)} \label{134}$

## Discrete Variable Representation (DVR)

Numerical solutions to the wavefunctions of a bound potential in the position representation require truncating and discretizing a region of space that is normally spanned by an infinite dimensional Hilbert space. The DVR approach uses a real space basis set whose eigenstates $$\varphi _ {i} (x)$$ we know and that span the space of interest—for instance harmonic oscillator wavefunctions—to express the eigenstates of a Hamiltonian in a grid basis ($$\theta _ {j}$$) that is meant to approximate the real space continuous basis $$\delta (x)$$. The two basis sets, which we term the eigenbasis ($$\varphi$$) and grid basis ($$\theta$$), will be connected through a unitary transformation

$$\Phi^{\dagger} \varphi (x) = \theta (x) \label{135}$$ $$\Phi \theta (x) = \varphi (x)$$

For $$N$$ discrete points in the grid basis, there will be $$N$$ eigenvectors in the eigenbasis, allowing the properties of projection and completeness will hold in both bases. Wavefunctions can be obtained by constructing the Hamiltonian in the eigenbasis,

$$H = T ( \hat {p} ) + V ( \hat {x} ),$$ transforming to the DVR basis, $$H^{D V R} = \Phi H \Phi,$$ and then diagonalizing.

Here we will discuss a version of DVR in which the grid basis is set up to mirror the continuous $$| \mathcal {X} \rangle$$ eigenbasis. We begin by choosing the range of $$x$$ that contain the bound states of interest and discretizing these into $$N$$ points ($$x_i$$) equally spaced by $$δx$$. We assume that the DVR basis functions $$\theta _ {j} \left( x _ {i} \right)$$ resemble the infinite dimensional position basis

$\theta _ {j} \left( x _ {i} \right) = \sqrt {\Delta x} \delta _ {i j} \label{136}$

Our truncation is enabled using a projection operator in the reduced space

$P _ {N} = \sum _ {i = 1}^{N} | \theta _ {i} \rangle \langle \theta _ {i} | \approx 1 \label{137}$

which is valid for appropriately high $$N$$. The complete Hamiltonian can be expressed in the DVR basis DVR

$H^{D V R} = T^{D V R} + V^{D V R}.$

For the potential energy, since $$\left\{\theta _ {i} \right\}$$ is localized with $$\left\langle \theta _ {i} | \theta _ {j} \right\rangle = \delta _ {i j}$$, we make the DVR approximation, which casts $$V^{DVR}$$ into a diagonal form that is equal to the potential energy evaluated at the grid point:

$V _ {i j}^{D V R} = \left\langle \theta _ {i} | V ( \hat {x} ) | \theta _ {j} \right\rangle \approx V \left( x _ {i} \right) \delta _ {i j} \label{138}$

This comes from approximating the transformation as $$\Phi V ( \hat {x} ) \Phi^{\dagger} \approx V \left( \Phi \hat {x} \Phi^{\dagger} \right) .$$

For the kinetic energy matrix elements $$\left\langle \theta _ {i} | T ( \hat {p} ) | \theta _ {j} \right\rangle$$, we need to evaluate second derivatives between different grid points. Fortunately, Colbert and Miller have simplified this process by finding an analytical form for the $$T^{DVR}$$ matrix for a uniformly gridded box with a grid spacing of $$∆x$$.

$T _ {i j}^{\mathrm {DVR}} = \dfrac {\hbar^{2} ( - 1 )^{i - j}} {2 m \Delta x^{2}} \left\{\begin{array} {c c} {\pi^{2} / 3} & {i = j} \\ {2 / ( i - j )^{2}} & {i \neq j} \end{array} \right\} \label{139}$

This comes from a Fourier expansion in a uniformly gridded box. Naturally this looks oscillatory in $$x$$ at period of $$δx$$. Expression becomes exact in the limit of $$N \rightarrow \infty$$ or $$\Delta x \rightarrow 0$$. The numerical routine becomes simple and efficient. We construct a Hamiltonian filling with matrix elements whose potential and kinetic energy contributions are given by Equations \ref{138} and \ref{139}. Then we diagonalize $$H^{DVR}$$, from which we obtain $$N$$ eigenvalues and the $$N$$ corresponding eigenfunctions.

## References

1. I. N. Levine, Quantum Chemistry, 5th ed. (Prentice Hall, Englewood Cliffs, NJ, 2000).
2. J. C. Light and T. Carrington, "Discrete-Variable Representations and Their Utilization" in Advances in Chemical Physics (John Wiley & Sons, Inc., 2007), pp. 263-310;
3. D. J. Tannor, Introduction to Quantum Mechanics: A Time-Dependent Perspective. (University Science Books, Sausilito, CA, 2007).
4. D. T. Colbert and W. H. Miller, A Novel Discrete Variable Representation for Quantum Mechanical Reactive Scattering Via the S‐Matrix Kohn Method, J. Chem. Phys. 96, 1982-1991 (1992).

1.5: Numerically Solving the Schrödinger Equation is shared under a CC BY-NC-SA 4.0 license and was authored, remixed, and/or curated by Andrei Tokmakoff via source content that was edited to conform to the style and standards of the LibreTexts platform; a detailed edit history is available upon request.