19
$\begingroup$

I'm in the desire to plot the Poincare Section of a differential equation defined by a Hamiltonian system. The Hamiltonian is as follows:

H[x_,a_,y_,b_] = (a^2+b^2)/2+f*y+(1-f-Sqrt[x^2+(1-y)^2])/2

Here x and y and generalized coordinates and a, b their respective momentum. In case you are curious, this is the Hamiltonian of a spring pendulum.

In my case, to evaluate numeric, I set f = 0.2 and H = 0.03. And after applying Hamilton equations, I have:

x''[t] = x[t]/(2*Sqrt[x^2 + (1 - y)^2])
y''[t] = -0.2 - (1 - y[t])/Sqrt[x^2 + (1 - y)^2]

while x'[t] = a[t] and y'[t] = b[t].

Now, I want to plot the Poincare Map/Section for y[t] = 0 and b[t] > 0. How am I supposed to do it? I have tried using this and this with no success. An expected graph can be seen at these 2 links: one, two.


EDIT: Following this help to plot Poincaré section for double pendulum question, I basically tried changing this to my problem. Here is my code:

equations = {x'[t] == px[t],
    y'[t] == py[t],
    px'[t] == x[t]/(2*Sqrt[x[t]^2 + (1 - y[t])^2]),
    py'[t] == -0.2 - (1 - y[t])/Sqrt[x[t]^2 + (1 - y[t])^2]};

psect[{x0_, px0_, y0_, py0_}] := 
  Reap[NDSolve[{equations, x[0] == x0, px[0] == px0, y[0] == y0, 
  py[0] == py0, WhenEvent[y[t] + 0*x[t] + 0*px[t] + 0*py[t] == 0,
  If[py[t] > 0,
    Sow[{px[t], x[t]}]]]}, {x, px, y, py}, {t, 0, 200}, 
       MaxSteps -> \[Infinity], MaxStepSize -> .0005]][[2, 1]];

ps = psect /@ Table[{0.01, 0.005, 0.01, i}, {i, -0.01, .01, 0.00125}]

ListLinePlot[#[[FindCurvePath[#][[1]]]] & /@ ps,
  Mesh -> All, PlotRange -> All]

I got some errors at the ps and plot part:

ps = psect /@ Table[{0.01, 0.005, 0.01, i}, {i, -0.01, .01, 0.00125}]
Part::partw: Part 1 of {} does not exist. >>
Part::partw: Part 1 of {} does not exist. >>
Part::partw: Part 1 of {} does not exist. >>
General::stop: Further output of Part::partw will be suppressed during this calculation. >>
ListLinePlot[#[[FindCurvePath[#][[1]]]] & /@ ps, Mesh -> All, 
  PlotRange -> All]
Part::pkspec1: The expression <<1>> cannot be used as a part specification. >>
Part::pkspec1: The expression <<1>> cannot be used as a part specification. >>
Part::pkspec1: The expression <<1>> cannot be used as a part specification. >>
General::stop: Further output of Part::pkspec1 will be suppressed during this calculation. >>

How can I work this around?

$\endgroup$
3
  • 1
    $\begingroup$ Have you seen for example this $\endgroup$ Commented Oct 9, 2014 at 14:14
  • $\begingroup$ i haven't found it before. thanks for the tip, i will try this out $\endgroup$
    – Geo
    Commented Oct 9, 2014 at 15:47
  • $\begingroup$ Cross-posted: community.wolfram.com/groups/-/m/t/366589 $\endgroup$
    – Michael E2
    Commented Jun 20, 2023 at 18:17

1 Answer 1

30
$\begingroup$

Try with the following code:

H[px_, pz_, x_, z_] := 1/2 (px^2 + pz^2) + z + 2 (3/4 - Sqrt[x^2 + (z - 1)^2])^2 - 1/8
xpp = D[x[t], t, t] - x[t] (-4 + 3/Sqrt[x[t]^2 + (-1 + z[t])^2]);
zpp = D[z[t], t, t] - 3 + 3/Sqrt[x[t]^2 + (-1 + z[t])^2] - (-4 + 3/Sqrt[x[t]^2 + (-1 + z[t])^2]) z[t];
n = 140;
list = Range[n];
For[i = 1;, i < n + 1, i++,
    ics = {Random[Real, {-0.55, 0.55}], Random[Real, {-0.45, 0.45}], Random[Real, {-0.44, 0.9}], Random[Real, {-0.17, 0.9}]};
inisol = Solve[H[ini, ics[[4]], ics[[1]], ics[[2]]] == 0.866, ini];
inival = ini /. inisol;
ics = {ics[[1]], ics[[2]], inival[[2]], ics[[4]]};
TimeMin = 0;
TimeMax = 12000;
Section = Reap[NDSolve[{xpp == 0, zpp == 0, x[0] == ics[[1]], 
 z[0] == ics[[2]], x'[0] == ics[[3]], z'[0] == ics[[4]]}, {x, 
 z}, {t, TimeMin, TimeMax}, MaxSteps -> ∞, Method -> {"EventLocator","Event" -> x[t], "EventAction" :> Sow[{z[t], z'[t]}]}]][[2]];
list[[i]] = ListPlot[Transpose[Table[Section, {j, Length[Section]}]][[1]][[1]], PlotStyle -> {Black, Opacity[0.25], AbsolutePointSize[0.5]}]]
Show[list, PlotRange -> All, Frame -> True, Axes -> False, LabelStyle -> Directive[Black, Small], ImageSize -> Large, Background -> Lighter[Gray, 0.95]]

The resulting Poincaré section is as follows:

enter image description here

Also, if you want to paint each orbit with a different color, you just have to place the following:

PlotStyle -> {RGBColor[Random[], Random[], Random[]], AbsolutePointSize[0.5]}

and

LabelStyle -> Directive[White, Small], Background -> GrayLevel[0.1]

enter image description here

It is even possible to paint Poincaré sections interactively:

enter image description here

$\endgroup$
2
  • 2
    $\begingroup$ Very nice indeed! $\endgroup$
    – chris
    Commented Mar 12, 2021 at 9:17
  • $\begingroup$ How do you do this last interactive bit? Also, is it possible to easily extract the last values of $(x,z,x',z')$ to verify energy conservation? $\endgroup$ Commented Nov 10, 2022 at 19:03

Not the answer you're looking for? Browse other questions tagged or ask your own question.