Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

(Mathematica) Poincarè map of restricted three body problem

  1. Mar 18, 2012 #1
    Hi everybody,

    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 queston 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 !!
     
  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted