4
$\begingroup$

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[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$
2
  • $\begingroup$ Have you seen for example this $\endgroup$ – Dr. belisarius Oct 9 '14 at 14:14
  • $\begingroup$ i haven't found it before. thanks for the tip, i will try this out $\endgroup$ – Geo Oct 9 '14 at 15:47
7
$\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]}]}]];
Section2 = Delete[Section, 1];
Section3 = Section2[[1]];
Section4 = Section3[[1]];
list[[i]] = ListPlot[Transpose[Table[Section3, {j, Length[Section3]}]][[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

$\endgroup$
1
  • $\begingroup$ Very nice indeed! $\endgroup$ – chris Mar 12 at 9:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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