Numerov's method

From Wikipedia, the free encyclopedia
Jump to: navigation, search

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[edit]

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[edit]

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[edit]

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[edit]

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.

References[edit]

External links[edit]