% lyapunov exponents in a flow in R^3, eg Lorenz attractor
% Needs the function lorenz_time1map.m which returns the time-1 map
% barnett 11/19/07

s = 10;
b = 8/3;
r = 25.0;
F = @(t, y) [ -s*y(1,:)+s*y(2,:); r*y(1,:)-y(2,:)-y(1,:)*y(3,:); -b*y(3,:)+y(1,:)*y(2,:)];   % vec flow field
yo = [0.1;0.1;0.1];                              % IC
[ts, ys]  = ode45(F, [0 300], yo);    % numerically solve in t domain
figure;plot3(ys(:,1), ys(:,2), ys(:,3), '-'); axis vis3d; % show the attractor

% test lorenz time-1 map
[x, J] = lorenz_time1map(yo, s,b,r)



% Measure Lyapunov exponents...
% ----- Re-orthogonalizing version, repeated averaging
M = 20;                              % how many averaging loops
N = 50;                                  % how many its per meas step
x = [0.1;0.1;0.1];
h = zeros(3,1);                    % place to store averaged lyap exps
for m=1:M
  J = eye(3);                   % Id is where Jacobean starts
  for n=1:N
    [x Jx] = lorenz_time1map(x,s,b,r);
    J = Jx*J;                        % update Jacobean
    [Q,R] = qr(J);                   % re-orthogonalize
    J = Q*diag(diag(R));             % but keep them correct lengths
  end
  rN = abs(diag(R))                  % print out progress
  h = h + log(rN)/N;                 % sum up lyap exps from each run
end
h = h/M                          % final answer: mean (over runs) of lyap exps

% you can see the middle lyap exp is v close to 0 - in fact all flows have
% a zero exponent if they are bounded and hit no equilibrium points.