The usual test for chaos is calculation of the largest 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 experimental data record, such a calculation is difficult to impossible, and that case will not be considered here.  The general idea is to follow two nearby orbits and to calculate their average logarithmic rate of separation.  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 conservative procedure is to do this at each iteration.  The complete procedure is as follows:

1. Start with any initial condition in the basin of attraction.

Even better would be to start with a point known to be on the attractor, in which case step 2 can be omitted.


2. 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 any 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.


3. Select (almost any) nearby point (separated by d₀).

An appropriate choice of d₀ is one that is on the order of the square root of the precision of the floating point numbers that are being used.  For example, in (8-byte) double-precision (minimum recommended for such calculations), variables have a 52-bit mantissa, and the precision is thus 2^-52 = 2.22 x 10^-16.  Therefore a value of d₀ = 10^-8 will usually suffice.


4. Advance both orbits one iteration and calculate the new separation d₁.

The separation is calculated from the sum of the squares of the differences in each variable.  So for a 2-dimensional system with variables x and y, the separation would be d = [(x_a - x_b)² + (y_a - y_b)²]^1/2, where the subscripts (a and b) denote the two orbits respectively.


5. Evaluate log |d₁/d₀| in any convenient base.

By convention, the natural logarithm (base-e) is usually used, but for maps, the Lyapunov exponent is often quoted in bits per iteration, in which case you would need to use base-2.  (Note that log₂x = 1.4427 log_e x).  You may get run-time errors when evaluating the logarithm if d₁ becomes so small as to be indistinguishable from zero.  In such a case, try using a larger value of d₀.  If this doesn't suffice, you may have to ignore values where this happens, but in doing so, your calculation of the Lyapunov exponent will be somewhat in error.


6. Readjust one orbit so its separation is d₀ and is in the same direction as d₁.

This is probably the most difficult and error-prone step.  As an example (in 2-dimensions), suppose orbit b is the one to be adjusted and its value after one iteration is (x_b1, y_b1).  It would then be reinitialized to x_b0 = x_a1 + d₀(x_b1 - x_a1) / d₁ and y_b0 = y_a1 + d₀(y_b1 - y_a1) / d₁.



7. Repeat steps 4-6 many times and calculate the average of step 5.

You might wish to discard the first few values you obtain to be sure the orbits have oriented themselves along the direction of maximum expansion.  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 an indication of the reliability of the calculation.  Sometimes, the result converges rather slowly, but a few thousand iterates of a map usually suffices to obtain an estimate accurate to about two significant digits.  It is a good idea to verify that your result is independent of initial conditions, the value of d₀, and the number of iterations included in the average.  You may also want to test for unbounded orbits, since you will probably get numerical errors and the Lyapunov exponent will not be meaningful in such a case.

If the system consists of ordinary differential equations (a flow) instead 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.  You will typically need millions of iterations of the differential equations to get a result good to better than about two significant digits. An example for the Lorenz attractor is available. See also the code for calculating the whole spectrum of Lyapunov exponents.

Sometimes you can get the whole spectrum of exponents using the above method, for example when the system is a two dimensional chaotic map or a three dimensional chaotic flow if you know the rate of state space contraction averaged along the orbit (the dissipation) which is the sum of the Lyapunov exponents (easily calculated from the trace of the Jacobian matrix averaged along the orbit for a flow or from the average determinant of the Jacobian matrix for a map) and using the fact that one exponent must be zero for a continuous flow.

To estimate the uncertainty in your calculated Lyapunov exponent, you can repeat the calculation for many different initial conditions (within the basin of attraction) and perturbation directions. For a chaotic system, the initial condition need only be changed slightly since orbits quickly become uncorrelated due to the sensitive dependence on initial conditions. You can then calculate a mean and standard deviation of the calculated values so as to avoid the all too common mistake of quoting more digits than are significant.

Ref: J. C. Sprott, Chaos and Time-Series Analysis (Oxford University Press, 2003), pp.116-117.

Back to Sprott's Technical Notes