This article focuses on deriving the equations of motion of the coupled elastic pendulums depicted in Figure 1. If you would like to see their solution, they can be found here.

We will derive our own equations of motion based on the Euler-Lagrange equation with dissipative forces

$\begin{aligned} \dfrac{d}{dt} \dfrac{\partial \mathcal{L}}{\partial \dot{q}_i} - \dfrac{\partial \mathcal{L}}{\partial q_i} &= \sum_j \vec{F}_{D,j} \cdot \dfrac{\partial \vec{r}_j}{\partial q_i}. \end{aligned}$

Where

• $j$ refers to the component of the system we are analysing.

• $(q_i)$ are the generalized coordinates of the system.

• $(\dot{q}_i)$ are the first time derivatives of the generalized coordinates of the system.

• $\vec{r}_j$ is the position vector of component $j$ of the system.

• $\mathcal{L}$ is the Lagrangian — the difference between the kinetic and potential energy — of the system.

• $p_i = \dfrac{\partial \mathcal{L}}{\partial \dot{q}_i}$ is the generalized momentum canonical to $q_i$.

• $F_i = \dfrac{\partial \mathcal{L}}{\partial q_i}$ is the generalized force canonical to $q_i$.

• $\vec{F}_{D,j}$ is the dissipative force vector for component $j$.

• $\hat{e}_{j,i} = \dfrac{\partial \vec{r}_j}{\partial q_i}$ is the generalized basis vector canonical to $q_i$ for component $j$ of the system.

• The left-hand side of Equation Equation (1) can also be represented as $-\dfrac{\delta \mathcal{L}}{\delta q_i}$, where $\dfrac{\delta \mathcal{L}}{\delta q_i}$ is the functional derivative of the Lagrangian with respect to $q_i$. To simplify things, we will call $-\dfrac{\delta \mathcal{L}}{\delta q_i} = \dfrac{\delta' \mathcal{L}}{\delta' q_i}$

• The right-hand side of Equation Equation (1) is also called the generalized dissipative force and can be represented as $Q_i$.

Caveats

The masses of the pendulum rods (or springs) and the friction they experience are ignored as including them into the calculation for rigid double pendulums does not make things more interesting and merely complicates the calculation. I did try deriving the equations of motion while including the rods using SymPy, specifically with the code

from sympy import symbols, Function, diff, cos, sin, simplify, sqrt, Abs, Eq, solve, latex
from sympy.vector import CoordSys3D
from multiprocessing import Pool, cpu_count
N = CoordSys3D('N');
# Set up symbols
t = symbols('t')
m1b = symbols('m1b'); 
m2b = symbols('m2b'); 
m1r = symbols('m1r'); 
m2r = symbols('m2r'); 
l1 = symbols('11');
l2 = symbols('l2');
k1 = symbols('k1'); 
k2 = symbols('k2'); 
g = symbols('g'); 
b1r = symbols('b1r'); 
b2r = symbols('b2r'); 
b1b = symbols('b1b'); 
b2b = symbols('b2b'); 
c1r = symbols("c1r"); 
c1b = symbols("c1b"); 
c2r = symbols("c2r"); 
c2b = symbols("c2b");
# Functions
r1 = Function('r1')(t); 
theta1=Function('theta1')(t); 
r2 = Function('r2')(t); 
theta2=Function('theta2')(t);

# Bob 1
x1b = r1*cos(theta1);
y1b = r1*sin(theta1);
x1bdot = diff(x1b, t);
y1bdot = diff(y1b, t);
v1b_mag = simplify(sqrt(x1bdot**2 + y1bdot**2));
v1b_sq = simplify(x1bdot**2 + y1bdot**2);
r1b = x1b * N.i + y1b * N.j;
v1b = x1bdot * N.i + y1bdot * N.j;
e1b_r1 = diff(r1b, r1);
e1b_r2 = diff(r1b, r2);
e1b_th1 = diff(r1b, theta1);
e1b_th2 = diff(r1b, theta2);

# Bob 2
x2b = x1b+r2*cos(theta2);
y2b = y1b + r2*sin(theta2);
x2bdot = diff(x2b, t);
y2bdot = diff(y2b, t);
v2b_mag = simplify(sqrt(x2bdot**2 + y2bdot**2));
v2b_sq = simplify(x2bdot**2 + y2bdot**2);
r2b = x2b * N.i + y2b * N.j;
v2b = x2bdot * N.i + y2bdot * N.j;
e2b_r1 = diff(r2b, r1);
e2b_r2 = diff(r2b, r2);
e2b_th1 = diff(r2b, theta1);
e2b_th2 = diff(r2b, theta2);


# Rod 1
x1r = x1b/2;
x1rdot = diff(x1r, t);
y1r = y1b/2;
y1rdot = diff(y1r, t);
v1r_mag = simplify(sqrt(x1rdot**2 + y1rdot**2));
v1r_sq = simplify(x1rdot**2 + y1rdot**2);
v1r = x1rdot * N.i + y1rdot * N.j;
r1r = x1r * N.i + y1r * N.j;
v1r = x1rdot * N.i + y1rdot * N.j;
e1r_r1 = diff(r1r, r1);
e1r_r2 = diff(r1r, r2);
e1r_th1 = diff(r1r, theta1);
e1r_th2 = diff(r1r, theta2);

# Rod 2
x2r = x1b +r2*cos(theta2)/2;
x2rdot = diff(x2r, t);
y2r = y1b + r2*sin(theta2)/2;
y2rdot = diff(y2r, t);
v2r_mag = simplify(sqrt(x2rdot**2 + y2rdot**2));
v2r_sq = simplify(x2rdot**2 + y2rdot**2);
v2r = x2rdot * N.i + y2rdot * N.j;
r2r = x2r * N.i + y2r * N.j;
v2r = x2rdot * N.i + y2rdot * N.j;
e2r_r1 = diff(r2r, r1);
e2r_r2 = diff(r2r, r2);
e2r_th1 = diff(r2r, theta1);
e2r_th2 = diff(r2r, theta2);

Frod1 = -(b1r+c1r*Abs(v1r_mag))*v1r;
Frod2 = -(b2r+c2r*Abs(v2r_mag))*v2r;
Fbob1 = -(b1b+c1b*Abs(v1b_mag))*v1b;
Fbob2 = -(b2b+c2b*Abs(v2b_mag))*v2b;
Qr1 = Frod1.dot(e1r_r1) + Frod2.dot(e2r_r1) + Fbob1.dot(e1b_r1) + Fbob2.dot(e2b_r1);
Qr1 = simplify(Qr1);
Qr2 = Frod1.dot(e1r_r2) + Frod2.dot(e2r_r2) + Fbob1.dot(e1b_r2) + Fbob2.dot(e2b_r2);
Qr2 = simplify(Qr2);
Qth1 = Frod1.dot(e1r_th1) + Frod2.dot(e2r_th1) + Fbob1.dot(e1b_th1) + Fbob2.dot(e2b_th1);
Qth1 = simplify(Qth1);
Qth2 = Frod1.dot(e1r_th2) + Frod2.dot(e2r_th2) + Fbob1.dot(e1b_th2) + Fbob2.dot(e2b_th2);
Qth2 = simplify(Qth2);
T = m1b/2 * v1b_sq + m2b/2 * v2b_sq + m1r/2 * v1r_sq + m2r/2 * v2r_sq + m1r/24*r1**2*diff(theta1, t)**2 + m2r/24 * r2**2*diff(theta2,t)**2
V = m1b * g * y1b + m1r * g * r1 * y1r + m2b * g * y2b + m2r * g * y2r + k1*(r1-l1)**2/2 + k2*(r2-l2)**2/2;
L = T - V; 
def compute_eq_of_motion(args):
    L, Q, coord = args
    t = symbols('t')
    lhs = diff(diff(L, diff(coord, t)), t) - diff(L, coord)
    return Eq(lhs, Q)
	
# Below line was commented out because it caused individual equations to be 
# incredibly long, so I decided to treat generalized dissipation forces as 
# variables. 
#Qs = [Qr1, Qr2, Qth1, Qth2];
Qs = [symbols("Qr1"), symbols("Qr2"), symbols("Qtheta1"), symbols("Qtheta2")]
coords = [r1, r2, theta1, theta2]
d2r1 = diff(r1, t, 2);
d2r2 = diff(r2, t, 2);
d2th1 = diff(theta1, t, 2);
d2th2 = diff(theta2, t, 2);

with Pool(cpu_count()) as pool:
    equations = pool.map(compute_eq_of_motion, [(L, Qs[i], coords[i]) for i in range(4)])

d2 = [d2r1, d2r2, d2th1, d2th2];
# Solve the system
sols = solve(equations, d2, simplify=True)

secdernames = ["d2r1", "d2r2", "d2theta1", "d2theta2"]
for i in range(4):
	print(secdernames[i] + " = \n" + latex(sols[d2[i]]))

And it caused IPython to crash.

Coordinates, velocities and generalized basis vectors

As can be seen, we have four degrees of freedom in this system. The angles the two pendulums make with the positive $x$-axis — $\theta_1$ and $\theta_2$, respectively — are among our degrees of freedom. We will also need degrees of freedom corresponding to the lengths of the pendulum rods. These degrees of freedom could either be the extent to which they are extended beyond their rest length or their total length. For the sake of simplicity, we will opt to use their total lengths — $r_1$ and $r_2$, respectively. Hence

This means that the velocity of the first pendulum bob is

Hence

As for the velocity of the second pendulum bob, it is

Let $\Delta = \theta_2-\theta_1$, then square of the velocity is

Let us define $|\Delta \vec{v}_{21}|^2 = |\vec{v}_2|^2 - |\vec{v}_1|^2$, as this will simplify our Lagrangian later.

Dissipative forces

We will assume that the dissipative forces are proportion to the velocity and velocity squared of the pendulum bobs. Meaning they will have the form

Where $j$ is the pendulum bob of interest, $b_j$ and $c_j$ are constants.

Kinetic energy

The kinetic energy of the system is given by

Potential energy

The potential energy of the system is given by

Lagrangian

Hence the Lagrangian of the system is

We will not expand this Lagrangian, as doing so just adds to its complexity. Instead, we will calculate the derivatives of each of its components.

Derivative of components of the Lagrangian

Square of the velocity of the first pendulum's bob

The relevant partial and standard derivatives are:

Hence the negative functional derivatives are

Difference in the square of each bob's velocity

Hence the partial and standard derivatives of the difference in the square of each bob's velocity is

Let us define $\dot{\Delta}_1 = 2\dot{\theta}_1 - \dot{\theta}_2$ and $\dot{\Delta}_2 = 2\dot{\theta}_2 - \dot{\theta}_1$.

Hence the negative functional derivative for $r_1$ is

Where $\dot{\Delta} + \dot{\theta}_1 = \dot{\theta}_2 - \dot{\theta}_1 + \dot{\theta}_1 = \dot{\theta}_2$ and $\dot{\Delta}_2 + \dot{\theta}_1 = 2\dot{\theta}_2 - \dot{\theta}_1 + \dot{\theta}_1 = 2\dot{\theta}_2$.

As for $r_2$

Hence $\dot{\Delta}-\dot{\theta}_2 = \dot{\theta}_2-\dot{\theta}_1-\dot{\theta}_2 = -\dot{\theta}_1$ and $\dot{\Delta}_1 + \dot{\theta}_2 = 2\dot{\theta}_1 - \dot{\theta}_2 + \dot{\theta}_2 = 2\dot{\theta}_1$.

As for $\theta_1$

As for $\theta_2$

Euler-Lagrange equations with dissipation

$r_1$

It is important to note that $\dfrac{\delta' f(q_i)}{\delta' q_i} = -\dfrac{\partial f}{\partial q_i}$ and of course if a term does not depend on $q_i$ or $\dot{q}_i$ its functional derivative with respect to $q_i$ is zero. Hence

We have deliberately ignored the $m_2gr_2\sin{\theta_2}$ and $-\dfrac{k_2(r_2-l_2)^2}{2}$ as they are independent of $r_1$.

The generalized dissipation force canonical to $r_1$ is hence

Hence the Euler-Lagrange equation for $r_1$ with dissipative forces is

Dividing by $m_1+m_2$

Expanding out second derivative terms and placing them first on the left-hand side yields

Moving all terms that do not involve second derivatives to the right-hand side yields

$r_2$

As for $r_2$

Hence the Euler-Lagrange equation for $r_2$ with dissipative forces is

Dividing by $m_2$ yields

Next we will expand out second time derivatives and moving everything else to the right-hand side

$\theta_1$

As for $\theta_1$

Hence Equation (1) is

Dividing by $(m_1+m_2)r_1^2$ yields

Expanding out all second time derivatives and moving all other terms to the right-hand side yields

$\theta_2$

As for $\theta_2$

Hence Equation (1) for $\theta_2$ is

Dividing by $m_2r_2^2$ yields

Next we will expand out second time derivatives and moving everything else to the right-hand side

Analysis

There are three ways we could solve this problem; each of which involves numerical integration to obtain the final solution. Firstly, we could algebraically manipulate this ordinary differential equation (ODE) system until only one second time derivative appears in each equation. The final ODE system after this manipulation could, in turn, be numerically integrated using any standard scheme (e.g. the Runge-Kutta-Fehlberg method). Secondly, we could numerically integrate it as is using a differential-algebraic equation (DAE) solver, but these tend to be more prone to give convergence errors in my experience. Finally, we could convert the system to matrix form and invert it to obtain numerical approximations of each of our second time derivatives and use these to numerically integrate the system with a standard ODE solver. We will opt for this last method, as the first approach is almost guaranteed to introduce errors and the second gives convergence errors, at least in Julia.

Essentially, we will write our differential equation system as

where $\mathbf{\ddot{q}}$ is a vector containing the second time derivatives of our generalized coordinates. Hence $\mathbf{\ddot{q}}=\mathbf{A}^{-1} \mathbf{b}$. Or, in full form, this above equation is

So the solution is

Ordinarily, I would use a symplectic method to integrate a classical mechanical system, as it minimizes cumulative errors in the Hamiltonian over long-term integration. That being said, symplectic methods typically are either implicit (which significantly increases their computational complexity) or require a separable Hamiltonian. The Hamiltonian of this system would not separable as the kinetic energy cannot be written entirely in terms of momenta. Additionally, rewriting the kinetic energy in terms of the momenta would be tedious and prone to error.