# Numerov's method

Numerov'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 Boris Vasil'evich Numerov.

## The method

The Numerov method can be used to solve differential equations of the form

$\left ( \frac{d^2}{dx^2} + a(x) \right ) y(x) = 0$

The function $a(x)$ is sampled in the interval [a..b] at equidistant positions $x_n$. Starting from function values at two consecutive samples $x_{n-1}$ and $x_n$ the remaining function values can be calculated as

$y_{n+1} = \frac {\left( 2-\frac{5 h^2}{6} a_n \right) y_n - \left( 1+\frac{h^2}{12}a_{n-1} \right)y_{n-1}}{1+\frac{h^2}{12}a_{n+1}}$

where $a_n=a(x_n)$ and $y_n=y(x_n)$ are the function values at the positions $x_n$ and $h=x_n-x_{n-1}$ is the distance between two consecutive samples.

### Nonlinear equations

For nonlinear equations of the form

$\frac{d^2}{dx^2} y = f(x,y)$

the method is given by

$y_{n+1} = 2y_n - y_{n-1} + \tfrac{1}{12} h^2 (f_{n+1} + 10f_n + f_{n-1}).$

This is an implicit linear multistep method, which reduces to the explicit method given above if f is linear in y by setting $f(x,y) = -a(x) y(x)$. It achieves order 4 (Hairer, Nørsett & Wanner 1993, §III.10).

## Application

In numerical physics the method is used to find solutions of the radial Schrödinger equation for arbitrary potentials.

$\left [ -{\hbar^2 \over 2\mu} \left ( \frac{1}{r} {\partial^2 \over \partial r^2} r- {l(l+1) \over r^2} \right ) + V(r) \right ] R(r) = E R(r)$

The above equation can be rewritten in the form

$\left [ {\partial^2 \over \partial r^2} - {l(l+1) \over r^2} + { 2\mu \over \hbar^2} \left( E - V(r)\right) \right ] u(r) = 0$

with $u(r) = r R(r)$. If we compare this equation with the defining equation of the Numerov method we see

$a(x) = \frac{2\mu}{\hbar^2} \left(E - V(x) \right) - \frac{l(l+1)}{x^2}$

and thus can numerically solve the radial Schrödinger equation.

## Derivation

Start with the Taylor expansion of $y(x)$ about a point $x_0$:

$y(x) = y(x_0) + (x-x_0)y'(x_0) + \frac{(x-x_0)^2}{2!}y''(x_0) + \frac{(x-x_0)^3}{3!}y'''(x_0) + \frac{(x-x_0)^4}{4!}y''''(x_0) + \frac{(x-x_0)^5}{5!}y'''''(x_0) + \mathcal{O} (h^6)$

Denote the distance from $x$ to $x_0$ by $h=x-x_0$ and, noting that this means $x=x_0+h$, we can write the above equation as

$y(x_0+h) = y(x_0) + hy'(x_0) + \frac{h^2}{2!}y''(x_0) + \frac{h^3}{3!}y'''(x_0) + \frac{h^4}{4!}y''''(x_0) + \frac{h^5}{5!}y'''''(x_0) + \mathcal{O} (h^6)$

Computationally, this amounts taking a step forward by an amount h. If we want to take a step backwards, replace every h with -h for the equation of $y(x_0-h)$:

$y(x_0-h) = y(x_0) - hy'(x_0) + \frac{h^2}{2!}y''(x_0) - \frac{h^3}{3!}y'''(x_0) + \frac{h^4}{4!}y''''(x_0) - \frac{h^5}{5!}y'''''(x_0) + \mathcal{O} (h^6)$

Note that only the odd powers of h experienced a sign change. On an evenly spaced grid, the nth site on a computational grid corresponds to position $x_n$ if the step-size between grid points are of length $h$ (hence h should be small for the computation to be accurate). This means we have sampling points $(x_{n-1},y_{n-1})$ and $(x_{n+1},y_{n+1})$. Taking the equations for $y_{n-1}$ and $y_{n+1}$ from continuous space to discrete space, we see that

$y_{n+1} = y(x_n+h) = y(x_n) + hy'(x_n) + \frac{h^2}{2!}y''(x_n) + \frac{h^3}{3!}y'''(x_n) + \frac{h^4}{4!}y''''(x_n) + \frac{h^5}{5!}y'''''(x_n) + \mathcal{O} (h^6)$
$y_{n-1} = y(x_n-h) = y(x_n) - hy'(x_n) + \frac{h^2}{2!}y''(x_n) - \frac{h^3}{3!}y'''(x_n) + \frac{h^4}{4!}y''''(x_n) - \frac{h^5}{5!}y'''''(x_n) + \mathcal{O} (h^6)$

The sum of those two equations gives

$y_{n-1} + y_{n+1} = 2y_n + {h^2}y''_n + \frac{h^4}{12}y''''_n + \mathcal{O} (h^6)$

We solve this equation for $y''_n$ and replace it by the expression $y''_n = -a_n y_n$ which we get from the defining differential equation.

$h^2 a_n y_n = 2y_n-y_{n-1} - y_{n+1} + \frac{h^4}{12}y''''_n + \mathcal{O} (h^6)$

We take the second derivative of our defining differential equation and get

$y''''(x) = - \frac{d^2}{d x^2} \left[ a(x) y(x) \right]$

We replace the second derivative $\frac{d^2}{d x^2}$ with the second order difference quotient and insert this into our equation for $a_n y_n$ (note that we take the mixed forward and backward finite difference, not the double forward difference or the double backward difference)

$h^2 a_n y_n = 2y_n-y_{n-1} - y_{n+1} - \frac{h^4}{12} \frac{a_{n-1} y_{n-1} -2 a_{n} y_{n} + a_{n+1} y_{n+1}}{h^2} + \mathcal{O} (h^6)$

We solve for $y_{n+1}$ to get

$y_{n+1} = \frac {\left( 2-\frac{5 h^2}{6} a_n \right) y_n - \left( 1+\frac{h^2}{12}a_{n-1} \right)y_{n-1}}{1+\frac{h^2}{12}a_{n+1}} + \mathcal{O} (h^6).$

This yields Numerov's method if we ignore the term of order $h^6$. It follows that the order of convergence (assuming stability) is 4.