Calculate Lyapunov spectrum for Lorenz system

52 vues (au cours des 30 derniers jours)
F.O
F.O le 30 Août 2020
Commenté : jiacheng feng le 9 Fév 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 commentaires
Anderanik Rafi
Anderanik Rafi le 12 Sep 2021
Modifié(e) : Anderanik Rafi le 12 Sep 2021
Hi.
You can use the attached code. I tested it, it works
jiacheng feng
jiacheng feng le 9 Fév 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?

Connectez-vous pour commenter.

Réponses (3)

Alan Stevens
Alan Stevens le 31 Août 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 commentaires
F.O
F.O le 31 Août 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 le 31 Août 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.

Connectez-vous pour commenter.


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

Connectez-vous pour commenter.


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

Catégories

En savoir plus sur Matrix Computations dans Help Center et File Exchange

Produits


Version

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by