lyapunovexponentbyslope.m

用斜率的方法来求时间序列的lyapunov exponent,在非线性动力学和时间序列处理中经常用到。很实用的小程序。

% Below the power method will be applied to the jacobian matrix of the
% 2-D henon map to approximate the first Lyapunov exponent by creating
% a graph of ln|yn| vs. n, where n is the number of iterations of the
% power method and yn = 1/n*ln|DG^n(xo)*yo|. The slope will be an 
% approximation to the largest Lyapunov exponent.

% get rid of transients

NN = 100;
xo = 1;
yo = .52;

for i = 1:NN,
	  if (i == 1)
         	  [x,y] = henon(xo,yo);
	  else
                  [x,y] = henon(x,y);
          end %if statement
end  %for loop



% choose arbitrary yo vector. Good practice to make sure it contains
% all components since nothing special is known about matrix itself.

%N = number of iterations
N = 1500;

%assign the starting values to the last iterate from above
xo = x;
yo = y;

%initialize your vectors for Henon Map values of x and y
x = zeros(1,N);
y = zeros(1,N);

%initialize the Jacobian matrix to the starting values
DGo = [(-2*xo),0.3;1,0];

%choose an arbitrary vector, yno
yno = [1;.2];

%initialize the yn vector for storage of ln(yn) data and n for graphing
yn = zeros(N+1,2);
n = zeros(1,N+1);

%loop which iterates the jacobian matrix using Henon values calculated above

for i = 1:N,
	  if (i == 1)
         	  [x(i),y(i)] = henon(xo,yo);
	  else
                  [x(i+1),y(i+1)] = henon(x(i),y(i));
          end %if statement
end  %for loop

% power method being applied to matrix
for j = 1:N,
    n(j) = j;
    if (j == 1)
     temp = DGo*yno;
     yn(1,1) = log(abs(temp(1,1)));
     yn(1,2) = log(abs(temp(2,1)));
   
    else
      % calculate jacobian matrix
         DG = [(-2*x(j)),0.3;1,0];
         temp = DG*temp;
         yn(j+1,1) = log(abs(temp(1,1)));
         yn(j+1,2) = log(abs(temp(2,1)));
    end % if statement
     
end %for loop
n(N+1) = N+1;

% graph both elements of the vector, both slopes will be dominate or largest
% lyapunov exponent.

for k = 1:2,
  figure
  hold
  grid on
  plot(n,yn(:,k),'-.b');
        
  title('Plot of yn vs. n for n = 1500')
  xlabel('n integer values')  
  ylabel('ln(yn(1)) (arbitrary units)')
     for j = 1:N-1,
       ynew(j) = yn(j+1,k);
     end

  nsquare = 0;
  ytotal = 0;
  nyntotal = 0;
  ntotal = 0;

    for i = 1:N,
	      ntotal = ntotal + i;
	      nsquare = nsquare + i*i;
              ytotal = ytotal + yn(i,k);
              nyntotal = nyntotal + i*yn(i,k);
     
    end

	      % using least squares, calculate slope
	      slope = (N*nyntotal - ntotal*ytotal)/(N*nsquare - ntotal*ntotal)
end