The logistic map revisited
Jerzy Ombach, Cracow, Poland ombach@im.uj.edu.pl October 8, 1999
This worksheet explores the period-doubling bifurcation sequence and ther phenomena associated with the discrete logistic map f(x) =a*x*(1-x). I refer to Bob Corless' worksheet logmap.mws in the share library, where you can see the power of algebraic utilities in Maple. Here I will show more geometric and numerical approach than Bob did, which is now possible due to new features of Maple, release 5.
> restart:
> libname := "C:/Program Files/Maple Release 5.1/Ombach",libname:
> with(onedd):
Introduction
Define the logistic function
> eval(f);
Plot the graph for a particular value of the parameter:
> a := 3.8:
> plot(f,0..1);
Animate the above plot.
> plots[animate]('a'*x*(1 - x), x = 0..1,'a' = 0..4);
Use the context menu to see the animation. Note that for parameters a : 0 <= a <= 4 the unit interval [0,1] is maped onto itself by the map .
Next we will iterate points from the interval [0,1] by the map . Thus, for a given in [0,1] we will look for the points defined inductively as , for n = 0, 1, 2, .... Identifying points x from the interval with points (x,x) on the diagonal we may interpret successive 's geometrically.
We find 10 iterations, i.e. the initial 11 elements of the orbit starting at = 0.2324 under map . We also request for the value of .
> orbit(20,3.5,0.02324,xn);
> xn;
The above animation may suggest that our orbit tends to some 4-periodic orbit. Make sure by continuing from . In this case we do not require for the value of the last point, and thus the fourth argument is omitted.
> orbit(10,3.5,xn);
Define the fourth iteration g(x) = f(f(f(f(x)))) and plot its graph along with the graph of f and the diagonal.
> g := f@@4;
> a := 3.5: plot([f(x),g(x),x],x=0..1, axes = BOXED, color =[BLACK,RED,BLUE]);
The plot shows two fixed points of (or equilibrium states), they are the points of intersection of the graph of with the diagonal. They are also fixed points of = . Still, has six additional fixed points. Previous examination of iterations indicates that four of them form an attracting 4-periodic orbit. A basic theorem of dynamical system says that attractivity of fixed points and periodic orbits depends on the absolute value of the derivative of the corresponding iteration. In particular, if p is an n-periodic point of and the derivative of at p has absolute value < 1, then the periodic orbit of p is attracting, and is not if this value is > 1. In our case we then have by inspection of the plot one 4-peridic attracting orbit. Alternatively, we can numerically find fixed points of and and compute the appropriate derivatives.
> [solve(f(x) = x)];
> map(unapply(diff(f(x),x),x),%);
> [solve(g(x) = x)];
> select(x ->(Im(x)= 0),%);
> map(unapply(diff(g(x),x),x),%);
Let us look for the limit set of a given point..
> x[0] := 0.6343:
> for n from 1 to 200 do x[n] := evalf(f(x[n-1])) od:
> plots[pointplot]([seq([n,x[n]],n = 101..200)], view = [100..200,0..1]);
Another way to make sure that for a = 3.5 we have a 4-periodic attractor is to examine graphs of successive noting the change of their shape. We will see the animation of 20 successive iterations of f when the parameter a = 3.5.
> iteration(20,3.5);
I suggest you to repeat the above steps with different parameters. In particular:
use 0 < a <= 3 to identify an attracting fixed point. use 3 < a <= 3.569 to identify -periodic attracting periodic orbits. use a = 3.9 or a = 4 to identify chaotic behavior.
Try other a <= 4.
Bifurcation of periodic orbits and transition to chaos
Going further we will see how the limit set of an orbit changes while the parameter a increases. First we plot the limit set (or more formally the segment , ... , of the orbit starting at x0 = 0.23) versus the parameter a = 4..
> limset(4,0.23,200,100,DIAMOND);
In the following procedure we make bifurcation diagram or animation (if whatever as the third argument is present) for the values of parameter a ranging over interval [0, 4].
> bif(0,4,232);
> bif(0,4);
Decrease the interval for parameter a.
> bif(2.5,4);
We want to see the reasons for which fixed point zero is an attractor for 0 <= a < 1, then a non-zero fixed point appears which is an attractor when a passes by 1, then a 2- periodic attracting periodic orbit is born when a passes by 3, which in turn bifurcates into 4-perioddic attractor and so on.
How about bifurcation at ? Remember the meaning of the derivative for the attractivity of the fixed point. We animate the graph of f for a in the interval [0.8, 2.4] with the step = 0.2.
> aiter(1,0.8,2.4,0.2);
How about bifurcation at ? Remember that fixed points of that are not fixed point of f form a 2-periodic orbit.
> aiter(2,2.8,3.2,0.02);
How long is this 2-periodic orbit an attractor?
> aiter(2,3,3.5,0.01);
It seems that for some a below 3.5 the 2-periodic orbit is no longer an attractor. We find an appropriate value of a .
> a := 'a': x := 'x':
> r1 := expand(eval((f@@2)(x)));
> r2 := diff(r1,x);
> fsolve({r1 = x,r2 = -1},{a,x},{a = 3.4..3.6,x = 0.3..0.5});
This bifurcation diagram suggests that a 4-periodic attractor appears at and losses its attractivity at some a above 3.5. Then, animate the fourth iteration of f .
> aiter(4,3.4,3.6,0.01);
> r1 := (f@@4)(x):
> r2 := diff(r1,x):
> fsolve({r1 = x,r2 = -1},{a,x},{a = 3.5..3.6,x = 0.3..0.5});
Thus the 4-peridodic orbit loses its attractivity at the parameter . Does 8-peridic orbit then appear? The concept of renormalization developed by Feigenbaum may help with understanding the whole process.
> afeig(3,3.6,0.05);
We see that the family of essential part of - deciding in fact of the dynamics of - marked with the small rectangle - left, and expanded by the linear change of variables to the feig(f) defined on the interval [0,1] - right, behaves for a ranging over certain interval [3,b] like the family of f does while a ranges over the interval [1,4]. First, we find this b and then check by inspection of the animation that our method of finding b is correct.
> a := 'a': x :='x':
> fsolve({f(x) = x,(f@@3)(0.5) = x},{a,x},{a = 3.6..3.7,x = 0.5..1});
> afeig(3,3.678573510,0.02);
In particular, one can prove that any k-periodic (resp. attracting) orbit of feig(f) produces a 2k-periodic orbit of f (resp. attracting). Moreover, bifuration of non-zero fixed point into 2-periodic orbit for feig(f) implies bifurcation of 2-periodic into 4-periodic orbit for f, and so on .
Now, for a > we can define feig ( feig(f)), etc. Thus, attracting -periodic orbits bifurcate into attracting -periodic orbits more and more quickly. There exists (approx. =3.570) such that f has periodic orbits with all -periods, k = 1,2,3, ..., none of them attracting. For a > periodic orbits with another periods appear. In particular, for a > b = periodic orbits with odd periods appear.
We look for a 3-peridic orbit.
> aiter(3,3.5,4,0.02);
We can see that 3-periodic orbit appears for a .> 3.8, when the graph of crosses the diagonal. In fact, they are two such orbits, one attracting for a while and the other repelling. We find the values of a for which 3-periodic orbits are born and one of them becomes an attractor.
> a := 'a':
> r1 := (f@@3)(x):
> fsolve({r1 = x,r2 = 1}, {a,x},{a = 3.8..3.9, x = 0.4..0.6});
When 3-periodic orbit is no longer an attractor (and bifurcates into 6-periodic attractor)?
> fsolve({r1 = x,r2 = -1},{a,x},{a = 3.8..3.9, x = 0.4..0.6});
By the Sharkovskii Theorem the existence of period three implies the existence of all other periods! So, for example, 5-periodic orbit should have appeared for some a < . Check it now.
> aiter(5,3.6,4,0.02);
In fact, 5-peridic orbits (an attractor including) appear for a close to 3.74. Incidentally, we can see two other values of a , when extra 5-periodic orbits appear. I reccommend examining the rising of periodic orbits with period 7, 9, 11, ... . It would be a good test of your computer power!
And what about a 6-periodic attractor? From this what we have already have learned, there are at least two values of a for which a 6-periodic attractor is born. One is a = mentioned above, when 3-periodic orbit bifurcates. We should also anticipate another 6-periodic orbit, which corresponds to a 3-periodic orbit of the feig(f).
> bif(3.6,3.7);
> r1 := (f@@6)(x):
Below we use Newton procedure from the share library instead of fsolve routine to find a value of the parameter a and a corresponding 6-periodic point.
> with(share): with(linalg): with(Newton):
See ?share and ?share,contents for information about the share library
Warning, new definition for diag
Warning, new definition for norm
Warning, new definition for trace
Share Library: Newton
Author: Taylor, Ross.
Description: An implementation of a Newton iteration for solvingnon-linear systems of equations
> Newton([r1 = x,r2 = 1],[a =3.62,x =0.3],output={functions,variables},iterations = 10,tolerance = 1e-6);
We check the result graphically.
> assign(%);
> orbit(6,a,x,xn);
> xn, x - xn;
> x := 'x':
Chaos and statistical behavior
They are parameter values, such that there is no periodic attractor. One can prove, that is just the case. Let us look for iterations of some point.
> x0 := stats[random, uniform](1);
> orbit(100,4,x0);
However, there are many orbits whith very simple behavior. In particular, there are infinite many periodic orbits and they are dense in the interval [0,1]. Look for 5-periodic orbits for example.
> a := 4:
> plots[display]([onedd[diag],plot(f@@5,0..1)]);
Another simple orbits are the eventually periodic orbits. For example the orbit starting at .
> orbit(10,4,0.5,xn);
> xn ;
Going farther, we can consider orbits starting at preimages/onedd od under the iterations of .
> pre := array(0..4);
> pre[0] := [1/2]: for k from 1 to 4 do pre[k] := [seq(solve(f(x) = pre[k-1][i]), i = 1..nops(pre[k-1]))]; od;
Any orbit starting from one of the above points collapse at zero. For example:
> orbit(10,4,1/2+1/4*sqrt(2-sqrt(2+sqrt(2+sqrt(2)))));
Have a look on the all computed preaimages/onedd:
> preimages/onedd := 0,1,seq(op(pre[k]), k = 0..4):
> nops({preimages/onedd});
> plots[pointplot]([seq([x,x],x = preimages/onedd)]);
Points which orbits are eventually zero are dense in [0,1].
Similarly, points which orbits are eventually the nonzero fixed point are dense in [0,1]. The same is true for eventually periodic points with any given period.
> solve(f(x) = x,x);
> pre[0] := [3/4]: for k from 1 to 4 do pre[k] := [seq(solve(f(x) = pre[k-1][i]), i = 1..nops(pre[k-1]))]; od;
> preimages/onedd := 3/4,seq(op(pre[k]), k = 0..4):
On the other hand, there are also points not eventually periodic. For example, points , k, n = 1, 2, 3, ..., with < form a dense set in in [0,1] but are not eventually periodic.
We can thus agree that for a = 4 the system is chaotic. However, we will see that even in this case the system behaves nicely from a statistical point of view.
We will count how many times the orbit starting at x0 and having 5 000 points visits one of each of the 100 disjoint intervals of lengths .
> emp_distr(4,x0,100,5000);
One can change the starting point x0, still the distribution will be more less the same. It appears that there is a density distribution function corresponding to the empirical distribution with the following density function.
> phi := x -> 1/(Pi*sqrt(x*(1-x)));
> plot(phi, 0..1);
> int(phi, 0..1), 5000/100;
> plots[display]({emp_distr(4,x0,100,5000), plot(50*phi(x),x = 0..1,numpoints = 400)});
A word of caution
For some values of parameter a , especially for these responsible for chaos, the system is sensitive to initial conditions and the roundoff error. Perform 25 iterations for two very close points.
> orbit(25,4,0.2324,xn):
> orbit(25,4,0.2324001,xn):
Increase the accuracy and repeat:
> Digits := 50:
> Digits := 10:
We see that the roundoff error is perhaps of less importance, and the differences are caused by the dynamical system itself. We thus may ask, if the 'orbits' produced by the computer give information of any value about the 'true orbits'. This problem has not been completely solved so far. Still for the case a = 4, it is known that the system has the shadowing property, i.e. any 'orbit' produced by computer is very close shadowed (traced) by some 'true one' .
