# Mathematica Understanding this Computer program

Tags:
1. Feb 26, 2015

### Slightly

Below is a computer program for Mathematica that details how to calculate the Lyapunov exponent of a differential equation. What is the purpose of the code colored in red?

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]

The link to the code: http://mathematica.stackexchange.com/questions/17593/lyapunov-exponent/17608#17608

2. Feb 26, 2015

### Staff: Mentor

This is how you numerically integrate a function computationally.

As an example, say you want to integrate dy=2xdx over some interval of 0 to 5

Then you'd do something like this:

Code (Text):

dx = 0.01
x = 0
y = 0
For i from 0 to 500
x = x + dx
dy = 2 * x * dx
y = y + dy
Print x,y
Next

3. Feb 26, 2015

### Slightly

I don't think that answers my question.