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 the Russian astronomer 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 {5h^{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 unidimensional Schrödinger equation for arbitrary potentials. An example of which is solving the radial equation for a spheically symmetric potential. In this example, after seperating the variables and analytically solving the angular equation, we are left with the following equation of the radial :$R(r)$ function:

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

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

${\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}}R(r)$

$-{\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 potential that is analogous to the classical effective potential.

$V_{eff}(r)=V(r)+{\frac {\hbar ^{2}}{2m}}{\frac {l(l+1)}{r^{2}}}$

$L^{2}=l(l+1)\hbar ^{2}$

$V_{eff}(r)=V(r)+{\frac {L^{2}}{2mr^{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.

$\left({\frac {d^{2}}{dr^{2}}}+{\frac {2m}{\hbar ^{2}}}(E-V_{eff}(r))\right)u(r)=0$

$a(x)={\frac {2m}{\hbar ^{2}}}(E-V_{eff}(r))$

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}}{dx^{2}}}\left[a(x)y(x)\right]$

We replace the second derivative ${\frac {d^{2}}{dx^{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}-2a_{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 {5h^{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.

References

• Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard (1993), Solving ordinary differential equations I: Nonstiff problems, Berlin, New York: Springer-Verlag, ISBN 978-3-540-56670-0.
  This book includes the following references:
• Numerov, Boris Vasil'evich (1924), "A method of extrapolation of perturbations", Monthly Notices of the Royal Astronomical Society, 84: 592–601, Bibcode:1924MNRAS..84..592N, doi:10.1093/mnras/84.8.592{{citation}}: CS1 maint: unflagged free DOI (link).
• Numerov, Boris Vasil'evich (1927), "Note on the numerical integration of d²x/dt² = f(x,t)", Astronomische Nachrichten, 230: 359–364, Bibcode:1927AN....230..359N, doi:10.1002/asna.19272301903.