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