Calculate Lyapunov spectrum for Lorenz system

Hello,

I am trying to use the following code https://www.math.tamu.edu/~mpilant/math614/Matlab/Lyapunov/LorenzSpectrum.pdf to calculate the spectrum of Lyapunov for a system of 3 ODE but the code gives wrong results. However in the paper the results seems to be reasonable. What is going wrong in the following code:

function lorenz_demo(time)
 [t,x] = ode45('g',[0:0.01:time],[1;2;3]);
 save x
 disp('press any key to continue ...')
 pause
 plot3(x(:,1),x(:,2),x(:,3))

 function xdot = g(t,x)
 xdot = zeros(3,1);
 sig = 10.0;
 rho = 28.0;
 bet = 8.0/3.0;
 xdot(1) = sig*(x(2)-x(1));
 xdot(2) = rho*x(1)-x(2)-x(1)*x(3);
 xdot(3) = x(1)*x(2)-bet*x(3);
 

 function lorenz_spectra(T,dt)
% Usage: lorenz_spectra(T,dt)
% T is the total time and dt is the time step
% parameters defining canonical Lorenz attractor
sig=10.0;
rho=28;
bet=8/3;
%T=100; dt=0.01; %time step
N=T/dt; %number of time intervals
% calculate orbit at regular time steps on [0,T]
% using matlab's built-in ode45 runke kutta integration routine
% begin with initial conditions (1,2,3)
x1=1; x2=2; x3=3;
% integrate forwards 10 units
[t,x] = ode45('g',[0:1:10],[x1;x2;x3]);
n=length(t);
% begin at this point, hopefully near attractor!
x1=x(n,1); x2=x(n,2); x3=x(n,3);
[t,x] = ode45('g',[0:dt:T],[x1;x2;x3]);
e1=0;
e2=0;
e3=0;
% show trajectory being analyzed
plot3(x(:,1),x(:,2),x(:,3),'.','MarkerSize',2);
JN = eye(3);
w = eye(3)
J = eye(3);
for k=1:N
% calculate next point on trajectory
x1 = x(k,1);
x2 = x(k,2);
x3 = x(k,3);
% calculate value of flow matrix at orbital point
% remember it is I+Df(v0)*dt not Df(v0)
J = (eye(3)+[-sig,sig,0;-x3+rho,-1,-x1;x2,x1,-bet]*dt);
% calculate image of unit ball under J
% remember, w is orthonormal ...
w = orth(J*w);
% calculate stretching
% should be e1=e1+log(norm(w(:,1)))/dt; but scale after summing
e1=e1+log(norm(w(:,1)));
e2=e2+log(norm(w(:,2)));
e3=e3+log(norm(w(:,3)));
% e1=e1+norm(w(:,1))-1;
% e2=e2+norm(w(:,2))-1;
% e3=e3+norm(w(:,3))-1;
% renormalize into orthogonal vectors
w(:,1) = w(:,1)/norm(w(:,1));
w(:,2) = w(:,2)/norm(w(:,2));
w(:,3) = w(:,3)/norm(w(:,3));
end
% exponent is given as average e1/(N*dt)=e1/T
e1=e1/T; % Lyapunov exponents
e2=e2/T;
e3=e3/T;
l1=exp(e1); % Lyapunov numbers
l2=exp(e2);
l3=exp(e3);
[e1,e2,e3]
[l1,l2,l3]
trace=e1+e2+e3

What I get is the following:

>> lorenz_spectra(1000,0.001)

w =

     1     0     0
     0     1     0
     0     0     1


ans =

   1.0e-14 *

   -0.1992   -0.2569   -0.3375


ans =

    1.0000    1.0000    1.0000


trace =

  -7.9360e-15

The results are not reasonable and not the same s in the paper. Could anyone help out with dicovering the source of the this?

[Anderanik Rafi]

Hi.

You can use the attached code. I tested it, it works

[jiacheng feng]

Hello!What method does your code use to calculate the Lyapunov?I'm sorry for my poor foundation. Your code is too complicated for me. Could you please give me a more detailed explanation?

Answer 1

[Alan Stevens]

It seems to be the

w = orth(J*w);

command that creates the problem! if you replace it with

w = J*w;

you get more reasonable looking values. Hwever, they don't match those in the paper, and are almost certainly still not correct!

[F.O]

In doing so all the values became the same and this is not possibe. They should be approximatly (0.9, 0, -14) which are the refenece value in many papers.

[Alan Stevens]

I got answers that weren't all the same (though they were similar), and I noted that they were probably not correct (I used T = 10 and dt = 0.001). Unfortunately I couldn't see any other way of getting away from e1, e2 and e3 being of the order of 10^-14 or so.

Answer 2

[Mujeeb Oyeleye]

function lorenz_spectra(T,dt) % Usage: lorenz_spectra(T,dt) % T is the total time and dt is the time step % parameters defining canonical Lorenz attractor sig=10.0; rho=28; bet=8/3; %T=100; dt=0.01; %time step N=T/dt; %number of time intervals % calculate orbit at regular time steps on [0,T] % using matlab's built-in ode45 runke kutta integration routine % begin with initial conditions (1,2,3) x1=1; x2=2; x3=3; % integrate forwards 10 units [t,x] = ode45('g',[0:1:10],[x1;x2;x3]); n=length(t); % begin at this point, hopefully near attractor! x1=x(n,1); x2=x(n,2); x3=x(n,3); [t,x] = ode45('g',[0:dt:T],[x1;x2;x3]); e1=0; e2=0; e3=0; % show trajectory being analyzed plot3(x(:,1),x(:,2),x(:,3),'.','MarkerSize',2); JN = eye(3); w = eye(3) J = eye(3); for k=1:N % calculate next point on trajectory x1 = x(k,1); x2 = x(k,2); x3 = x(k,3); % calculate value of flow matrix at orbital point % remember it is I+Df(v0)*dt not Df(v0) J = (eye(3)+[-sig,sig,0;-x3+rho,-1,-x1;x2,x1,-bet]*dt); % calculate image of unit ball under J % remember, w is orthonormal ... w = J*w; % calculate stretching % should be e1=e1+log(norm(w(:,1)))/dt; but scale after summing e1=e1+log(norm(w(:,1))); e2=e2+log(norm(w(:,2))); e3=e3+log(norm(w(:,3))); % e1=e1+norm(w(:,1))-1; % e2=e2+norm(w(:,2))-1; % e3=e3+norm(w(:,3))-1; % renormalize into orthogonal vectors w(:,1) = w(:,1)/norm(w(:,1)); w(:,2) = w(:,2)/norm(w(:,2)); w(:,3) = w(:,3)/norm(w(:,3)); end

[Lazaros Moysis]

So what is the correction you propose? I cannot point it out. I am having the same issue.

Answer 3

[Lazaros Moysis]

I am having the same trouble with this code. DId anyone eventually debug it?