# Poincare Section of an Hamiltonian

I'm in desire to plot the Poincare Section of a differential equation defined by a hamiltonian system.

The hamiltonian is as follow:

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 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.

A expected graph can be seen at those 2 links: one, two.

EDIT: Following this help to plot Poincaré section for double pendulum question, I tried basically 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 == x0, px == px0, y == y0,
py == 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[#][]]] & /@ 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[#][]]] & /@ 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?

• Have you seen for example this – Dr. belisarius Oct 9 '14 at 14:14
• i haven't found it before. thanks for the tip, i will try this out – Geo Oct 9 '14 at 15:47