- #1
Slightly
- 29
- 0
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
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