Mathematica Understanding this Computer program

  • Thread starter Thread starter Slightly
  • Start date Start date
  • Tags Tags
    Computer Program
AI Thread Summary
The discussion centers on a Mathematica program designed to calculate the Lyapunov exponent for a system of differential equations. The code includes two sets of differential equations, each representing a different dynamic system. The purpose of the code highlighted in red involves initializing the variables and setting up the initial conditions for the integration process. It establishes the starting points for the variables and introduces a small perturbation to one of the initial conditions, which is crucial for calculating the Lyapunov exponent. The program then numerically integrates the equations over a specified time interval, updating the state variables and calculating the distance between trajectories to derive the Lyapunov exponent. The output includes a plot of the logarithm of the Lyapunov exponent over time, illustrating the system's stability characteristics. The discussion also touches on a simpler example of numerical integration, demonstrating the fundamental concept of updating variables iteratively.
Slightly
Messages
29
Reaction score
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
 
Physics news on Phys.org
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:
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
 
I don't think that answers my question.
 

Similar threads

Back
Top