Numerov's method

From HandWiki

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

[math]\displaystyle{ \frac{d^2 y}{dx^2} = - g(x) y(x) + s(x). }[/math]

In it, three values of [math]\displaystyle{ y_{n-1}, y_n, y_{n+1} }[/math] taken at three equidistant points [math]\displaystyle{ x_{n-1}, x_n, x_{n+1} }[/math] are related as follows:

[math]\displaystyle{ y_{n+1} \left(1 + \frac{h^2}{12} g_{n+1}\right) = 2 y_n \left(1 - \frac{5 h^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} + 10 s_n + s_{n-1}) + \mathcal{O}(h^6), }[/math]

where [math]\displaystyle{ y_n = y(x_n) }[/math], [math]\displaystyle{ g_n = g(x_n) }[/math], [math]\displaystyle{ s_n = s(x_n) }[/math], and [math]\displaystyle{ h = x_{n+1} - x_n }[/math].

Nonlinear equations

For nonlinear equations of the form

[math]\displaystyle{ \frac{d^2 y}{dx^2} = f(x,y), }[/math]

the method gives

[math]\displaystyle{ y_{n+1} - 2 y_n + y_{n-1} = \frac{h^2}{12} (f_{n+1} + 10 f_n + f_{n-1}) + \mathcal{O} (h^6). }[/math]

This is an implicit linear multistep method, which reduces to the explicit method given above if [math]\displaystyle{ f }[/math] is linear in [math]\displaystyle{ y }[/math] by setting [math]\displaystyle{ f(x,y) = - g(x) y(x) + s(x) }[/math]. It achieves order-4 accuracy (Hairer Nørsett).

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 [math]\displaystyle{ R(r) }[/math]:

[math]\displaystyle{ \frac{d}{dr} \left(r^2 \frac{dR}{dr} \right) - \frac{2 m r^2}{\hbar^2} (V(r) - E) R(r) = l (l + 1) R(r). }[/math]

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

[math]\displaystyle{ u(r) = r R(r) \Rightarrow R(r) = \frac{u(r)}{r}, }[/math]
[math]\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) \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}. }[/math]

And when we make the substitution, the radial equation becomes

[math]\displaystyle{ r \frac{d^2 u}{dr^2} - \frac{2 m r}{\hbar^2} (V(r) - E) u(r) = \frac{l (l + 1)}{r} u(r), }[/math]

or

[math]\displaystyle{ -\frac{\hbar^2}{2 m} \frac{d^2 u}{dr^2} + \left(V(r) + \frac{\hbar^2}{2 m} \frac{l (l + 1)}{r^2} \right) u(r) = E u(r), }[/math]

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

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

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:

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

Derivation

We are given the differential equation

[math]\displaystyle{ y''(x) = - g(x) y(x) + s(x). }[/math]

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

[math]\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). }[/math]

Denoting the distance from [math]\displaystyle{ x }[/math] to [math]\displaystyle{ x_0 }[/math] by [math]\displaystyle{ h = x - x_0 }[/math], we can write the above equation as

[math]\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). }[/math]

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

[math]\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). }[/math]

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

[math]\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). }[/math]

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

[math]\displaystyle{ y_{n+1} - 2 y_n + y_{n-1} = h^2 y''_n + \frac{h^4}{12}y''''_n + \mathcal{O}(h^6). }[/math]

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

[math]\displaystyle{ y''''_n = \frac{d^2}{d x^2} (-g_n y_n + s_n), }[/math]
[math]\displaystyle{ h^2 y''''_n = -g_{n+1} y_{n+1} + s_{n+1} + 2 g_n y_n - 2 s_n - g_{n-1} y_{n-1} + s_{n-1} + \mathcal{O}(h^4). }[/math]

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

[math]\displaystyle{ y_{n+1} - 2 y_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} + 2 g_n y_n - 2 s_n - g_{n-1} y_{n-1} + s_{n-1}) + \mathcal{O}(h^6), }[/math]

or

[math]\displaystyle{ y_{n+1} \left(1 + \frac{h^2}{12} g_{n+1} \right) - 2 y_n \left(1 - \frac{5 h^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} + 10 s_n + s_{n-1}) + \mathcal{O}(h^6). }[/math]

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

References



Retrieved from "https://handwiki.org/wiki/index.php?title=Numerov%27s_method&oldid=2641923"