< >
 Ch. 4 Sec. 2

## 4.2 Initial Value Problems of Second Order

$\begin{eqnarray} \frac{\textstyle d^{2}y}{\textstyle dt^{2}}=b[y,dy/dt] \end{eqnarray}$

Example: generic equation of motion,
$\frac{\textstyle d^{2} r}{\textstyle dt^{2}} = \frac{\textstyle 1}{\textstyle m} K[r(t)]$

### 4.2.1 Verlet Method

Apply DDST,
$\frac{\textstyle d^{2}y}{\textstyle dt^{2}} |_{n} = \frac{\textstyle \delta^{2}y_{n}}{\textstyle (\Delta t)^{2}}+O[(\Delta t)^{2}]$

to find
 Verlet, classic: $y_{n+1}=2y_{n}-y_{n-1}+b_{n}(\Delta t)^{2}+O[(\Delta t)^{4}]$

Note 1: velocity $v \equiv \dot{y}$ does not appear explicitly. If $v$ is needed use the crude estimate

$\begin{eqnarray} v_{n} = \frac{\textstyle 1}{\textstyle 2\Delta t}[y_{n+1}-y_{n-1}]+O[(\Delta t)^{2}] \end{eqnarray}$

Note 2: A similar algorithm was devised in the early 1900s by Norwegian astrophysicist Carl Størmer; when investigating the polar aurora he analysed the motion of charged particles in the earth magnetic field. We may honor his contribution by calling the method "Størmer-Verlet algorithm".

Stability of the Verlet scheme:

Harmonic oscillator:
$\begin{eqnarray} y_{n+1}=2y_{n}-y_{n-1}-\omega_{0}^{2}y_{n}(\Delta t)^{2} \end{eqnarray}$

$\begin{eqnarray} g^{2}-(2-\alpha^{2})g+1=0 \end{eqnarray}$

with $\alpha\equiv \omega_{0} \Delta t$. The root
$\begin{eqnarray} g=(1-\frac{\textstyle \alpha^{2}}{\textstyle 2}) \pm \sqrt{\frac{\textstyle \alpha^{4}}{\textstyle 4}-\alpha^{2}} \end{eqnarray}$

is imaginary for $\alpha <2$, with $\vert g\vert^{2}=1$.

Equivalent formulations: Verlet leapfrog and Velocity Verlet.

 Verlet leapfrog: $\begin{eqnarray} v_{n+1/2} & = & v_{n-1/2} + b_{n} \Delta t \\ v_{n} & = & \frac{\textstyle 1}{\textstyle 2}(v_{n+1/2} + v_{n-1/2}) \;\;\; {\rm (if\;\; desired)} \\ y_{n+1} & = & y_{n}+ v_{n+1/2} \Delta t + O[(\Delta t)^{4}] \end{eqnarray}$

 Velocity Verlet (Swope): $\begin{eqnarray} y_{n+1} & = & y_{n}+ v_{n} \Delta t + b_{n} \frac{\textstyle (\Delta t)^{2}}{\textstyle 2} + O[(\Delta t)^{4}] \\ v_{n+1/2} & = & v_{n} + b_{n} \frac{\textstyle \Delta t}{\textstyle 2} \\ {\rm Evaluation \; step\;\;} && y_{n+1} \rightarrow b_{n+1} \\ v_{n+1} & = & v_{n+1/2} + b_{n+1} \frac{\textstyle \Delta t}{\textstyle 2} \end{eqnarray}$

### 4.2.2 Predictor-Corrector Method for 2nd order ODE

Predictor step:

In $d^{2}y/dt^{2}=b(t)$ replace the function $b(t)$ by a NGB polynomial and integrate twice.

$\begin{eqnarray} \dot{y}^{P}_{n+1}\Delta t-\dot{y}_{n}\Delta t &=& (\Delta t)^{2} \left[ b_{n} + \frac{\textstyle 1}{\textstyle 2} \nabla b_{n} + \frac{\textstyle 5}{\textstyle 12} \nabla^{2} b_{n} + \frac{\textstyle 3}{\textstyle 8} \nabla^{3} b_{n} + \dots \right] \\ && \\ y^{P}_{n+1}\!-\!y_{n}\!-\!\dot{y}_{n} \Delta t &=& \frac{\textstyle (\Delta t)^{2}}{\textstyle 2} \left[ b_{n} + \frac{\textstyle 1}{\textstyle 3} \nabla b_{n} + \frac{\textstyle 1}{\textstyle 4} \nabla^{2} b_{n} + \frac{\textstyle 19}{\textstyle 90} \nabla^{3} b_{n} + \dots \right] \end{eqnarray}$

A specific predictor of order $k$ is found by using terms up to order $\nabla^{k-2}b_{n}$. Thus the predictor of third order reads

$\begin{eqnarray} \dot{y}^{P}_{n+1}\Delta t-\dot{y}_{n}\Delta t&=&(\Delta t)^{2} \left[ \frac{\textstyle 3}{\textstyle 2} b_{n} - \frac{\textstyle 1}{\textstyle 2} b_{n-1} \right] + O[(\Delta t)^{4}] \\ y^{P}_{n+1}-y_{n}-\dot{y}_{n}\Delta t &=& \frac{\textstyle (\Delta t)^{2}}{\textstyle 2} \left[ \frac{\textstyle 4}{\textstyle 3} b_{n} - \frac{\textstyle 1}{\textstyle 3} b_{n-1} \right] + O[(\Delta t)^{4}] \end{eqnarray}$

For a compact notation we define the vector
$\begin{eqnarray} b_{k} &\equiv& \{ b_{n},b_{n-1},\dots b_{n-k+2} \}^{T} \end{eqnarray}$

and the coefficient vectors $c_{k}$ and $d_{k}$. Then

 Predictor of order k for second order DE: $\begin{eqnarray} \dot{y}^{P}_{n+1}\Delta t-\dot{y}_{n}\Delta t & = & (\Delta t)^{2} c_{k} \cdot b_{k}+O[(\Delta t)^{k+1}] \\ y^{P}_{n+1}-y_{n}-\dot{y}_{n}\Delta t & = & \frac{\textstyle (\Delta t)^{2}}{\textstyle 2} d_{k} \cdot b_{k} +O[(\Delta t)^{k+1}] \end{eqnarray}$

The first few vectors $c_{k}$, $d_{k}$ are given by
$\begin{eqnarray} c_{2}=1 \;\;\;\;\; && \;\; d_{2}=1 \\ && \\ c_{3}=\left(\begin{array}{r} 3/2 \\ -1/2 \end{array} \right) \;\;\; && \;\; d_{3}=\left(\begin{array}{r} 4/3 \\-1/3 \end{array} \right) \\ && \\ c_{4}=\left(\begin{array}{r} 23/12 \\-16/12 \\5/12 \end{array} \right) \; && \;\; d_{4}=\left(\begin{array}{r} 19/12 \\-10/12 \\3/12 \end{array} \right) \\ && \\ c_{5}=\left(\begin{array}{r} 55/24 \\-59/24 \\37/24 \\-9/24 \end{array}\right) \; && \;\; d_{5}=\left(\begin{array}{r} 323/180 \\-264/180 \\159/180 \\-38/180 \end{array} \right) \end{eqnarray}$

Evaluation step:

Insert the preliminary result $y_{n+1}^{P},\dot{y}_{n+1}^{P}$ in the physical law for $b[y,\dot{y}]$:

$\begin{eqnarray} b_{n+1}^{P} &\equiv& b\left[ y_{n+1}^{P},\dot{y}_{n+1}^{P}\right] \end{eqnarray}$

Corrector step:

Insert $b_{n+1}^{P}$ in a NGB formula centered on $t_{n+1}$ and re-integrate twice:

$\begin{eqnarray} \dot{y}_{n+1}\Delta t - \dot{y}_{n}\Delta t& = &(\Delta t)^{2} \left[ b_{n+1}^{P} - \frac{\textstyle 1}{\textstyle 2} \nabla b_{n+1} - \frac{\textstyle 1}{ \textstyle12} \nabla^{2} b_{n+1} - \frac{\textstyle 1}{\textstyle 24} \nabla^{3} b_{n+1} - \dots \right] \\ && \\ y_{n+1} - y_{n} - \dot{y}_{n}\Delta t & = & \frac{\textstyle (\Delta t)^{2}}{\textstyle 2} \left[ b_{n+1}^{P} - \frac{\textstyle 2}{\textstyle 3} \nabla b_{n+1} - \frac{\textstyle 1}{\textstyle 12} \nabla^{2} b_{n+1} - \frac{\textstyle 7}{\textstyle 180}\nabla^{3} b_{n+1} - \dots \right] \end{eqnarray}$

Defining the vector
$\begin{eqnarray} b_{k}^{P} &\equiv& \{ b_{n+1}^{P},b_{n},\dots b_{n-k+3} \}^{T} \end{eqnarray}$

and coefficient vectors $e_{k}$, $f_{k}$, we write

 Corrector of order k for second-order DE: $\begin{eqnarray} \dot{y}_{n+1}\Delta t-\dot{y}_{n}\Delta t & = & (\Delta t)^{2} e_{k} \cdot b_{k}^{P}+O[(\Delta t)^{k+1}] \\ y_{n+1}-y_{n}-\dot{y}_{n}\Delta t & = & \frac{\textstyle (\Delta t)^{2}}{\textstyle 2} f_{k} \cdot b_{k}^{P}+O[(\Delta t)^{k+1}] \end{eqnarray}$

The first few coefficient vectors are
$\begin{eqnarray} e_{2}=1 \;\;\;\;\; && \;\; f_{2}=1 \\ && \\ e_{3}=\left(\begin{array}{r} 1/2 \\1/2 \end{array}\right) && \;\; f_{3}=\left(\begin{array}{r} 1/3 \\2/3 \end{array}\right) \\ && \\ e_{4}=\left(\begin{array}{r} 5/12 \\8/12 \\-1/12 \end{array}\right) \; && \;\; f_{4}=\left(\begin{array}{r} 3/12 \\10/12 \\-1/12 \end{array}\right) \\ && \\ e_{5}=\left(\begin{array}{r} 9/24 \\19/24 \\-5/24 \\1/24 \end{array}\right) \; && \;\; f_{5}=\left(\begin{array}{r} 38/180 \\171/180 \\-36/180 \\7/180 \end{array}\right) \end{eqnarray}$

### 4.3.3 Nordsieck Formulation of the PC Method

Instead of threading a NGB polynomial through preceding points, expand $y(t)$ about $t_{n}$. The resulting "Taylor predictor" may be written as

$\begin{eqnarray} y_{n+1}^{P} & = & y_{n} + y\,'_{n}\Delta t + y\,''_{n} \frac{\textstyle (\Delta t)^{2}}{\textstyle 2!} + y\,'''_{n} \frac{\textstyle (\Delta t)^{3}}{\textstyle 3!} + O[(\Delta t)^{4}] \\ {y\,'}^{P}_{n+1} \Delta t & = & \;\;\;\;\;\;\;\;\; y\,'_{n} \Delta t + y\,''_{n} (\Delta t)^{2} + y\,'''_{n} \frac{\textstyle (\Delta t)^{3}}{\textstyle 2!} + O[(\Delta t)^{4}] \\ {y\,''}^{P}_{n+1} \frac{\textstyle (\Delta t)^{2}}{\textstyle 2!} & = & \;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\; \;\;\;\;\; y\,''_{n} \frac{\textstyle (\Delta t)^{2}}{\textstyle 2!} + y\,'''_{n} \frac{\textstyle (\Delta t)^{3}}{\textstyle 2!} + O[(\Delta t)^{4}] \\ {y\,'''}^{P}_{n+1} \frac{\textstyle (\Delta t)^{3}}{\textstyle 3!} & = & \;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\; \;\;\;\;\;\;\;\; \;\;\;\;\;\;\; y\,'''_{n}\frac{\textstyle (\Delta t)^{3}}{\textstyle 3!} + O[(\Delta t)^{4}] \end{eqnarray}$

Defining the vector

$\begin{eqnarray} z_{n} \equiv \left( \begin{array}{c} y_{n} \\ y\,'_{n} \Delta t \\ y\,''_{n} \frac{\textstyle (\Delta t)^{2}}{\textstyle 2!} \\ \vdots \end{array} \right) \end{eqnarray}$

and the (Pascal triangle) matrix

$\begin{eqnarray} A & \equiv & \left( \begin{array}{ccccc} 1 & 1 & 1 & 1 & \dots \\ 0 & 1 & 2 & 3 & \dots \\ 0 & 0 & 1 & 3 & \dots \\ & \ddots & \ddots & 1 & \ddots \\ & & & & \ddots \end{array} \right) \end{eqnarray}$

we have

$\begin{eqnarray} z_{n+1}^{P} = A \cdot z_{n} \end{eqnarray}$

Now evaluate: $\Longrightarrow$

Evaluation step:

Insert $z_{n+1}^{P}$ in the force law: $b_{n+1}^{P} \equiv b[y_{n+1}^{P}, {y\,'}^{P}_{n+1}]$

Corrector step:

Define the deviation $\gamma \equiv [b_{n+1}^{P} - {y\,''}^{P}_{n+1}] \frac{\textstyle (\Delta t)^{2}}{\textstyle 2}$ and write the corrector as

$\begin{eqnarray} z_{n+1} = z_{n+1}^{P} + \gamma c \end{eqnarray}$

with an optimized coefficient vector $c$. The first few vectors are

$\begin{eqnarray} c = \left( \begin{array}{c} 1/6 \\5/6 \\1 \\1/3 \end{array} \right) \,,\;\;\; \left( \begin{array}{c} 19/120 \\3/4 \\1 \\1/2 \\1/12 \end{array} \right)\,,\;\;\; \left( \begin{array}{c} 3/20 \\251/360 \\1 \\11/18 \\1/6 \\1/60 \end{array} \right)\,,\dots \end{eqnarray}$

- Self-starting

Stability: Again, between explicit (bad) and implicit (good).

### 4.3.4 Runge-Kutta Method for 2nd order ODE

 4th order RK for velocity independent forces: Let $b=b(y)$; then $\begin{eqnarray} b_{1} & = & b \left[ y_{n} \right] \\ && \\ b_{2} & = & b \left[y_{n}+\dot{y}_{n} \frac{\textstyle \Delta t}{\textstyle 2} \right] \\ && \\ b_{3} & = & b \left[ y_{n}+\dot{y}_{n}\frac{\textstyle \Delta t}{\textstyle 2} +b_{1}\frac{\textstyle (\Delta t)^{2}}{\textstyle 4} \right] \\ && \\ b_{4} & = & b \left[ y_{n}+\dot{y}_{n}\Delta t +b_{2}\frac{\textstyle (\Delta t)^{2}}{\textstyle 2} \right] \\ && \\ \dot{y}_{n+1}=\dot{y}_{n}\!&\!+\!&\!\frac{\textstyle \Delta t}{\textstyle 6}[b_{1} +2b_{2}+2b_{3}+b_{4}] +O[(\Delta t)^{5}] \\ y_{n+1}=y_{n}\!&\!+\!&\!\dot{y}_{n}\Delta t + \frac{\textstyle (\Delta t)^{2}}{\textstyle 6} [b_{1}+b_{2}+b_{3}]+O[(\Delta t)^{5}] \end{eqnarray}$

 4th order RK for velocity dependent forces: If $b=b(y,\dot{y})$, use $\begin{eqnarray} b_{1} &=& b [ y_{n},\dot{y}_{n}] \\ && \\ b_{2} &=& b \left[ y_{n}+\dot{y}_{n}\frac{\textstyle \Delta t}{\textstyle 2} +b_{1}\frac{\textstyle (\Delta t)^{2}}{\textstyle 8},\dot{y}_{n} +b_{1}\frac{\textstyle \Delta t}{\textstyle 2} \right] \\ && \\ b_{3} &=& b \left[ y_{n}+\dot{y}_{n}\frac{\textstyle \Delta t}{\textstyle 2} +b_{1} \frac{\textstyle (\Delta t)^{2}}{\textstyle 8},\dot{y}_{n} +b_{2}\frac{\textstyle \Delta t}{\textstyle 2} \right] \\ && \\ b_{4} &=& b \left[ y_{n}+\dot{y}_{n} \Delta t +b_{3}\frac{\textstyle (\Delta t)^{2}}{\textstyle 2}, \dot{y}_{n}+b_{3} \Delta t \right] \\ && \\ && \\ \dot{y}_{n+1}=\dot{y}_{n} & + & \frac{\textstyle \Delta t}{\textstyle 6} [b_{1}+2b_{2}+2b_{3}+b_{4}] +O[(\Delta t)^{5}] \\ y_{n+1}=y_{n} & + & \dot{y}_{n}\Delta t + \frac{\textstyle (\Delta t)^{2}}{\textstyle 6} [b_{1}+b_{2}+b_{3}]+O[(\Delta t)^{5}] \end{eqnarray}$

- Self-starting

But: Repeated evaluation of the acceleration $b(y)$ in a single time step

### 4.2.5 Symplectic Algorithms

[For some background information on Lagrange's and Hamilton's formulation of mechanics, see
here .]

Hamiltonian dynamics is an important field of application for ODE algorithms. One simple performance test is the conservation of total mechanical energy.

Another quantity that should be conserved in Hamiltonian dynamics is the symplectic form. A number of algorithms have been devised which by construction conserve this quantity.

In a classical system with $M$ degrees of freedom, let the complete set of (generalized) coordinates be denoted by $q$, and the conjugate momenta by $p$. Then

$\frac{\textstyle d q}{\textstyle dt} = \nabla_{p} H( q, p) \;\;\;\;\;\;\; \frac{\textstyle d p}{\textstyle dt} = -\nabla_{q} H( q, p)$

where $H(q, p)$ is the Hamiltonian.

Rewrite this by defining the phase space vector $z \equiv \{ q, p \}$:

$\frac{\textstyle d z}{\textstyle dt} = J \cdot \nabla_{z} H(z)$

where
$J \equiv \left(\begin{array}{cc} 0 & I \\ - I & 0 \end{array} \right)$

is the "symplectic matrix". ("symplectic" means "intertwined").

The canonical transformation

$z(t_{0}) \; \Longrightarrow \; z(t)$

conserves not only the energy (= numerical value of the Hamiltonian), but also the symplectic form

$s(z_{1}, z_{2}) \equiv z_{1}^{T} \cdot J \cdot z_{2}$

EXAMPLE: Harmonic oscillator:
$H( z) \equiv H(q,p) = \frac{\textstyle k}{\textstyle 2} q^{2} + \frac{\textstyle p^{2}}{\textstyle 2m}$

Canonical transformation:
$z(t) = \left(\begin{array}{r} q \\ \\p \end{array} \right) = \left( \begin{array}{cc} \cos \omega t & \frac{\textstyle 1}{\textstyle m \omega } \sin \omega t \;\; \\ \\ -m \omega \sin \omega t&\cos \omega t \end{array} \right) \cdot \left( \begin{array}{r} q(0) \\ \\ p(0) \end{array} \right) \equiv A \cdot z(0)$

Symplectic structure: Let $\{q_{1}(0), p_{1}(0) \}$ and $\{q_{2}(0), p_{2}(0) \}$ be two different initial conditions; then

$\begin{eqnarray} s(q_{1}(0) \dots p_{2}(0)) & \equiv & (q_{1}(0), p_{1}(0)) \cdot \left(\begin{array}{cc} 0 & 1 \\ - 1 & 0 \\ \end{array} \right) \cdot \left(\begin{array}{r} q_{2}(0) \\ p_{2}(0) \\ \end{array} \right) \\ && \\ & = & q_{1}(0) p_{2}(0) - p_{1}(0) q_{2}(0) \end{eqnarray}$

Geometric interpretation: $s$ is the area of a parallelogram defined by the two initial state vectors $z_{1,2}$.
$s$ is indeed constant:
$\begin{eqnarray} s(z_{1}(t),z_{2}(t)) & = & z_{1}^{T}(t) \cdot J \cdot z_{2}(t) \\ && \\ & = & z_{1}^{T}(0) \cdot A^{T} \cdot J \cdot A \cdot z_{2}(0) \\ && \\ & = & z_{1}^{T}(0) \cdot J \cdot z_{2}(0) \end{eqnarray}$

This means that $A^{T} \cdot J \cdot A = J$.

Let us check the performance of the Euler-Cauchy scheme:

$z(t) = \left(\begin{array}{r} q \\ \\ p \end{array} \right) = \left(\begin{array}{cc} 1 & \frac{\textstyle \Delta t}{\textstyle m} \; \\ \\-m \omega^{2} \Delta t & 1 \; \end{array} \right) \cdot \left(\begin{array}{r} q(0) \\ \\ p(0) \end{array} \right) \equiv E \cdot z(0)$

The EC algorithm enhances the symplectic form by a factor $1+ (\omega \Delta t)^{2}$ at each time step. It is therefore not a symplectic algorithm.

Symplectic integrator of fourth order (Neri and Candy-Rozmus):

 Symplectic algorithm of fourth order: Let the Hamiltonian be separable in terms of coordinates and momenta: $H( q, p)=U(q)+T(p)$. For the derivatives of $H$ we use the notation $F(q) \equiv - \nabla_{q} U(q) , \;\;\;\; P(q) \equiv \nabla_{p} T(p)$ The state at time $t$ is given by $\{ q_{0},p_{0}\}$. For $i=1$ to $4$ do $p_{i}=p_{i-1}+b_{i} F(q_{i-1})\Delta t , \;\;\;\;\; q_{i}=q_{i-1}+a_{i} P(p_{i})\Delta t$ where $\begin{eqnarray} a_{1}=a_{4} & = & (2+2^{1/3}+2^{-1/3})/6 \\ a_{2}=a_{3} & = & (1-2^{1/3}-2^{-1/3})/6 \\ b_{1}&=&0 \\ b_{2}=b_{4} & = & 1/(2-2^{1/3}) \\ b_{3} & = & 1/(1-2^{2/3}) \end{eqnarray}$ The state at time $t_{n+1}$ is $\{ q_{4},p_{4}\}$.

Third-order scheme (Ruth): same structure as Candy's algorithm, but with coefficients

$\begin{eqnarray} (a_{1},a_{2},a_{3}) & = & (2/3,-2/3,1) \\ (b_{1},b_{2},b_{3}) & = & (7/24, 3/4, -1/24) \end{eqnarray}$

Of the algorithms discussed previously only the Størmer-Verlet formula is symplectic.

A little theory:
For non-integrable Hamiltonians there can be no algorithm that conserves both energy and symplectic structure. But symplectic integrators conserve a Hamiltonian function that is different from, but close to, the given Hamiltonian (Yoshida.) Therefore they display no long-time energy trend. (Contrast: for example, Runge-Kutta integrators have good short-time accuracy but show a regularly increasing deviation in energy.)

Example: the first-order symplectic algorithm

$q_{n+1}=q_{n}+P(p_{n})\Delta t , \;\;\;\; p_{n+1}=p_{n}+F(q_{n+1})\Delta t$

exactly conserves
$\tilde{H} \equiv H + H_{1}\Delta t + H_{2}(\Delta t)^{2}+ H_{3}(\Delta t)^{3}+ \dots$

where
$H_{1}=\frac{\textstyle 1}{\textstyle 2}H_{p}H_{q} , \;\;\; H_{2}=\frac{\textstyle 1}{\textstyle 12}(H_{pp}H^{2}_{q}+H_{qq}H^{2}_{p}) , \;\;\; H_{3}=\frac{\textstyle 1}{\textstyle 12}H_{pp}H_{qq}H_{p}H_{q}\;\dots$

($H_{q} =\nabla_{q}H$ etc.) To take the harmonic oscillator, the perturbed Hamiltonian

$\tilde{H}=H_{ho} + \frac{\textstyle \omega^{2} \Delta t}{\textstyle 2}pq$

is conserved exactly.

A similar first order symplectic algorithm, known as the Euler-Cromer formula,

$p_{n+1}=p_{n}+F(q_{n})\Delta t , \;\;\;\; q_{n+1}=q_{n}+P(p_{n+1})\Delta t$

exactly conserves the perturbed Hamiltonian $\tilde{H}=H_{ho} - \frac{\textstyle \omega^{2} \Delta t}{\textstyle 2}pq$
When applied to oscillator-like equations of motion it is a definite improvement over the (unstable) Euler-Cauchy method.

Apply the (non-symplectic) RK method and the (symplectic) Størmer-Verlet algorithm (or the Candy procedure) to the one-body Kepler problem with elliptic orbit. Perform long runs to assess the long-time performance of the integrators. (For RK the orbit should eventually spiral down towards the central mass, while the symplectic procedures should only give rise to a gradual precession of the perihelion.)

### 4.2.6 Numerov's Method

Mostly used for a certain kind of initial value problems (IVPs) that arise when boundary value problems (BVPs) are rewritten as IVPs. Many BVP have the general form

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

with given boundary values $y(x_{1})$ and $y(x_{2})$.

Example: one-dimensional Poisson equation

$\frac{\textstyle d^{2}\phi}{\textstyle dx^{2}}= -\rho(x)$

with $\phi$ given at $x_{1}$ and $x_{2}$. Thus $g(x)=0$ and $s(x)=-\rho(x)$.

Let us assume we have, instead of $y(x_{1})$ and $y(x_{2})$, the full set of initial values $y(x_{1})$ and $y'(x_{1})$.

Divide the interval $[x_{1},x_{2}]$ into sub-intervals of length $\Delta x$ and at each intermediate point $x_{n}$ expand $y(x)$ into a power series. Adding the Taylor formulae for $y_{n+1}$ and $y_{n-1}$ we find

$y_{n+1}=2y_{n}-y_{n-1}+{y_{n}}''(\Delta x)^{2}+y_{n}^{(4)} \frac{\textstyle (\Delta x)^{4}}{\textstyle 12}+O[(\Delta x)^{6}]$

Inserting ${y_{n}}'' = -g_{n} y_{n}+s_{n}$ we have
$y_{n+1}=2y_{n}-y_{n-1}+(\Delta x)^{2}[-g_{n}y_{n}+s_{n}] +\frac{\textstyle (\Delta x)^{4}}{\textstyle 12}y_{n}^{(4)} +O[(\Delta x)^{6}]$

The fourth derivative $y^{(4)}$ is approximated by

$\begin{eqnarray} y_{n}^{(4)} & = & \frac{\textstyle d^{2}y''}{\textstyle dx^{2}} |_{n} = \frac{\textstyle d^{2}(-gy+s)}{\textstyle dx^{2}} |_{n} \approx \frac{\textstyle 1}{\textstyle (\Delta x)^{2}} \delta_{n}^{2}(-gy+s) = \\ && \\ &=& \frac{\textstyle 1}{\textstyle (\Delta x)^{2}} [-g_{n+1}y_{n+1}+2g_{n}y_{n}-g_{n-1}y_{n-1} + s_{n+1}-2s_{n}+s_{n-1}] \end{eqnarray}$

In this way we find Numerov's formula

$\begin{eqnarray} y_{n+1}[1+\frac{\textstyle (\Delta x)^{2}}{\textstyle 12}g_{n+1}] &=& 2y_{n}[1-\frac{\textstyle 5}{\textstyle 12}(\Delta x)^{2}g_{n}]- y_{n-1}[1+\frac{\textstyle (\Delta x)^{2}}{\textstyle 12}g_{n-1}]+ \\ && \\ && +\frac{\textstyle (\Delta x)^{2}}{\textstyle 12} [s_{n+1}+10s_{n}+s_{n-1}]+O[(\Delta x)^{6}] \end{eqnarray}$