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).

Error creating thumbnail: File missing

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 negligible
• no damping 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}B=m_2{l}_1{l}_{2}C=(m_1+m_2)l_{1}gD=m_2{l}_1^{2}E=m_2{l}_{2}g$

Hence,

$A{\theta }_1^{{\text{'}}^{\prime }}+B{\theta }_2^{{\text{'}}^{\prime }}+C{\theta }_1=0$

$D{\theta }_2^{{\text{'}}^{\prime }}+B{\theta }_1^{{\text{'}}^{\prime }}+E{\theta }_2=0$

Solving for ${\theta }_1^{{\text{'}}^{\prime }}$ and ${\theta }_2^{{\text{'}}^{\prime }}$ we obtain,

${\theta }_1^{{\text{'}}^{\prime }}=-\left(\frac{B}{A}\right){\theta }_2^{{\text{'}}^{\prime }}-\left(\frac{C}{A}\right){\theta }_1$

${\theta }_2^{{\text{'}}^{\prime }}=-\left(\frac{B}{D}\right){\theta }_1^{{\text{'}}^{\prime }}-\left(\frac{E}{D}\right){\theta }_2$

Therefore,

${\theta }_1^{{\text{'}}^{\prime }}=-\left(\frac{CD}{AD+B^2}\right){\theta }_1-\left(\frac{BE}{AD+B^2}\right){\theta }_2$

${\theta }_2^{{\text{'}}^{\prime }}=\left(\frac{BC}{AD+B^2}\right){\theta }_1-\left(\frac{AE}{AD+B^2}\right){\theta }_2$

State Space

$\left[\begin{array}{c}{\theta }_1^{\text{'}}\\ {\theta }_1^{{\text{'}}^{\prime }}\\ {\theta }_2^{\text{'}}\\ {\theta }_2^{{\text{'}}^{\prime }}\end{array}\right]=\hat{A}\underset{\_}{x}(t)+\hat{B}\underset{\_}{u}(t)=\left[\begin{array}{cccc}0& 1& 0& 0\\ & & & \\ \frac{-CD}{AD-B^2}& 0& \frac{BE}{AD-B^2}& 0\\ & & & \\ 0& 0& 0& 1\\ & & & \\ \frac{BC}{AD-B^2}& 0& \frac{-AE}{AD-B^2}& 0\\ \end{array}\right]\left\{\begin{array}{c}{\theta }_1\\ {\theta }_1^{\text{'}}\\ {\theta }_2\\ {\theta }_2^{\text{'}}\end{array}\right\}+\hat{0}$

Let's plug some numbers. It's known that $g=32$. In addition, we assume that $m_1=3$, $m_2=1$, and $l_1=l_2=16$, so the constants defined previously become,

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

Hence, the state space matrix is,

$\left[\begin{array}{c}{\theta }_1^{\text{'}}\\ {\theta }_1^{{\text{'}}^{\prime }}\\ {\theta }_2^{\text{'}}\\ {\theta }_2^{{\text{'}}^{\prime }}\end{array}\right]=\left[\begin{array}{cccc}0& 1& 0& 0\\ -\frac{8}{3}& 0& \frac{2}{3}& 0\\ 0& 0& 0& 1\\ \frac{8}{3}& 0& -\frac{8}{3}& 0\\ \end{array}\right]\left\{\begin{array}{c}{\theta }_1\\ {\theta }_1^{\text{'}}\\ {\theta }_2\\ {\theta }_2^{\text{'}}\end{array}\right\}$

Eigenvalues & Eigenvectors

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

$[\lambda I-A]=\left[\begin{array}{cccc}\lambda & 1& 0& 0\\ -\frac{8}{3}& \lambda & \frac{2}{3}& 0\\ 0& 0& \lambda & 1\\ \frac{8}{3}& 0& -\frac{8}{3}& \lambda \\ \end{array}\right]$

Once we define the $\hat{A}$ matrix, the eigenvalues are determined by using the eigVi() 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=\left[\begin{array}{c}-0.2\mathbf{i}\\ 0.4\\ 0.4\mathbf{i}\\ -0.8\\ \end{array}\right]k_2=\left[\begin{array}{c}0.2\mathbf{i}\\ 0.4\\ -0.4\mathbf{i}\\ -0.8\\ \end{array}\right]k_3=\left[\begin{array}{c}-0.29277\mathbf{i}\\ 0.33806\\ 0.58554\mathbf{i}\\ -0.67621\\ \end{array}\right]k_4=\left[\begin{array}{c}0.29277\mathbf{i}\\ 0.33806\\ 0.58554\mathbf{i}\\ 0.67621\\ \end{array}\right]$

Standard Equation

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

$\underset{\_}{x}=c_1\underset{\_}{k_1}{e}^{{\lambda }_{1}t}+c_2{\underset{\_}{k}}_2{e}^{{\lambda }_{2}t}+c_3\underset{\_}{k_3}{e}^{{\lambda }_{3}t}+c_4\underset{\_}{k_4}{e}^{{\lambda }_{4}t}$

$\underset{\_}{x}=c_1\left[\begin{array}{c}-0.2\mathbf{i}\\ 0.4\\ 0.4\mathbf{i}\\ -0.8\\ \end{array}\right]e^{2\mathbf{i}}+c_2\left[\begin{array}{c}0.2\mathbf{i}\\ 0.4\\ -0.4\mathbf{i}\\ -0.8\\ \end{array}\right]e^{-2\mathbf{i}}+c_3\left[\begin{array}{c}-0.29277\mathbf{i}\\ 0.33806\\ 0.58554\mathbf{i}\\ -0.67621\\ \end{array}\right]e^{1.1547\mathbf{i}}+c_4\left[\begin{array}{c}0.29277\mathbf{i}\\ 0.33806\\ 0.58554\mathbf{i}\\ 0.67621\\ \end{array}\right]e^{-1.1547\mathbf{i}}$

Matrix Exponential

The matrix exponential is,

$\dot{\overline{z}}=\hat{A}\overline{z}=TA{T}^{-1}\overline{z}$

where

$A=\left[\begin{array}{cccc}0& 1& 0& 0\\ -\frac{8}{3}& 0& \frac{2}{3}& 0\\ 0& 0& 0& 1\\ \frac{8}{3}& 0& -\frac{8}{3}& 0\\ \end{array}\right]$,

and

$T^{-1}=[k_1|k_2|k_3|k_4]=\left[\begin{array}{cccc}-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{array}\right]$,

so

$T=(T^{-1}{)}^{-1}=[k_1|k_2|k_3|k_4{]}^{-1}=\left[\begin{array}{cccc}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{array}\right]$

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 $\hat{A}$ matrix. Thus,

$\dot{\overline{z}}=\hat{A}\overline{z}=TA{T}^{-1}\overline{z}=\left[\begin{array}{cccc}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{array}\right]\overline{z}$

So, the exponential matrix becomes,

$\dot{\underset{\_}{x}}=\hat{A}\underset{\_}{x}=T^{-1}{e}^{\hat{A}t}T\underset{\_}{x}=e^{At}\underset{\_}{x}$

where

$e^{\hat{A}t}=\left[\begin{array}{cccc}{e}^{2\mathbf{i}t}& 0& 0& 0\\ 0& e^{-2\mathbf{i}t}& 0& 0\\ 0& 0& e^{1.1547\mathbf{i}t}& 0\\ 0& 0& 0& e^{-1.1547\mathbf{i}t}\\ \end{array}\right]$

Hence,

$\dot{\underset{\_}{x}}=\left[\begin{array}{cccc}-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{array}\right]\left[\begin{array}{cccc}{e}^{2\mathbf{i}t}& 0& 0& 0\\ 0& e^{-2\mathbf{i}t}& 0& 0\\ 0& 0& e^{1.1547\mathbf{i}t}& 0\\ 0& 0& 0& e^{-1.1547\mathbf{i}t}\\ \end{array}\right]\left[\begin{array}{cccc}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{array}\right]$