3
$\begingroup$

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?

$\endgroup$

2 Answers 2

7
$\begingroup$

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

enter image description here

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

enter image description here enter image description here

$\endgroup$
1
  • $\begingroup$ Thanks a lot for this illuminating answer! $\endgroup$
    – codebpr
    Commented Apr 23, 2023 at 6:28
3
$\begingroup$

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

enter image description here

$\endgroup$
1
  • $\begingroup$ Thank you for the details! $\endgroup$
    – codebpr
    Commented Apr 23, 2023 at 6:29

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.