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.

  • 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. Sep 1 '19 at 19:09
  • 1
    $\begingroup$ Hey, I've edited my question in order to add more details. $\endgroup$ – Herr Schrödinger Sep 1 '19 at 19:41
  • $\begingroup$ What exactly is $s$? $\endgroup$ – Lukas Lang Sep 1 '19 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 Sep 1 '19 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 Sep 1 '19 at 21:58

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 | |
  • $\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 Sep 1 '19 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 Sep 1 '19 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 Sep 1 '19 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 Sep 1 '19 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 Sep 1 '19 at 21:55

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.