Evaluating the Poincaré section for Hénon-Heiles potential through Hénon Method

I must find a poincaré section of a Hénon-Heiles system as described in Hénon-Heiles 1964 paper.

The Hénon-Heiles Hamiltonian is the following, $$H = \frac{1}{2}(p_{1}^{2}+p_{2}^{2}+q_{1}^{2}+q_{2}^{2})+q_{1}^{2} q_{2} - \frac{1}{3}q_{2}^{3} \,\, ,$$ from which one can find the motion equations, $$\dot{q_{1}} = p_{1} \,\, ,$$ $$\dot{p}_{1} = -(q_{1}+2q_{1}q_{2}) \,\, ,$$ $$\dot{q}_{2} = p_{2} \,\, ,$$ $$\dot{p}_{2} = -(q_{2} + q_{1}^{2}-q_{2}^{2}) \,\, ,$$

I want to build a Poincaré section in the plane defined by $$q_{1}=0$$ and $$p_{1} = 0$$, this gives me the system, $$\dot{q}_{2} = p_{2} \,\, ,$$ $$\dot{p}_{2} = -(q_{2}-q_{2}^{2}) \,\, ,$$ Setting $$q_{2}$$ as independent variable (to get the values of crossing) we find the system, $$\frac{dt}{dq_{2}} = \frac{1}{p_{2}} \,\, ,$$ $$\frac{d{p}_{2}}{dq_{2}} = -\frac{(q_{2}-q_{2}^{2})}{p_{2}} \,\, .$$

Our symplectic numerical integrator assume the equations, $$q^{(n+1)} = q^{(n)}+\Delta t\ \dot{q}|_{q^{(n)},p^{(n)}}$$ $$p^{(n+1)} = p^{(n)}+\Delta t\ \dot{p}|_{q^{(n+1)},p^{(n)}}$$

A sketch of the algorithm can be found here in P. Palaniyandi paper and it is summarized in three steps,

Step 1. Integrate Eq. (1) with a fixed integration step size $$h$$.

Step 2. Stop the integration immediately after the trajectory crosses the Poincaré surface ($$\Sigma$$). The $$x_N$$ component of the distance between the Poincaré surface and the first integration point (after crossover) is noted as $$\Delta x_N$$.

Step 3. Compute the next integration point by integrating Eq. (4) with a fixed integration step size $$-\Delta x_N$$.

Step 4. Return to Step 1.

My first attempt to the problem is the following code,

Clear[states, times]

tfinal = 10;
t0 = 0;
q20 = .5;
p20 = .1;
Δt = 10^-3;
s = .6;

times = {0};
states = {{q20, p20}};
cruz = {};

While[t0 < tfinal, t1 = t0 + Δt;
q2 = q20 + Δt p20;
p2 = p20 - Δt (q2 q2 - q2);
If[(q2 - s) (q20 - s) < 0, Break[];
]; t0 = t1;
q20 = q2; p20 = p2; AppendTo[times, t0];
AppendTo[states, {q20, p20}]]


I know that when the "If" condition is satisfied I have to integrate once backwards to get the crossing point and then keep the forward integrations again, but I cannot how to implement this figure this. My main goal is to get the figure 4 of the paper. PS: The Poicaré section must be evaluated by the Hénon algorithm, and the final plot that I need is the following

• It is unclear to me how your code corresponds to the Hénon-Heiles paper. Is it (12) that you are attempting to solve? I might suggest stating what the differential equations are in the question. Sep 1 '19 at 19:09
• Hey, I've edited my question in order to add more details. Sep 1 '19 at 19:41
• What exactly is $s$? Sep 1 '19 at 21:45
• $s$ is the value which defines the crossing with Poincaré section through the expresion $(x_{N} - s)(x_{N+1}-s) \leq 0$. Sep 1 '19 at 21:52
• So you want to have a 1-dimensional Poincaré section (considering that you require $p_1=0,q_1=0,p_2=s$)? Because in the current example, this would yield exactly two points... Sep 1 '19 at 21:58

2 Answers

With the following code, you can plot the Poincaré sections of the Hénon-Heiles system:

With[{icv = {0, 0.36169437164930385, 0.20100851639176504, 0.029106357137938632}},
psection = Reap[NDSolve[{x'[t] == px[t], px'[t] == -(x[t] + 2 x[t]* y[t]),
y'[t] == py[t], py'[t] == -(y[t] + x[t]^2 - y[t]^2),
x[0] == icv[[1]], px[0] == icv[[2]], y[0] == icv[[3]],
py[0] == icv[[4]]}, {x, px, y, py}, {t, 0, 1000},
MaxSteps -> ∞,
Method -> {"EventLocator", "Event" -> x[t],
"EventAction" :> Sow[{y[t], py[t]}]}]][[2]];]


We load the MaTeX package for labels with LaTeX:

Needs["MaTeX"];


Finally, we plot de Poincaré section:

ListPlot[section, PlotRange -> All, AspectRatio -> 1, PlotStyle -> Black,
Frame -> True, FrameStyle -> Black, Axes -> False, LabelStyle -> Directive[Black, Small],
FrameLabel -> {{MaTeX["p_{y}", Magnification -> 1.5], None},
{MaTeX["y", Magnification -> 1.5],
MaTeX["\\text{Hénon-Heiles system}", Magnification -> 1.3]}},
RotateLabel -> False, Epilog -> Inset[MaTeX["E=0.08333", Magnification -> 1],
{0.21, 0.15}, Automatic, 1], ImageSize -> Medium]


The Poincaré section:

You can use a For loop to compute with more initial conditions subject to the energy constraint of the Hénon-Heiles system. For more details, see my answer Poincaré Sections for spring pendulum.

With many initial conditions:

With[{icv = {{0, 0.4082401254164024, 0, 0},
{0, 0.24008038662195985, -0.15316368171208566, 0.28838672839135326},
{0, 0.24352018790129148, 0.09848498536170303, 0.3135210500944904},
{0, 0.20788819727584018, 0.18687849877360435, 0.3047456328865395},
{0, 0.11124773284074371, 0.24234454959281115, 0.32410152798705066},
{0, 0.09612790055948518, 0.2562151235263743, 0.32091473672760257},
{0, 0.14504873906640348, 0.2433184249171519, 0.3098719090788399},
{0, 0.21466692008332666, 0.21341797098009158, 0.2855018063099425},
{0, 0.2242200636897092, 0.21057030870037158, 0.27976766448196655},
{0, 0.3683386516212832, 0.18778901046261137, 0.011696534024095057},
{0, 0.3451153373835715, 0.23762310035771178, 0.005962392196119196},
{0, 0.33504538978074816, 0.256132905175892, 0.0016617858251372995},
{0, 0.37097564355830726, -0.03575247849541091, 0.1665183633794433},
{0, 0.37507484064199426, -0.06280527015275122, 0.14788240243852177},
{0, 0.3774580162329546, -0.08416273725065143, 0.12924644149760023},
{0, 0.3779249463635693, -0.09982487978911159, 0.11491108692766056},
{0, 0.3818158505457779, -0.10979169776813169, 0.08910744870176919},
{0, 0.3850735819245613, -0.11406319118771174, 0.06617088138986574},
{0, 0.3571135001050445, -0.02151416709681081, 0.1966226079763166},
{0, 0.3537329768656702, -0.0015805311387706023, 0.20379028526128642}}},
psection = Reap[
Table[NDSolve[{x'[t] == px[t], px'[t] == -(x[t] + 2 x[t]* y[t]),
y'[t] == py[t], py'[t] == -(y[t] + x[t]^2 - y[t]^2),
x[0] == Part[icv[[i]], 1], px[0] == Part[icv[[i]], 2],
y[0] == Part[icv[[i]], 3], py[0] == Part[icv[[i]], 4]}, {x,
px, y, py}, {t, 0, 5000}, MaxSteps -> \[Infinity],
Method -> {"EventLocator", "Event" -> x[t],
"EventAction" :> Sow[{y[t], py[t]}]}], {i, 1, Length[icv]}]][[2]];]


Trajectories:

With[{icv = {{0, 0.4082401254164024, 0, 0},
{0, 0.24008038662195985, -0.15316368171208566, 0.28838672839135326},
{0, 0.24352018790129148, 0.09848498536170303, 0.3135210500944904},
{0, 0.20788819727584018, 0.18687849877360435, 0.3047456328865395},
{0, 0.11124773284074371, 0.24234454959281115, 0.32410152798705066},
{0, 0.09612790055948518, 0.2562151235263743, 0.32091473672760257},
{0, 0.14504873906640348, 0.2433184249171519, 0.3098719090788399},
{0, 0.21466692008332666, 0.21341797098009158, 0.2855018063099425},
{0, 0.2242200636897092, 0.21057030870037158, 0.27976766448196655},
{0, 0.3683386516212832, 0.18778901046261137, 0.011696534024095057},
{0, 0.3451153373835715, 0.23762310035771178, 0.005962392196119196},
{0, 0.33504538978074816, 0.256132905175892, 0.0016617858251372995},
{0, 0.37097564355830726, -0.03575247849541091, 0.1665183633794433},
{0, 0.37507484064199426, -0.06280527015275122, 0.14788240243852177},
{0, 0.3774580162329546, -0.08416273725065143, 0.12924644149760023},
{0, 0.3779249463635693, -0.09982487978911159, 0.11491108692766056},
{0, 0.3818158505457779, -0.10979169776813169, 0.08910744870176919},
{0, 0.3850735819245613, -0.11406319118771174, 0.06617088138986574},
{0, 0.3571135001050445, -0.02151416709681081, 0.1966226079763166},
{0, 0.3537329768656702, -0.0015805311387706023, 0.20379028526128642}}},
trajectories = Reap[
Table[NDSolve[{x'[t] == px[t], px'[t] == -(x[t] + 2 x[t]* y[t]),
y'[t] == py[t], py'[t] == -(y[t] + x[t]^2 - y[t]^2),
x[0] == Part[icv[[i]], 1], px[0] == Part[icv[[i]], 2],
y[0] == Part[icv[[i]], 3], py[0] == Part[icv[[i]], 4]}, {x,
px, y, py}, {t, 0, 5000}, MaxSteps -> \[Infinity],
Method -> {"EventLocator", "Event" -> x[t],
"EventAction" :> Sow[{y[t], py[t]}]}], {i, 1, Length[icv]}]][[1]];]
(*A trajectory*)
tr = ParametricPlot[Evaluate[{x[t], y[t]} /. trajectory[[2]]], {t, 0, 500}, PlotPoints -> 500, PlotStyle -> {Blue,Thickness[0.002]}, AxesStyle -> Black, LabelStyle -> Directive[Black, Small]];


A trajectory:

Chaotic regime:

Interactive way:

• It would be nice if you share this last interactive code with us, can you add this to your answer? :) Jun 30 at 15:47
• The interactive code is the result of a collaboration with a couple of friends, and I cannot share it yet because it is included in an article that is in review. However, I will gladly share it as soon as the article is published. Greetings! Jun 30 at 20:38

Unless I'm missing something, you can simply integrate your equations using NDSolve, and plot the cut using ParametricPlot:

tfinal = 10;
t0 = 0;
q20 = .5;
p20 = .1;

sol = NDSolve[
{
q2'[t] == p2[t],
p2'[t] == -(q2[t] - q2[t]^2),
q2[t0] == q20,
p2[t0] == p20
},
{p2[t], q2[t]},
{t, t0, tfinal}
]

ParametricPlot[{p2[t], q2[t]} /. sol // Evaluate, {t, t0, tfinal}]


Of course, your example is rather trivial since the solution never leaves the $$q_1=0,p_1=0$$ plane. For non-trivial cases, I would try to use WhenEvent to get the crossings. If you have a non-trivial example, please add it to the question.

• As I said in the OP, I need to evaluate this Poincaré section by the Hénon method. A non-trivial case would be with $p_{1} = 1$. Sep 1 '19 at 21:37
• Can you add the equations for the $p_1=1$ case to the question? I don't see how you can apply the same method without some serious changes. Sep 1 '19 at 21:42
• The point is, I must implement the Hénon algorithm as described in the references, I cannot use such routines as NDSolve... Sep 1 '19 at 21:45
• Why do you need to use that method if there are easier ones? Also, please mention that requirement in the question. It was not clear to me from reading the question that you need to use a specific method. Sep 1 '19 at 21:53
• That is an exercise which I'm having serious trouble to finish. I'm going to edit my post to make it clear. Sep 1 '19 at 21:55