BACK

Lab 4: Lyapunov exponent

Section 1: Strange attractors in two-dimensional maps

In the lectures, you have seen pictures of chaotic dynamics of various two-dimensional maps.  You were told that for some parameter values, the dynamics of such maps does not settle into a fixed point or a periodic orbit, but instead wanders around sets with very peculiar shapes and containing infinitely many points.  Such sets are called strange attractors.  The existence of such sets is a beautiful manifestation of chaotic dynamics in nonlinear systems.  In this Section of Lab 4, you will have a chance to construct such sets using Matlab for two maps: the Henon map and the Ikeda map.

The Henon map is given by the equations:

xn+1= a - xn2 + byn,
yn+1= xn,

where (a,b) are the map parameters, with 0 < b <= 1.  A strange attractor can be observed for parameter many values, e.g.  (a,b) = (1.4,0.3),  (1.0, 0.54), etc.

The Ikeda map is, in fact, a model for a string of self-interacting light pulses emitted by a laser and travelling through a nonlinear medium.  The map equations are:

xn+1= a + b (xncos phi n -  yn sin phi n ),
yn+1= b (xnsin phi n +  yn cos phi n ),

where phi n= k - eta /(1 + xn2 + yn2 ) and (a, b, k, eta ) are the map parameters.  The traditional parameter values for the Ikeda map exhibiting chaotic dynamics on a strange attractor are  (a, b, k, eta ) =  (1.0, 0.9, 0.4, 6.0), but other parameter values also produce strange attractors of different sizes and shapes.


Below is a Matlab function that computes an orbit of the Henon map


Program 1: the Henon map (Matlab file: henon.m)

function [x,y] = henon(a, b, x0, y0, N)
%  This function computes N points of the Henon map
%  with parameter values (a,b) starting from initial point (x0,y0)

x = zeros(1,N); y = zeros(1,N);
x(1) = x0;  y(1) = y0;

for n = 1:N-1,
  x(n+1) = a - x(n).^2 + b*y(n);
  y(n+1) = x(n);
end,

To view the strange attractor for the Henon map, execute the following commands at the Matlab prompt:

>> [x,y] = henon(1.0, 0.54, 0.0, 0.0, 10000);
>> plot(x,y,'.','markersize',2);

Now you can explore the dynamics of the map for different values of the parameters.

Q1.1 Using Program 1 as a guide, write Matlab function ikeda.m that computes an orbit of the Ikeda map.  Plot strange attractors for parameter values:
a) (a, b, k, eta ) =  (1.0, 0.9, 0.4, 6.0) and
b) (a, b, k, eta ) =  (1.0, 0.8, 0.5, 7.0).


Section 2: Lyapunov exponent of one-dimensional maps

1. Sensitive dependence on initial conditions.  A defining attribute of a chaotic orbit is that it displays exponentially sensitive dependence on initial conditions.   To see what this means, do the following exercise.

Q2.1 For the logistic map with r = 3.8  compute two orbits,  x1 and x2, originating from the initial conditions x1(1) = 0.3 and x2(1) = 0.3 + 1e-10  (that is, 0.3000000001).   Plot the two orbits on the same graph.  You will see that the orbits stay close to one another for a number of iterates, but eventually drift apart.  After how many iterates does the separation between the two orbits first exceed 0.2?  To see that the separation grows exponentially fast with the number of iterates, plot the natural logarithm of the separation between the orbits (that is,  log(abs(x2-x1))  ) versus the number of iterates.  What you should see is that this quantity grows approximately linearly until it reaches a plateau just below 0.  Estimate the slope h of the linear growth region.
The separations stops growing because both orbits are bounded withing the interval [0, 1], so the separation between them cannot be larger than 1.  Since log(1) = 0, the logarithm of the separation cannot be larger than zero.  In the linear growth region we see that
log | x2 - x1 |  a + h n ,
where a = log(10-10) is the logarithm of the initial separation between the orbits and n is the number of iterates.  Alternatively, we can write
| x2 - x1 |  A exp(h n),
where A = exp(a) = 10-10 , which shows that the separation indeed grows exponentially with n.   The same result can be obtained for any two nearby initial conditions.  Moreover, for a given parameter r, the value of the coefficient h is always the same.  This coefficient is called the Lyapunov exponent of the dynamical system.  It characterises the rate at which nearby trajectories diverge away from each other.

Obviously, for chaotic orbits h > 0, the Lyapunov exponent is positive, because the separation increases.  However, if we consider parameter values where the orbits converge to a fixed point or a periodic orbit, then the separation between the orbits will decrease.  The rate of decrease will also be exponential.

Q2.2  Compute two orbits of the logistic map with r = 2.9 starting at x1(1) = 0.3 and x2(1) = 0.4.  After a few iterates, the two orbits converge to the same period-2 orbit.  Plot the logarithm of the separation between the two orbits.  Again we see approximately linear dependence of this quantity on the number of iterates, but now the slope of the line is negative.  Estimate the value of the slope.
It is clear that the Lyapunov exponent allows us to distinguish between different types of orbits in dynamical systems.  Therefore, it is useful to compute the Lyapunov exponent whenever exploring the behaviour of a dynamical system.
 

2.1. Computing Lyapunov exponent.  Obviously, estimating the Lyapunov exponent with the method described in Q2.1 is far from efficient.  Besides, the exponential growth of the separation between two orbits continues only while the separation is smaller than the system size, so it is hard to accurately estimate the Lyapunov exponent from a short segment of exponential growth.  In order to have a more efficient method for computing Lyapunov exponents it is necessary to use some analysis.

Consider two orbits starting at x0 and  x0 +  0 , where  0 is infinitesimally small. (Roughly speaking, an infinitesimal number is so small that we can multiply it by an arbitrarily large number and the result will still be infinitesimally small. Also it means that in equations, we can neglect all powers of infinitesimals higher than one, that is, if  is infinitesimal and a, b, c are some non-zero numbers, then a + b 2 + c 3 = a .)  The next points of the two orbits are

x1 = f(x0) 
 (1)
and
x1 +  1 = f(x0 +  0).
 (2)
Note that, since  0 is infinitesimal, we can use the Taylor series for f(x0 +  0) at x0 to write equation (2) as
x1 + delta 1 = f(x0) + f'(x0)  0,
where f'(x0) is the derivative of f evaluated at x0.  Therefore, using equation (1), we have
1 = f '(x0)  0 .
Similarly, for the n-th point along the orbits,  xn = f(xn-1) ,  xn +  n = f(xn-1 +  n-1), so that
 
n = f'(xn-1)  n-1 .
 (3)

Equation (3) will be useful in deriving the formula for computing the Lyapunov exponent of the map f(x).  Endeed, as we have seen in exercises Q2.1 and Q2.2, the logarithm of the separation between the two orbits starting at x0 and  x0 +  0 changes approximately linearly with n, that is,

log|  n|    log|  0| + h n ,
Where h is the Lyapunov exponent.  From this equation, we can estimate the exponent as

h = (1/n) log |  n/ 0| .

for some large n.  Now, we use equation (3) to rewrite the ratio  n/ 0 as follows

n/ 0 = ( n/ n-1)( n-1/ n-2)...( 2/ 1)( 1/ 0) = f'(xn-1) f'(xn-2) ... f'(x1) f'(x0) .
Finally, the formula for computing the Lyapunov exponent is given by
 
h = (1/n) {log | f'(xn-1)| + log | f'(xn-2)| + ... + log | f'(x1)| + log | f'(x0)|} .
 (4)

Now, all we have to do is to write a Matlab program that will compute the Lyapunov exponent using equation (4).  To this end,

Q2.3  Write a Matlab function starting with the line
function h = lyapexp(r, x0, Npre, n);
which computes the Lyapunov exponent of the logistic map with the parameter r.  Since we want the Lyapunov exponent for the eventual state of the system (same as for the bifurcation diagram), the orbit is started at x0, but we compute Npre pre-iterates before starting the computation of the Lyapunov exponent according to equation (4) during the next n iterates.
Lab Report No. 4 is due on Thursday, May 15th at 11:30am.

BACK