\section{ORDINARY DIFFERENTIAL EQUATIONS}
Consider the initial value problem $$\eqalign{y'(x)&=f(x,y(x)),\quad
a\le x\le b,\cr y(a)&=A,\cr}\eqno(1)$$ where $f(x,y)$ is continuous for
$a\le x\le b$ and all $y$. There does not necessarily exist a unique
solution $y(x)$ throughout $[a,b]$. The problem $y'=\sqrt{|y|}$, $y(0)=0$,
$0\le x\le b$ has infinitely many solutions $\displaystyle y(x)=\cases{0,
&$0\le x\le c$,\cr {1\over4}(x-c)^2, &$c<x\le b$\cr}$. The problem
$y'=y^2, y(0)=1, 0\le x\le5$ has the unique solution $\displaystyle 
y(x)={1\over 1-x}$, which does not exist at $x=1$. Precise conditions
for the existence of a unique solution are given by

\medskip\noindent{\bf Theorem (Picard)}. Let $f(x,y)$ be continuous for
$a\le x\le b$ and all $y$, and satisfy the Lipschitz condition 
$\left|f(x,u)-f(x,v)\right|\le L|u-v|$ for $a\le x\le b$, all $u,v$.
Then there exists a unique solution $y(x)$ to the initial value problem
$$y'(x)=f(x,y(x)) \hbox{ for }a\le x\le b,\quad y(a)=A.$$

\noindent Picard's Theorem also holds for systems of equations like
$$\displaylines{{y'}_1(x)=f_1(x,y_1,y_2,\ldots, y_n),\qquad\qquad\qquad\qquad
y_1(a)=A_1,\cr \vdots\qquad\qquad \hbox{ with initial conditions } 
\qquad\qquad\vdots\cr
{y'}_n(x)=f_n(x,y_1,y_2,\ldots,y_n),\qquad\qquad\qquad\qquad y_n(a)=A_n,\cr}$$
where each $f_i$ satisfies a Lipschitz condition in each $y_j$. Note
that any higher order initial value problem
$y^{(n)}=g(x,y,y',\ldots$, $y^{(n-1)})$, $a\le x\le b$, $y(a)=A_1$,
$y'(a)=A_2,\ldots,y^{(n-1)}(a)=A_n$ can always be converted (in many
different ways) to a first order initial value problem.

\noindent{\bf Example 1}. $$\eqalign{y''(x)+xy'(x)+x^2y(x)=0\Leftrightarrow 
\quad&Y_1'=Y_2,\qquad
Y_2'=-x^2Y_1-xY_2 \hbox{ or }\cr &Y_1'=e^{-x^2/2}Y_2,\qquad
Y_2'=-x^2e^{x^2/2}Y_1 \hbox{ or }\cr &Y_1'={-x\over2}Y_1+Y_2,\qquad
Y_2'=\left({1\over2}-{3\over4}x^2\right)Y_1-{x\over2}Y_2.\cr}$$

\noindent{\bf Example 2}. $y''(x)+f(y(x))y'(x)+y(x)=0$. Let
$G'(t)=f(t)$, $Y_1(x)=y(x)$, $Y_2(x)=y'(x)+G(y(x))$. Then an equivalent
system is $Y_1'(x)=Y_2(x)-G(Y_1(x))$, $Y_2'(x)=-Y_1(x)$. Note that $G$
is smoother then $f$. 

\vfil\eject

\noindent{\it In the following it will be tacitly understood that $y$ and $f$ can
be vector functions.}

\centerline{\bf One-step methods.}

Euler's method is the simplest example of a one-step method. Assume that
Picard's Theorem applies and that the solution $y(x)$ to (1) has as many
derivatives as needed. Let $\displaystyle h={b-a\over N}$, $x_i=a+ih$, and $y_i$
denote an approximation to the exact solution $y(x_i)$. Then
$$y(x_{k+1})=y(x_k)+hf(x_k,y(x_k))+{h^2\over2!}y''(\xi)\approx
y(x_k)+hf(x_k,y(x_k))$$ suggests the algorithm
$$\eqalign{y_{k+1}&=y_k+hf(x_k,y_k)_,\qquad k=0,1,\ldots,N-1,\cr
y_0&=A,\qquad\qquad\qquad\qquad \hbox{(Euler's method)}.\cr}$$ Let $e_k=y(x_k)-y_k$,
$M_2=\max\limits_{[a,b]}|y''(x)|$, and $L$ be the Lipschitz constant for
$f$. Then an induction argument (Henrici, p. 271) proves that
$$\eqalign{|y(x_n)-y_n|&\le e^{L(x_n-a)}|y(x_0)-y_0|+{M_2h\over 2L}
(e^{L(x_n-a)}-1)\cr &\le e^{L(b-a)}|y(x_0)-y_0|+{M_2h\over 2L}
(e^{L(b-a)}-1)\cr}$$ for $n=1,\ldots,N$. Since normally $y_0=y(x_0)=A$,
the error in Euler's method $|y(x_n)-y_n|={\cal O}(h)\rightarrow 0$ as
$h\rightarrow 0$. 

The general form of a one-step method is
$$y_{n+1}=y_n+h\Phi(x_n,y_n),\quad y_0=A,$$ where $\Phi$, in addition to
$f$, is required to satisfy a Lipschitz condition:
$$|\Phi(x,u)-\Phi(x,v)|\le L|u-v| \hbox{ for }a\le x\le b, \hbox{ all }u,v, 
\hbox{ and }
0<h\le h_0.$$

\medskip\noindent{\bf Definition}. The one-step method $y_0=A$,
$y_{n+1}=y_n+h\Phi(x_n,y_n)$ is of order $p$ if
$y(x_{n+1})=y(x_n)+h\Phi(x_n,y(x_n))+h\tau_n$, where $\tau_n={\cal
O}(h^p)$. $\tau_n$ is called the {\it local truncation error}.

\medskip\noindent{\bf Theorem}. Let the one-step method of order $p$:
$y_0=A$, $y_{n+1}=y_n+h\Phi(x_n,y_n)$ be applied to the problem
$y'=f(x,y)$, $y(a)=A$, let $f(x,y)$ satisfy the hypotheses of Picard's
Theorem, and let $\Phi(x,y)$ satisfy a Lipschitz condition in $y$ with
Lipschitz constant $L$. Then for any $x_n=a+nh\in [a,b]$,
$$\left|y(x_n)-y_n\right|\le e^{L(x_n-a)}\left|y(x_0)-y_0\right|+{Ch^p\over
L}\left(e^{L(x_n-a)}-1\right).$$

\noindent Note that usually $y_0=y(x_0)=A$, so the error $|y(x_n)-y_n|={\cal
O}(h^p)$, whence the method is called {\it $p$th order}.
\medskip\hrule\medskip

\centerline{\bf Examples of one-step algorithms.}

\noindent{\bf Euler's method} (order $p=1$): $\Phi(x,y)=f(x,y)$.

\noindent{\bf Taylor algorithm} (order $p$):
$\displaystyle\Phi(x,y)=f(x,y)+{h\over2}f^{(1)}(x,y)+\cdots+{h^{p-1}\over
p!}f^{(p-1)}(x,y)$.

\noindent{\bf Heun's method} (order $p=2$):
$\displaystyle\Phi(x,y)_={1\over2}f(x,y)+{1\over2}f(x+h,y+hf(x,y))$.

\noindent{\bf Classical Runge-Kutta} (order $p=4$): $k_0=hf(x_n,y_n)$,
$\displaystyle k_1=hf(x_n+{h\over2},y_n+{1\over2}k_0)$,\hfil\break
$\displaystyle k_2=hf(x_n+{h\over2},y_n+{1\over2}k_1)$, $k_3=hf(x_n+h,y_n+k_2)$,
$\displaystyle y_{n+1}=y_n+{1\over6}[k_0+2k_1+2k_2+k3]$.

\noindent The derivation of Runge-Kutta formulas is very complicated, but the
concept is simple. Since
$$\eqalign{y(x_{n+1})&=y(x_n)+hf(x_n,y(x_n))+{h^2\over
2!}f^{(1)}(x_n,y(x_n))+\cdots+{h^p\over
p!}f^{(p-1)}(x_n,y(x_n))\cr &\qquad+{h^{p+1}\over (p+1)!}y^{(p+1)}(\xi_n),\cr
y(x_{n+1})&=y(x_n)+h\Phi(x_n,y(x_n))+h\tau_n,\quad \tau_n={\cal
O}(h^p),\cr}$$ the function $\Phi(x,y)$ must satisfy
$$\Phi(x,y)=f(x,y)+{h\over2}f^{(1)}(x,y)+\cdots+{h^{p-1}\over
p!}f^{(p-1)}(x,y)+{\cal O}(h^p).$$

\medskip\hrule\medskip

\centerline{\bf Error and stepsize control.}

Let $u(x)$ be the exact solution of $u'=f(x,u),u(x_n)=y_n$. The
{\it local error\/} from $x_n$ to $x_{n+1}$ is $u(x_{n+1})-y_{n+1}$,
i.e., the error in $y_{n+1}$ if $(x_n,y_n)$ had been exactly on the
curve $y(x)$. However, the point $(x_n,y_n)$ is probably not exactly
correct, so the {\it global error\/} $y(x_{n+1})-y_{n+1}$ may be very
different from the local error. Still, the global error can be
controlled by controlling the local error.

\medskip\noindent{\bf Theorem (Shampine and Allen)}. With the notation
and hypotheses of the previous Theorem, suppose that the local error is
kept down to $\epsilon$ per unit step at each step. Then
$\displaystyle |y(x_n)-y_n|\le {\epsilon\over L}\left(e^{
L(x_n-a)}-1\right)$ for any $x_n\in [a,b]$.
\medskip
The local error must be estimated before it can be
controlled.\hfil\break local error $=u(x_{n+1})-y_{n+1}
=u(x_n)+h\Phi(x_n,y_n) +h\tau_n -y_{n+1} =h\tau_n.$ Suppose that $\Phi$
and $\tilde\Phi$ correspond to methods of order $p$ and $q$
respectively, $p<q$. Let $$\eqalign{&y_{n+1}=y_n+h\Phi(x_n,y_n), \quad
\tilde y_{n+1}=y_n+h\tilde\Phi(x_n,y_n),\cr &u(x_{n+1})
-y_{n+1}=h\tau_n, \quad u(x_{n+1})-\tilde y_{n+1}=h\tilde\tau_n.\cr}$$
Then $\tilde y_{n+1}-y_{n+1}=h\tau_n-h\tilde \tau_n\approx h\tau_n$,
since $\tau_n={\cal O}(h^p)\gg{\cal O}(h^q)=\tilde\tau_n$. This
provides an estimate of the local error incurred by $\Phi$. The catch
is that a higher order method with $\tilde\Phi$ must be used
simultaneously with $\Phi$. If $\Phi$ and $\tilde\Phi$ could
``overlap'', i.e., use the same function values, this would be fairly
efficient. Such formulas are due to Fehlberg $(p=4,q=5)$ (see Shampine
and Allen):  $$\eqalign{&k_1=f(x_n,y_n),\quad
k_2=f(x_n+{h\over4},y_n+{h\over4}k_1),\cr &k_3=f\left(x_n+{3h\over8},
y_n+h\left({3\over32}k_1+{9\over32}k_2\right)\right),\cr
&k_4=f\left(x_n+{12h\over13},y_n+h\left({1932\over2197}k_1-{7200\over
2197}k_2+ {7296\over2197}k_3\right)\right),\cr &k_5=f\left(
x_n+h,y_n+h\left({439\over216}k_1-8k_2+
{3680\over513}k_3 -{845\over4104}k_4\right)\right),\cr
&k_6= f\left(x_n+ {h\over2},y_n+ h\left(-{8\over27}k_1+
2k_2-{3544\over2565}k_3 +
{1859\over4104}k_4- {11\over40}k_5\right)\right),\cr
&y_{n+1}= y_n+ h\left({25\over216}k_1+ {1408\over2565}k_3+
{2197\over4104}k_4- {1\over5} k_5\right),\cr &\tilde y_{n+1}=y_n+
h\left({16\over135}k_1+ {6656\over12825}k_3
 +{28561\over56430}k_4- {9\over50}k_5+ {2\over55}k_6\right).\cr}$$

\medskip\hrule\medskip
\centerline{\bf Algorithm for changing the step size.}

Let est$\displaystyle ={\tilde y_{n+1}-y_{n+1}\over h}
\approx{\hbox{local error}\over h}$ be an estimate of the local error
relative to $h$, the step size. If $|\hbox{est}|\le\epsilon$ is
satisfied, then accept $y_{n+1}$ and continue from $(x_{n+1},y_{n+1})$
with a new stepsize yet to be determined. If $|\hbox{est}|>\epsilon$,
then $h$ must be reduced to $\gamma h$ and a step of length $\gamma h$
taken from $(x_n,y_n)$. If a step of length $\gamma h$ were used, the
local error would be given by $$\hbox{local error }=(\gamma
h)^5\tau_0+{\cal O}(h^6)=\gamma^5h\hbox{ est }+{\cal O}(h^6)$$ for some
constant $\tau_0$. Now the step should be as large as possible, so
local error $\approx \gamma h\epsilon$. Therefore $\displaystyle \gamma
h\epsilon\approx \gamma^5h\hbox{ est }\Rightarrow
\gamma=\left({\epsilon\over\hbox{ est}}\right)^{1/4}$. This is the
``optimal'' $\gamma$.  Rejecting a step is costly, so in practice the
conservative step $.8\gamma h$ is used. In summary, starting from
$(x_n,y_n)$:  \item{1)} Take a step of length $h$, and compute the
estimate est of the local error relative to $h$.  \item{2)} Compute
$\displaystyle \gamma=\left({\epsilon\over est}\right)^{1/4}$, replace
$h$ by $.8\gamma h$.  \item{3)} If $|\hbox{est}|\le\epsilon$, accept
$y_{n+1}$ and continue from $(x_{n+1},y_{n+1})$.  \item{4)} If
$|\hbox{est}|>\epsilon$, reject $y_{n+1}$ and continue from
$(x_n,y_n)$.

Other practical considerations are that $h$ should never exceed some
maximum step size, nor be smaller than the floating point granularity at
$x_n$, nor be permitted to increase or decrease by more than some fixed
factors. The maximum step size limit prevents stepping over sudden
changes in the solution, and the increase/decrease factors mollify
``chattering.''