Calculate Lyapunov spectrum for Lorenz system

48 次查看（过去 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 中查找有关 State-Space Control Design 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by