# Numerical Errors in Logistic Map Calculation

J. C. Sprott
Department of Physics, University of Wisconsin, Madison, WI 53706, USA
September 15, 1997
(Modified April 20, 2000)

When numerically iterating the logistic equation

Xn+1 = AXn(1 - Xn)

round-off errors can cause spurious results.  Consider the case A = 4 which maps the unit interval (0, 1) back onto itself exactly twice.  Because of the shape of the parabolic logistic function, many values of Xn near 0.5 map into values of Xn+1 near 1.0 which in turn map into values of Xn+2 close to zero.  If these values of Xn+2 are indistinguishable from zero to within machine precision, they will thereafter remain at zero because zero is an unstable fixed point solution of the logistic map.  As a result, initial conditions X0 will iterate to zero on average after a relatively small number of iterations N whose value we estimate below.

Suppose that the logistic equation with A = 4 is iterated using (4-byte) single-precision arithmetic.  Four bytes is 32 bits, of which 23 are used for the mantissa, 8 for the exponent, and 1 for the sign according to the IEEE standard for floating point arithmetic.  Thus the computer can represent 223 = 8,388,608 unique values in the range (0, 1) for a given exponent.  Values near zero can be represented more accurately by changing the exponent.  Thus the problem occurs for values of Xn+1 that are indistinguishable from 1.  Values larger than 1 - 2-24 will round to 1.0.  With this value of Xn+1, we can calculate the range of Xn values whose next iterates round to 1.0 by solving the quadratic equation

X2 - X + 1/4 - 2-26 = 0

which has solutions X = 0.5 ± 2-13.  Thus the region of the interval (0, 1) over which  Xn is an eventually fixed point is of width 2-12 = 2.44 x 10-4 centered on 0.5.  Since the probability density of finding an iterate of the logistic map with A = 4 is

P = 1 / pi[X(1 - X)]1/2

which has a value of P = 2 / pi = 0.637 for X = 0.5, the probability that an iterate will fall within the region above is 1.55 x 10-4.  Thus on average, we would expect the reciprocal of this number of iterates before the fixed point at zero is reached, or N = 6400.

To test this prediction, the logistic map was iterated in single precision many times using a uniform random initial condition in (0,1), and the number of iterations required for zero to be reached was determined.  The source and object code used for this test are available.  Typical results follow:

Solution reached 0 after 3334 iterations - Average = 1981 - Cases = 989
Solution reached 0 after 1082 iterations - Average = 1980 - Cases = 990
Solution reached 0 after 3153 iterations - Average = 1981 - Cases = 991
Solution reached 0 after 4017 iterations - Average = 1984 - Cases = 992
Solution reached 0 after 2941 iterations - Average = 1985 - Cases = 993
Solution reached 0 after 2581 iterations - Average = 1985 - Cases = 994
Solution reached 0 after 1183 iterations - Average = 1984 - Cases = 995
Solution reached 0 after 3890 iterations - Average = 1986 - Cases = 996
Solution reached 0 after 2604 iterations - Average = 1987 - Cases = 997
Solution caught in a periodic cycle
Solution reached 0 after 3136 iterations - Average = 1988 - Cases = 998
Solution reached 0 after 542 iterations - Average = 1987 - Cases = 999
Solution reached 0 after 2082 iterations - Average = 1987 - Cases = 1000
Solution reached 0 after 927 iterations - Average = 1986 - Cases = 1001
Solution reached 0 after 1984 iterations - Average = 1986 - Cases = 1002
Solution reached 0 after 538 iterations - Average = 1984 - Cases = 1003
Solution caught in a periodic cycle
Solution reached 0 after 297 iterations - Average = 1982 - Cases = 1004
Solution reached 0 after 3455 iterations - Average = 1984 - Cases = 1005
Solution reached 0 after 2529 iterations - Average = 1984 - Cases = 1006
Solution reached 0 after 2051 iterations - Average = 1985 - Cases = 1007
Solution reached 0 after 1509 iterations - Average = 1984 - Cases = 1008
Solution reached 0 after 1056 iterations - Average = 1983 - Cases = 1009
Solution reached 0 after 2723 iterations - Average = 1984 - Cases = 1010
Occasionally the orbit gets stuck in some very high periodic cycle and never reaches zero.  Excluding such cases, the average over 32,000 runs is N = 2031.  This result is about 3 ( or pi?) times worse than expected for reasons that are not clear.  Furthermore, the expectation is that the number N is inversely proportional to the square root of the precision.  Thus for (8-byte) double-precision with a mantissa of 52 bits, the expectation is for N = 1.5 x 108, which is in reasonable agreement with numerical results.  This particular problem can be solved by using a value of A slightly less than 4.0, such as 3.999999, although the orbit still occasionally gets caught in a cycle of some presumably high period.