- #1

federes

- 1

- 0

sorry for the inconvenience.

I try to plot the poincarè map of the restricted three body problem. I find in this forum the follow script that do this for the Lorenz system:

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

And this for the Hénon-Heiles (from mathematica " tutorial/NDSolveEventLocator " ):

Needs["DifferentialEquations`NDSolveProblems`"];

system = GetNDSolveProblem["HenonHeiles"];

vars = system["DependentVariables"]; eqns = {system["System"],

system["InitialConditions"]}

Off[NDSolve::noout];

sdata =

Reap[

NDSolve[eqns, {}, {T, 3000},

Method -> {"EventLocator", "Event" -> Subscript[Y, 1][T],

"EventAction" :> Sow[{Subscript[Y, 2][T], Subscript[Y, 4][T]}],

"EventLocationMethod" -> "LinearInterpolation",

"Method" -> {"SymplecticPartitionedRungeKutta",

"DifferenceOrder" -> 4,

"PositionVariables" -> {Subscript[Y, 1][T],

Subscript[Y, 2][T]}}},

StartingStepSize -> 0.25, MaxSteps -> \[Infinity]];

];

psdata = sdata[[-1, 1]];

ListPlot[psdata, Axes -> False, Frame -> True, AspectRatio -> 1]

Now the question is:

Could everyone help me to find a method to show the Poincarè surface for the restricted three body problem?

Thanks a lot to everybody !