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.
