% measure lyapunov exponents in map in R^2 using Henon as example
% Barnett 10/29/15

a = 1.4; b = 0.3;
f = @(x) [a-x(1,:).^2+b*x(2,:);  x(1,:)]; % Henon map (vec func) w/ params a,b
Df = @(x) [-2*x(1,:) b; 1 0];              % Jacobean matrix as func of x

% ----- quick and dirty: meas largest lyap exp h_1 only 
N = 100;                                  % how many its (start at 100 & grow)
                                          % at 2000, overflows!
x = zeros(2, N+1);
x(:,1) = [0;0];
u = [1;0];                  % arbitrary starting tangent vector
for n=1:N
  x(:,n+1) = f(x(:,n));
  u = Df(x(:,n))*u;
end
plot(x(1,:), x(2,:), '.', 'markersize', 10);            % show orbit
axis equal; axis([-2 2 -2 2]);

r1 = norm(u)             % exp large but make sure doesn't exceed 1e+308 !
h1 = log(r1)/N           % the answer



% ----- Try to get all lyapunov exponents: Crude version
N = 40;                                   % how many its
x = zeros(2, N+1);
x(:,1) = [0;0];
J = eye(2);
for n=1:N
  x(:,n+1) = f(x(:,n));
  J = Df(x(:,n))*J;
end
rN = svd(J)                   % extract semiaxes lengths of J = J_n
h = log(rN)/N
% need to be able to run for longer so averaged over the chaotic attractor
% BUT have overflow giving bad accuracy on h_2 especially.
% Garbage for N>50



% ----- Re-orthogonalizing version, Alligood-Sauer-Yorke Sec 5.2
N = 400;                                  % how many its
x = zeros(2, N+1);
x(:,1) = [0;-1];
J = eye(2);                   % orthogonal set as column vectors
for n=1:N
  x(:,n+1) = f(x(:,n));
  J = Df(x(:,n))*J;
  [Q,R] = qr(J);                   % re-orthogonalize
  J = Q*diag(diag(R));             % but keep them correct lengths
end
rN = abs(diag(R))
h = log(rN)/N
% better but still breaks down for N>400 since the h_2 direction causes
% to underflow (smaller than 1e-308) in the diagonal elements of R

% ----- optional: Eckmann-Ruelle 1985 Sec V.C, re-orthogonalizing version
N = 400;                                  % how many its
x = zeros(2, N+1);
x(:,1) = [0;-1];
rN = ones(2,1); Q = eye(2);      % initial diag(R) entries, and initial Q
for n=1:N
  x(:,n+1) = f(x(:,n));
  [Q,R] = qr( Df(x(:,n))*Q );
  rN = abs(diag(R)).*rN;
end
h = log(rN)/N
% Seem to be *exactly* same as previous re-orthog version! Not sure why



% ----- Best: Re-orthogonalizing version with averaging over long traj
M = 100;            % number of repeats of the N iterations
N = 400;                                  % how many its (NM total)
x = [0;1];                     % note: don't keep whole history any more
h = zeros(2,1);                 % cumulative lyap exp storage
for m=1:M
  J = eye(2);                   % orthogonal set as column vectors
  for n=1:N
    J = Df(x)*J;
    x = f(x);
    [Q,R] = qr(J);                   % re-orthogonalize
    J = Q*diag(diag(R));             % but keep them correct lengths
  end
  rN = abs(diag(R));
  h = h + log(rN)/N;                 % add in the contribs from this run
end
h = h/M
% final estimate is within 1% of values in the book.

  
