% HW3 code solns: Matlab code to measure lyapunov exponent for logistic
% map.    barnett 10/17/07

% ---------------------    METHOD 1- using 2 different orbits.
eps = 1e-8;   % small but not too small (geometric mean of 1 and 1e-16)
as = 0:0.001:4;
hs = zeros(size(as));
for j=1:numel(as)            % note this approach lets you index hs and as
  a = as(j);                 % choose the current `a' value.
  f = @(x) a*x.*(1-x);       % logistic map
  
  N = 30;                   % number of iterations (keep small or get log of 0)
  x = zeros(1, N+1);         % 1-by-(N+1) sized list to store the history of x
  y = zeros(1, N+1);
  x(1) = 0.3;               % IC shouldn't matter much, but actually seems to!
  y(1) = x(1) + eps;        % slightly different IC
  
  for n=1:N
    x(n+1) = f(x(n)); y(n+1) = f(y(n)); 
  end
  
  hs(j) = log(abs(x(end)-y(end)) / eps) / N;     % h is slope of log ratio
end
figure; plot(as, hs, '+-'); xlabel('a'); ylabel('lyap exp h');

% notice that for some values of a, x_k-y_k became too small, so difference
% was rounded off to zero, and the datapoints are missing. It's still ok.
% However, only -0.6<h<0.6 are correctly computed (due to reaching 1e-16 or 1)

% ----------------------- METHOD 2- using derivative along 1 orbit (better)
as = 0:0.001:4;    % note this takes about 10 sec to run!
hs = zeros(size(as));
for j=1:numel(as)
  a = as(j);
  f = @(x) a*x.*(1-x);
  df = @(x) a*(1-x) - a*x;          % derivative of logistic map
  
  N = 100;                   % number of iterations (can use more than above)
  x = zeros(1, N+1);
  x(1) = 0.1;                % IC
  deriv = 1;                 % start counting the 
  
  for n=1:N
    x(n+1) = f(x(n));
    deriv = deriv * df(x(n));        % update product of derivatives
  end

  hs(j) = log(abs(deriv)) / N;     % implement Defn 3.1 in book.
end
figure; plot(as, hs, '+-'); xlabel('a'); ylabel('lyap exp h'); % as before
% Note that all lyapunov exponents are well estimated here