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}]