# 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.

## The method

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

${\displaystyle {\frac {d^{2}y}{dx^{2}}}=-g(x)y(x)+s(x).}$

In it, three values of ${\displaystyle y_{n-1},y_{n},y_{n+1}}$ taken at three equidistant points ${\displaystyle x_{n-1},x_{n},x_{n+1}}$ are related as follows:

${\displaystyle y_{n+1}\left(1+{\frac {h^{2}}{12}}g_{n+1}\right)=2y_{n}\left(1-{\frac {5h^{2}}{12}}g_{n}\right)-y_{n-1}\left(1+{\frac {h^{2}}{12}}g_{n-1}\right)+{\frac {h^{2}}{12}}(s_{n+1}+10s_{n}+s_{n-1})+{\mathcal {O}}(h^{6}),}$

where ${\displaystyle y_{n}=y(x_{n})}$, ${\displaystyle g_{n}=g(x_{n})}$, ${\displaystyle s_{n}=s(x_{n})}$, and ${\displaystyle h=x_{n+1}-x_{n}}$.

### Nonlinear equations

For nonlinear equations of the form

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

the method gives

${\displaystyle y_{n+1}-2y_{n}+y_{n-1}={\frac {h^{2}}{12}}(f_{n+1}+10f_{n}+f_{n-1})+{\mathcal {O}}(h^{6}).}$

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

## Application

In numerical physics the method is used to find solutions of the unidimensional Schrödinger equation for arbitrary potentials. An example of which is solving the radial equation for a spherically symmetric potential. In this example, after separating the variables and analytically solving the angular equation, we are left with the following equation of the radial function ${\displaystyle R(r)}$:

${\displaystyle {\frac {d}{dr}}\left(r^{2}{\frac {dR}{dr}}\right)-{\frac {2mr^{2}}{\hbar ^{2}}}(V(r)-E)R(r)=l(l+1)R(r).}$

This equation can be reduced to the form necessary for the application of Numerov's method with the following substitution:

${\displaystyle u(r)=rR(r),}$
${\displaystyle R(r)={\frac {u(r)}{r}},}$
${\displaystyle {\frac {dR}{dr}}={\frac {1}{r}}{\frac {du}{dr}}-{\frac {u(r)}{r^{2}}}={\frac {1}{r^{2}}}\left(r{\frac {du}{dr}}-u(r)\right),}$
${\displaystyle {\frac {d}{dr}}\left(r^{2}{\frac {dR}{dr}}\right)={\frac {du}{dr}}+r{\frac {d^{2}u}{dr^{2}}}-{\frac {du}{dr}}=r{\frac {d^{2}u}{dr^{2}}}.}$

And when we make the substitution, the radial equation becomes

${\displaystyle r{\frac {d^{2}u}{dr^{2}}}-{\frac {2mr}{\hbar ^{2}}}(V(r)-E)u(r)={\frac {l(l+1)}{r}}u(r),}$

or

${\displaystyle -{\frac {\hbar ^{2}}{2m}}{\frac {d^{2}u}{dr^{2}}}+\left(V(r)+{\frac {\hbar ^{2}}{2m}}{\frac {l(l+1)}{r^{2}}}\right)u(r)=Eu(r),}$

which is equivalent to the one-dimensional Schrödinger equation, but with the modified effective potential

${\displaystyle V_{\text{eff}}(r)=V(r)+{\frac {\hbar ^{2}}{2m}}{\frac {l(l+1)}{r^{2}}}=V(r)+{\frac {L^{2}}{2mr^{2}}},\quad L^{2}=l(l+1)\hbar ^{2}.}$

This equation we can proceed to solve the same way we would have solved the one-dimensional Schrödinger equation. We can rewrite the equation a little bit differently and thus see the possible application of Numerov's method more clearly:

${\displaystyle {\frac {d^{2}u}{dr^{2}}}=-{\frac {2m}{\hbar ^{2}}}(E-V_{\text{eff}}(r))u(r),}$
${\displaystyle g(r)={\frac {2m}{\hbar ^{2}}}(E-V_{\text{eff}}(r)),}$
${\displaystyle s(r)=0.}$

## Derivation

We are given the differential equation

${\displaystyle y''(x)=-g(x)y(x)+s(x).}$

To derive the Numerov's method for solving this equation, we begin with the Taylor expansion of the function we want to solve, ${\displaystyle y(x)}$, around the point ${\displaystyle x_{0}}$:

${\displaystyle 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}).}$

Denoting the distance from ${\displaystyle x}$ to ${\displaystyle x_{0}}$ by ${\displaystyle h=x-x_{0}}$, we can write the above equation as

${\displaystyle 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}).}$

If we evenly discretize the space, we get a grid of ${\displaystyle x}$ points, where ${\displaystyle h=x_{n+1}-x_{n}}$. By applying the above equations to this discreet space, we get a relation between the ${\displaystyle y_{n}}$ and ${\displaystyle y_{n+1}}$:

${\displaystyle y_{n+1}=y_{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}).}$

Computationally, this amounts to taking a step forward by an amount ${\displaystyle h}$. If we want to take a step backwards, we replace every ${\displaystyle h}$ with ${\displaystyle -h}$ and get the expression for ${\displaystyle y_{n-1}}$:

${\displaystyle y_{n-1}=y_{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}).}$

Note that only the odd powers of ${\displaystyle h}$ experienced a sign change. By summing the two equations, we derive that

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

We can solve this equation for ${\displaystyle y_{n+1}}$ by substituting the expression given at the beginning, that is ${\displaystyle y''_{n}=-g_{n}y_{n}+s_{n}}$. To get an expression for the ${\displaystyle y''''_{n}}$ factor, we simply have to differentiate ${\displaystyle y''_{n}=-g_{n}y_{n}+s_{n}}$ twice and approximate it again in the same way we did this above:

${\displaystyle y''''_{n}={\frac {d^{2}}{dx^{2}}}(-g_{n}y_{n}+s_{n}),}$
${\displaystyle h^{2}y''''_{n}=-g_{n+1}y_{n+1}+s_{n+1}+2g_{n}y_{n}-2s_{n}-g_{n-1}y_{n-1}+s_{n-1}+{\mathcal {O}}(h^{4}).}$

If we now substitute this to the preceding equation, we get

${\displaystyle y_{n+1}-2y_{n}+y_{n-1}={h^{2}}(-g_{n}y_{n}+s_{n})+{\frac {h^{2}}{12}}(-g_{n+1}y_{n+1}+s_{n+1}+2g_{n}y_{n}-2s_{n}-g_{n-1}y_{n-1}+s_{n-1})+{\mathcal {O}}(h^{6}),}$

or

${\displaystyle y_{n+1}\left(1+{\frac {h^{2}}{12}}g_{n+1}\right)-2y_{n}\left(1-{\frac {5h^{2}}{12}}g_{n}\right)+y_{n-1}\left(1+{\frac {h^{2}}{12}}g_{n-1}\right)={\frac {h^{2}}{12}}(s_{n+1}+10s_{n}+s_{n-1})+{\mathcal {O}}(h^{6}).}$

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