
% Computation of Lyapunov Exponents for a discrete
% dyn. syst. using QR factorization
% RCG Nov'14

format compact;

%params
Ntrans=1000; % transients
Nit=1000; % number of iterations
disp=100; % points todisplay for convergence
fs=20;

%Define map (column vector)
Henon_map = inline('[a-P(1).^2+b.*P(2);P(1)]','a','b','P');
a = 1.4; b = 0.3;

%Separation for contraction
epsilon = 1e-8;

%ICs
x0 = 0.1; y0 = 0.1;
H=[x0;y0];
sumLog = 0;
allLE=[];
allt=[];
figure(1);clf;set(gca,'FontSize',fs)

%Transients
for ii = 1:Ntrans
    H = Henon_map(a,b,H);
end

%Initial basis
p0 = H;
p1 = H + epsilon.*[1;0];
p2 = H + epsilon.*[0;1];

%Mainloop
for jj = 1:Nit
    H_p0 = Henon_map(a,b,p0);
    H_p1 = Henon_map(a,b,p1);
    H_p2 = Henon_map(a,b,p2);
    
    M = [H_p1-H_p0, H_p2-H_p0];
    [Q,R] = qr(M);
    sumLog = sumLog + log(abs(diag(R))./epsilon); 
    
    if(mod(jj,Nit/disp)==0)
     allt=[allt jj];
     allLE=[allLE sumLog/jj];
     plot(allt,allLE,'o',[0 allt(end)],kron(allLE(:,end),ones(2,2)),'k-')
     title(['LE: ',num2str(allLE(:,end)',8)],'FontSize',fs)
     xlabel('n','FontSize',fs);ylabel('LE','FontSize',fs);drawnow
    end
    
    p0 = H_p0;
    p1 = H_p0 + epsilon.*Q(:,1);
    p2 = H_p0 + epsilon.*Q(:,2);
end

LE = sumLog/Nit
plot(allt,allLE,'o',[0 allt(end)],kron(allLE(:,end),ones(2,2)),'k-')
title(['LE: ',num2str(LE',8)],'FontSize',fs)
xlabel('n','FontSize',fs);ylabel('LE','FontSize',fs)


