The Logistic Equation:
Computable Chaos

Steven Jaffe
Mathematics Department
Vassar College
Poughkeepsie NY 12601
E-mail: stjaffe@vassar.edu

The following is designed to examine the strange dynamics of the logistic equation when r=4:  x(t+1)=f(x(t)); f(x)=4x(1-x).
see James V. Whittaker, An Analytical Description of Some Simple Cases of Chaotic Behaviour, American Mathematical Monthly 98: 489--504 (June--July 1991).

The result is that if we make the change of variables x=Sin^2(Pi/2*v),
v between 0 and 1, then the dynamics can be described very simply in terms of the binary digits of v:
	
	Let the initial value of v be
			v0=0.a1a2a3...  in binary: each ai is either 0 or 1.
	  		(for v=1 take all the ai=1).

Then the next iterate, v1, will be given by:
	v1= 0.a2a3...   if a1=0
	v1= 0.b2b3...   if a1=1, where bi is the opposite digit to ai: 
			1 if ai is 0, 0 if ai is 1.
			
For example, if v0=.01001010... then v1=0.1001010...
If v0=.1010010001... then v1=.101101110...

This description leads to some immediate, but nevertheless remarkable conclusions.  For instance, if v has a finite binary expansion, i.e. only a finite number n of nonzero binary digits, then after at most n+1 steps the solution will reach v=0, and of course stay there, i.e. it will be absorbed by the fixed point at v=0.  (Note that v=0 corresponds to x=sin^2(0)=0.)  Thus there is a dense set of points whose trajectories eventually hit, and stick at, 0.

Define the logistic function:

f[x_]:=4x(1-x)

Define the change-of-variables taking v to x:

x[v_]:=Sin[(Pi/2)*v]^2

Start at v=1/2; this corresponds to x=

x[1/2]

NestList[f,x[1/2],10]

(iterates the difference equation 10 times)

Try v=1/4 which is .01 in binary: 
(to avoid very lengthy exact calculations, we evaluate the starting point numerically; but note that there is a  danger that this will affect the long-term behavior because the solution depends very sensitively on the precise starting value.)

NestList[f,N[x[1/4]],10]

Now with v=1/8 which is .001 in binary; the trajectory should  follow, in binary v coordinates, the sequence
.00100...,  .0100...,  .100...,  .111... =1,  .000..., and then remain at 0  
		
		(note the way .111... came from .1000...:  
		drop the first digit, leaving all 0's; 
		but since the dropped digit was a 1, we invert all the rest, 
		giving us all 1's; this number is the same as 1.0.)

NestList[f,N[x[1/8]],10]

Note: here we see the problem we warned of earlier: because we have approximated the starting value numerically, the solution does not reach 0 exactly, but reaches a very small number. 

This seems understandable, and we are probably willing to treat a number like 2.1684x10^-19 as “as good as” 0.

However, if we continue this approximate trajectory, it may eventually get large, whereas the true trajectory should stick at exactly 0 forever.  

Let's see:

NestList[f,N[x[1/8]],100]

ListPlot[%,PlotJoined->True]

Thus is demonstrated sensitivity to initial conditions!

Mathematica normally does arithmetic to 20 decimal places, but we can get it to use more by specifying the initial condition to more places as follows:

NestList[f,N[x[1/8],40],100]

By working with 40 decimal places, we have achieved the true trajectory.

N[Last[Out[3]],100]

Another starting value, this time given directly in base 2
notation:
(note that there are 26 digits, so this should reach 0 after 26 or 27 steps)

v0=2^^.10010001000010000010000001
NestList[f,N[x[v0]],35]

Count them!  (Note that the first number in the list is the starting value x0.)

Now let's make a periodic orbit; say of period 3.
If v0=.100100100100... repeating, in binary, then the trajectory will go

.100100100100,  .11011011011..., .0100100100..., .100100100...
 
 in other, since the digits of v0 repeat every 3, the trajectory will have period 3:  
 
 Now in order to use this v0 we need a way to express it numerically in a finite form: 
 note that 8*v0= 100.100100100... (in binary!), so
 		8*v0 = 4 + v0
 which gives v0 = 4/7.

v0=4/7
n=NestList[f,N[x[v0]],35]

The graph shows the period:

p1=ListPlot[n,
	PlotStyle->{PointSize[.02],RGBColor[1,0,0]},
	DisplayFunction->Identity];
p2=ListPlot[n,PlotJoined->True,DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction]

One last example: an orbit of period 7:
let's take v = .10100101010010...
Then 2^7*v = v+(1010010)(base 2),
 or 128v = v+82, 
 or v = 82/127

v0=82/127
n=NestList[f,N[x[v0]],35]

The graph shows the period more clearly:

p1=ListPlot[n,
	PlotStyle->{PointSize[.02],RGBColor[1,0,0]},
	DisplayFunction->Identity];
p2=ListPlot[n,PlotJoined->True,DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction]

Try for period 5: 
Note that in general to get period p we need to start with 
v0 = a / (2^p-1), where a is between 1 and 2^p-1.

x0=N[x[6/31],20]

n=NestList[f,x0,20]

p1=ListPlot[n,
	PlotStyle->{PointSize[.02],RGBColor[1,0,0]},
	DisplayFunction->Identity];
p2=ListPlot[n,PlotJoined->True,DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction]

Period 17: note that we need to work with many more decimal places in order to keep an accurate solution over such a long period:

x0=N[x[1/(2^17-1)],100];

n=NestList[f,x0,70];

p1=ListPlot[n,
	PlotStyle->{PointSize[.02],RGBColor[1,0,0]},
	DisplayFunction->Identity];
p2=ListPlot[n,PlotJoined->True,DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction,
		PlotRange->All]

n[[1]]

n[[2]]

n[[18]]

n[[19]]

A fine point: when we make a period p orbit this way, the starting point x0 may or may not itself be on the orbit: sometimes the orbit starts instead with x1.  The reason is that when we follow through the iterations on the digits of v0, there will be reversals when we hit  a leading digit of 1; 
if, after p steps, there have been an even number of reversals, we are back to v0;
but, if after p steps there are an odd number of reversals we come back, not to v0, but to v0 with the digits all reversed: call this v0'. Considering carefully how the reversals come about, we can show that v0' with the the digits reversed will now repeat with period p.  Furthermore the next iterate that follows v0' will be exactly v1, the iterate that follows v0, i.e. F(v0')=F(v0), so the period orbit will be 
v0', v1, v2, v3, ..., vp=v0', etc.

For example, earlier we got a period 5 orbit by using v0=6/31:

v0=6/31; x0=N[x[v0]]; n=NestList[f,x0,12]

If we started with v0=7/31 instead, we get

v0=7/31; x0=N[x[v0]]; n=NestList[f,x0,12]

Note that x5 is not equal to x0, but x6=x1, x7=x2, etc.
The binary iterates in the v coordinates are
v0 = .001110011100111...
v1 = .01110011100111...
v2 = .1110011100111...
v3 = .001100011000...  (reversal!)
v4 = .01100011000...
v5 = .1100011000...  (this is v0 reversed!)
v6 = .011100111...   (same as v1 - so now it all repeats)

For further play: define "try" to do everything once v0 and the number of iterates nits are chosen. We can also specify the number of digits in the numerical calculations by ndigits:

ndigits=19

try:=
(
x0=N[x[v0],ndigits];
n=NestList[f,x0,nits];
p1=ListPlot[n,
	PlotStyle->{PointSize[.02],RGBColor[1,0,0]},
	DisplayFunction->Identity];
p2=ListPlot[n,PlotJoined->True,DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction,
		PlotRange->All]
)

Try it!

v0=1/8;nits=10;try

v0=2^^.01001000100001;nits=20;try

v0=1/(2^35-1);nits=100;try

Need more digits!

v0=1/(2^35-1);nits=100;ndigits=100;try

That's better! An orbit of period 35.

If we start with a value for x0 chosen “at random”, we see:

n=NestList[f,.3275,100];

p1=ListPlot[n,
	PlotStyle->{PointSize[.02],RGBColor[1,0,0]},
	DisplayFunction->Identity];
p2=ListPlot[n,PlotJoined->True,DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction]

a “chaotic” trajectory; this is “typical” behavior for this equation.

Another periodic trajectory: p=5 (so 2^p-1=31)

v0=16/31; nits=35; ndigits=20; try

A trajectory that hits 0 after 10 or 11 steps:

v0=2^^.1010010011; nits=35; ndigits=20; try

This was the starting value x0:

N[x[v0]]

To make your own solutions:

First method: use the “try” command: first set the number of digits in the computation, number of iterations, and starting value v0.
Note that v0 can be given directly in binary form using the notation 2^^.a0a1a2...:

ndigits=20; nits=35; v0=2^^.10010001; try

Second method: copy and edit the following sequence of cells; this allows you more control of the parameters; for example you can specify x0 directly.

x0=N[x[1/(2^17-1)],100];

n=NestList[f,x0,70];

p1=ListPlot[n,
	PlotStyle->{PointSize[.02],RGBColor[1,0,0]},
	DisplayFunction->Identity];
p2=ListPlot[n,PlotJoined->True,DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction,
		PlotRange->All]