1. Bifurcations. In Lab 2 we have seen that the behaviour of a dynamical system can change quite dramatically with the change of system parameters.  It is important to understand that these changes are not only quantitative, such as, for example, change in the location of a fixed point, but also qualitative: fixed points can be created or destroyed, their stability can change, the system behaviour can change from regular (stationary or periodic) to irregular - chaotic.  It is the qualitative changes in the system dynamics that are the subject of investigation in the theory of dynamical systems.

Definition: Qualitative changes in the system dynamics are called bifurcations, and the parameter values at which they occur are called bifurcation points.

2. Constructing bifurcation diagrams. In order to study bifurcations in dynamical systems, it is convenient to visualise the bifurcations that happen at different parameter values.  You might have realised from Lab 2 that, even though it is possible to find bifurcation points using the plots of orbits or the cobweb diagram, it is not the most convenient and efficient way of identifying the bifurcation points.  The reason is that with such plots you can only visualise the system behaviour at one parameter value at a time.  A better way to see the general behaviour of the system at different parameter values is to plot the orbits as a function of the parameter.  That is, we will plot the orbit points x_n along the vertical axis against the values of parameter r along the horizontal axis.   Such a plot is called the bifurcation diagram.

Note that we will not plot all the points of an orbit, just the ones that represent the eventual state of the system (asymptotic behaviour) at each parameter value.   Therefore, the first hundred or so orbit points will be discarded allowing the orbit to settle down into its asymptotic behaviour.

Below is the program that constructs a bifurcation diagram for the logistic map with parameter r in the range from 2.5 to 4.  As before, save the program text into a Matlab script file and run it by simply typing the name of the file at the Matlab prompt.

Program 1

Npre = 200; Nplot = 100;
x = zeros(Nplot,1);
for r = 2.5:0.005:4.0,
  x(1) = 0.5;
  for n = 1:Npre,
    x(1) = r*x(1)*(1 - x(1));
  end,
  for n = 1:Nplot-1,
    x(n+1) = r*x(n)*(1 - x(n));
  end,
  plot(r*ones(Nplot,1), x, '.', 'markersize', 2);
  hold on;
end,
title('Bifurcation diagram of the logistic map');
xlabel('r');  ylabel('x_n');
set(gca, 'xlim', [2.5 4.0]);
hold off;

Comments to the program:

• In the 1st line we specify: Npre, the number of orbit points we want to discard (number of pre-iterates),  and Nplot, the number of orbit points we want to plot for each value of parameter r (number of iterates).
• Line 2 defines x as a Matlab array with Nplot rows and 1 column.
• In the 3rd line we start a loop over the values of parameter r, starting from 2.5 and ending at 4.0 with the increment of 0.005.
• Line 4 specifies the starting point for the orbit.  The structure of the bifurcation diagram should not depend of a particular choice of the starting point, as long as the value is in the range from 0 to 1.  Initial points outside this range result in orbits that diverge to  -  .
• The next three lines are the loop where the pre-iterates are computed.  Note that, since we do not want to store the values of the orbit points in the pre-iteration loop,  we use x(1) on the left-hand-side, so that the old value is discarded and replaced by the new value of the orbit point.
• Next is the loop that computes the points that will be plotted.  This loop is identical to the one in Program 1 of Lab 2.
• Next line plots the orbit stored in x for a given value of r.  Recall that in the Matlab function plot(x,y), the sizes of arrays x and y must be identical.  Since r is a 1-by-1 array (a scalar) while x is an  Nplot-by-1 array (a vector), we use expression r*ones(Nplot,1) to create a vector whose size is identical to that of x and whose elements are all equal to r.   Alternatively, we can use the Matlab function repmat.  That is, repmat(r,Nplot,1) produces the same result by replicating the scalar r  into Nplot rows and 1 column, thus creating the Nplot-by-1 array.
• The plot function also specifies the type of symbols to be used (dots), and contains a pair of arguments ('markersize',2) that change the size of the plotted dots.  The default markersize is 6, which is too large and causes many dots to overlap.
• We use hold on to retain all the plotted points while plotting orbits for subsequent values of parameter r.
• Matlab function set can be used to modify the appearance of the figure and its components (axes, lines, dots, etc.).  Here it is used to ensure that the limits of the horizontal axis coinside with the range of r values.  Similarly, the statement  set(gca, 'ylim', [0.5 1.0]); can be used to 'zoom in' on the upper half of the bifurcation diagram.  Try it.

 
Q1  Explore the influence of parameter Npre on the appearance of the bifurcation diagram.   Describe the apparent differences in the diagrams produced with different values of Npre (say 10, 50, 200) and explain these differences.

Program 1 can be easily modified to create a bifurcation diagram for any range of parameter values with any increment and any number of points plotted for each r.  However, specifying the increment directly is not very convenient, since it is difficult to estimate how many different values of r are necessary to obtain good-looking diagram.  Large increment values result in a very sparse diagram while very small increment values increase computation time and may overload computer memory (if this ever happens, use Control-C to interrupt program execution).   It is more convenient to specify the number of different r values to be used in the given range and then tell the program to calculate the necessary increment.  To this end,

Q2  Transform Program 1 into a Matlab function bifdiag(Npre, Nplot, rfrom, rto, Nr), where rfrom and rto give the range of r, and Nr is the number of different r values in this range.
Hint: bifdiag(200, 100, 2.5, 4.0, 301) should plot a diagram identical to that produced by Program 1.

The diagram produced by Program 1 combines information about the asymptotic behaviour of the logistic map for different parameter values in the range from 2.5 to 4.0.  Remember that each vertical line contains Nplot points.  So, if we see only a small number of points appearing on a given vertical line, it means that these points are repeatedly visited by the orbit and, therefore, the orbit is periodic.  It is thus clear that for 2.5 < r < 3.0 the orbits converge to a fixed point, while for 3.0 < r < 3.45 (approx.) they converge to a period-2 orbit.   At approximately r = 3.45 there is another period-doubling bifurcation, where a period-2 orbit is replaced by a period-4 orbit.  In fact, a closer look reveals that further increasing r results in a whole cascade of period-doubling bifurcations occuring closer and closer to each other and producing orbits of periods 8, 16, 32, 64, and so on.

Q3 (The Feigenbaum constant)  Let r_n denote a bifurcation point where a period-2ⁿ orbit is created.  For example, r₁ = 3.0 and r₂  3.45.  We see from the bifurcation diagram that the intervals between successive period-doubling bifurcations become smaller and smaller.  To quantify this decrease, we introduce ratios of the neighbouring intervals as follows: _n = (r_n+1 - r_n)/(r_n+2 - r_n+1) .

It was shown in 1975 by Feigenbaum that, in the limit of large n, the ratios  _n converge to a single number  called the Feigenbaum constant.  The remarkable thing about this constant that it has the same value for all maps where a period-doubling bifurcation cascade can be observed (and it is indeed quite common).  Because of this, it is often referred to as a universal constant and enjoys the same level of respect in the dynamical systems theory as  = 3.1415... in geometry or e = 2.718... in calculus.  The subject of this exercise is to estimate the value of the Feigenbaum constant. To this end, modify Program 1, or use function bifdiag to construct a bifurcation diagram showing orbits of period 4, 8, 16 and 32 with enough detail to allow you to estimate the bifurcation points r₃ , r₄ and r₅.  (To get better estimates of these values, it may be helpful to learn how to use the 'zoom in' feature available in Matlab figures and put a grid on the plot by typing grid on at the Matlab prompt.) Use these values to calculate  ₁,  ₂ and  ₃ .

3. Self-similarity.  Closer inspection of different regions of the bifurcation diagram reveals remarkable structure and patterns.   The period-doubling cascade we have just discussed ends at about r = 3.57.  Then we see many values of r where the orbits do not appear to be periodic.  Indeed, these orbits are not periodic and, no matter how many orbit points we compute, we never get back to the same point.  Such orbits are called chaotic.  Note, however, that in the range 3.57 < r < 4.0, there are also regions where the orbits are again periodic.  Such regions are called periodic windows. The largest such window begins near r  3.83 and contains a stable period-3 orbit.   A closer look at the period-3 window (as well as any other periodic window) reveals that the structure of the window repeats the structure of the overall bifurcation diagram!  That is, the stable orbit suffers the same fate as the original fixed point of the logistic map, i.e., it undergoes a cascade of period-doubling bifurcations, creating orbits of period 6, 12, 24, and so on.  What is most remarkable, is that each branch of the diagram looks exactly like the whole diagram.

Q4 To see this, construct a bifurcation diagram in the range 3.83 < r < 3.856 and zoom in on the upper branch (for example, by typing set(gca, 'ylim', [.953 .964]); at the Matlab prompt).

Lab Report No. 3 is due on Thursday, May 8th at 11:30am.

BACK