Numerov Algorithm

Numerov’s method

Numerov’s method¹ (also called Cowell’s method) is a numerical method to solve ordinary differential equations of second order in which the first-order term does not appear. It is a fourth-order linear multistep method. The method is implicit, but can be made explicit if the differential equation is linear. Numerov’s method was developed by the Russian astronomer Boris Vasil’evich Numerov.

Given the differential equation of this form

\begin{matrix} (1) & \frac{d^{2} y}{d x^{2}} = - g (x) y (x) + s (x); x \in [a, b] \end{matrix}

where $g (x)$ and $s (x)$ are continuous function on the domain $[a, b]$ . To derive the Numerov’s method, we first discretize the interval $[a, b]$ using equally spaced grid, where the grid space is $h = x_{n + 1} - x_{n}$ . We then proceed with the Taylor expansion of the function $y (x)$ around the point $x_{n}$ :

\begin{aligned} y (x_{n + 1}) & = y (x_{n}) + y^{'} (x_{n}) h + \frac{y^{″} (x_{n})}{2!} h^{2} + \frac{y^{‴} (x_{n})}{3!} h^{3} \\ (2) & + \frac{y^{⁗} (x_{n})}{4!} h^{4} + \frac{y^{'''''} (x_{n})}{5!} h^{5} + \frac{y^{''''''} (x_{n})}{6!} h^{6} + \dots \end{aligned}

Similarly, we have

\begin{aligned} y (x_{n - 1}) & = y (x_{n}) - y^{'} (x_{n}) h + \frac{y^{″} (x_{n})}{2!} h^{2} - \frac{y^{‴} (x_{n})}{3!} h^{3} \\ (3) & + \frac{y^{⁗} (x_{n})}{4!} h^{4} - \frac{y^{'''''} (x_{n})}{5!} h^{5} + \frac{y^{''''''} (x_{n})}{6!} h^{6} + \dots \end{aligned}

By summing the two equations, we have

\begin{matrix} (4) & y_{n + 1} - 2 y_{n} + y_{n - 1} = h^{2} y_{n}^{″} + \frac{h^{4}}{12} y_{n}^{⁗} + O (h^{6}) \end{matrix}

To get an expression for $y_{n}^{⁗}$ , remember that $y_{n}^{″} = - g_{n} y_{n} + s_{n} \equiv z_{n}$ . By repeating the same procedure we did above

\begin{aligned} h^{2} y_{n}^{⁗} & = z_{n + 1} - 2 z_{n} + z_{n - 1} + O (h^{4}) \\ = - g_{n + 1} y_{n + 1} + s_{n + 1} \\ + 2 g_{n} y_{n} - 2 s_{n} \\ (5) & - g_{n - 1} y_{n - 1} + s_{n - 1} + O (h^{4}) \end{aligned}

If we substitute this into the proceeding equation, we get

\begin{aligned} y_{n + 1} - 2 y_{n} + y_{n - 1} & = h^{2} (- g_{n} y_{n} + s_{n}) \\ (6) & + \frac{h^{2}}{12} (- g_{n + 1} y_{n + 1} + s_{n + 1} + 2 g_{n} y_{n} - 2 s_{n} - g_{n - 1} y_{n - 1} + s_{n - 1}) + O (h^{6}) \end{aligned}

or equivalently

\begin{matrix} (7) & y_{n + 1} (1 + \frac{h^{2}}{12} g_{n + 1}) - 2 y_{n} (1 - \frac{5 h^{2}}{12} g_{n}) + y_{n - 1} (1 + \frac{h^{2}}{12} g_{n - 1}) = \frac{h^{2}}{12} (s_{n + 1} + 10 s_{n} + s_{n - 1}) + O (h^{6}) \end{matrix}

This yields the Numerov’s method if we ignore the term of order $h^{6}$ .

\begin{aligned} y_{n + 1} & = \frac{2 y_{n} (1 - \frac{5 h^{2}}{12} g_{n}) - y_{n - 1} (1 + \frac{h^{2}}{12} g_{n - 1}) + \frac{h^{2}}{12} (s_{n + 1} + 10 s_{n} + s_{n - 1})}{1 + \frac{h^{2}}{12} g_{n + 1}} \\ (8) & = \frac{(12 - 10 f_{n}) y_{n} - y_{n - 1} f_{n - 1} + \frac{h^{2}}{12} (s_{n + 1} + 10 s_{n} + s_{n - 1})}{f_{n + 1}} \end{aligned}

where $f_{n} = 1 + \frac{h^{2}}{12} g_{n}$ . It is easy to show that

\begin{matrix} (9) & y_{n - 1} = \frac{(12 - 10 f_{n}) y_{n} - y_{n + 1} f_{n + 1} + \frac{h^{2}}{12} (s_{n + 1} + 10 s_{n} + s_{n - 1})}{f_{n - 1}} \end{matrix}

With this relation, we can now find the value of $y (x_{n + 1})$ given two previous values $y (x_{n})$ and $y (x_{n - 1})$ . All we have to do to find an approximate solution is to find the values $y (x_{1})$ and $y (x_{2})$ at the first two points $x_{1}$ and $x_{2}$ , then we use the Numerov’s method to find the value of all following points in the interval $[a, b]$ .

Case 1. Hydrogen-like atom

The Schrödinger equation for hydrogen-like atom is ²

\begin{matrix} (10) & [- \frac{ℏ^{2}}{2 m} \nabla^{2} - \frac{Z e^{2}}{4 π ε_{0} r}] ψ (r) = E ψ (r) \end{matrix}

Here, we assume that there is only one electron and thus neglecting the Coulomb interaction between the electrons. By variable separation, the wavefunction can be written as $ψ (r) = \frac{u (r)}{r} Y_{l m} (\hat{r})$ , where $Y_{l m} (\hat{r})$ is the spherical harmonics and $u (r)$ follows the equation

\begin{matrix} (11) & - \frac{ℏ^{2}}{2 m} u^{″} + [\frac{ℏ^{2}}{2 m} \frac{l (l + 1)}{r^{2}} - \frac{Z e^{2}}{4 π ε_{0} r} - E] u (r) = 0 \end{matrix}

Introducing the dimensionless variables

\begin{aligned} (12) & x = \frac{r}{a_{0}}; & a_{0} = \frac{4 π ε_{0} ℏ^{2}}{m e^{2}} \approx 0.529 Å \\ (13) & ϵ = \frac{E}{E_{0}}; & E_{0} = \frac{ℏ^{2}}{2 m a_{0}^{2}} = 1 Ry \approx 13.60 eV \end{aligned}

where $a_{0}$ is the Bohr radius and $E_{0}$ is the Rydberg energy.

The equation for $u (r)$ then becomes

\begin{matrix} (14) & u^{″} (x) - [\frac{l (l + 1)}{x^{2}} - \frac{2 Z}{x} - ϵ] u (x) = 0 \end{matrix}

The boundary conditions are

\begin{matrix} (15) & u (0) = 0; u (\infty) = 0 \end{matrix}

1. From known energy eigenvalue to the wavefunctions

As is well known that the energy levels of hydrogen-like atoms are

\begin{matrix} (16) & E_{n} = - \frac{Z^{2}}{n^{2}}; n = 1, 2, 3, \dots \end{matrix}

Let’s first show that given the correct energy levels, we can obtain the corresponding wavefunctions.

CLICK TO SHOW THE CODE!

#!/usr/bin/env python

import numpy as np
from scipy.integrate import odeint, simps

def radial_schrod_deriv(u, r, l, E, Z=1):
    '''
    The radial Schrödinger equation for hydrogen-like atom is

        u''(x) - [(l(l+1)/x^2) - (2Z/x) - E]u(x) = 0 

    Since it is a second-order equation, i.e. 
    
        y'' + g(x)y(x) = 0

    To use scipy.integrate.odeint, we can turn this equation into two
    first-order equation by defining a new dependent variables
    
        y'(x) = z(x)
        z'(x) = -g(x)y(x)

    Then we can solve this system of ODEs using "odeint" wit list. 
    '''
    
    y, z = u

    return np.array([
        z, 
        ((l*(l+1) / r**2) - (2*Z / r) - E) * y
    ])


def radial_wfc_scipy(r0, n=1, l=0, Z=1, direction='F', du=0.1):
    '''
    Get the radial wavefunction by integrating the equation with
    scipy.integrate.odeint.
    '''

    assert direction.upper() in ['F', 'B']

    E = -Z**2 / n**2

    # forward integration 
    if direction.upper() == 'F':
        ur = odeint(radial_schrod_deriv, [0.0, du], r0, args=(l, E, Z))[:,0]

    # back integration 
    else:
        ur = odeint(radial_schrod_deriv, [0.0, -du], r0[::-1], args=(l, E, Z))[:,0][::-1]

    # normalization
    ur /= np.sqrt(simps(ur**2, x=r0))

    return ur


def radial_wfc_numerov(r0, n=1, l=0, Z=1, du=0.001):
    '''
    Numerov algorithm
    
                  [12 - 10f(n)]*y(n) - y(n-1)*f(n-1)
        y(n+1) = ------------------------------------
                               f(n+1)

    where
        
        f(n) = 1 + (h**2 / 12)*g(n)

        g(n) = [E + (2*Z / x) - l*(l+1) / x**2]
    '''

    ur = np.zeros(r0.size)
    fn = np.zeros(r0.size)

    E      = -Z**2 / n**2
    ur[-1] = 0.0
    ur[-2] = du

    dr  = r0[1] - r0[0]
    h12 = dr**2 / 12.

    gn = (E + 2*Z / r0 - l*(l+1) / r0**2)
    fn = 1. + h12 * gn

    for ii in range(r0.size - 3, -1, -1):
        ur[ii] = (12 - 10*fn[ii+1]) * ur[ii+1] - \
                 ur[ii+2] * fn[ii+2]
        ur[ii] /= fn[ii]

    # normalization
    ur /= np.sqrt(simps(ur**2, x=r0))

    return ur

In the code above, we have defined two integration method to integrate the radial Schrödinger equation

radial_wfc_scipy uses odeint implemented in scipy.³
radial_wfc_numerov uses the Numerov’s method.
We implemented forward/backward integration scheme in radial_wfc_scipy function, where forward integration means integrating from $u (0)$ and backward integraton starts from $u (\infty)$ .
Of course we can not integrate from $r = \infty$ , we used a large enough cutoff $r_{c}$ and assumes $u (r_{c}) = 0$ .

CLICK TO SHOW THE CODE!

################################################################################
# Define the radial grid
r0 = np.linspace(1E-10, 20, 500)

Z = 1
l = 0
n = 1

# forward integration from u(0)
ur1 = radial_wfc_scipy(r0, n, l, direction='F')
# backward integration from u(\infty)
ur2 = radial_wfc_scipy(r0, n, l, direction='B')
# backward integration from u(\infty) with Numerov method
ur3 = radial_wfc_numerov(r0, n, l)

################################################################################
import matplotlib as mpl
import matplotlib.pyplot as plt
plt.style.use('ggplot')

fig = plt.figure(
    figsize=(7.2, 3.0),
    dpi=300
)
axes = [plt.subplot(1,2,ii+1) for ii in range(2)]

################################################################################
ax = axes[0]

ax.plot(
    r0, ur1,
    ls='none',
    ms=4, marker='o', mfc='white', mew=1.0, 
    zorder=1,
    label=r'Forward integration',
)

ax.plot(
    r0, ur2,
    ls='none',
    ms=4, marker='o', mfc='white', mew=1.0, 
    zorder=1,
    label=r'Backward integration',
)
ax.plot(r0, 2*r0*np.exp(-r0), lw=1.0, color='cyan', zorder=2, label=r'Exact: $2r\cdot e^{-r}$')

ax.set_xlabel(r'$r$ [$a_0$]',   labelpad=5)
ax.set_ylabel(r'$u(r)$ [a.u.]', labelpad=5)

ylim = list(ax.get_ylim())
ylim[1] = 1.0
ax.set_ylim(ylim)
ax.legend(loc='best', fontsize='small')

#################################################################################
ax = axes[1]

ax.semilogy(
    r0, 2*r0*np.exp(-r0) - ur2,
    ls='-', alpha=0.8,
    zorder=1,
    label=r'scipy.integrate.odeint',
)

ax.semilogy(
    r0, 2*r0*np.exp(-r0) - ur3,
    ls=':', alpha=0.8,
    zorder=1,
    label=r'Numerov',
)
ax.legend(loc='lower right', fontsize='small')

ax.set_xlabel(r'$r$ [$a_0$]',   labelpad=5)
ax.set_ylabel(r'Error', labelpad=5)

plt.tight_layout()
plt.savefig('fig1.png')
plt.show()

**Figure 1.** Left: forward and backward integration of the radial Schrödinger using scipy.integrate.odeint. The solid cyan line shows the exact solution. Right: the difference between the exact solution and the numerical results. Two different integration schemes are used, i.e. odeint in scipy and Numerov method.

We then try to get the radial wavefunction of hydrogen ( $Z = 1$ ) with quantum number $n = 1, l = 0$ . Some remarks here:

It turns out that the forward integration of the equation is not stable, as can be seen in the left panel of Figure 1.
It is better to use logarithmic mesh for radial variable rather than linear. Radial functions need smaller number of points in logarithmic mesh. For example, use np.logspace(1E-6, 2, 500) to generate the radial grid and pass to the radial_wfc_scipy function. Note that Numerov’s method only support linear mesh, as can be seen in the method section.

One can show that the Numerov’s method also works fine for higher excited states. For example, below we compare the radial wavefunction $u_{n l} (r)$ for $n = 1, 2, 3$ obtained from Numerov method with the exact solutions. The script can be found here.

**Figure 2.** The lowest few radial wavefunctions $u_{n l} (r)$ with $n = 1, 2, 3$ .

2. Unknown energy levels

If we don NOT know the correct energy levels, we can utilize the boundary condition to get the energy. Recall that the boundary condition is

\begin{matrix} (17) & u (0) = 0; u (r_{c}) = 0 \end{matrix}

where $r_{c}$ is large enough so that $u (r_{c}) = 0$ . If $E_{0}$ corresponds to the correct energy level, when we integrate the radial wavefunction $u (r)$ from $r_{c}$ backward using Numerov’s method, we should get $u (0) = 0$ . A deviation of the energy from $E_{0}$ will result in $u (0) \neq 0$ . Therefore we can use the so called “shooting” method to get the correct energy. The basic procedure is

Start wite a guess energy, e.g. $E_{1}$ . In practice, the guess energy should be smaller than the smallest potential energy. In our case, $E_{1}$ should be smaller than $- Z^{2}$ . With the guess energy $E_{1}$ , integrate the equation and get the value of the function at $r = 0$ , which we will denote as $u_{1}$ . Meanwhile, Set another energy $E_{2} = E_{1}$ .
Increase $E_{2}$ by an amount $δ E$ and get a new energy $E_{2} = E_{2} + δ E$ . Integrate the equation to get the value at $r = 0$ , which we will denote as $u_{2}$ .
Go back to step 2 if $u_{1} * u_{2} > 0$ .
At this step, we should have the correct energy enclosed in $[E_{1}, E_{2}]$ . Use root finding method, e.g. scipy.optimize.brentq to get the correct energy.

CLICK TO SHOW THE CODE!

#!/usr/bin/env python

import numpy as np
from scipy.integrate import odeint, simps
from scipy.optimize import brentq


def radial_wfc_at0(E, r0, n=1, l=0, Z=1, du=0.01):
    '''
    Numerov algorithm
    
                  [12 - 10f(n)]*y(n) - y(n-1)*f(n-1)
        y(n+1) = ------------------------------------
                               f(n+1)

    where
        
        f(n) = 1 + (h**2 / 12)*g(n)

        g(n) = [E + (2*Z / x) - l*(l+1) / x**2]

    here, we use reverse integration from the other end.
    '''
    ur = np.zeros_like(r0)
    ur[-1] = 0.0
    ur[-2] = du

    dr  = r0[1] - r0[0]
    h12 = dr**2 / 12.
    gn  = E + 2*Z / r0 - l*(l+1) / r0**2
    fn  = 1. + h12 * gn

    for ii in range(gn.size - 3, -1, -1):
        ur[ii] = (12 - 10*fn[ii+1]) * ur[ii+1] - \
                 ur[ii+2] * fn[ii+2]
        ur[ii] /= fn[ii]

    # normalization
    ur /= np.sqrt(simps(ur**2, x=r0))

    # now extrapolate the wavefunction to u(0)
    # the first derivative equals at ur[0]
    # (u0 - ur[0]) / (0 - r0[0]) = (ur[1] - ur[0]) / (r0[1] - r0[0])
    u0 = ur[0] + (ur[1] - ur[0]) * (0 - r0[0]) / dr

    return u0


if __name__ == "__main__":
    r0 = np.linspace(1E-6, 30, 3000)

    E_lower = E_upper = -1.1
    dE      = 0.10
    u1      = radial_wfc_at0(E_lower, r0, n=1)

    while True:
        E_upper += dE
        u2 = radial_wfc_at0(E_upper, r0, n=1)
        if u1 * u2 < 0:
            break

    E = brentq(radial_wfc_at0, E_lower, E_upper, args=(r0, 1, 0, 1, 0.001)) 
    print(E)

Numerov Algorithm

Numerov’s method

Case 1. Hydrogen-like atom

1. From known energy eigenvalue to the wavefunctions

2. Unknown energy levels

References

Trending Tags

Numerov Algorithm

Numerov’s method

Case 1. Hydrogen-like atom

1. From known energy eigenvalue to the wavefunctions

2. Unknown energy levels

References

Further Reading

Numerical Solutions to the Time-dependent Schrödinger Equation

Phonopy: Rutile TiO2

Brillouin Zone Integration: Linear Tetrahedron Method

Trending Tags

Phonopy: Rutile TiO₂