Numerov's method
# 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

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

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

    $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 $y_{n}=y(x_{n})$, $g_{n}=g(x_{n})$, $s_{n}=s(x_{n})$, and $h=x_{n+1}-x_{n}$.

### Nonlinear equations

For nonlinear equations of the form

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

the method gives

    $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 $f$ is linear in $y$ by setting $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 $R(r)$:

    ${\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:

    $u(r)=rR(r)\Rightarrow R(r)={\frac {u(r)}{r}},$
    ${\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)\Rightarrow {\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

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

or

    $-{\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

    $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:

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

## Derivation

We are given the differential equation

    $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, $y(x)$, around the 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}).$

Denoting the distance from $x$ to $x_{0}$ by $h=x-x_{0}$, 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}).$

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

    $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 $h$. If we want to take a step backwards, we replace every $h$ with $-h$ and get the expression for $y_{n-1}$:

    $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 $h$ experienced a sign change. By summing the two equations, we derive that

    $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 $y_{n+1}$ by substituting the expression given at the beginning, that is $y''_{n}=-g_{n}y_{n}+s_{n}$. To get an expression for the $y''''_{n}$ factor, we simply have to differentiate $y''_{n}=-g_{n}y_{n}+s_{n}$ twice and approximate it again in the same way we did this above:

    $y''''_{n}={\frac {d^{2}}{dx^{2}}}(-g_{n}y_{n}+s_{n}),$
    $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

    $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

    $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 $h^{6}$. It follows that the order of convergence (assuming stability) is 4.

