Archive for October, 2013

Liapunov exponent computation

October 4, 2013 Leave a comment

Liapunov exponent method is one of the clever ways to investigate a certain dynamical system to be chaotic. In this post, I shall try to compute Liapunov exponent of one famous discrete dynamical system, i.e. the logistic map. After this, I shall try to compute the Liapunov exponent of a continuous dynamical system.

Actually, I write this post because I want to know how to numerically compute the Liapunov exponent. I often read a book or a website such as wolfram.alpha that shows the Liapunov exponent but it didn’t explain how to get this numerically. After I am able to post this post, and when I forget how to compute this exponent, I’d know where to look.

Suppose we have an initial condition x_0 and a nearby initial point x_0 + \delta_0, where \delta_0 is an initial difference between initial points, which by the way will be extremely small. These two nearby points will be iterated by some map. Let \delta_n be the difference between these two points after n iterations. Consider the following equation

(1) …  |\delta_n| =|\delta_0 |e^{n \lambda}.

The coefficient \lambda is called the Liapunov exponent. We want to compute this coefficient. Note that \delta_n = f^n(x_0)-f^n(x_0+\delta_0), where f is the given map. Taking logarithm of (1), we get

(2) … \lambda = \frac{1}{n} \ln \left|\frac{\delta_n}{\delta_0}\right|=\frac{1}{n} \ln \left|\frac{f^n(x_0)-f^n(x_0+\delta_0)}{\delta_0}\right|=\frac{1}{n} \ln \left|{(f^n)'(x_0)}\right|

as \delta_0 \to 0. The term inside the logarithm can be simplified by using the chain rule: (f^n)'(x_0) = f'(x_{n-1}) f'(x_{n-2}) \dots f'(x_0). Thus (2) becomes

(3) … \lambda = \frac{1}{n} \sum_{i=0}^{n-1} \ln |f'(x_i)|.

When n \rightarrow \infty, \lambda is defined as the Liapunov exponent. When this exponent is positive, it is a sign of a chaotic attractor.

Let us consider the logistic map that was in my previous post here. Let us fix the parameter value r. Starting at a random initial point, we will iterate the initial point until it reaches the steady state (which will be a fixed point, periodic or aperiodic), say 300 iterates. Now we compute a large number iterates, say 10000 and compute (3) starting from the 301st point. Repeat this procedure for the next r. Then we plot the Liapunov exponent with respect to r. We shall find a figure as follows.


As we can see from the above figure. The map is chaotic in some values of parameter between 3.5 to 4, in which the Liapunov exponent is positive.

The code to produce the above figure is the following

tic; %to compute how long is the computation
clear;clc; %clear screem and clear all the previous variables
a = 3:0.001:4; %parameter of the logistic map
N = 300; %number of iterations until it is stable
N1 = 10000; %number of iterations to compute the Liapunov exponent
y = zeros(length(a),1); %declares variable for Liapunov exponents
for i=1:length(a)
x(1) = rand; %random initial conditions%
for j = 2:N
x(j) = a(i)*x(j-1)*(1-x(j-1));
sum = 0;
for j = 1:N1
x(N+j) = a(i)*x(N+j-1)*(1-x(N+j-1));
sum = sum + log(abs(a(i)-2*a(i)*x(N+j)));
y(i) = sum/10000;
fprintf('%d %d\n',i,length(a));
plot(a,0); hold on;
plot(a,y); %plot Liapunov exponent against parameter
toc; %display how long computation is