Off[General::spell1];

beginAssignment;

Iteration of real functions, 3:
from order to chaos

In this Notebook we consider a 1-parameter family of simple and very similar quadratic functions quad[a].  
We will look at what happens to a number when you repeatedly apply quad[a] to it.  
When a is 2 the behaviour turns out to be very simple, and when a is 4 the behaviour can only be described as chaos.  What happens in between?

In this exercise you will discover something that you have probably never encountered before.
The aim is not to just evaluate each cell, but to appreciate the result of the calculations. 
Answer all the questions asked in a way that is clear and unambiguous. 
You answers should display your understanding. 
Use appropriate words as sentences as needed to convey this — provide more rather than less.

Here is the definition of our family of quadratic functions:

quad[a_][x_]:= a*x(1 - x)

Here are the graphs of the functions quad[2], quad[3] and quad[4] (and of the identity function).

Plot[{x,quad[2][x],quad[3][x],quad[4][x]},{x,-0.2`,1.2`},AspectRatio->Automatic]

You can see that each of the three quadratic functions maps the interval [0,1] into itself and has two fixed points: the x-coordinates of the points where the graph of the function cuts the graph of the identity function.

The main Mathematica functions used in this Notebook are Nest and NestList.

Clear[f,c]

Nest[f,c,3]

NestList[f,c,3]

Order

We are going to look at quad[2].

quad[2][x]

Here is what happens when you repeatedly apply quad[2] to the number 0:

NestList[quad[2],0,20]

Here is what happens when you repeatedly apply quad[2] to a number very close to 0:

NestList[quad[2],0.001,20]

And here is what happens when you repeatedly apply quad[2] to a number close to 1:

NestList[quad[2],.9,20]

There are two fixed points: 0 is repelling 0.5 attracts

Solve[quad[2][x]==x,x]

Try a few other starting values from the interval [0,1] to get a good feeling for what happens.

ListPlot[NestList[quad[2],0.9`,20],PlotRange->{0,1},PlotStyle->Hue[0.`,1,1]]

Here are some cobwebs that display the fact that the number 0.5 is an attracting fixed point of quad2 in a different way.

cobweb[quad[2],0.05,10];

cobweb[quad[2],0.95,10];

Starting at the large dot the cobweb goes repeatedly from a point on the graph of quad[2] to a point on the graph of the identity function.  Think about it and you will see that the x-coordinates of the vertical segments of this cobweb are the successive iterates of the starting number 0.95 under the function quad[2] (and so are the y-coordinates of the horizontal segments).

Chaos

Let's generate 101 0f the numbers you get by repeatedly applying the function quad[4], starting with 0.8.

NestList[quad[4],0.8,100]

Can you see any pattern here?  Perhaps a picture will help.

ListPlot[NestList[quad[4],Nest[quad[4],0.8`,100],20]
,PlotRange->{0,1},PlotStyle->Hue[0.`,1,1]]

We can make this picture of the sequence of iterates clearer by joining the successive points, as below:

plotA=ListPlot[NestList[quad[4],0.8`,40]
,PlotRange->{0,1},Joined->True,PlotStyle->Hue[0.`,1,1]]

There is no discernible regularity.  Worse still, consider the orbit of a starting number that is extremely close to 0.8 but not equal to it.

plotB=ListPlot[NestList[quad[4],0.80000001`,40]
,PlotRange->{0,1},Joined->True,PlotStyle->Hue[0.6`,1,1]]

Compare the two graphs:

Show[plotA,plotB]

As you would expect from the continuity of our function quad[4], initially the orbits of the two points 0.8 and 0.80000001 are so close together that one curve hides the other; what is surprising is that after less than 30 iterations the orbits differ wildly from one another.  
This is chaos.

Imagine a physical system that was modelled by iteration of the function quad[4].  
Typically we might be interested in predicting its state some time in the future on the basis of the current state.  
Our example shows that a minute and perhaps unavoidable experimental error (of 0.00000001) in the initial state would render such predictions completely useless after as little as 30 iterations.  
It is the same phenomenon that makes detailed weather prediction over more than a few days impossible.

We are looking at a rough mathematical justification of the insight that the loss of a nail in the shoe of a horse may affect the course of history, though not of course in a predictable way.

Let's consider our example from the perspective of cobwebs.

cobweb[quad[4],0.8,10];

After 10 iterations there is no evident regularity.  
Here's a smart way of looking at 10 iterations of the orbit after 200 have already taken place:

cobweb[quad[4],Nest[quad[4],0.8,200],10];

Again there is no evident regularity.

The emergence of chaos from order

What happens to points under iteration of quad[2] is very regular.  What happens to points under iteration of quad[4] is very chaotic.  As a moves from 2 to 4 the behaviour of quad[a] changes between these extremes.  The way that the transition takes place is full of fascinating detail which you can explore using Mathematica.

Here we shall make a beginning by looking at quad[3.1]. 
You can find the fixed points of quad[3.1] by solving the quadratic equation quad[3.1][x]== x. 
They turn out to be 0 and 21/31.

You should be able to anticipate what the next picture will look like.

cobweb[quad[3.1],21/31,10];

Here is what happens if we start near but not at the fixed point 21/31.

cobweb[quad[3.1],2/3,10,0.5,0.8];

It appears that both fixed points 21/31 and 0 are repelling.  
So the orbit of the number 2/3, which is trapped in the interval [0,1], does not converge to either.  
What does happen?

cobweb[quad[3.1],2/3,40,0.5,0.8];

If you wait 40 iterations and then draw the next 10 iterations of the orbit we see something new:

cobweb[quad[3.1],Nest[quad[3.1],2/3,40],10,0.5,0.8];

If you replace the 40 by a larger number you will get the same picture.  
Can you see what has happened?

Suppose we try the same thing but with another starting point:

cobweb[quad[3.1],Nest[quad[3.1],.1,40],10,0.5,0.8];

We have discovered an attracting 2-cycle.

Let’s look directly at the sequence of iterates.  You can plainly see convergence to the 2-cycle.

list2=NestList[quad[3.1],.8,40]

ListPlot[NestList[quad[3.1`],0.8`,20],PlotRange->{0,1},PlotStyle->Hue[0.`,1,1]]

ListPlot[NestList[quad[3.1`],0.8`,20]
,PlotRange->{0,1},Joined->True,PlotStyle->Hue[0.`,1,1]]

Epilogue

The function we needed for drawing cobwebs is defined under here.  There's no need to look at the definition unless you really want to.

The heart of the following is a function F that transforms lists in a certain way:
F[{a,b,c,d,e}]->{{a,b},{b,b},{b,c},{c,c},{c,d},{d,d},{d,e},{e,e}}

cobweb::usage = "
cobweb[function,start,numberOfSteps,xMin,xMax] \n
generates a picture of numberOfSteps iterations \n
of function, starting at the number start; \n
the graph goes from xMin to xMax. \n
\n
cobweb[function,start,numberOfSteps] \n
generates a picture of numberOfSteps iterations \n
of function, starting at start; the graph goes \n
from 0 to 1.";

cobweb[function_,start_,numberOfSteps_,xMin_:0.0,xMax_:1.0]:=
Block[{list,quadruple},
list=NestList[function,start,numberOfSteps];
quadruple[x_]:=Transpose[{x,x,x,x}];
points=
Partition[Flatten[Prepend[quadruple[Rest[list]],First[list]]],2];
Plot[{x,Evaluate[function[x]]},{x,xMin,xMax},
PlotRange->{xMin,xMax},
AspectRatio->Automatic,
Epilog->{Hue[.0,1,1],Line[points],PointSize[.03],
         Point[{start,function[start]}]}] ]

?cobweb