Skip to main content

Triple Pendulum

Acrobatics of a simple physical model

Original code and equations by Peter Tenisheff (tenishefff@yandex.ru)

This Demonstration simulates the motion of a free moving triple pendulum using basic numerics of Wolfram.

Download original notebook

Equations of Motion​

The equations of motion for a triple pendulum system can be more compactly written using notations Δθij(t)=θi(t)−θj(t)\Delta\theta_{ij}(t) = \theta_i(t) - \theta_j(t):

A1θ1′′(t)+B12cos⁡(Δθ12(t))θ2′′(t)+B13cos⁡(Δθ13(t))θ3′′(t)+C12sin⁡(Δθ12(t))+C13sin⁡(Δθ13(t))+gD1sin⁡(θ1(t))=0,A2θ2′′(t)+B21cos⁡(Δθ12(t))θ1′′(t)+B23cos⁡(Δθ23(t))θ3′′(t)−C21sin⁡(Δθ12(t))+C23sin⁡(Δθ23(t))+gD2sin⁡(θ2(t))=0,A3θ3′′(t)+B31cos⁡(Δθ13(t))θ1′′(t)+B32cos⁡(Δθ23(t))θ2′′(t)−C31sin⁡(Δθ13(t))−C32sin⁡(Δθ23(t))+gD3sin⁡(θ3(t))=0\begin{align*} & A_1 \theta_1''(t) + B_{12} \cos(\Delta\theta_{12}(t)) \theta_2''(t) + B_{13} \cos(\Delta\theta_{13}(t)) \theta_3''(t) \\ & + C_{12} \sin(\Delta\theta_{12}(t)) + C_{13} \sin(\Delta\theta_{13}(t)) + g D_1 \sin(\theta_1(t)) = 0, \\ & A_2 \theta_2''(t) + B_{21} \cos(\Delta\theta_{12}(t)) \theta_1''(t) + B_{23} \cos(\Delta\theta_{23}(t)) \theta_3''(t) \\ & - C_{21} \sin(\Delta\theta_{12}(t)) + C_{23} \sin(\Delta\theta_{23}(t)) + g D_2 \sin(\theta_2(t)) = 0, \\ & A_3 \theta_3''(t) + B_{31} \cos(\Delta\theta_{13}(t)) \theta_1''(t) + B_{32} \cos(\Delta\theta_{23}(t)) \theta_2''(t) \\ & - C_{31} \sin(\Delta\theta_{13}(t)) - C_{32} \sin(\Delta\theta_{23}(t)) + g D_3 \sin(\theta_3(t)) = 0 \end{align*}

where:

  • A1=(m1+m2+m3)l12A_1 = (m_1 + m_2 + m_3) l_1^2
  • A2=(m2+m3)l22A_2 = (m_2 + m_3) l_2^2
  • A3=m3l32A_3 = m_3 l_3^2
  • Bij=mjliljB_{ij} = m_j l_i l_j
  • Cij=mjliljθj′(t)2C_{ij} = m_j l_i l_j \theta_j'(t)^2
  • Di=(m1+m2+m3)liD_i = (m_1 + m_2 + m_3) l_i

and gg is the acceleration due to gravity, m1,m2,m3m_1, m_2, m_3 are the masses, l1,l2,l3l_1, l_2, l_3 are the lengths of the pendulums, and θ1,θ2,θ3\theta_1, \theta_2, \theta_3 are the angles with the vertical.

Here is corresponding system rewritten using WL:

ClearAll[g, m, l, timeMax, \[Theta]1, \[Theta]2, \[Theta]3];

g = 9.81;
(*SbB[*)Subscript[m(*|*),(*|*)1](*]SbB*) = 1; (*SbB[*)Subscript[m(*|*),(*|*)2](*]SbB*) = 1; (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*) = 1;
(*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*) = 1; (*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*) = 1; (*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*) = 1;
timeMax = 10;


eqns = {
  ((*SbB[*)Subscript[m(*|*),(*|*)1](*]SbB*) + (*SbB[*)Subscript[m(*|*),(*|*)2](*]SbB*) + (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*)) (*SpB[*)Power[(*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*)(*|*),(*|*)2](*]SpB*) \[Theta]1''[t] + ((*SbB[*)Subscript[m(*|*),(*|*)2](*]SbB*) + (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*)) (*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*) \[Theta]2''[t] Cos[\[Theta]1[t] - \[Theta]2[t]] +
  (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*) \[Theta]3''[t] Cos[\[Theta]1[t] - \[Theta]3[t]] +
  ((*SbB[*)Subscript[m(*|*),(*|*)2](*]SbB*) + (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*)) (*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*) (*SpB[*)Power[\[Theta]2'[t](*|*),(*|*)2](*]SpB*) Sin[\[Theta]1[t] - \[Theta]2[t]] +
  (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*) (*SpB[*)Power[\[Theta]3'[t](*|*),(*|*)2](*]SpB*) Sin[\[Theta]1[t] - \[Theta]3[t]] +
  g ((*SbB[*)Subscript[m(*|*),(*|*)1](*]SbB*) + (*SbB[*)Subscript[m(*|*),(*|*)2](*]SbB*) + (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*)) (*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*) Sin[\[Theta]1[t]] == 0,
  
  ((*SbB[*)Subscript[m(*|*),(*|*)2](*]SbB*) + (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*)) (*SpB[*)Power[(*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*)(*|*),(*|*)2](*]SpB*) \[Theta]2''[t] + ((*SbB[*)Subscript[m(*|*),(*|*)2](*]SbB*) + (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*)) (*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*) \[Theta]1''[t] Cos[\[Theta]1[t] - \[Theta]2[t]] +
  (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*) \[Theta]3''[t] Cos[\[Theta]2[t] - \[Theta]3[t]] -
  ((*SbB[*)Subscript[m(*|*),(*|*)2](*]SbB*) + (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*)) (*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*) \[Theta]1'[t]^2 Sin[\[Theta]1[t] - \[Theta]2[t]] +
  (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*) (*SpB[*)Power[\[Theta]3'[t](*|*),(*|*)2](*]SpB*) Sin[\[Theta]2[t] - \[Theta]3[t]] +
  g ((*SbB[*)Subscript[m(*|*),(*|*)2](*]SbB*) + (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*)) (*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*) Sin[\[Theta]2[t]] == 0,

  (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*) (*SpB[*)Power[(*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*)(*|*),(*|*)2](*]SpB*) \[Theta]3''[t] + (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*) \[Theta]1''[t] Cos[\[Theta]1[t] - \[Theta]3[t]] +
  (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*) \[Theta]2''[t] Cos[\[Theta]2[t] - \[Theta]3[t]] -
  (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*) (*SpB[*)Power[\[Theta]1'[t](*|*),(*|*)2](*]SpB*) Sin[\[Theta]1[t] - \[Theta]3[t]] -
  (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*) (*SpB[*)Power[\[Theta]2'[t](*|*),(*|*)2](*]SpB*) Sin[\[Theta]2[t] - \[Theta]3[t]] +
  g (*SbB[*)Subscript[m(*|*),(*|*)3](*]SbB*) (*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*) Sin[\[Theta]3[t]] == 0
};

initCond = {
  \[Theta]1[0] == (*FB[*)((\[Pi])(*,*)/(*,*)(2))(*]FB*) - 0.1,
  \[Theta]2[0] == \[Pi],
  \[Theta]3[0] == \[Pi],
  \[Theta]1'[0] == 0,
  \[Theta]2'[0] == 0,
  \[Theta]3'[0] == 0
};

Now let's solve it using numerical methods:

sol = NDSolve[{eqns, initCond}, {\[Theta]1, \[Theta]2, \[Theta]3}, {t, 0, timeMax}] // First

Some shorthand notations:

ClearAll[xC, yC];

With[{
  \[Phi]1 = \[Theta]1 /. sol,
  \[Phi]2 = \[Theta]2 /. sol,
  \[Phi]3 = \[Theta]3 /. sol
},

  xC[t_, 1] := (*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*) Sin[\[Phi]1[t]];
  yC[t_, 1] := -(*SbB[*)Subscript[l(*|*),(*|*)1](*]SbB*) Cos[\[Phi]1[t]];
  xC[t_, 2] := xC[t, 1] + (*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*) Sin[\[Phi]2[t]];
  yC[t_, 2] := yC[t, 1] - (*SbB[*)Subscript[l(*|*),(*|*)2](*]SbB*) Cos[\[Phi]2[t]];
  xC[t_, 3] := xC[t, 2] + (*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*) Sin[\[Phi]3[t]];
  yC[t_, 3] := yC[t, 2] - (*SbB[*)Subscript[l(*|*),(*|*)3](*]SbB*) Cos[\[Phi]3[t]];

];

trajectories = Table[Table[{xC[t, i], yC[t, i]}, {i, 3}], {t, 0, timeMax, 0.05}];
trajectories = trajectories // Transpose;

Animation​

Let's animate the motion of our triple pendulum. Initially, we will use points for the joints of the pendulum connected by a black line. Here, we use a custom UpdateFunction for Animate to avoid costly reevaluation of the entire plot:

Module[{points, tracks}, Animate[
  Graphics[
    {
      {
        Directive[TransitionType->None],
        {Gray, Opacity[0.25], Line[tracks[[3]] // Offload ]},
        {Red, Opacity[0.15], Line[tracks[[2]] // Offload ]},
        {Blue, Opacity[0.15], Line[tracks[[1]] // Offload ]}
      },
      
      Black, 
      Line[{{0, 0}, points[[1]], points[[2]], points[[3]]} // Offload],
      Blue, PointSize[0.03], Point[points[[1]] // Offload],
      Red, PointSize[0.03], Point[points[[2]] // Offload],
      Green, PointSize[0.04], Point[points[[3]] // Offload] 
    },
    
    PlotRange -> {{-3, 3}, {-3, 3}},
    Prolog->{Line[{{-3,0}, {3,0}}]},
    Axes -> True,  Frame->True,
    Background -> White, 
    ImageSizeRaw -> {450,450}
  ], 
  {t, 0, timeMax, 0.05}, 
  
  AnimationRate -> 30,
  
  "UpdateFunction" -> Function[t,

    tracks = trajectories[[All, ;; Floor[t / 0.05] + 1]];
    points = Table[{xC[t, i], yC[t, i]}, {i, 3}];
    
    False
  ]
] ]
(*VB[*)(FrontEndRef["9b7a3593-b78b-4d80-b84f-b4d42c41192c"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKWyaZJxqbWhrrJplbJOmapFgY6CZZmKTpJpmkmBglmxgaWholAwCCLxV4"*)(*]VB*)

By default, Graphics smooths the motion of lines and points, which might lead to some artifacts for continuous lines representing the full trajectory of each joint. Therefore, we explicitly disable it in Directive for them.


Extras​

Here is another version of it with pictures placed on a graph instead of basic points:

Module[{points}, Animate[
  Graphics[
    {

      Brown,
      {AbsoluteThickness[4], Line[{{0, 0}, points[[1]], points[[2]], points[[3]]} // Offload]},
      Inset[ColorNegate[(*VB[*)(FrontEndRef["8b049d7d-d9ab-4a67-a50c-19a72ab7f93c"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKWyQZmFimmKfoplgmJumaJJqZ6yaaGiTrGlommhslJpmnWRonAwCK1hYH"*)(*]VB*)], {0,0}],
      Blue, PointSize[0.03], Inset[(*VB[*)(FrontEndRef["8b049d7d-d9ab-4a67-a50c-19a72ab7f93c"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKWyQZmFimmKfoplgmJumaJJqZ6yaaGiTrGlommhslJpmnWRonAwCK1hYH"*)(*]VB*), points[[1]] // Offload],
      Red, PointSize[0.03], Inset[(*VB[*)(FrontEndRef["8b049d7d-d9ab-4a67-a50c-19a72ab7f93c"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKWyQZmFimmKfoplgmJumaJJqZ6yaaGiTrGlommhslJpmnWRonAwCK1hYH"*)(*]VB*), points[[2]] // Offload],
      Green, PointSize[0.04], Inset[(*VB[*)(FrontEndRef["8b049d7d-d9ab-4a67-a50c-19a72ab7f93c"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKWyQZmFimmKfoplgmJumaJJqZ6yaaGiTrGlommhslJpmnWRonAwCK1hYH"*)(*]VB*), points[[3]] // Offload]  
      
    },
    PlotRange -> {{-3, 3}, {-3, 3}},
    Prolog->{Line[{{-3,0}, {3,0}}]},
    Axes -> True,  Frame->True,
    Background -> White, 
    TransitionType->None,
    ImageSizeRaw -> {450,450}
  ], 
  
  {t, 0, timeMax, 0.05}, 
  AnimationRate -> 30,
  
  "UpdateFunction" -> Function[t,
    points = Table[{xC[t, i], yC[t, i]}, {i, 3}];
    
    False
  ]
] ]
(*VB[*)(FrontEndRef["8564d026-49fc-4c29-b787-b26363dff908"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKW5iamaQYGJnpmlimJeuaJBtZ6iaZW5jrJhmZGZsZp6SlWRpYAAB4WhUt"*)(*]VB*)
info

Click on a play above button to run the animation