% Except otherwise noted,
% this work is © 2009 Dmitri Fedorov, under CC-BY-SA license, version 3+.

\chapter{Ordinary differential equations}
\section{Introduction}
Ordinary differential equations (ODE) are generally defined as
differential equations in one variable where the highest order derivative
enters linearly.  Such equations invariably arise in many different
contexts throughout mathematics and science as soon as changes in the
phenomena at hand are considered, usually with respect to variations of
certain parameters.

Systems of ordinary differential equations can be generally reformulated
as systems of first-order ordinary differential equations,
	\begin{equation}\label{ode:eq:o}
\ve{y}'(x)=\ve{f}(x,\ve{y}) \;,
	\end{equation}
where $\ve{y}^\prime\doteq{d\ve{y}}/{dx}$, and the variables $\ve{y}$
and the {\em right-hand side function} $\ve{f}(x,\ve{y})$ are understood
as column-vectors.  For example, a second order differential equation
in the form
	\begin{equation}
u'' = g(x,u,u')
	\end{equation}
can be rewritten as a system of two first-order equations,
	\begin{equation}
\left\{\begin{array}{l} y'_1=y_2 \\ y'_2=g(x,y_1,y_2)
\end{array}\right. \;,
	\end{equation}
using the variable substitution $y_1=u$, $y_2=u'$.

In practice ODEs are usually supplemented with boundary conditions
which pick out a certain class or a unique solution of the ODE.
In the following we shall mostly consider {\em initial value
problems}\index{initial value problem}: ODE with the boundary condition
in the form of an initial condition at a given point $a$,
	\begin{equation}
\ve{y}(a) = \ve{y}_{0} \:.
	\end{equation}
The problem then is to find the value of the solution $\ve{y}$ at some
other point $b$.  Finding a solution to an ODE is often referred to as
{\em integrating} the ODE.

An ODE integration algorighm typically advances the solution from the
initial point $a$ to the final point $b$ in a number of discrete steps
	\begin{equation}
\{x_0\doteq a,x_1,\dots,x_{n-1},x_n\doteq b\}.
	\end{equation}
An efficient algorithm tries to integrate an ODE using as few steps
as possible under the constraint of the given accuracy goal.  For this
purpose the algorthm should continuously adjust the step-size during the
integration, using few larger steps in the regions where the solution
is smooth and perhaps many smaller steps in more treacherous regions.

Typically, an adaptive step-size ODE integrator is implemented as two
routines.  One of them---called {\em driver}---monitors the local errors
and tolerances and adjusts the step-sizes.  To actually perform a step
the driver calls a separate routine---the {\em stepper}---which advances
the solution by one step, using one of the many available algorithms, and
estimates the local error.  The GNU Scientific Library, GSL, implements
about a dozen of different steppers and a tunable adaptive driver.

In the following the chapter describes several of the popular driving
algorithms and stepping methods for solving initial value ODE problems.

\section{Error estimate}

In an adaptive step-size algorithm the stepping routine must provide
an estimate of the integration error, upon which the driver bases
its strategy to determine the optimal step-size for a user-specified
accuracy goal.

A stepping method is generally characterized by its {\em order\:}:
a method has order $p$ if it can integrate exactly an ODE where the
solution is a polynomial of order $p$. In other words, for small $h$
the error of the order-$p$ method is $O(h^{p+1})$.

For sufficiently small steps the error $\delta y$ of an integration
step for a method of a given order $p$ can be estimated by
comparing the solution $\ve{y}_\mathrm{full\_step}$ obtained with
one full-step integration with a potentially more precise solution,
$\ve{y}_\mathrm{two\_half\_steps}$, obtained with two consecutive
half-step integrations,
	\begin{equation}
\delta \ve{y} = 
\frac{\ve{y}_\mathrm{full\_step}-\ve{y}_\mathrm{two\_half\_steps}}
{2^{p+1}-1} \:.
	\end{equation}
where $p$ is the order of the algorithm used.

Indeed, if the step-size $h$ is small, we can assume
	\begin{eqnarray}
\delta \ve{y}_\mathrm{full\_step} &=& Ch^{p+1}	\:,\\ \delta
\ve{y}_\mathrm{two\_half\_steps} &=& C\left(\frac{h}{2}\right)^{p+1} \:,
	\end{eqnarray}
where $\delta \ve{y}_\mathrm{full\_step}$ and $\delta
\ve{y}_\mathrm{two\_half\_steps}$ are the errors of the full-step and
two half-steps integrations, and $C$ is an unknown constant.  The two
can be combined as
	\begin{eqnarray}
\ve{y}_\mathrm{full\_step}-\ve{y}_\mathrm{two\_half\_steps} &=& \delta
\ve{y}_\mathrm{full\_step} -\delta \ve{y}_\mathrm{two\_half\_steps}
\nonumber\\ &=& C\left(\frac{h}{2}\right)^{p+1}(2^{p+1}-1) \:,
	\end{eqnarray}
from which it immedliately follows that
	\begin{equation}
C\left(\frac{h}{2}\right)^{p+1} =
\frac{\ve{y}_\mathrm{full\_step}-\ve{y}_\mathrm{two\_half\_steps}}
{2^{p+1}-1} \:.
	\end{equation}
One has, of course, to take the potentially more precise
$\ve{y}_\mathrm{two\_half\_steps}$ as the solution $\ve{y}$.  Its error
is then given as
	\begin{equation}
\delta \ve{y} \doteq \delta \ve{y}_\mathrm{two\_half\_steps}
= C\left(\frac{h}{2}\right)^{p+1} =
\frac{\ve{y}_\mathrm{full\_step}-\ve{y}_\mathrm{two\_half\_steps}}
{2^{p+1}-1} \:,
	\end{equation}
which had to be demonstrated.  This prescription is often referred to
as the {\em Runge's principle}.

One drawback of the Runge's principle is that the full-step and the
two half-step calculations generally do not share evaluations of the
right-hand side function $\ve{f}(x,\ve{y})$, and therefore many extra
evaluations are needed to estimate the error.

An alternative prescription for error estimation is to make the same
step-size integration using two methods of {\em different orders},
with the difference between the two solutions providing the estimate of
the error.  If the lower order method mostly uses the same evaluations of
the right-hand side function---in which case it is called {\em embedded}
in the higher order method---the error estimate does not need additional
evaluations.

Predictor-corrector methods are naturally of embedded type: the
correction --- which generally increases the order of the method --- itself
can serve as the estimate of the error.


\section{Runge-Kutta methods}

Runge-Kutta methods are one-step methods which advance the solution over
the current step using only the information gathered from withing the
step itself.  The solution $\ve{y}$ is advanced from the point $x_i$
to $x_{i+1}=x_i+h$, where $h$ is the step-size, using a one-step
formula,
	\begin{equation}
\ve{y}_{i+1}=\ve{y}_i+h\ve{k},
	\end{equation}
where $\ve{y}_{i+1}$ is the approximation to $\ve{y}(x_{i+1})$,
and the value $\ve{k}$ is chosen such that the method integrates exactly
an ODE whose solution is a polynomial of the highest possible order.

The Runge-Kutta methods are distinguished by their {\em order\:}: a method
has order $p$ if it can integrate exactly an ODE where the solution is
a polynomial of order $p$. In other words, for small $h$ the error of
the order-$p$ method is $O(h^{p+1})$.

The first order Runge-Kutta method is the {\em Euler's
method}\index{Euler's method},
    \begin{equation}
\ve{k}=\ve{f}(x_{0},\ve{y}_{0})\:.
    \end{equation}

Second order Runge-Kutta methods advance the solution by an auxiliary
evaluation of the derivative, e.g. the {\em mid-point method},
	\begin{eqnarray}
\ve{k}_0&=&\ve{f}(x_0,\ve{y}_0)\;,\nonumber\\
\ve{k}_{1/2}&=&\ve{f}(x_{0}+\tfrac{1}{2}h,
\ve{y}_0+\tfrac{1}{2}h\ve{k}_0 )\;,\nonumber\\
\ve{k}&=&\ve{k}_{1/2} \;,
	\end{eqnarray}
or the {\em two-point method}, also called the {\em Heun's method}
	\begin{eqnarray}\label{ode:eq:2p}
\ve{k}_0 &=& \ve{f}(x_0,\ve{y}_0),\;\nonumber\\ \ve{k}_1
&=& \ve{f}(x_{0}+h, \ve{y}_0+h\ve{k}_0),\nonumber\\ \ve{k}
&=& \frac{1}{2}(\ve{k}_{0}+\ve{k}_{1}) \;.
	\end{eqnarray}

These two methods can be combined into a third order method,
\begin{equation}
\ve{k}=\frac{1}{6}\ve{k}_0+\frac{4}{6}\ve{k}_{1/2}
+\frac{1}{6}\ve{k}_1 \;.
\end{equation}

The most commont is the fourth-order method, which is called {\em RK4} or
simply {\em the Runge-Kutta method},
	\begin{eqnarray} 
\ve{k}_0 &=& \ve{f}(x_0, \ve{y}_0) \:,\nonumber\\
\ve{k}_1 &=& \ve{f}(x_0 + \tfrac{1}{2}h,
	\ve{y}_0 + \tfrac{1}{2}h\ve{k}_0) \:,\nonumber\\
\ve{k}_2 &=& \ve{f}(x_0 + \tfrac{1}{2}h,
	\ve{y}_0 + \tfrac{1}{2}h\ve{k}_1) \:,\nonumber\\
\ve{k}_3 &=& \ve{f}(x_0 + h,
	\ve{y}_0 + h\ve{k}_2 )\:,\nonumber\\
\ve{k}&=&
\tfrac{1}{6}\ve{k}_0+\tfrac{1}{3}\ve{k}_1
+\tfrac{1}{3}\ve{k}_2+\tfrac{1}{6}\ve{k}_3\:.
	\end{eqnarray}

A general Runge-Kutta method can be written as
	\begin{equation}
\ve{y}_{n+1} = \ve{y}_n + \sum_{i=1}^s b_i \ve{K}_i \:,
	\end{equation}
where
	\begin{eqnarray}
\ve{K}_1 &=& h\ve{f}(x_n, \ve{y}_n) \:,\nonumber\\
\ve{K}_2 &=& h\ve{f}(x_n+c_2h,\ve{y}_n+a_{21}\ve{K}_1) \:,\nonumber\\
\ve{K}_3 &=& h\ve{f}(x_n+c_3h,\ve{y}_n+a_{31}\ve{K}_1+a_{32}\ve{K}_2) \:,\\
&\vdots& \nonumber \\
\ve{K}_s &=& h\ve{f}(x_n+c_sh,\ve{y}_n+a_{s1}\ve{K}_1+a_{s2}\ve{K}_2+
\cdots+a_{s,s-1}\ve{K}_{s-1}) \:. \nonumber
	\end{eqnarray}
The upper case $\ve{K}_i$ are simply the lower case $\ve{k}_i$ multiplied
by the step-size $h$.

To specify a particular Runge-Kutta method one needs to provide the
coefficients $\{a_{ij} | 1\le j < i\le s\}$, $\{b_i | i=1..s\}$ and
$\{c_i | i=1..s\}$. The matrix $[a_{ij}]$ is called the Runge-Kutta
matrix, while the coefficients $b_i$ and $c_i$ are known as the weights
and the nodes.  These data are usually arranged in the so called
{\em Butcher's tableau}\index{Butcher's tableau},
	\begin{equation}
\begin{array}{c|ccccc}
0 & & & & & \\
c_2 & a_{21} & & & & \\
c_3 & a_{31} & a_{32} & & & \\
\vdots & \vdots & & \ddots & & \\
c_s & a_{s1} & a_{s2} & \cdots & a_{s,s-1} \\
\hline
 & b_1 & b_2 & \cdots & b_{s-1} & b_s
\end{array} \;.
	\end{equation}

For example, the Butcher's tableau for the RK4 method is
	\begin{equation}
\begin{array}{c|cccc}
0 & & & & \\
1/2 & 1/2 & & & \\
1/2 & 0 & 1/2 & & \\
1 & 0 & 0 & 1 &  \\
\hline
 & 1/6 & 1/3 & 1/3 & 1/6
\end{array} \;.
	\end{equation}

\subsection{Embeded methods with error estimates}
The embedded Runge-Kutta methods---in addition to advancing the solution
by one step---produce an estimate of the local error of the step.  This is
done by having two methods in the tableau, one with a certain order $p$
and another one with order $p-1$.  The difference bitween the two method
gives the the estimate of the local error.  The lower order method is
{\em embedded} in the higher order method, that is, it uses the same
$\ve{K}$-values. This allows a very effective estimate of the error.

The embedded lower order method is written as
	\begin{equation}
\ve{y}^*_{n+1} = \ve{y}_n + \sum_{i=1}^s b^*_i \ve{K}_i \:,
	\end{equation}
where $\ve{K}_i$ are the same as for the higher order method. The error
estimate is then given as
	\begin{equation}
\ve{e}_{n} = \ve{y}_{n+1} - \ve{y}^*_{n+1} =
\sum_{i=1}^s (b_i - b^*_i) \ve{K}_i \:.
	\end{equation}

The Butcher's tableau for this kind of method is extended by one row to
give the values of $b^*_i$.

The simplest embedded methods are Heun-Euler method,
	\begin{equation}
	\begin{array}{c|cc}
0 & & \\
1 & 1 & \\
\hline
 & 1/2 & 1/2 \\
 & 1 & 0 \\
\end{array} \;,
	\end{equation}
and midpoint-Euler method,
	\begin{equation}
	\begin{array}{c|cc}
0 & & \\
1/2 & 1/2 & \\
\hline
 & 0 & 1 \\
 & 1 & 0 \\
\end{array} \;,
	\end{equation}
which both combine methods of orders 2 and 1. An implementation of the
midpoint-Euler method in C is shown in Table~\ref{ode:t:me}.
	\begin{table}[htb]
\caption{Embedded Runge-Kutta midpoint-Euler stepper with error
estimate.\label{ode:t:me}} \lstinputlisting{c/ode/rkstep12.c}
	\end{table}

The {\em Bogacki-Shampine method}\index{Bogacki-Shampine
method}~\cite{bogacki-shampine} combines methods of orders 3 and 2,
	\begin{equation} \begin{array}{c|cccc}
0 & & & &\\ 1/2 & 1/2 &  &\\ 3/4 & 0 & 3/4 & &\\ 1 & 2/9 & 1/3 & 4/9
&\\ \hline
 & 2/9 & 1/3 & 4/9 & 0\\ & 7/24 & 1/4 & 1/3  & 1/8\\
\end{array} \;.
	\end{equation}
Bogacki and Shampine argue that their method has better stability
properties and actually outperforms higher order methods at lower
accuracy goal calculations.  This method has the FSAL---first same as
last---property: the  value $\ve{k}_4$ at one step equals $\ve{k}_1$ at
the next step; thus only three function evaluations are needed per step.
Table~\ref{ode:t:bs} shows a simple implementation which does not utilise
this property for the sake of presentational clarity.
	\begin{table}[htb]
\caption{Embedded Runge-Kutta Bogacki-Shampine stepper} \label{ode:t:bs}
\lstinputlisting{c/ode/rkstep23.c}
	\end{table}

The Runge-Kutta-Fehlberg method~\cite{rkf} --- called {\em
RKF45} --- implemented in the renowned {\tt rkf45} Fortran routine, has two
methods of orders 5 and 4,
	\begin{equation}
	\begin{array}{c|cccccc}
0 & & & & & &\\
1/4 & 1/4 & & & & &\\
3/8 & 3/32 & 9/32 & & & &\\
12/13 & 1932/2197 & -7200/2197 & 7296/2197 & & &\\
1 & 439/216 & -8 & 3680/513 & -845/4104 & \\
1/2 & -8/27 & 2  & -3544/2565 & 1859/4104 & -11/40 & \\
\hline
 & 16/135 & 0 & 6656/12825 & 28561/56430 & -9/50 & 2/55 \\
 & 25/216 & 0 & 1408/2565 & 2197/4104 & -1/5 & 0 \\
	\end{array} \nonumber
	\end{equation}

\section{Multistep methods}
Multistep methods try to use the information about the function gathered
at the previous steps. They are generally not {\em self-starting} as
there are no previous points at the start of the integration. The first
step must be done with a one-step method like Runge-Kutta.

A number of multistep have been devised (and named after different
mathematicians); we shall only consider a few simple ones here to get
the idea of how it works.

\subsection{Two-step method}
Given the previous point, $(x_{i-1},\ve{y}_{i-1})$, in addition to the
current point $(x_i,\ve{y}_i$), the
sought function $\ve{y}$ can be approximated in the vicinity of the
point $x_i$ as
	\begin{equation}\label{eq:2}
\bar{\ve{y}}(x)=\ve{y}_i+\ve{y}'_i\cdot(x-x_i)
+\bar{\ve{c}}\cdot(x-x_i)^2,
	\end{equation}
where $\ve{y}'_i=\ve{f}(x_i,\ve{y}_i)$ and
the coefficient $\bar{\ve{c}}$ can be found from the condition
$\bar{\ve{y}}(x_{i-1})=\ve{y}_{i-1}$, which gives
	\begin{equation}
\bar{\ve{c}} =
\frac{\ve{y}_{i-1}-\ve{y}_{i}+\ve{y}'_{i}\cdot(x_{i}-x_{i-1})}
{(x_{i}-x_{i-1})^{2}} \:.
	\end{equation}
The value of the function at the next point, $x_{i+1}\doteq x_i+h$,
can now be estimated as $\bar{\ve{y}}(x_{i+1})$ from~(\ref{eq:2}).

The error of this second-order two-step stepper can be estimated by a
comparison with the first-order Euler's step, which is given by the
linear part of~(\ref{eq:2}). The correction term $\bar{\ve{c}}h^2$
can serve as the error estimate,
	\begin{equation}\label{}
\delta\ve{y} = \bar{\ve{c}}h^2 \:.
	\end{equation}

%\subsection{A three-step method}
%Including yet another point, $(x_{-2},\ve{y}_{-2})$, allows to
%further increase the order of the approximation~(\ref{eq:2}),
%	\begin{equation}\label{eq:-2}
%\bar{\bar{\ve{y}}}(x)=\bar{\ve{y}}(x)
%+\bar{\bar{\ve{d}}}\cdot(x-x_0)^2(x-x_{-1}) \:.
%	\end{equation}
%The coefficient $\bar{\bar{\ve{d}}}$ can be found from the condition
%$\bar{\bar{\ve{y}}}(x_{-2})=\ve{y}_{-2}$, which gives
%	\begin{equation}\label{}
%\bar{\bar{\ve{d}}} = \frac{\ve{y}_{-2}-\bar{\ve{y}}(x_{-2})}
%{(x_{-2}-x_0)^2(x_{-2}-x_{-1})} \:.
%	\end{equation}
%The error estimate at the point $x_1\doteq x_0+h$ is again estimated as
%the difference between the higher and the lower order steppers,
%	\begin{equation}\label{}
%\mathrm{err} \approx \|\bar{\bar{\ve{y}}}(x_1)-\bar{\ve{y}}(x_1)\|
%\:.
%	\end{equation}

\subsection{Two-step method with extra evaluation}
One can further increase the order of the
approximation~(\ref{eq:2}) by adding a third order term,
	\begin{equation}\label{eq:-2}
\bar{\bar{\ve{y}}}(x)=\bar{\ve{y}}(x)
+\bar{\bar{\ve{d}}}\cdot(x-x_i)^2(x-x_{i-1}) \:.
	\end{equation}
The coefficient $\bar{\bar{\ve{d}}}$ can be found from the matching
condition at a certain point inside the invervarl, for example at half-step,
	\begin{equation}\label{}
\bar{\bar{\ve{y}}}^\prime(x_{1/2})=
\ve{f}(x_{1/2},\bar{\ve{y}}(x_{1/2})) \doteq \bar{\ve f}_{1/2}\;,
	\end{equation}
where $x_{1/2}\doteq x_i+h/2$. This gives
	\begin{equation}\label{}
\bar{\bar{\ve d}}= \frac{\bar{\ve{f}}_{1/2}-\ve{y}'_i
-2\bar{\ve{c}}\cdot(x_{1/2}-x_i)}
{2(x_{1/2}-x_i)(x_{1/2}-x_{i-1})+(x_{1/2}-x_i)^{2}}.
	\end{equation}
The error estimate at the point $x_{i+1}\doteq x_0+h$ is again given as
the difference between the higher and the lower order methods,
	\begin{equation}\label{}
\delta\ve{y} = \bar{\bar{\ve{y}}}(x_{i+1})-\bar{\ve{y}}(x_{i+1}) \:.
	\end{equation}

\section{Predictor-corrector methods}
A predictor-corrector method uses extra iterations to improve the
solution.  It is an algorithm that proceeds in two steps.  First, the
predictor step calculates a rough approximation of $\ve{y}(x+h)$. Second,
the corrector step refines the initial approximation.  Aditionally
the corrector step can be repeated in the hope that this achieves an
even better approximation to the true solution.

For example, the two-point Runge-Kutta method~(\ref{ode:eq:2p})
is as actually a predictor-corrector method, as it first calculates
the {\em prediction} $\tilde{\ve{y}}_{i+1}$ for $\ve{y}(x_{i+1})$,
	\begin{equation}
\tilde{\ve{y}}_{i+1}=\ve{y}_{i}+h\ve{f}(x_{i},\ve{y}_{i})\:,
	\end{equation}
and then uses this prediction in a {\em correction} step,
	\begin{equation}
\check{\tilde{\ve{y}}}_{i+1}=
\ve{y}_{i}+h\frac12
\left(
\ve{f}(x_{i},\ve{y}_{i})
+\ve{f}(x_{i+1},\tilde{\ve{y}}_{i+1})
\right) \:.
	\end{equation}

\subsection{Two-step method with correction}
Similarly, one can use the two-step approximation (\ref{eq:2}) as a
predictor, and then improve it by one order with a correction step, namely
	\begin{equation}\label{eq:3}
\check{\bar{\ve{y}}}(x)=
\bar{\ve{y}}(x)+\check{\ve{d}}\cdot(x-x_i)^{2}(x-x_{i-1}).
	\end{equation}
The coefficient $\check{\ve{d}}$ can be found from the condition
$\check{\bar{\ve{y}}}'(x_{i+1})=\bar{\ve{f}}_{i+1}$, where
$\bar{\ve{f}}_{i+1}\doteq\ve{f}(x_{i+1},\bar{\ve{y}}(x_{i+1}))$,
	\begin{equation}
\check{\ve{d}}=
\frac{\bar{\ve{f}}_{i+1}-\ve{y}'_i-2\bar{\ve{c}}\cdot(x_{i+1}-x_i)}
{2(x_{i+1}-x_i)(x_{i+1}-x_{i-1})+(x_{i+1}-x_i)^{2}}.
	\end{equation}

Equation (\ref{eq:3}) gives a better estimate,
$\ve{y}_{i+1}=\check{\bar{\ve{y}}}(x_{i+1})$, of the sought function at
the point $x_{i+1}$.  In this context the formula (\ref{eq:2}) serves
as {\em predictor}, and (\ref{eq:3}) as {\em corrector}.  The difference
between the two gives an estimate of the error.

This method is equivalent to the two-step method with an extra evaluation
where the extra evaluation is done at the full step.

\section{Adaptive step-size control}

Let {\em tolerance} $\tau$ be the maximal accepted error consistent
with the required accuracy to be achieved in the integration of an ODE.
Suppose the inegration is done in $n$ steps of size $h_i$ such that
$\sum_{i=1}^n h_i=b-a$.  Under assumption that the errors at the
integration steps are random and statistically uncorrelated, the local
tolerance $\tau_i$ for the step $i$ has to scale as the square root of
the step-size,
	\begin{equation}
\tau_i=\tau \sqrt{\frac{h_i}{b-a}}.
	\end{equation}
Indeed, if the local error $e_i$ on the step $i$ is less than the local
tolerance, $e_i \le \tau_i$, the total error $E$ will be consistent with
the total tolerance $\tau$,
	\begin{equation}
E\approx\sqrt{\sum_{i=1}^n e_i^2} \le \sqrt{\sum_{i=1}^n \tau_i^2} =
\tau\sqrt{\sum_{i=1}^n \frac{h_i}{b-a}} = \tau\;.
	\end{equation}

The current step $h_i$ is accepted if the local error $e_i$ is
smaller than the local tolerance $\tau_i$, after which the next step
is attempted with the step-size adjusted according to the following
empirical prescription~\cite{gsl},
	\begin{equation}\label{ode:eq:h}
h_{i+1}=h_{i}\times \left(\frac{\tau_i}{e_i}\right)^\mathrm{Power}
\times\mathrm{Safety},
	\end{equation}
where $\mathrm{Power}\approx0.25$ and $\mathrm{Safety}\approx0.95$.

If the local error is larger than the local tolerance the step is rejected
and a new step is attempted with the step-size adjusted according to the
same prescription~(\ref{ode:eq:h}).

One simple prescription for the local tolerance $\tau_i$ and the local
error $e_i$ to be used in~(\ref{ode:eq:h}) is
	\begin{equation}
\tau_i=(\epsilon\|\ve{y}_i\|+\delta)\sqrt{\frac{h_i}{b-a}} \:,\;
e_i=\|\delta\ve{y}_i\| \:,
	\end{equation}
where $\delta$ and $\epsilon$ are the required absolute and relative
precision and $\delta\ve{y}_i$ is the estimate of the integration error
at the step~$i$.

A more elaborate prescription considers components of the
solution separately,
	\begin{equation}
(\tau_i)_k =
\big(\epsilon |(\ve{y}_i)_k| + \delta\big)\sqrt{\frac{h_i}{b-a}} 
\:,\; (\ve{e}_i)_k = |(\delta\ve{y}_i)_k|
%\;\Bigg|_{k=1..n}
\;,
	\end{equation}
where the index $k$ runs over the components of the solution.  In this
case the step acceptence criterion also becomes component-wise: the step
is accepted, if
	\begin{equation}
\forall k:\:(\ve{e}_i)_k<(\tau_i)_k \:.
	\end{equation}
The factor $\tau_i/e_i$ in the step adjustment formula~(\ref{ode:eq:h})
is then replaced by
	\begin{equation}
\frac{\tau_i}{e_i} \rightarrow \min_k\frac{(\tau_i)_k}{(\ve{e}_i)_k} \:.
	\end{equation}

Yet another refinement is to include the derivatives $\ve{y}'$ of the
solution into the local tolerance estimate, either overally,
	\begin{equation}
\tau_i=\Big(
\epsilon\alpha\|\ve{y}_i\|+\epsilon\beta\|\ve{y}'_i\|
+\delta\Big)\sqrt{\frac{h_i}{b-a}} \:,\;
	\end{equation}
or commponent-wise,
	\begin{equation}
(\tau_i)_k = \Big(
\epsilon\alpha|(\ve{y}_i)_k|+\epsilon\beta|(\ve{y}'_i)_k|
+\delta \Big)\sqrt{\frac{h_i}{b-a}} \:.
	\end{equation}
The weights $\alpha$ and $\beta$ are chosen by the user.

Table~(\ref{ode:t:d}) shows a simple implementation of the described
algorithm.

	\begin{table}[thb]
\caption{An ODE driver with adaptive step-size control in C.
}\label{ode:t:d}
\lstinputlisting{c/ode/drive.c} 
	\end{table}

