Understanding this Computer program

  • Context: Mathematica 
  • Thread starter Thread starter Slightly
  • Start date Start date
  • Tags Tags
    Computer Program
Click For Summary
SUMMARY

This discussion focuses on a Mathematica program designed to calculate the Lyapunov exponent of a differential equation. The code utilizes numerical integration through the NDSolve function, employing the Adams method with specified precision goals. Key variables include initial conditions for the differential equations and a loop to compute the Lyapunov exponent over a defined time interval. The output is visualized using ListPlot to display the logarithm of the Lyapunov exponent over time.

PREREQUISITES
  • Familiarity with Mathematica programming language
  • Understanding of differential equations and their numerical solutions
  • Knowledge of Lyapunov exponents and their significance in dynamical systems
  • Experience with data visualization techniques in Mathematica
NEXT STEPS
  • Explore Mathematica's NDSolve function for advanced numerical integration techniques
  • Study the Adams method for solving ordinary differential equations
  • Learn about Lyapunov stability theory and its applications in chaos theory
  • Investigate data visualization options in Mathematica, focusing on ListPlot and its customization
USEFUL FOR

This discussion is beneficial for mathematicians, physicists, and engineers interested in dynamical systems, particularly those working with numerical methods in Mathematica to analyze stability and chaos.

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

  • · Replies 10 ·
Replies
10
Views
6K
  • · Replies 1 ·
Replies
1
Views
2K
  • · Replies 1 ·
Replies
1
Views
4K
  • · Replies 14 ·
Replies
14
Views
4K
  • · Replies 1 ·
Replies
1
Views
4K
  • · Replies 11 ·
Replies
11
Views
7K
Replies
1
Views
3K
  • · Replies 14 ·
Replies
14
Views
4K
  • · Replies 15 ·
Replies
15
Views
3K
  • · Replies 6 ·
Replies
6
Views
3K