Mathematica for Poincaré sections - or maybe a different tool?

[keenPenguin]

Hi,

I have a system of ODE's for which I want to compute points for a Poincaré section. I used NDSolve and EventHandler in Mathematica but EventHandler appears rather messy/limited for my condition: Suppose I want to save 4 iterations of a Poincaré map and then abort integration imediately. It's an ugly task in Mathematica because one EventHandler condition is, say, psi mod 2pi = 0 while the other is "abort integration when the first condition has been met 4 times". The latter condition is hard to do using EventHandler because it refers to the first condition (so I have to introduce a counter in the first condition...)

Do you know any solution?

Maybe it's a good idea to leave Mathematica at this point and write the whole thing in a language capable of handling loops ;-) I would use Fortran probably.

EDIT: Sorry, by EventHandler I always mean EventLocator

 

[jackmell]

You good with Mathematica programming? Here's code to generate the Lorenz attractor and then construct the Poincare map representing where the flow passes through the manifold x[t]==0 shown as the light-red plane in the plot below. It's kinda messy. You may wish to cut and paste the code into your Mathematica and then run it to see what it looks like. In the first plot below, the Poincare section is superimposed onto the attractor. The second plot is the actual Poincare section but you should run it in Mathematica so you can rotate it interactively (ver 6 or better).

mysol = NDSolve[{x'[t] == -3 (x[t] - y[t]), 
   y'[t] == -x[t] z[t] + 26.5 x[t] - y[t], z'[t] == x[t] y[t] - z[t], 
   x[0] == z[0] == 0, y[0] == 1}, {x, y, z}, {t, 0, 200}, 
  MaxSteps -> 100000]

myx[t_] := Evaluate[x[t] /. mysol];

myplot = Plot[myx[t], {t, 0, 200}]

pp1 = ParametricPlot3D[
  Evaluate[{x[t], y[t], z[t]} /. mysol], {t, 0, 200}, 
  AxesLabel -> {"x", "y", "z"}, PlotPoints -> 150]
cp1 = ContourPlot3D[x == 0, {x, -20, 20}, {y, -20, 20}, {z, 0, 50}, 
  Mesh -> None, ContourStyle -> {LightPurple, Opacity[0.4]}]

(* get line segments from plot *)
pts = Cases[myplot, Line[{t__}] -> t, \[Infinity]];

(* from set above, select all pairs which exhibit a sign change to \
identify a root *)
myselection = 
  Select[Split[pts, Sign[Last[#2]] == -Sign[Last[#1]] &], 
   Length[#1] == 2 &];

(* get first element of the array *)
mymap = Map[First, myselection, {2}];

(* find in line segments where each changes sign *)
mymap2 = Map[FindRoot[myx[t] == 0, {t, #[[1]], #[[2]]}] &, mymap];
mypoints = t /. mymap2;

(* plot Poincare map  where x[t]=0 *)
myplots = 
  Point@ {First[Evaluate[x[#] /. mysol]], 
      First[Evaluate[y[#] /. mysol]], 
      First[Evaluate[z[#] /. mysol]]} & /@ mypoints;
sg1 = Show[Graphics3D[{Red, PointSize[0.008], myplots}]]

(* show everything *)
Show[{pp1, cp1, sg1}]

 

[Simon_Tyler]

Here's a couple of Demonstrations that generate Poincare Sections:
http://demonstrations.wolfram.com/DuffingOscillator/
http://demonstrations.wolfram.com/KAMToriReforming/
You can view and download the code for the demonstrations.

 

[kwabbles]

The below code works great. However, I can't seem to change it, to look at the intersection in the xy plane. I changed myx, myplot and cp1 (to look where z=27) and those work there is something down the line that i am not changing correctly. Any pointers

You good with Mathematica programming? Here's code to generate the Lorenz attractor and then construct the Poincare map representing where the flow passes through the manifold x[t]==0 shown as the light-red plane in the plot below. It's kinda messy. You may wish to cut and paste the code into your Mathematica and then run it to see what it looks like. In the first plot below, the Poincare section is superimposed onto the attractor. The second plot is the actual Poincare section but you should run it in Mathematica so you can rotate it interactively (ver 6 or better).

mysol = NDSolve[{x'[t] == -3 (x[t] - y[t]), 
   y'[t] == -x[t] z[t] + 26.5 x[t] - y[t], z'[t] == x[t] y[t] - z[t], 
   x[0] == z[0] == 0, y[0] == 1}, {x, y, z}, {t, 0, 200}, 
  MaxSteps -> 100000]

myx[t_] := Evaluate[x[t] /. mysol];

myplot = Plot[myx[t], {t, 0, 200}]

pp1 = ParametricPlot3D[
  Evaluate[{x[t], y[t], z[t]} /. mysol], {t, 0, 200}, 
  AxesLabel -> {"x", "y", "z"}, PlotPoints -> 150]
cp1 = ContourPlot3D[x == 0, {x, -20, 20}, {y, -20, 20}, {z, 0, 50}, 
  Mesh -> None, ContourStyle -> {LightPurple, Opacity[0.4]}]

(* get line segments from plot *)
pts = Cases[myplot, Line[{t__}] -> t, \[Infinity]];

(* from set above, select all pairs which exhibit a sign change to \
identify a root *)
myselection = 
  Select[Split[pts, Sign[Last[#2]] == -Sign[Last[#1]] &], 
   Length[#1] == 2 &];

(* get first element of the array *)
mymap = Map[First, myselection, {2}];

(* find in line segments where each changes sign *)
mymap2 = Map[FindRoot[myx[t] == 0, {t, #[[1]], #[[2]]}] &, mymap];
mypoints = t /. mymap2;

(* plot Poincare map  where x[t]=0 *)
myplots = 
  Point@ {First[Evaluate[x[#] /. mysol]], 
      First[Evaluate[y[#] /. mysol]], 
      First[Evaluate[z[#] /. mysol]]} & /@ mypoints;
sg1 = Show[Graphics3D[{Red, PointSize[0.008], myplots}]]

(* show everything *)
Show[{pp1, cp1, sg1}]

 

[blue_raver22]

.I understand the frustration of encountering limitations in tools that we use for our research. In this case, it seems like Mathematica's EventHandler function may not be the best option for your specific task of computing points for a Poincaré section. While Mathematica is a powerful tool for many scientific applications, it may not always be the best choice for every task.

I would suggest exploring other software options that may better suit your needs. For example, there are specific programs designed for analyzing dynamical systems and generating Poincaré sections, such as XPPAUT or AUTO. These programs are specifically designed for handling the types of calculations and conditions you are describing, and may be a better fit for your task than trying to work around the limitations of Mathematica.

Alternatively, if you are comfortable with Fortran, you could certainly write your own code to handle the task. This may require more effort and time, but it would allow you to have complete control over the calculations and conditions without any limitations imposed by a specific software program.

In conclusion, while Mathematica may be a useful tool for many scientific applications, it may not always be the best choice for every task. It is important to explore different options and choose the best tool for the specific task at hand. I hope this helps and wish you success in your research.