%{
calculation of the largest lyapunov exponent (lle) by the method suggested by J.C.Sprott,
Department of Physics, University of Wisconsin, Madison, WI 53706, USA
%}
% http://sprott.physics.wisc.edu/chaos/lyapexp.htm
%{
used to find the lyapunov exponent spectrum of a system of nonlinearly
coupled oscillators as the coupling parameter is varied
%}
global be;
global gm;
%'gm' and 'be' are the parameters of the system
be=1;
c1=1;
lya3=zeros(101,2); %stores the lle and the corresponding value of gm
for gm=10:0.05:40 %any range may be chosen by the user
c=1;
t0=0;
tf=8;
N=80;
E=0.2; px=0.2;
x1=0;
x2=px;
x3=0;
x4=sqrt(2*E-(px^2));
%(x1,x2,x3,x4) and (x11,x22,x33,x44) are arbitrarily close initial conditions
x11=10^-8;
x22=px+10^-8;
x33=10^-8;
x44=x4+10^-8;
d0=sqrt((x11-x1)^2+(x22-x2)^2+(x33-x3)^2+(x44-x4)^2); % initial separation between the trajectories
lam=zeros(401,1); % stores the exponent calculated at each iteration
for i=1:400
%solving the equations of the system of interest for each of the trajectories
[t,x] = ode_RK4_Yang([t0 tf],[x1 x2 x3 x4],N);
[tt,xx] = ode_RK4_Yang([t0 tf],[x11 x22 x33 x44],N);
d1=sqrt((xx(81,1)-x(81,1))^2+(xx(81,2)-x(81,2))^2+(xx(81,3)-x(81,3))^2+(xx(81,4)-x(81,4))^2); %separation at the end of each iteration
lam(c,1)=log(d1/d0);
%reinitialization of the trajectories after each iteration
x11=x(81,1)+(xx(81,1)-x(81,1))*(d0/d1);
x22=x(81,2)+(xx(81,2)-x(81,2))*(d0/d1);
x33=x(81,3)+(xx(81,3)-x(81,3))*(d0/d1);
x44=x(81,4)+(xx(81,4)-x(81,4))*(d0/d1);
x1=x(81,1);
x2=x(81,2);
x3=x(81,3);
x4=x(81,4);
t0=tf;
tf=tf+8;
c=c+1;
end
lya3(c1,1)=mean(lam(100:400)); % finally lle is the mean of all the calculated exponents
lya3(c1,2)=gm;
c1=c1+1;
end
