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
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

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



Similar Discussions: (Mathematica) Poincarè map of restricted three body problem
  1. Mathematica problem (Replies: 2)

Loading...