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*.

For example, consider the so-called: Qualitative changes in the system dynamics are calledDefinitionbifurcations, and the parameter values at which they occur are calledbifurcation 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 parameterNpreon the appearance of the bifurcation diagram. Describe the apparent differences in the diagrams produced with different values ofNpre(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 functionbifdiag(Npre, Nplot, rfrom, rto, Nr), whererfromandrtogive the range ofr,andNris the number of differentrvalues 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) Letr_{n}denote a bifurcation point where a period-2^{n}orbit is created. For example,r_{1}= 3.0 and r_{2 }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 theFeigenbaum 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 auniversalconstant and enjoys the same level of respect in the dynamical systems theory as = 3.1415... in geometry ore= 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 functionbifdiagto 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_{3 }, r_{4}and r_{5}. (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 typinggrid onat the Matlab prompt.) Use these values to calculate_{1},_{2 }and_{3 }.

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 <The result is an almost exact copy of the whole diagram. Actually, the diagram contains infinitely many such copies of itself. Thisr< 3.856 and zoom in on the upper branch (for example, by typingset(gca, 'ylim', [.953 .964]);at the Matlab prompt).

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