Calculate Lyapunov spectrum for Lorenz system

61 views (last 30 days)
F.O
F.O on 30 Aug 2020
Commented: jiacheng feng on 9 Feb 2023
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?
  2 Comments
Anderanik Rafi
Anderanik Rafi on 12 Sep 2021
Edited: Anderanik Rafi on 12 Sep 2021
Hi.
You can use the attached code. I tested it, it works
jiacheng feng
jiacheng feng on 9 Feb 2023
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?

Sign in to comment.

Answers (3)

Alan Stevens
Alan Stevens on 31 Aug 2020
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!
  2 Comments
F.O
F.O on 31 Aug 2020
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
Alan Stevens on 31 Aug 2020
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.

Sign in to comment.


Mujeeb Oyeleye
Mujeeb Oyeleye on 10 Oct 2021
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
  1 Comment
Lazaros Moysis
Lazaros Moysis on 4 Nov 2021
So what is the correction you propose? I cannot point it out. I am having the same issue.

Sign in to comment.


Lazaros Moysis
Lazaros Moysis on 4 Nov 2021
I am having the same trouble with this code. DId anyone eventually debug it?

Categories

Find more on Matrix Computations in Help Center and File Exchange

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!