Coupled Oscillator: Double Pendulum

By Jimmy Apablaza

This problem is described in Page 321-322, Section 7.6 of the A first Course in Differential Equations textbook, 8ED (ISBN 0-534-41878-3).

Figure 1. Coupled Pendulum.‎

Problem Statement

Consider the double-pendulum system consisting of a pendulum attached to another pendulum shown in Figure 1.

Assumptions:

• the system oscillates vertically under the influence of gravity.
• the mass of both rod are neligible
• no dumpung forces act on the system
• positive direction to the right.

The system of differential equations describing the motion is nonlinear

$(m_{1}+m_{2})l_{1}^{2}\theta _{1}^{\prime \prime }+m_{2}l_{1}l_{2}\theta _{2}^{\prime \prime }cos(\theta _{1}-\theta _{2})+m_{2}l_{1}l_{2}(\theta _{2}^{\prime })^{2}sin(\theta _{1}-\theta _{2})+(m_{1}+m_{2})l_{1}gsin\theta _{1}=0$

$m_{2}l_{1}^{2}\theta _{2}^{\prime \prime }+m_{2}l_{1}l_{2}\theta _{1}^{\prime \prime }cos(\theta _{1}-\theta _{2})-m_{2}l_{1}l_{2}(\theta _{1}^{\prime })^{2}sin(\theta _{1}-\theta _{2})+m_{2}l_{2}gsin\theta _{2}=0$

In order to linearize these equations, we assume that the displacements $\theta _{1}$ and $\theta _{2}$ are small enough so that $cos(\theta _{1}-\theta _{2})\approx 1$ and $sin(\theta _{1}-\theta _{2})\approx 0$. Thus,

$(m_{1}+m_{2})l_{1}^{2}\theta _{1}^{\prime \prime }+m_{2}l_{1}l_{2}\theta _{2}^{\prime \prime }+(m_{1}+m_{2})l_{1}g\theta _{1}=0$

$m_{2}l_{1}^{2}\theta _{2}^{\prime \prime }+m_{2}l_{1}l_{2}\theta _{1}^{\prime \prime }+m_{2}l_{2}g\theta _{2}=0$

Solution

Since our concern is about the motion functions, we will assign the masses $m_{1}$ and $m_{2}$, the rod lenghts $l_{1}$ and $l_{1}$, and gravitational force $g$ constants to different variables as follows,

$A=(m_{1}+m_{2})l_{1}^{2}\quad B=m_{2}l_{1}l_{2}\quad C=(m_{1}+m_{2})l_{1}g\quad D=m_{2}l_{1}^{2}\quad E=m_{2}l_{2}g$

Hence,

$A\theta _{1}^{''}+B\theta _{2}^{''}+C\theta _{1}=0$

$D\theta _{2}^{''}+B\theta _{1}^{''}+E\theta _{2}=0$

Solving for $\theta _{1}^{''}$ and $\theta _{2}^{''}$ we obtain,

$\theta _{1}^{''}=-\left({\dfrac {B}{A}}\right)\theta _{2}^{''}-\left({\dfrac {C}{A}}\right)\theta _{1}$

$\theta _{2}^{''}=-\left({\dfrac {B}{D}}\right)\theta _{1}^{''}-\left({\dfrac {E}{D}}\right)\theta _{2}$

Therefore,

$\theta _{1}^{''}=-\left({\dfrac {CD}{AD+B^{2}}}\right)\theta _{1}-\left({\dfrac {BE}{AD+B^{2}}}\right)\theta _{2}$

$\theta _{2}^{''}=\left({\dfrac {BC}{AD+B^{2}}}\right)\theta _{1}-\left({\dfrac {AE}{AD+B^{2}}}\right)\theta _{2}$

State Space

${\begin{bmatrix}\theta _{1}^{'}\\\theta _{1}^{''}\\\theta _{2}^{'}\\\theta _{2}^{''}\end{bmatrix}}={\widehat {A}}\,{\underline {x}}(t)+{\widehat {B}}\,{\underline {u}}(t)={\begin{bmatrix}0&1&0&0\\&&&\\{\dfrac {-CD}{AD-B^{2}}}&0&{\dfrac {BE}{AD-B^{2}}}&0\\&&&\\0&0&0&1\\&&&\\{\dfrac {BC}{AD-B^{2}}}&0&{\dfrac {-AE}{AD-B^{2}}}&0\\\end{bmatrix}}{\begin{Bmatrix}\theta _{1}\\\theta _{1}^{'}\\\theta _{2}\\\theta _{2}^{'}\end{Bmatrix}}+{\widehat {0}}$

Let's plug some numbers. Knowing $g=32$, and assuming that $m_{1}=3$, $m_{2}=1$, and $l_{1}=l_{2}=16$, the constants defined previously become,

$A=1024\quad B=256\quad C=2048\quad D=256\quad E=512$

Hence, the state space matrix is,

${\begin{bmatrix}\theta _{1}^{'}\\\theta _{1}^{''}\\\theta _{2}^{'}\\\theta _{2}^{''}\end{bmatrix}}={\begin{bmatrix}0&1&0&0\\-{\dfrac {8}{3}}&0&{\dfrac {2}{3}}&0\\0&0&0&1\\{\dfrac {8}{3}}&0&-{\dfrac {8}{3}}&0\\\end{bmatrix}}{\begin{Bmatrix}\theta _{1}\\\theta _{1}^{'}\\\theta _{2}\\\theta _{2}^{'}\end{Bmatrix}}$

Eigenvalues & Eigenvectors

The eigenvalues and eigenvectors are easily obtained with the help of a TI-89 calculator. First, we consider the ${\widehat {A}}$'s identity matrix,

$[\lambda I-A]={\begin{bmatrix}\lambda &1&0&0\\-{\dfrac {8}{3}}&\lambda &{\dfrac {2}{3}}&0\\0&0&\lambda &1\\{\dfrac {8}{3}}&0&-{\dfrac {8}{3}}&\lambda \\\end{bmatrix}}$

Once we define the ${\widehat {A}}$ matrix, the eigenvalues are determined by using the eigVl() function,

$\lambda _{1}=2\mathbf {i}$

$\lambda _{2}=-2\mathbf {i}$

$\lambda _{3}=1.1547\mathbf {i}$

$\lambda _{4}=-1.1547\mathbf {i}$

On the other hand, we use the eigVc() function to find the eigenvectors,

$k_{1}={\begin{bmatrix}-0.2\mathbf {i} \\0.4\\0.4\mathbf {i} \\-0.8\\\end{bmatrix}}\quad k_{2}={\begin{bmatrix}0.2\mathbf {i} \\0.4\\-0.4\mathbf {i} \\-0.8\\\end{bmatrix}}\quad k_{3}={\begin{bmatrix}-0.29277\mathbf {i} \\0.33806\\0.58554\mathbf {i} \\-0.67621\\\end{bmatrix}}\quad k_{4}={\begin{bmatrix}0.29277\mathbf {i} \\0.33806\\0.58554\mathbf {i} \\0.67621\\\end{bmatrix}}$

Standard Equation

Now, we plug the eigenvalues and eigenvectors to produce the standar equation,

$x=c_{1}k_{1}e^{\lambda _{1}t}+c_{2}k_{2}e^{\lambda _{2}t}+c_{3}k_{3}e^{\lambda _{3}t}+c_{4}k_{4}e^{\lambda _{4}t}$

$x=c_{1}{\begin{bmatrix}-0.2\mathbf {i} \\0.4\\0.4\mathbf {i} \\-0.8\\\end{bmatrix}}e^{2\mathbf {i} }+c_{2}{\begin{bmatrix}0.2\mathbf {i} \\0.4\\-0.4\mathbf {i} \\-0.8\\\end{bmatrix}}e^{-2\mathbf {i} }+c_{3}{\begin{bmatrix}-0.29277\mathbf {i} \\0.33806\\0.58554\mathbf {i} \\-0.67621\\\end{bmatrix}}e^{1.1547\mathbf {i} }+c_{4}{\begin{bmatrix}0.29277\mathbf {i} \\0.33806\\0.58554\mathbf {i} \\0.67621\\\end{bmatrix}}e^{-1.1547\mathbf {i} }$

Matrix Exponential

The matrix exponential is,

${\dot {\bar {z}}}={\hat {A}}{\bar {z}}=TAT^{-1}{\bar {z}}$

where

${\hat {A}}={\begin{bmatrix}0&1&0&0\\-{\dfrac {8}{3}}&0&{\dfrac {2}{3}}&0\\0&0&0&1\\{\dfrac {8}{3}}&0&-{\dfrac {8}{3}}&0\\\end{bmatrix}}$,

and

$T^{-1}=[k_{1}|k_{2}|k_{3}|k_{4}]={\begin{bmatrix}-0.2\mathbf {i} &0.2\mathbf {i} &-0.29277\mathbf {i} &0.29277\mathbf {i} \\0.4&0.4&0.33806&0.33806\\0.4\mathbf {i} &-0.4\mathbf {i} &0.58554\mathbf {i} &0.58554\mathbf {i} \\-0.8&-0.8&-0.67621&0.67621\\\end{bmatrix}}$,

so

$T=(T^{-1})^{-1}=[k_{1}|k_{2}|k_{3}|k_{4}]^{-1}={\begin{bmatrix}1.25\mathbf {i} &0.625&-0.625\mathbf {i} &-0.3125\\-1.25\mathbf {i} &0.625&0.625\mathbf {i} &-0.3125\\0.853913\mathbf {i} &0.73951&0.426956\mathbf {i} &0.369755\\-0.853913\mathbf {i} &0.73951&-0.426956\mathbf {i} &0.369755\\\end{bmatrix}}$

Again, we can resort to the TI-89 calculator. As it is mentioned above, the matrix exponential is obtained by typing eigVc(a)^-1*a*eigVc(a), where a is the ${\widehat {A}}$ matrix. Thus,

${\dot {\bar {z}}}={\hat {A}}{\bar {z}}=TAT^{-1}{\bar {z}}={\begin{bmatrix}2\mathbf {i} &0&0&0\\0&-2\mathbf {i} &0&0\\0&0&1.1547\mathbf {i} &0\\0&0&0&-1.1547\mathbf {i} \\\end{bmatrix}}{\bar {z}}$