Lab 1R: Maps and Orbits

Section 1: Introduction

1.1 Flows.  Broadly speaking, a dynamical system is any system whose evolution in time is described by deterministic laws.  An example could be Newton's laws of motion.  Given the velocities and positions of particles at some instant of time, these laws enable the prediction of the velocities and positions at all future times.  In mathematical terms, such laws can be best represented by a set of ordinary differential equations:

dx/dt = F(x,t) ,
where x is a vector whose components are the positions and velocities of the particles. Since time t is a continuous quantity (a real number), such dynamical systems are called continuous, or flows. The solutions of the differential equations can be though of as lines in the space of  the vector x  (called the phase space) and are called orbits. Orbits starting from different points in phase space indeed resemble the flow of a fluid (the orbits cannot intersect each other), hence the name.

1.2. Maps.  Another type of dynamical systems can be imagined, where the evolution law is given by a function f(x) that maps the system from one state, say x0, to another state, say x1, as follows:  x1= f(x0).  Next, x1 is mapped into x2 according to the same rule: x2= f(x1), and so on, so that, in general

xn+1 = f(xn) .
In this case, the states of the system are numbered by an integer n which can be though of as a "time" variable.  Such dynamical system is called discrete, or simply a map.  An orbit of the discrete dynamical system is an infinite sequence of points {x0,  x1,  x2,  x3, ...}.

1.3. The Logistic map.  Just like the continuous dynamical systems model real processes occuring in nature, a map can also model a real process.  Some of the simplest, yet interesting, models of this type can be found in biology.  Imagine insects which, during one year of their life, hatch out of eggs, grow, mature, lay eggs and die.  Next year a new generation of insects hatches out of the eggs layed by the previous generation.  In the simplest possible model, we can assume that each mature insect lays, on average, k eggs, so that if, in year 1, there were z1 insects, in year 2 there will be z2 = k z1, in year 3 there will be z3 = k z2 insects, and so on, assuming that the number k does not change from year to year (is a constant).  In general, we have a map:

zn+1 = k zn .
Note, however, that this model is not very realistic since, if k < 1 the insects will eventually die out: zn -> 0 as n -> $\infty$ , while for k > 1, the number of insects will be multiplied by k every year, so that zn -> $\infty$   as n -> $\infty$ .  An obvious reason preventing the number of insects from becoming infinitely large is the availability of food.  When the food supply is limited, some of the insects will die out before reaching maturity.  In the simplest possible model of this kind, we can assume that k is not a constant, but decreases (linearly) with increasing number of insects:
kn = r (1 - zn/zmax) ,
where r is another constant and zmax is the population number at which insects eat all available food and die before reaching maturity. For this more realistic model, the map is given by zn+1 = r zn (1 - zn/zmax) , or, introducing the relative number of insects xn = zn/zmax ,

xn+1 = r xn (1 - xn) .

This simple model of the insect population dynamics is called the logistic map.  Despite its apparent simplicity, this map is capable of producing extremely complicated orbits.
 
 

Section 2: Plotting orbits

During this lab session, we shall be examining the orbits generated by the logistic map, with a particular focus on how the properties of the orbits depend on the parameter r.  To get started, we will write a simple Matlab program to study the behaviour of orbits generated by the logistic map.

Below is a program which computes an orbit of length N for a logistic map with a given parameter r starting from a specified initial condition.  Copy this program into a file and then run it in the Matlab command window (by typing the name of the file on the command line).  You should be able to copy and paste the lines from the Explorer window into the Matlab editor without re-typing them.

Program 1

r = 0.8;  N = 50;
x = zeros(N,1);
x(1) = 0.9;
for n = 1:N-1,
  x(n+1) = r*x(n)*(1 - x(n));
end,
plot(1:N, x, 'o-');
title('Orbit of the logistic map');
xlabel('n'); ylabel('x_n');

Comments to the program:

 
Q2.1.  Use the program above to explore the behaviour of the orbits for different values of r and starting from different initial points in the range from 0 to 1.  (Remember to save the file each time you change it, before running the program).  Describe the orbit behavior for r = 0.8,  2, 3, 3.2, 3.5, 3.83, 4, 4.5.  Some of the possible behaviours are: orbits from different initial points all tend to 0; ... all tend to some number other than zero;   ... oscillate between two values;  ... oscillate among ? values; no pattern - different initial conditions result in different orbits; orbits escape from the interval [0, 1].  You should also classify each orbit according to its eventual (asymptotic) state: fixed point, period-? orbit, chaotic orbit.


As a consequence of this experiment, we see that even a very simple dynamical system - as simple as a quadratic function - can exhibit fairly complicated behaviour.  Moreover, the behaviour changes as the parameter r changes.

Next, we shall try to identify values of the parameter r for which the behaviour of orbits changes from one type to another.  Before doing that, note that it is not very convenient to have to save the program each time you make changes to the parameters.  It would be much more convenient if the parameters could be entered directly into the program from the command line.  To this end,
 

Q2.2.a.  Transform the program into a Matlab function, so that the same results can be achieved by typing, for example,

>> orbit(0.8,50,0.9);

where 'orbit' is the name of the file containing the function.  The result should be the same as that of the program above.
 

The study of dynamical systems concerns with understanding the eventual state of the system or, in other words, its asymptotic behaviour.  Therefore, it is often desirable to look only at the last N points of the orbit.  To this end,
 
Q2.2.b.  Further modify the program so that it computes first M points of the orbit without displaying them, and then displays the next N orbit points.  Note that there will be M + N orbit points computed altogether.  The value of M can be either specified in the function body (say,  M = 100;) or entered as an additional function argument, namely

>> orbit(0.8,100,50,0.9);

Q2.3.   Use the function from Q2.2. or the program above to find (approximately) the values of the parameter r for which all orbits
 a) stop tending to 0 and begin tending to a different fixed point;
 b) stop tending to a fixed point and start tending to a period-2 orbit;
 c) stop tending to a period-2 orbit and start tending to a period-4 orbit.


Section 3: Graphical analysis

3.1. Graph of a function.  The goal of this section is to understand the behaviour of the orbits of a given map without having to compute the orbits, but simply by looking at the graph of the function f(x).   To plot the graph of a function in Matlab, we can use the command 'fplot'.  For instance, to plot the function f(x) = 2.5 x (1-x) in the range of x values from 0 to 1, we type

>> fplot('2.5*x*(1-x)',[0 1],'r-');

(The string 'r-' specifies that the graph should be a red line, the default is a blue line).

Alternatively, we can plot the function using 'plot' command: first specify the range of x values:

>> x = 0:0.01:1;

Next, compute the values of the function and plot the result:

>> f = 2.5.*x.*(1-x);
>> plot(x, f, 'r-');

(Note that we use '.*' for multiplication instead of just '*'.  The reason is that '*' is interpreted by Matlab as a vector multiplication.   The '.*' symbol tells Matlab to multiply vectors element by element.)

3.2. The cobweb diagram.  There is a simple geometric procedure for describing the behaviour of orbits using only the graph of f(x).  It is sometimes referred to as the cobweb diagram.  Let's say we want to show on the function graph an orbit starting from x0.  To help us do that, we first draw the diagonal line y = x, which makes a 45o angle with x- and y-axes.

>> hold on; plot([0 1],[0 1],'k-');

(The command 'hold on' holds the current plot of the function, so that the black line from (0, 0) to (1, 1) is added to it).

The next point on the orbit is the number f(x0).  The graph y = f(x) allows us to read off this point, since (x0, f(x0)) is the point on the function graph directly over x0.  To visualise this,  we can draw a vertical line from the point (x0, x0) on the diagonal to the point (x0, f(x0)) on the graph of the function f(x).

>> x0 = 0.1;
>> x = x0; f = 2.5.*x.*(1-x);
>> plot([x x], [x f]);

Next, we need to assign the number f(x0) to the next value of x:  x1 = f(x0) .  In order to do it on the graph, we draw a horizontal line from the point  (x0, f(x0))  till it meets the diagonal line x = y precisely at the point (f(x0),  f(x0)) .

>> plot([x f], [f f]);

We can continue graphing the orbit by assigning the value of f to x:

>> x = f; f = 2.5.*x.*(1-x);

and repeating the above steps again.  All this can be arranged in a loop and written as a Matlab program.  As before, copy the program into a file and run it in Matlab.  Explore different values of r, N and x0.

Program 2

% Plotting the graph of the function and the diagonal
r = 2.5;  N = 10;  x0 = 0.1;
x = 0:0.01:1;
f = r.*x.*(1-x);
clf;
plot(x,f,'r-',[0 1],[0 1],'k-');
hold on;
% Plotting the 'cobweb' diagram
f = x0;
for i = 1:N,
  x = f;
  f = r.*x.*(1-x);
  plot([x x], [x f]);
  pause;
  plot([x f], [f f]);
  pause;
end,

Comments to the program:

In the diagram, each successive point of the orbit appears on the diagonal y = x after one 'reflection' off of the function graph.
 
 

Section 4: Fixed points and periodic orbits

It is immediately clear from the diagram, that for smaller values of r, the orbits converge to the fixed points, which are given by the intersection of the diagonal y = x and the function y = f(x).  Thus, we recognise that the points where

x* = r x* (1-x*)
are the fixed points of the logistic map.  One such point is x* = 0, and orbits converge to it for small values of r (try running Program 2 at r = 0.8).  For larger r, the orbits converge to the other fixed point, which can be easily determined from the equation above for any r:
 x* = 1 - 1/r

In fact, for any map xn+1 = f(xn) it is possible to determine its fixed points simply by plotting the function y = f(x) together with the line y = x  and looking for the intersection points.

Q4.1.  Consider a map with f(x) = 2.5 sin(x2).  Use Matlab to plot the function and the diagonal line and determine the number and approximate values of all the fixed points of this map.


Returning back to the logistic map, we see that it always has two fixed points.  Only one of them can attract orbits, and is therefore stable, while the other one is unstable.  Note however, that for large enough values of r both fixed points become unstable (Plot the cobweb diagrams for r = 3.2 starting at x0 = 0.01 and at x0 = 0.68).  In fact, for r = 3.2 the orbits converge to a period-2 orbit, which is given by two points x1 and x2 such that  f(x1) = x2 and  f(x2) = x1.  Alternatively, we can write that x1= f(f(x1)) or x2= f(f(x2)) and we see that x1 and x2 are fixed points of the map f(f(x)).  To see that this is indeed so, we can plot the function y = f(f(x))  on the cobweb diagram:

>> x = 0:0.01:1; f = r.*x.*(1-x); y = r.*f.*(1-f);
>> plot(x, y, 'g-');

So, for r = 3.2, period-2 orbit is the stable orbit of the logistic map.
 

Q4.2.  What determines the stability of fixed points of a map?  Plot the cobweb diagrams for r = 3.2 starting at x0 = 0.01 and at  x0 = 0.68.  Analyse the diagrams and try to understand what makes the diagram drift away from the fixed points.  Formulate the stability condition for a fixed point in terms of the properties of the function f(x).
Hint: Just like in the case of one-dimensional flows considered in the lectures, the answer has something to do with the derivative (slope) of the function at the fixed point.


For the Lab Report No. 1R, give answers to all the questions above, providing texts of the programs and figures wherever requested.