# Recreating Poincaré section for a double pendulum

I am trying to recreate the Poincaré Sections for the double pendulum as described in this paper. I am unable to figure out a way to include energy values in the code. The code that I used until now is given below:

\$Assumptions = {l1 ∈ Reals, l1 > 0, l2 ∈ Reals,
l2 > 0, m1 ∈ Reals, m1 > 0, m2 ∈ Reals, m2 > 0};

lagrangian =
1/2 (l1^2 (m1 + m2) θ1'[t]^2 + l2^2 m2 θ2'[t]^2 +
2 l1 l2 m2 Cos[θ1[t] - θ2[t]] θ1'[
t] θ2'[t]) + g l1 (m1 + m2) Cos[θ1[t]] +
g l2 m2 Cos[θ2[t]];

<< VariationalMethods

eom = EulerEquations[lagrangian, {θ1[t], θ2[t]}, {t}] // FullSimplify

e1 = eom[[1]]
e2 = eom[[2]]

m2 = 1; l2 = 1; g = 1; m1 = 3; l1 = 2;

momentum1 = D[lagrangian, θ1'[t]];
momentum2 = D[lagrangian, θ2'[t]];

hamiltonian = {momentum1 θ1'[t] + momentum2 θ2'[t]} -
lagrangian // Simplify

sol = Reap@
NDSolve[{e1,
e2, θ1[0] == 0, θ2[0] == 0.15, θ1'[0] ==
0, θ2'[0] == 0,
WhenEvent[θ1[t] == 0 && θ1'[t] > 0,
Sow[{θ2[t], θ2'[
t]}]]}, {θ1, θ2}, {t, 0, 1000}];

ListPlot[sol[[2]], PlotTheme -> "Scientific"]


How do I recreate Figures 2 and 3 from the paper?

The Poincaré section in Fig. 2 is shown for $$E_1 = E_{01} + 0.01 = -1-l_1(1+m_1) + 0.01.$$

To sample different parts of the section, you have to change two initial parameters, for example $$\dot {\theta_1}(0)$$ and $$\dot{\theta_2}(0)$$, but in a way that the total energy is $$E_1$$. You can use Mathematica to calculate the relationship between $$\dot {\theta_1}(0)$$ and $$\dot{\theta_2}(0)$$ if the other two initial conditions are $$\theta_1(0) = \theta_2(0) = 0$$:

Solve[(hamiltonian == -1 - l1 (1 + m1) + 1/100) /. {θ1[t] -> 0, θ2[t] -> 0},
{θ1'[t]}]


$$\dot{\theta_1}(0) = \frac{1}{40} \left(\sqrt{2-75 \dot{\theta_2}(0)^2}-5 \dot{\theta_2}(0)\right)$$

Now you have to rerun NDSolve multiple times with different initial conditions, sampling along $$\dot{\theta_2}(0)$$ (which I denote with ω2):

ω1[ω2_] := 1/40 (-5 ω2 + Sqrt[2 - 75 ω2^2]);

sol = Quiet@Table[Reap@NDSolve[{e1, e2, θ1[0] == 0, θ2[0] ==
0.0, θ1'[0] == ω1[ω2], θ2'[0] == ω2,
WhenEvent[θ1[t] == 0 && θ1'[t] > 0,
Sow[{θ2[t], θ2'[t]}]]}, {θ1, θ2}, {t, 0, 1000}],
{ω2, -.2, .02, .01}];

ListPlot[Catenate@sol[[All, 2]], PlotTheme -> "Scientific",
AspectRatio -> Automatic, PlotStyle -> PointSize[.006]]


To get Fig. 3, you have to change the total energy condition in the equation above, and calculate new relationship between $$\dot{\theta_1}(0)$$ and $$\dot{\theta_2}(0)$$.

Of course you can further automatize this with a general solver for $$\dot{\theta_1}(0)$$:

ω1[ω2_, ℰ_] := Derivative[1][θ1][t] /.
First@Solve[(hamiltonian == ℰ) /. {θ1[t] -> 0, θ2[t] -> 0,
Derivative[1][θ2][t] -> ω2}, {θ1'[t]}]

plotPoincare[ℰ_, {ω2min_, ω2max_, δω2_}] := Module[{sol},
sol = Quiet@Table[Reap@NDSolve[{e1, e2, θ1[0] == 0, θ2[0] ==
0.0, θ1'[0] == ω1[ω2, ℰ], θ2'[0] == ω2,
WhenEvent[θ1[t] == 0 && θ1'[t] > 0,
Sow[{Mod[θ2[t] + Pi, 2 Pi] - Pi, θ2'[t]}]]}, {θ1, θ2}, {t, 0,
1000}], {ω2, ω2min, ω2max, δω2}];
ListPlot[Catenate@sol[[All, 2]], PlotTheme -> "Scientific",
AspectRatio -> Automatic, PlotStyle -> PointSize[.01]]
]

ℰ2 = +1 - l1 (1 + m1) - 1/100;
ℰ3 = -1 + l1 (1 + m1) - 1/100;

plotPoincare[ℰ2, {-2, 2, .1}]
plotPoincare[ℰ3, {-8, 8, .5}]


• Thanks a lot for this illuminating answer! Commented Apr 23, 2023 at 6:28

We can get a 2D Poincare section for a double pendulum by first creating the Lagragian equations for the masses m1 and m2. We call the angles t1 and t2, g is the gravity constant (for an example, we set all the parameters to 1):

v1 = l1 t1'[t]; (*velocity of m1*)
v21 = l2 t2' [t]; (*velocity of m2 relative to m1*)
v2 = Sqrt[v1^2 + v21^2 + 2 v1 v21 Cos[t2[t] - t1[t]]];(*velocity of m2*)
pot = - m1 g Cos[t1[t]] - m2 g Cos[t2[t]]; (*potential*)
lagrange =
1/2  m1 v1^2 + 1/2 m2 v2^2 - pot /. {l1 -> 1, m1 -> 1, l2 -> 1,
m2 -> 1, g -> 1};
eq = {D[D[lagrange, t1'[t]], t] - D[lagrange, t1[t]] == 0,
D[D[lagrange, t2'[t]], t] - D[lagrange, t2[t]] == 0};


We will solve the equations inside a manipulate, that allows to change the initial values easily. For an example we choose the section by requiring x1==0 and we plot t2 against t2'. To get the point of the Poincare section we use a "WhenEvent":

tmax = 1000;
Manipulate[
ini = {t1[0] == t10, t2[0] == t20, t1'[0] == v10, t2'[0] == v20};
pts = Reap[
NDSolve[{eq, ini,
WhenEvent[t1[t] == 0, Sow[{t2[t], t2'[t]}]]}, {t1, t2}, {t, 0,
tmax}];
][[2, 1]];
Graphics[Point[pts], Axes -> True, PlotRange -> {{-2, 2}, {-2, 2}}]
, {t10, 0, Pi/2}, {{t20, Pi/2}, 0, Pi/2}, {v10, 0, 1/2}, {v20, 0,
1/2}]
`

• Thank you for the details! Commented Apr 23, 2023 at 6:29