NUMERICAL solution of ordinary differential equations
is the most important technique in continuous time dynamics. Since most
ordinary differential equations are
not soluble analytically, numerical integration
is the only way to obtain
information about the trajectory. Many different methods have been proposed and
used in an attempt to solve accurately various types of ordinary differential
equations. However there are a handful of methods known and used universally
(i.e., * Runge--Kutta*,
* Adams--Bashforth--Moulton*
and * Backward Differentiation Formulae* methods).
All these * discretize* the differential system to produce a
* difference equation* or * map*. The methods obtain different maps
from
the same differential equation, but they have the same aim; that the dynamics
of the map should correspond closely to the dynamics of the differential
equation. From the Runge--Kutta family of algorithms
come arguably the most well-known and used methods for numerical
integration (see, for example, Henrici henrici, Gear
gear, Lambert lam73, Stetter stetter,
Chua & Lin chua, Hall & Watt hall, Butcher
butcher, Press * et al.* press, Parker & Chua
parker, or Lambert lambert).
Thus we choose to look at Runge--Kutta methods to investigate what
pitfalls there may be in the integration of nonlinear and chaotic systems.

We examine here the initial-value problem; the conditions on the solution of the differential equation are all specified at the start of the trajectory---they are initial conditions. This is in contrast to other problems where conditions are specified both at the start and at the end of the trajectory, when we would have a (two-point) boundary-value problem.

Problems involving ordinary differential equations can always be reduced to a system of first-order ordinary differential equations by introducing new variables which are usually made to be derivatives of the original variables. Thus for generality, we consider the non-autonomous initial value problem

where represents . This can be either a single equation or, more generally, a coupled system of equations (often, we will have ):

The variable **x** often represents time. Almost all numerical methods
for the initial value problem are addressed to solving it in the above
form.

A non-autonomous system can always be transformed into an
autonomous system of dimension one higher by letting
,
so that we add the equation to the system and
to the initial conditions. In this case, however, we will have unbounded
solutions since as . This can be prevented
for non-autonomous systems that
are periodic in **x** by identifying planes of constant separated
by one period, so that the system is put onto a cylinder. We are usually
interested in one of the two cases above: either an autonomous system, or
a non-autonomous system that is periodic in **x**. In these cases, we can
define the concepts of the limit sets
of the system and their associated basins
of attraction which are so useful in dynamics.

It is known that sufficient conditions for a unique, continuous, differentiable
function to exist as a solution to this problem are that be
defined and continuous and satisfy a
* Lipschitz condition* in **y** in
. The Lipschitz condition is that

Here **L** is the Lipschitz constant
which must exist for the condition to be
satisfied. We shall always assume that such a unique solution exists.

Our aim is to investigate how well Runge--Kutta methods do at modelling
ordinary differential equations by looking at the resulting maps as
dynamical systems. Chaos in numerical analysis has been investigated before:
the midpoint method
in the papers by Yamaguti & Ushiki yamaguti and Ushiki
ushiki,
the Euler method by Gardini * et al.*
gardini,
the Euler method and the Heun method
by Peitgen & Richter peitgen, and the
Adams--Bashforth--Moulton methods in a paper by Prüfer prufer.
These studies dealt with the chaotic dynamics of the maps produced in their
own right, without relating them to the original differential equations.

In recent papers by Iserles iserles
and Yee * et al.* yee, the
connection is examined between a map and the differential equation that it
models. Other studies by Kloeden & Lorenz kloeden, and Beyn
beyn1,beyn2, concentrate on showing how the limit
sets
of the map are related to those of the ordinary differential equations.
Sauer & Yorke sauer use shadowing theory to find orbits of the
map which are shadowed by trajectories of the differential equation.

We bring together here all the strands in these different papers, and extend the examination of the connection between the map and the differential equation from our viewpoint as dynamicists. This topic has begun to catch the awareness of the scientific community lately (see for example Stewart stewart), and several of the papers we discuss appeared after the initial submission of this work; we have included comments on them in this revised version.

Wed Sep 27 17:21:22 MET 1995