Program 4 - Lyapunov Exponents

In this program we first compute the Lyapunov exponent of an orbit of a one-dimensional system. We choose an initial condition, and we consider a significant number of states of the orbit starting at the selected initial condition. As we move along the orbit we compute the Lyapunov exponent up to the iteration where we are. We see that after an initial period of uncertainty, the Lyapunov exponent settles on a specific value. Since it is positive, the orbit is unstable.

f[x_,a_]:= a x(1-x)

g[x_,a_]:= a(1-2x)

x=.3;a=4;n=1;b=0;

q = Table[{ b=b+Log[Abs[g[x,a]]],n=n,b/n,n=n+1,x=f[x,a]}, {600}];

p=q[[Range[1,600],Range[2,3]]];

ListPlot[p, PlotStyle->PointSize[0.015],PlotRange->{Automatic, {0,1}}];

Now we compute the Lyapunov exponent of an orbit of a two-dimensional system. The function governing the system is
 F(x,y)=(3.67x(1–.3x–.7y),x).

f1[x_, y_] := 3.67  x (1 - .7  x - .3  y); f2[x_] := x;

x = .3; y = .2; n = 1;

m[x,y] is the derivative of the function F, q is the identity matrix of size 2.

m[x_, y_]
         := {{3.67 - 3.67 1.4  x - 3.67 .3  y, -3.67 .3  x}, {1, 0}}; 
    q = IdentityMatrix[2];

The program first compute the product of the derivatives of F at the states of the orbit in the appropriate order (namely the last derivative should go first), then computes the Lyapunov exponent up to the given state.

q1 = 
      Table[{q = m[x, y].q, n, Log[(Max[Abs[Eigenvalues[q]]])^(1/n)], 
          x = f1[u = x, y], y = f2[u], n = n + 1}, {400}];

q2 = q1[[Range[1, 400], Range[2, 3]]];

ListPlot[q2, PlotStyle -> PointSize[0.02],  
      PlotRange -> {Automatic, {-1, 1}}];

We see that the Lyapunov Exponent is slowly stabilizing at approximately -0.2. Consequently, the orbit is stable.