Numerical Calculation of Largest Lyapunov Exponent
Department of Physics, University of
Madison, WI 53706, USA
October 15, 1997
(Revised August 31, 2004)
The usual test for chaos is calculation of the
Lyapunov exponent. A positive largest Lyapunov exponent indicates
chaos. When one has access to the equations generating the chaos,
this is relatively easy to do. When one only has access to an
data record, such a calculation is difficult to impossible, and that
will not be considered here. The general idea is to follow two
orbits and to calculate their average logarithmic rate of
Whenever they get too far apart, one of the orbits has to be moved back
to the vicinity of the other along the line of separation. A
procedure is to do this at each iteration. The complete procedure
is as follows:
Sample software that implements this procedure while searching for
solutions in general 2-D quadratic maps is available in (DOS) BASIC source
and executable code. Sample software
calculates the Lyapunov exponent (-0.5 to 0.5) for the Hénon Map
= 1 - CXn2 + BXn-1
for B = 0.3 and 0 < C < 2 is also available in
BASIC source and executable
code. Output from that program is shown below:
- Start with any initial condition in the basin of attraction.
Even better would be to start with a point known to be on the
in which case step 2 can be omitted.
- Iterate until the orbit is on the attractor.
This requires some judgement or prior knowledge of the system under
study. For most systems, it is safe just to iterate a few hundred
times and assume that is sufficient. It usually will be, and in
case, the error incurred by being slightly off the attractor is usually
not large unless you happen to be very close to a bifurcation point.
- Select (almost any) nearby point (separated by d0).
An appropriate choice of d0 is one that is about
1000 times larger than the precision of the floating point numbers that
are being used. For example, in (8-byte) double-precision
recommended for such calculations), variables have a 52-bit mantissa,
the precision is thus 2-52 = 2.22 x 10-16.
Therefore a value of d0 = 10-12 will
- Advance both orbits one iteration and calculate the new
The separation is calculated from the sum of the squares of the
in each variable. So for a 2-dimensional system with variables x
and y, the separation would be d = [(xa
- xb)2 + (ya - yb)2]1/2,
where the subscripts (a and b) denote the two orbits
- Evaluate log |d1/d0| in
By convention, the natural logarithm (base-e) is usually used,
but for maps, the Lyapunov exponent is often quoted in bits per
in which case you would need to use base-2. (Note that log2x
= 1.4427 loge x). You may get run-time errors
evaluating the logarithm if d1 becomes so small as
be indistinguishable from zero. In such a case, try using a
value of d0. If this doesn't suffice, you may
to ignore values where this happens, but in doing so, your calculation
of the Lyapunov exponent will be somewhat in error.
- Readjust one orbit so its separation is d0
same direction as d1.
This is probably the most difficult and error-prone step. As
an example (in 2-dimensions), suppose orbit b is the one to be
and its value after one iteration is (xb1, yb1).
It would then be reinitialized to xb0 = xa1
+ d0(xb1 - xa1)
and yb0 = ya1 + d0(yb1
- ya1) / d1.
- Repeat steps 4-6 many times and calculate the average of
You might wish to discard the first few values you obtain to be sure
the orbits have oriented themselves along the direction of maximum
It is also a good idea to calculate a running average as an indication
of whether the values have settled down to a unique number and to get
indication of the reliability of the calculation. Sometimes, the
result converges rather slowly, but a few thousand iterates usually
to obtain an estimate accurate to about two significant digits.
is a good idea to verify that your result is independent of initial
the value of d0, and the number of iterations
in the average. You may also want to test for unbounded orbits,
you will probably get numerical errors and the Lyapunov exponent will
be meaningful in such a case.
If the system consists of ordinary differential equations (a flow)
of difference equations (a map), the procedure is the same except that
the resulting exponent is divided by the iteration step size so that it
has units of inverse seconds instead of inverse iterations. An
for the Lorenz attractor is available. See
also the code for calculating the whole
spectrum of Lyapunov exponents.
Ref: J. C. Sprott, Chaos and Time-Series
(Oxford University Press, 2003), pp.116-117.
Back to Sprott's Technical Notes