Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

}
share|improve this question
1  
user5267, please read this: mathematica.stackexchange.com/editing-help –  Mr.Wizard Jan 11 '13 at 6:45
    
@Mr.Wizard, Thank you for your help. But I really do not know what to do with this link. Can you please explain more about it? –  user5267 Jan 11 '13 at 9:52
    
It is a guide to properly formatting your posts here on StackExchange. Well formatted questions consistently get more attention and better answers, as well as reduce the burden on those who try to clean up poorly formatting posts. There is also a "tool bar" above each edit box with a series of icons that can be used for formatting (for example, select your code and click the { } icon). –  Mr.Wizard Jan 11 '13 at 10:07
    
Oh. Thanks. It would be helpful. –  user5267 Jan 11 '13 at 10:41
add comment

3 Answers

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

Mathematica graphics

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.

Mathematica graphics

share|improve this answer
    
It was awesome. Thanks man :) –  user5267 Jan 12 '13 at 4:35
add comment

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]
share|improve this answer
    
Thank you for your answer.By the way, I think that these two methods almost give the same results. Aren't they?! I appreciate it, if you share your idea about it. –  user5267 Jan 11 '13 at 9:44
add comment

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]
share|improve this answer
    
thank you so much. It was really great and helpful. –  user5267 Jan 12 '13 at 4:35
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.