Lyapunov Exponent

Question

Does anyone know a (simple) Mathematica code for computing the Lyuponov Exponent for the Rossler System?

Thank you

Rossler System:

rossler = {

 x'[t] == -(y[t] + z[t]),

 y'[t] == x[t] + 0.1 y[t],

 z'[t] == 0.2 + x[t] z[t] - 5.7 z[t],

 x[0] == 1, y[0] == 1, z[0] == 1

}

Answer 1

You can use the LCE package by Dr Sandri Marco. It is updated for version 7 and I just tried in on your system for V9 and it worked.

Download lcm.zip and use the package as instructed

Here is the result of running your system on it on my PC

<< lce.m

?LCEsC

These are the 3 Lyapunov Exponents for the Rossler system:

rossler[{x_, y_, z_}] := {-y - z, x +0.1 y, 0.2 + z (x - 5.7)};
x0 = {1,1,1};
T = 0.2; K = 2000; TR = 1;  stepsize = 0.001;
lcesrossler = LCEsC[rossler, x0, T, K, TR, stepsize ]
LyapunovDimension[First[lcesrossler]]
T = 100; TR = 20;
PhaseSpaceC[rossler, x0, T, TR, stepsize, {1, 2, 3}]

giving:

{0.0647984, 0.00535441, -5.23912}

These are close to the known values for this system as given at Prof Sprott site where he has done research in this. From the above page:

Lyapunov exponents (base-e): = 0.0714, 0, -5.3943

By changing the computation parameters and if you have more time to wait, the result can be improved more to become closer to the known values.

Answer 2

At first, let's solve the system of ODEs. Taking into account that you give specific initial conditions, then the solution of the system will correspond to a three-dimensional orbit. The following code solves the system of the ODEs and also plots the output 3D orbit. Now, about the Lyapunov Exponent. How exactly do you define this exponent. I mean, by using the variational equations or by monitoring the deviation between two initially nearby orbits? If it is the latter, then I could provide such a Mathematica code.

Clear["Global`*"];
deq1 = -(y[t] + z[t]);
deq2 = x[t] + 0.1 y[t];
deq3 = 0.2 + x[t] z[t] - 5.7 z[t];
x0 = y0 = z0 = 1;
tin = 0;
tfin = 50;
sol = NDSolve[{x'[t] == deq1, y'[t] == deq2, z'[t] == deq3, 
x[0] == x0, y[0] == y0, z[0] == z0}, {x[t], y[t], z[t]}, {t, tin, 
tfin}];
xt = x[t] /. sol[[1]];
yt = y[t] /. sol[[1]];
zt = z[t] /. sol[[1]];
P1 = ParametricPlot3D[{xt, yt, zt}, {t, tin, tfin}, 
AxesLabel -> {"x", "y", "z"}, BoxRatios -> {1, 1, 1}, 
PlotRange -> All]

Answer 3

Here is a sample code on how to compute the evolution of the Lyapunov Characteristic Exopnent (LCE) with Mathematica. Feel free to make any changes you like and let me know if this is what you wanted.

ClearAll["Global`*"];
deq1 = -(y1[t] + z1[t]);
deq2 = x1[t] + 0.1 y1[t];
deq3 = 0.2 + x1[t] z1[t] - 5.7 z1[t];

deq4 = -(y2[t] + z2[t]);
deq5 = x2[t] + 0.1 y2[t];
deq6 = 0.2 + x2[t] z2[t] - 5.7 z2[t];

x10 = 1; y10 = 1; z10 = 1;
dx0 = 10^-8;
x20 = x10 + dx0; y20 = y10; z20 = z10;
tin = 0; tfin = 10000;
tstep = 1;
acc = 12;

lcedata = {};
sum = 0;

d0 = Sqrt[(x10 - x20)^2 + (y10 - y20)^2 + (z10 - z20)^2 ];

For[i = 1, i < tfin/tstep, i++,

sdeq = {x1'[t] == deq1, y1'[t] == deq2, z1'[t] == deq3, 
x2'[t] == deq4, y2'[t] == deq5, z2'[t] == deq6, x1[0] == x10, 
y1[0] == y10, z1[0] == z10, x2[0] == x20, y2[0] == y20, 
z2[0] == z20};
sol = NDSolve[
sdeq, {x1[t], y1[t], z1[t], x2[t], y2[t], z2[t]}, {t, 0, tstep}, 
MaxSteps -> Infinity, Method -> "Adams", PrecisionGoal -> acc, 
AccuracyGoal -> acc];

xx1[t_]  = x1[t] /. sol[[1]];
yy1[t_]  = y1[t] /. sol[[1]];
zz1[t_] = z1[t] /. sol[[1]];

xx2[t_]  = x2[t] /. sol[[1]];
yy2[t_]  = y2[t] /. sol[[1]];
zz2[t_] = z2[t] /. sol[[1]];

d1 = Sqrt[(xx1[tstep] - xx2[tstep])^2 + (yy1[tstep] - yy2[tstep])^2 + 
          (zz1[tstep] - zz2[tstep])^2 ];

sum += Log[d1/d0];
dlce = sum/(tstep*i);
AppendTo[lcedata, {tstep*i, Log10[dlce]}];

w1 = (xx1[tstep] - xx2[tstep])*(d0/d1); 
w2 = (yy1[tstep] - yy2[tstep])*(d0/d1);
w3 = (zz1[tstep] - zz2[tstep])*(d0/d1); 

x10 = xx1[tstep];
y10 = yy1[tstep];
z10 = zz1[tstep];

x20 = x10 + w1;
y20 = y10 + w2;
z20 = z10 + w3;

i = i++;

If[Mod[tstep*i, 100] == 0, 
Print[" For t = ", tstep*i, " , ", " LCE = ", dlce]]

]

S0 = ListPlot[{lcedata}, Frame -> True, Axes -> False, 
PlotRange -> All, Joined -> True, 
          FrameLabel -> {"t", "log10(LCE)"}, 
FrameStyle -> Directive["Helvetica", 17], ImageSize -> 550]