33 views (last 30 days)

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?

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!

Alan Stevens
on 31 Aug 2020

Opportunities for recent engineering grads.

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

Start Hunting!
## 0 Comments

Sign in to comment.