Calculate Lyapunov spectrum for Lorenz system

75 次查看（过去 30 天）
F.O
F.O 2020-8-30
评论： jiacheng feng 2023-2-9
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 个评论
jiacheng feng
jiacheng feng 2023-2-9
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?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（3 个）

Alan Stevens
Alan Stevens 2020-8-31
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 个评论
F.O
F.O 2020-8-31
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 2020-8-31
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.

请先登录，再进行评论。

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

请先登录，再进行评论。

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

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Computations 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by