The Logistic Map

Introduction

One of the most challenging topics in science is the study of chaos.  As an example of chaos, consider fluid flowing round an object. If the velocity of the fluid is not very large the fluid flows in a smooth steady way, called "laminar flow", which can be calculated for simple geometries. However, if the velocity is increased the flow becomes more complicated,  and depends on time and space in a non-trivial way which is essentially impossible to calculate. This "turbulent" behavior is called "chaos".  

Note that chaos does not mean that the behavior is completely random; one does see recognizable structures, such as little whirlpools, called "vortices", but it is very complicated. The reason that chaotic behavior is essentially impossible to calculate is that the future behavior of the fluid depends in a very sensitive way on what the fluid is doing now.  If we consider two configurations of the velocity field, which differ only by a very small amount, then, in the chaotic regime, those differences grow exponentially time, (the coefficient in the exponential is called the "Lyapunov exponent" which we will calculate in this notebook for a particular model) so that eventually  the two velocity fields become completely different.  This "sensitivity to initial conditions" is generally taken as the definition of chaos.  Because the flow of air in the atmosphere is chaotic, it is impossible to predict the weather in detail for more than about a week because we cannot integrate the equations with sufficient accuracy, even assuming that we had sufficiently detailed data on the state of the weather now (which we don't).

While a lot of chaos theory was developed by Poincaré around the start of the 20th century,  there were few subsequent developments (because the field is so difficult) until, starting in the 1980s, the mass availability of computers meant that the equations could  easily be integrated numerically, and perhaps most importantly, the solutions visualized by computer graphics.  

As the velocity of the fluid increases, rather than there being a discontinuous change from laminar to turbulent flow at a critical velocity, the motion changes only gradually as it passes the instability where it first becomes non-laminar, and it is only at larger values of the velocity that fully developed turbulence appears.  An important question, then, is the nature of this "transition to chaos".  He we will use the graphical ability of Mathematica to consider  a very simple looking model which has a transtion to chaos. There does not seem to be a single route to chaos, but the model that we study here, shows one of them, known as "period doubling". It is, in fact, the world's simplest model with a transition to chaos.

The Model

Most realistic systems with chaotic behavior, such as fluid flow around an obstacle, are described by a non-linear partial differential equation, which determines how the velocity changes in space and time. We will see later in the course that chaos also occurs in many ordinary (non-linear) differential equations. Here we will consider something even simpler, an "iterative map". This means that starting from an initial value of a variable, (Subscript[x, 0]) say, we generate a sequence of values, (Subscript[x, 1]), (Subscript[x, 2]), etc. from the map (i.e. function),  (Subscript[x, n+1])  = f((Subscript[x, n])) where we here make a simple choice f(x) = 4 λ x (1 - x), i.e.

	(Subscript[x, n+1]) = 4λ (Subscript[x, n])( 1-(Subscript[x, n])), 
	
where λ is parameter. In other words, (Subscript[x, 1]) = 4λ (Subscript[x, 0])(1-(Subscript[x, 0])),   (Subscript[x, 2]) = 4λ (Subscript[x, 1])(1-(Subscript[x, 1])), etc.  We will be interested in the behavior of successive iterations of this map, as a function of the parameter λ. In particular we will study the behavior of the (Subscript[x, n]) for large n.

We consider λ in the range from 0 to 1, so, if (Subscript[x, 0]) is  between 0 and 1, it is easy to see that all subsequent values of x also lie in this range. In fact the largest value of (Subscript[x, n+1])(which occurs for (Subscript[x, n])= 1/2) is equal to λ. 

This so-called "logistic map" has been used as  model for population dynamics, but here we just treat it as a toy model which has a transition to chaos.

A first look at the properties of the model

We create the function f in Mathematica as follows:

Clear["Global`*"]

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

It is a parabola which vanishes at x = 0 and 1. We plot it for λ = 0.4:

Plot[f[x]/.λ->0.4,{x,0,1}]

Let us first see what happens to successive iterations for a few values of λ. If we start with λ < 0.25, then one can easily show analytically that (Subscript[x, n+1]) < (Subscript[x, n]) and (lim_( n -> ∞)) x(n) = 0.  To verify this in Mathematica we put λ = 0.1 and generate a list of successive values with the Mathematica function NestList[func, x0, n], in  which func(x) is the function to be iterated, x0 is the starting value of x, and n is the number of iterations

λ=0.1; NestList[f, 0.5, 10]

Here we have set λ = 0.1, and then iterated f(x) 10 times with the initial value  0.5. The initial value of x is unimportant provided it is not a 0 or 1 in which case all the subsequent values of x would be zero. In this example we see that successive values eventually tend to zero. We say that zero is a fixed point, where the initial and final values of x are equal, i.e. X is a fixed point of the function f if

f(X) = X

It is best to see successive iterations graphically, i.e.

λ=0.1;

pt = ListPlot [ {NestList[f, 0.5, 10],NestList[f, 0.5, 10]} , Joined -> {True, False},Axes -> False, Frame -> True, FrameLabel -> {"n","(Subscript[x, n])"}, PlotRange -> {0, 0.12}, PlotMarkers -> Automatic, RotateLabel -> False,PlotLabel -> "λ = 0.1"]

The options PlotMarkers -> Automatic prevents the points being ridiculously small. 

We already discussed fixed point iteration in an earlier part of the class, and there we showed that a fixed point is stable, i.e. you converge towards it if you start off not too far away, provided

|f'(X)| <1

Since f'(x) = 4 λ (1 - 2x), the fixed point at x = 0 is stable when 4λ < 1, i.e. λ < 0.25. For λ > 0.25 this fixed point is unstable, and we need to find out what happens instead.

If we  start from a value of λ greater than 0.25 but less than 0.75, successive points "flow" to a fixed point at a non-zero value of x. However, for values of λ a little larger than 0.75 the fixed  point "bifurcates" to a "limit cycle" of period 2. This  then bifurcates again (i.e. the period doubles) at a larger value of λ to a limit cycle with period  4. As λ increases the period continues to double at successively closer and closer values of λ until, at around λ = 0.89, the period becomes infinity and we have chaotic behavior. This is illustrated in the following figure with four values of λ  (λ = 0.6, fixed point), (λ = 0.8, period 2 limit cycle), (λ = 0.88, period 4 limit cycle), and ( λ = 0.91, chaos):

λ=0.6;g1 =  ListPlot [ {NestList[f, 0.5, 10],NestList[f, 0.5, 10]} ,Joined -> {True, False}, Axes -> False, Frame -> True, FrameLabel -> {"n","(Subscript[x, n])"}, PlotRange -> {0.5, 0.62}, PlotMarkers -> Automatic, RotateLabel -> False,PlotLabel -> "λ = 0.6" ];

λ=0.8; g2 = ListPlot [ {NestList[f, 0.5, 10],NestList[f, 0.5, 10]},Joined -> {True, False}, Axes -> False, Frame -> True, FrameLabel -> {"n","(Subscript[x, n])"}, PlotRange -> {0.4, 0.9}, PlotMarkers -> Automatic, RotateLabel -> False,PlotLabel -> "λ = 0.8" ];

λ=0.88;g3 = ListPlot [ {NestList[f, 0.5, 10],NestList[f, 0.5, 10]},Joined -> {True, False}, Axes -> False, Frame -> True, FrameLabel -> {"n","(Subscript[x, n])"}, PlotRange -> {0.3, 0.9}, PlotMarkers -> Automatic, RotateLabel -> False,PlotLabel -> "λ = 0.88" ];

λ=0.91; g4 = ListPlot [ {NestList[f, 0.5, 10],NestList[f, 0.5, 10]},Joined -> {True, False}, Axes -> False, Frame -> True, FrameLabel -> {"n","(Subscript[x, n])"}, PlotRange -> {0, 1}, PlotMarkers -> Automatic, RotateLabel -> False,PlotLabel -> "λ = 0.91" ];

Show[GraphicsGrid[{{g1,g2},{g3,g4}},Spacings->Scaled[0.05]]]

The sequence towards which the x values converge is called an attractor. We have seen that an attractor may be a fixed point, a limit cycle, or a chaotic attractor. We will be interested in the nature of the attractor as a function of λ.

Fixed Points

In previous section, we saw that, for certain values of λ, the sequence of (Subscript[x, n])converges to a fixed point value X where X = f(X).We illustrate a fixed point graphically for the case of λ = 0.4

λ = 0.4;

Plot[{f[x],x},{x,0,1},AxesLabel->{"x","y"}, PlotRange -> {0, 0.5}]

The fixed point value is the point of intersection of the line y = x with the curve y = f(x). 

Mathematica has a convenient function called FixedPoint[f, start] for computing fixed points of maps. The first argument is the function and the second is the starting value for x. For example, still considering λ = 0.4, we get

FixedPoint[f, 0.5]

FixedPoint gives a result when two successive values are exactly the same. If convergence is very slow, or if roundoff errors prevent exact agreement, it is best to put in a less stringent test for determining whether convergence has been reached. This is done with the option SameTest to the FixedPoint command. The following example shows how it works.

FixedPoint[f, 0.5 , SameTest -> (Abs[#1 - #2] < 10^-11 &)]

The notation is that of a pure function, in which #1 and #2 are the arguments and the function is terminated by &.

We can find the value of λ where the fixed point at non-zero X becomes unstable quite easily using pencil and paper. We need to simultaneously solve

x= f(x),       |f'(x)| = 1

The solution is

λ = 3/4, x= 2/3

(There is also a second solution λ = 1/4, x = 0 which is the other limit of stability of the x ≠ 0 fixed point; as we have seen, for λ < 0.25 the stable fixed point is at x = 0.) Mathematica  will also give us the above result

Clear[x, λ]

Solve[{1 ==  4 λ (1 - x), 1 == Abs[ 4 λ (1 - 2 x)]}, {x, λ}]

Let's plot f(x) and x at the borderline case λ = 0.75.

λ = 0.75;

Plot[{f[x],x},{x,0,1},AxesLabel->{"x","y"}, PlotRange -> {0, 0.8}]

We see that f '(X) < 0 (X is the intersection point). Since the magnitude of the derivative is 1, we have f '(X) = -1.

We have seen that for certain choices of λ bigger that 0.75 successive values converge to a limit cycle of length 2. This is a fixed point of the twice iterated function f(f(x)). Lets call this f2(x). First we Clear λ and x:

Clear[ λ, x]

f2[x_] = f[f[x]]

Let's take λ = 0.8, where we found above that the system goes into a  limit cycle  of length 2

λ = 0.8;

and plot y=f2 (as well as y=x) against x

Plot[{f2[x],x},{x,0,1},AxesLabel->{"x","y"}]

We see that there are now 3 fixed points. The middle one is also a fixed point of f (note that any fixed point of f is also a fixed point of f2). See the plot below which also includes f (the parabola).

Plot[{f2[x], x, f[x]},{x,0,1},AxesLabel->{"x","y"}]

We can locate this fixed point using FindRoot (not FixedPoint because it is untable)

unstableFP = FindRoot[ f[x] == x, {x, 0.7}]

This fixed point is unstable because | f2'| (and also | f ' |) is greater than 1 at this value of x. Indeed, one can see from the figure above that f2[x] has a greater slope at the middle fixed point than the slope of x (i.e. its slope is greater than 1). We can also find the precise  numerical values of the derivative of f and f2 at this fixed point,

{f'[x] /. unstableFP, f2'[x] /. unstableFP}

showing that they are greater in magnitude than 1, confirming that the fixed point at 0.6875 is unstable.

The other two are fixed points of f2 (the two values of x correspond to the values of x in the period-2 cycle of f) and these are stable because |f2'| < 1 at these values of x, as is clear from the above figure and which we shall confirm quantitatively.

We can find the stable fixed points of f2 using the FixedPoint command

stableX =FixedPoint[f2, 0.5,SameTest  -> (Abs[#1 - #2]  < 10^-11&)]

We verify that | f2' < 1 | at this value of X

f2'[stableX]

The value 0.513045 is the lower of the two stable fixed points of f2 in the above figure. Depending on the different initial value for (Subscript[x, 0])we would converge either to this one or the other one, which is close to x=0.8:

stableX2 = FixedPoint[f2, 0.8,SameTest  -> (Abs[#1 - #2]  < 10^-11&)]

f2'[stableX2]

These last results are for λ = 0.8, where the fixed point of f2 (2-cycle of f) is stable. We now ask Mathematica to tell us where this attractor becomes unstable. We will help Mathematica by specifying that the derivative is -1 (rather than +1) at the instability point (We can see this, for example from the figure below, which is for λ = 0.86, near the limit of stability. The stable fixed point values of f2 are the first and third fixed points (not counting the trivial one at 0), as discussed above, which are near 0.4 and 0.9. The derivative is clearly negative there. The intermediate fixed point, at about 0.7 is unstable since one can see from the figure that its derivative is greater than 1. (Remember that the slope of the straight line is one.)

Clear[λ]

Plot[{f2[x]/.λ->0.86,x},{x,0,1}]

Clear[λ]

Solve[{f2[x] == x, f2'[x] == -1}, {x, λ}]

There are several solutions. To see which we want it is useful to compute the numerical values:

λ /. N[%]

The only solution with λ real and positive is λ = (1 + (Sqrt[6]))/4 = 0.86237. At this value of λ, the 2-cycle  becomes unstable and gives way to a 4-cycle.

We can  test for 4-cycles by  looking at the fourth iterated function f4(x) = f2(f2(x)):

Clear[ λ, x]

f4[x_] = f2[f2[x]]

We take λ = 0.88 where above we found a length-4 cycle above:

λ = 0.88;

This length-4 cycle is confirmed because we find a fixed point of f4:

FixedPoint[f4, 0.5, SameTest  -> (Abs[#1 - #2 ] < 10^-11&) ]

This value for x is one of the 4 values on the 4-cycle of f to which the trajectory converges. Starting with different values for (Subscript[x, 0]),we could converge to the other values of x.

Using Mathematica, we can determine where the stable fixed point of f4 (4-cycle of f) becomes unstable. Since  this is too complicated to solve analytically, and since NSolve gives a large number of roots, I prefer to use FindRoot, giving a reasonable initial guess based on the results from the iterations.

λ /. FindRoot[{f4[x] == x, f4'[x] == -1}, {x, 0.88}, {λ, 0.9}]

For λ > 0.886023 the 4-cycle  gives way to an 8-cycle. We can also find where the higher order cycles become unstable, though the calculations become increasingly complicated and require quite accurate initial guesses for the parameters λ and x. For, example, to determine where the 8-cycle becomes unstable I find

Clear[λ]

f8[x_] = f4[f4[x]] ;

λ /. FindRoot[{f8[x] == x, f8'[x] == -1}, {x,0.891}, {λ,0.891}]

To conclude this section, we have seen that if the results of successive iterations of the map suggest that there is convergence to a fixed point, or limit cycle of a certain length, it is useful to confirm this hypothesis with the FixedPoint command. We can also use the stability criterion |fn'(X)| < 1 (where fn(x) is the n-th iterate of f(x)) to determine where cycles of length n become unstable, using  Mathematica's commands, Solve, NSolve and FindRoot.

General  Behavior

We now investigate systematically the behavior as a function of λ.  It is convenient to define a function which will iterate the map n times. Furthermore, the ultimate behavior, such as fixed point or limit cycle, only appears after an initial number of iterations during which the x values are in the process of "converging" to the limiting behavior. The length of this "initial transient" increases as one approaches a bifurcation point where the period doubles. We therefore want our function to drop the first m, say, iterations and only print a list of the remaining n-m values. The former  is accomplished by the Mathematica commnd Drop so we define

iterate [ m_,  n_]  := Drop [ NestList [f, 0.5, n], m]

The initial value of x has been set to 0.5. We will use Mathematica commands to display the points in the (x, λ) plane. We display a point by the command Point, so we define the following function:

drawpt[y_] := Point[{λ, y}]

The result of a the command Point, and related commands, is a "graphics element". Graphics elements can be combined together with the Graphics  command to  produce a "graphics object" which can then be plotted with the command Show.

Here, we want to plot a point for each of the values of x generated by our function iterate. For this we use the Mathematica function Map[drawpt, list] which applies (maps) the function drawpt to each element of list. The following function, graph, sets up the lists of points that will be plotted:

graph[λmin_,λmax_, nλ_, mdrop_, n_] := Graphics [ {PointSize[0.001],Table [ Map[drawpt, iterate [mdrop,n ] ], {λ, λmin, λmax, (λmax-λmin)/nλ} ]}]

This takes nλ values of λ between λmin and λmax and, for each of these, makes a list of the first n iterates of f, starting with x=0.5, but dropping the first mdrop values.  These points can now be plotted using the Show command below, in which we have also given some useful plot options.

Show[graph[0.72,1,400,300,700],Axes->False,Frame->True,FrameLabel->{"λ","x"},PlotRange->{{0.72,1.02},{0,1}},AspectRatio->1]

This figure just shows the region for λ > 0.72, which is the most interesting part. For smaller values of λ one always has a fixed point, which is at 0 for λ < 0.25 and at a non-zero value for 0.25 <  λ <  0.75. This figure just shows the attactors, the sets of values of x towards which the iterations converge, for different values of λ. There can also be unstable fixed points and limit cycles, which we don't see. For example the fixed point which is stable for λ < 0.75, continues smoothly for λ > 0.75, but it becomes unstable, so the iterations don't flow towards it and hence it does not appear on the plot for λ > 0.75.  Instead, the stable attractor, a length-2 cycle, appears in the figure for λ just greater than 0.75.

Notice how rich the behavior is for λ > 0.75.  We not only see the period doubling cascade, discussed above, starting at 0.75 and leading to chaos at about 0.89, but in addition there are "windows" of limit cycle behavior even in the region from 0.89 to 1. The largest of these, which starts at λ = ((Sqrt[8])+1)/4 =  0.9571...  has a period 3 limit cycle, which then undergoes period doubling with periods 6,  12,  24 etc. until one enters again into a chaotic regime. The plot below, which enlarges this region, shows the period doubling cascade of the middle branch of the period 3 limit cycle,

Show[graph[0.96,0.963,400,300,700],Axes->False,Frame->True,FrameLabel->{"λ","x"},PlotRange->{0.455,0.525}, AspectRatio -> 1]

Notice that, even in this greatly enlarged view, one sees the same period doubling structure on smaller and smaller scales. One can continue to expand the scale of the plot to show even smaller scales, limited in practice only by numerical precision.  

Note too that there is structure in the density of points even in the chaotic regime.

The separations between successive values of λ where the period doubles get closer and closer together as the period increases. Let us define (Subscript[λ, k]) to be the value of λ at which the (2^(k-1))period becomes unstable resulting in a (2^k) cycle, so, for example, (Subscript[λ, 0]) = 1/4, (Subscript[λ, 1]) = 3/4,  (Subscript[λ, 2]) = ((Sqrt[6])+1)/4 ≃  0.8624... , (Subscript[λ, 3]) ≃  0.8860. Successive differences (Subscript[λ, k]) - (Subscript[λ, k-1]) are then in a geometric ratio (at least for large k), and the inverse of this, denoted by δ,  is called the Feigenbaum constant,  i.e.

δ = (Underscript[Lim, k->∞]) (( (Subscript[λ, k]) - (Subscript[λ, k-1]))/((Subscript[λ, k+1]) - (Subscript[λ, k]) ))

Its value is about 4.6692... .  This number is important because it is "universal", i.e. the same value is obtained for a wide range of functions (actually any function with a quadratic maximum), not just f(x) = 4 λ x (1 - x). Furthermore, for any given function, it is the same for all the period doublings which occur at different parameter values (here the parameter varied is λ).

To conclude this section, who would have thought that the simple looking map in the first equation would have the amazingly rich behavior shown in the last two figures. This richness was only fully appreciated once computers and computer graphics had developed to the point that figures like these could readily be produced.

Lyapunov Exponent

A defining feature of a chaotic system is sensitivity to initial conditions. If two trajectories which start off close to each other deviate more and more with increasing time, the system is said to be chaotic. The rate at which nearby trajectories deviate from each other with time is characterized by a quantity called the Lyapunov exponent. Here we discuss the Lyapunov exponent for the logistic map.

Consider two iterations of the logistic map starting from two values of x which are close together.  Let the two starting values be (Subscript[x, 0]) and (Subscript[x, 0])+(Subscript[δx, 0]). These map to (Subscript[x, 1]) and (Subscript[x, 1])+(Subscript[δx, 1]), ..., (Subscript[x, n])+(Subscript[δx, n]).  Expanding f(x) about (Subscript[x, n])we have

(Subscript[δx, n])= f'((Subscript[x, n-1])) (Subscript[δx, n-1])

assuming that (Subscript[δx, n]) is sufficiently small. Hence the separation of two trajectories after n steps, (Subscript[δx, n]), is related to their initial separation, (Subscript[δx, 0]), by

|((Subscript[δx, n])/(Subscript[δx, 0]))| = (∏_(i=0)^(n-1))| f'((Subscript[x, i])) | .

We expect that this will vary exponentially at large n like

|((Subscript[δx, n])/(Subscript[δx, 0]))| =  (e^((Subscript[λ, L])n))   (large n)

and so we define the Lypunov exponent (Subscript[λ, L]) by

(Subscript[λ, L])= (Subscript[Lim,  n -> ∞ ])(1/n)(∑_(i=1)^n) ln |f'((Subscript[x, i]))| .

If (Subscript[λ, L])> 0 neighboring trajectories diverge from each other at large n, which corresponds to chaos. However if the trajectories converge to a fixed point or limit cycle they will get closer together, which corresponds to (Subscript[λ, L])< 0.

Hence we can determine whether or not the system is chaotic by the sign of the Lyapunov exponent.
    
Below we calculate the Lyapunov exponent for some values of the parameter λ (not to be confused with the Lyapunov exponent (Subscript[λ, L])).

The Lyapunov exponent, which we call lya, is obtained by generating a list of (Subscript[x, i]) using NestList. The average of the elements of the a list can be very elegantly calculated in Mathematica using Apply[Plus, list] and dividing by Length[list]. Hence the average of ln | f'(x)| is given by

lya [l_, xinit_, n_, ndrop_] := (λ = l;xlist =Drop[  NestList[f,xinit,  n] ,  ndrop+1 ]; 
Apply[ Plus, Log[ Abs[ f'[xlist] ] ] ] / Length[xlist])

The following deternines (Subscript[λ, L]) for λ=0.91, starting with (Subscript[x, 0])=0.7, and averaging over 50000 iterations of the map, except for dropping off the first 20 to allow for an initial transient.

lya[0.91, 0.7, 50000, 20]

The positive value indicates that λ=0.91 is in a region of chaos. By contrast if we specify λ = 0.78

lya[0.78, 0.7, 50000, 20]

we get a negative value, indicating that the trajectory of points (Subscript[x, i]), (i = 0, 1, 2,  .... , ) converges to an attractor, which, in this case we already know  is a  length 2 limit cycle.

Next lets plot the Lyapynov exponent for a range of values of λ (this is slow)

Plot[lya[λ,0.7,5000,2000],{λ,0.5,1.},AxesOrigin->{0.5,0},AxesLabel->{"λ","(Subscript[λ, L])"}]

The main significance of this figure is that one can easily distinguish the regions which are chaotic ((Subscript[λ, L])> 0) from the regions which tend to a fixed point or limit cycle ((Subscript[λ, L])< 0). You see several points (the first is at λ = 0.75) where the Lyapuov exponent hits 0 and then goes negative again. These are the period doubling bifurcations. Precisely at the period doubling point the system is at the limit of chaos, but then becomes non-chaotic when the period doubles. However, at the end of the period doubling regime, at λ about 0.8922, (Subscript[λ, L]) crosses the axis and the system enters a chaotic regime. 

Now let's speed things up with a compiled version in which the function f[x] and its derivative are hard wired in. In the Compile command, the first argument is a list of variables which are assumed to be numerical, and the second argument is the function. It may be necessary to declare certain variables to be of a particular type such as integer. In this case the variable name is replaced by a list of two elements, the first being the variable and the second element the data type such as _Real or _Integer, e.g. {n, _Integer}. In the following, it is necessary that n and ndrop are declared to be integer:

lyaC = Compile[{ λ, xinit, {n, _Integer},{ndrop, _Integer}},  xlist =Drop[  NestList[4 λ # (1 - #)& ,xinit,  n] ,  ndrop+1 ]; 
Apply[ Plus, Log[ Abs[ 4  λ (1 - 2 xlist)]  ] ] / Length[xlist] ];

Note that the function f(x) = 4 λ x(1 - x) is represented as a pure function (without a name) whose argument is represented by # and which is terminated by &.

lya[0.9, 0.7, 50000, 100]  //Timing

lyaC[0.9, 0.7, 50000, 100]  //Timing

We see that the compiled version is significantly faster (the time is the first element of the list in the output.) 

Next we blow up the chaotic region.

Plot[lyaC[λ,0.7,30000,5000],{λ,0.85,1.},AxesOrigin->{0.85,0},AxesLabel->{"λ","(Subscript[λ, L])"}]

Note that for ,λ in the range greater than the point where (Subscript[λ, L]) first goes positive, there are many regions where (Subscript[λ, L]) is negative, "islands of stability" where the behavior is fixed point or limit cycle.

Now we home in on the principle period doubling transition to chaos by enlarging the scale.

Plot[lyaC[λ,0.7,30000,5000],{λ,0.892,0.8926},AxesOrigin->{0.892,0},AxesLabel->{"λ","(Subscript[λ, L])"}]

In the above figure, the first value of  ,λ where (Subscript[λ, L]) hits zero is (Subscript[λ, 5])=0.89219... . To the left of this zero the stable attractor is a period 16 cycle, and between that zero and the next one at (Subscript[λ, 6])= 0.89242... we have a period 32 cycle.

We now enlarge the scale again,  (this figure may take some time).

Plot[lyaC[λ,0.7,50000,15000],{λ,0.892482,0.892488},AxesLabel->{"λ","(Subscript[λ, L])"},PlotPoints->200]

Notice that there is fine detail no matter how much one expands the scale. We see that chaos emerges (i.e. (Subscript[λ, L])> 0) for λ between 0.892486 and 0.892487.

Numerical Precision in the Chaotic Regime

Consider the logistic map in the chaotic regime. We know that two trajectories which start off close together get further and further apart on iterating the map. The rate of separation of the trajectories is determined by the Lyapunov exponent (Subscript[λ, L]), i.e.

|(Subscript[δx, n])|=  (e^((Subscript[λ, L])n))|(Subscript[δx, 0])|

Hence any roundoff error in the calculation will also grow like  (e^((Subscript[λ, L])n)). If (Subscript[ϵ, m])is machine precision, then, when  (Subscript[ϵ, m]) (e^((Subscript[λ, L])n))is about 1, all precision will have been lost and the results will be just noise. Hence it is not possible to predict the value of (Subscript[x, n]) for very large n in the chaotic regime. This is analogous to being unable to predict the weather arbitrarily far into the future.

To illustrate the loss or precision, we will take λ = 1, for which the Lyapunov exponent is known to be

(Subscript[λ, L]) =log(2),

so

|(Subscript[δx, n])|= (2^n)|(Subscript[δx, 0])|  (large n).

Let's verify this numerically:

{lyaC[1, 0.7, 200000, 15000], N[Log[2]]}

The numerical value of the Lyapunov exponent is very close to log 2, as claimed.

We define the map and set λ = 1:

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

λ = 1;

We calculate the error in the result, using machine precision, after n iterations as a function of n.

diff[n_] := N[ Nest[f, 0.3, n] - Nest[f, N[3/10, 70], n] ]

In the line above, the first term on the right hand side is the result using machine precision, and the second is the result using 70 digits of precision (which is essentially exact by comparison). The starting value is x = 3/10. We now plot how the error increases with n.

err=ListPlot[Table[{n,Log[10,Abs[diff[n]]]},{n,2,60}],AxesLabel->{"n","(Subscript[log, 10])(err)"},PlotStyle->{PointSize[0.02],Hue[0]},PlotLabel->"(Machine Precision)"]

The vertical axis is (Subscript[log, 10])(err) which is the negative of the number of digits of accuracy (we will call the number of digits of accuracy the precision.) After about 55 iterations all precision has been lost. The lost of accuracy is just what one expects with an initial precision of a little better than 16 decimal places and an error growing with the Lyapunov exponent log(2). The following is a plot of y = -16.5 + (Subscript[log, 10])(2) x (note that (Subscript[log, 10])(2) = 0.30103):

lin= Plot[-16.5 + Log[10, 2] x, {x, 1, 60}, PlotStyle -> {AbsoluteThickness[2], Hue[0.7]}];

Show[err,lin]

It fits the data well:

Alternatively, we can use Mathematica's arbitrary precision numbers. In this case Mathematica only keeps the digits that it is sure are correct.  We can start with much higher precision, essentially whatever we want, but the precision that Mathematica claims in subsequent results goes down faster with increasing n than what we found above using machine precision.  The difference occurs because, with arbitrary precision numbers, if Mathematica is not sure of the correctness of  higher order digits it does not give them.

Below we take the initial precision to be 17, hardly more than the default precision of 16, but still (since > 16) an arbitrary precision number.

prec=ListPlot[Table[{n,-Precision[Nest[f,N[(3/10),17],n]]},{n,2,60}],AxesLabel->{"n","-Precision"},PlotStyle->{PointSize[0.02],Hue[0]},PlotLabel->"(Initial Precision=17)"]

We see that the (guaranteed) precision is lost about twice as fast as the (not-guaranteed) precision found above with machine precision numbers. Here is a fit, with slope 0.6:

lin2= Plot[-17 + 0.6 x, {x, 1, 32}, PlotStyle -> {AbsoluteThickness[2], Hue[0.7]}];

Show[prec,lin2]

Finally we show the precision claimed by Mathematica after 200 iterations, as a function of the initial precision.

ListPlot[Table[Precision[Nest[f,N[(3/10),n],200]],{n,5,150}],PlotStyle->{PointSize[0.02],Hue[0.7]},AxesLabel->{"precision","final precision"},PlotLabel->"(200 iterations)"]

We see that at least 120 digits of initial precision are necessary to determine (Subscript[x, 200])from (Subscript[x, n+1])= 4 (Subscript[x, n ])(1 - (Subscript[x, n])) (with (Subscript[x, 1])= 3/10).  

Overall, we have found empirically that the number of digits of precision which must be specified is about 0.6 of the number of iterations. Hence, to determine is (Subscript[x, 10000]), we need about 6000 digits.  The answer, according to my computer, is 0.104373:

N[ Nest[f, N[3/10, 6100], 10000] ] // Timing

What about (Subscript[x, 100000])?  I found 0.625851 (it takes quite a long time) :

N[ Nest[f, N[3/10, 61000], 100000] ] // Timing

What is (Subscript[x, 1000000])? The calculation is extremely slow. A student in a previous year's class, Dustin Gilbert, ran his computer for a long time to get the answer, and here it is, 0.751139. (Because it takes so long to execute, this cell is currently set to not execute when the notebook is evaluated. This option can be changed by selecting the cell by clicking on its vertical bar on the right, going to the Cell menu, and selecting Cell Properties.)

N[Nest[f,N[3/10,610000],1000000]]//Timing

{47924.1,0.751139}

The calculation took 47,924 secs, a bit over 13 hours!