4
$\begingroup$

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

Poincaré Section of the system.

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ 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. $\endgroup$
    – C. E. ♦
    Commented Sep 1, 2019 at 19:09
  • 1
    $\begingroup$ Hey, I've edited my question in order to add more details. $\endgroup$
    – Herr Schrödinger
    Commented Sep 1, 2019 at 19:41
  • $\begingroup$ What exactly is $s$? $\endgroup$
    – Lukas Lang
    Commented Sep 1, 2019 at 21:45
  • $\begingroup$ $s$ is the value which defines the crossing with Poincaré section through the expresion $(x_{N} - s)(x_{N+1}-s) \leq 0$. $\endgroup$
    – Herr Schrödinger
    Commented Sep 1, 2019 at 21:52
  • $\begingroup$ 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... $\endgroup$
    – Lukas Lang
    Commented Sep 1, 2019 at 21:58

2 Answers 2

Reset to default
6
$\begingroup$

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:

enter image description here

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]];]

enter image description here

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:

enter image description here

Chaotic regime:

enter image description here

Interactive way:

hen1

hen2

Improve this answer
$\endgroup$
2
  • $\begingroup$ It would be nice if you share this last interactive code with us, can you add this to your answer? :) $\endgroup$
    – Herr Schrödinger
    Commented Jun 30, 2021 at 15:47
  • 2
    $\begingroup$ 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! $\endgroup$
    – E. Chan-López
    Commented Jun 30, 2021 at 20:38
1
$\begingroup$

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.

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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$. $\endgroup$
    – Herr Schrödinger
    Commented Sep 1, 2019 at 21:37
  • $\begingroup$ 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. $\endgroup$
    – Lukas Lang
    Commented Sep 1, 2019 at 21:42
  • $\begingroup$ The point is, I must implement the Hénon algorithm as described in the references, I cannot use such routines as NDSolve... $\endgroup$
    – Herr Schrödinger
    Commented Sep 1, 2019 at 21:45
  • 1
    $\begingroup$ 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. $\endgroup$
    – Lukas Lang
    Commented Sep 1, 2019 at 21:53
  • $\begingroup$ That is an exercise which I'm having serious trouble to finish. I'm going to edit my post to make it clear. $\endgroup$
    – Herr Schrödinger
    Commented Sep 1, 2019 at 21:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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