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.
Back to Sprott's Technical Notes