% 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