%
% Example program for invariant distribution and Lyapunov exponent of
% logistic map
%
% Jean-Luc Thiffeault (jeanluc@math.wisc.edu)
%

Npts = 20000; Nit = 50; nbins = 50; r = 4;
x = rand(1,Npts);       % random initial points
f = @(r,x) r.*x.*(1-x); % logistic map

figure(1)

% the initial points are uniformly distributed on the interval [0,1]
subplot(2,1,1)
[P,xb] = hist(x,nbins); P = P./trapz(xb,P);
bar(xb,P,1)
title('initial distribution')

% now iterate the map on the initial points
for i = 1:Nit, x = f(r,x); end

% the final points obey the invariant distribution
subplot(2,1,2)
[P,xb] = hist(x,nbins); P = P./trapz(xb,P);
bar(xb,P,1)
title(sprintf('distribution after %d iterations',Nit))
if r == 4
  hold on
  xi = linspace(0,1,300); xi = xi(2:end-1);
  plot(xi,1/pi*(xi.*(1-xi)).^-.5,'r','LineWidth',2) % analytic answer for r=4
  hold off
end

figure(2)

% two trajectories very close to each other initially
delta = 1e-10; x1 = [.1]; x2 = [x1 + delta];

% iterate the map on both trajectories
for i = 1:Nit, x1 = [x1 f(r,x1(end))]; x2 = [x2 f(r,x2(end))]; end

% plot the separation
semilogy(0:Nit,abs(x1-x2))
hold on
semilogy(0:Nit,delta*2.^(0:Nit),'r') % compare to 2^n separation rate
hold off
axis([0 Nit delta 10])
xlabel('n')
ylabel('distance between trajectories')

figure(3)

Nit = 2000; Nplot = 1000; r = linspace(0,4,400); x = .1*ones(size(r));

% iterate the map on the same initial condition, but different r
for i = 1:Nit, x = [x; f(r,x(end,:))]; end

clf, hold on
% plot the last Nplot iterates for each r
for i = 1:Nplot
  plot(r,x(Nit-i-1,:),'.','MarkerSize',2)
end
hold off
xlabel('r')
title('bifurcation diagram for logistic map')
