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

## Contents

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_2l_1l_2\theta_2^{\prime\prime}cos(\theta_1-\theta_2) + m_2l_1l_2(\theta_2^{\prime})^2sin(\theta_1-\theta_2) + (m_1+m_2)l_1gsin\theta_1 = 0$
$m_2l_1^2\theta_2^{\prime\prime} + m_2l_1l_2\theta_1^{\prime\prime}cos(\theta_1-\theta_2) - m_2l_1l_2(\theta_1^{\prime})^2sin(\theta_1-\theta_2) + m_2l_2gsin\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)\approx1$ and $sin(\theta_1-\theta_2)\approx0$. Thus,

$(m_1+m_2)l_1^2\theta_1^{\prime\prime} + m_2l_1l_2\theta_2^{\prime\prime} + (m_1+m_2)l_1g\theta_1 = 0$
$m_2l_1^2\theta_2^{\prime\prime} + m_2l_1l_2\theta_1^{\prime\prime} + m_2l_2g\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_2l_1l_2 \quad C=(m_1+m_2)l_1g \quad D=m_2l_1^2 \quad E=m_2l_2g$

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. 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=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 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 = \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,

$\underline{x} = c_1 \underline{k_1} e^{\lambda_1 t} + c_2 \underline{k}_2 e^{\lambda_2 t} + c_3 \underline{k_3} e^{\lambda_3 t} + c_4 \underline{k_4} e^{\lambda_4 t}$
$\underline{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

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

So, the exponential matrix becomes,

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

where

$e^{\hat{A}t} = \begin{bmatrix} 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{bmatrix}$

Hence,

$\dot{\underline{x}} = \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} \begin{bmatrix} 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{bmatrix} \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} \underline{x}$

so, the matrix exponential can be solved using matrix multiplication.