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:
x_{n+1}= a - x_{n}^{2} + by_{n},
y_{n+1}= x_{n},
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:
x_{n+1}= a + b (x_{n}cos _{n} - y_{n} sin _{n} ),
y_{n+1}= b (x_{n}sin _{n} + y_{n} cos _{n} ),
where _{n}= k - /(1 + x_{n}^{2} + y_{n}^{2} ) and (a, b, k, ) are the map parameters. The traditional parameter values for the Ikeda map exhibiting chaotic dynamics on a strange attractor are (a, b, k, ) = (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, ) = (1.0, 0.9, 0.4, 6.0) and
b) (a, b, k, ) = (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.
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 x_{0} and x_{0} + _{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
x_{1} = f(x_{0}) |
(1) |
x_{1} + _{1} = f(x_{0} + _{0}). |
(2) |
x_{1} + _{1} = f(x_{0}) + f'(x_{0}) _{0},where f'(x_{0}) is the derivative of f evaluated at x_{0}. Therefore, using equation (1), we have
_{1} = f '(x_{0}) _{0 }.Similarly, for the n-th point along the orbits, x_{n} = f(x_{n-1}) , x_{n} + _{n} = f(x_{n-1} + _{n-1}), so that
_{n} = f'(x_{n-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 x_{0} and x_{0} + _{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'(x_{n-1}) f'(x_{n-2}) ... f'(x_{1}) f'(x_{0}) .Finally, the formula for computing the Lyapunov exponent is given by
h = (1/n) {log | f'(x_{n-1})| + log | f'(x_{n-2})| + ... + log | f'(x_{1})| + log | f'(x_{0})|} . |
(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 lineLab Report No. 4 is due on Thursday, May 15th at 11:30am.
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.