0 down vote favorite
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
which gives us the equations:
x'[t]==px[t]
y'[t]==py[t]
px[t]==x[t]/(2*Sqrt[x^2+(1-y)^2])
py[t]==-0.2-(1-y[t])/Sqrt[x^2+(1-y)^2]
or
x''[t]=x[t]/(2*Sqrt[x^2+(1-y)^2])
y''[t]=-0.2-(1-y[t])/Sqrt[x^2+(1-y)^2]
x'[t]=px[t]
y'[t]=py[t]
With this system of equations, I want to plot the Poincare Section/Map for y[t]=0
and py[t]>0
.
Following this question, I tried to adapt to my problem, with no success after a while. Below you can find 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]
which gave me the following erros:
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 fix this problem and be successfully able to plot the desired Poincare Map?
A solution can be seen on thise two articles: one, two
Attachments: