Deriving the equations of motion for the double elastic pendulum

Figure 1: Diagram of the double elastic pendulum.

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

ddt∂L∂q˙i−∂L∂qi=∑jF⃗D,j⋅∂r⃗j∂qi.\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

  • jj refers to the component of the system we are analysing.

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

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

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

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

  • pi=∂L∂q˙ip_i = \dfrac{\partial \mathcal{L}}{\partial \dot{q}_i} is the generalized momentum canonical to qiq_i.

  • Fi=∂L∂qiF_i = \dfrac{\partial \mathcal{L}}{\partial q_i} is the generalized force canonical to qiq_i.

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

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

  • The left-hand side of Equation Equation (1) can also be represented as −δLδqi-\dfrac{\delta \mathcal{L}}{\delta q_i}, where δLδqi\dfrac{\delta \mathcal{L}}{\delta q_i} is the functional derivative of the Lagrangian with respect to qiq_i. To simplify things, we will call −δLδqi=δ′Lδ′qi-\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 QiQ_i.

  1. Coordinates, velocities and generalized basis vectors
    1. Bob 1
    2. Rod 1
      1. Parameterization
    3. Centre of mass
    4. Bob 2
    5. Rod 2
      1. Parameterization
    6. Centre of mass
  2. Dissipative forces
  3. Kinetic energy
    1. Rod 1
    2. Rod 2
    3. Final result
  4. Potential energy
  5. Lagrangian
  6. Derivative of components of the Lagrangian
    1. Square of the velocity of the first pendulum's bob
    2. Second bob's independent velocity squared
    3. Difference in independent velocities squared
  7. Euler-Lagrange equations with dissipation
    1. r1r_1
    2. r2r_2
    3. θ1\theta_1
    4. θ2\theta_2
    5. Matrix form

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 xx-axis — θ1\theta_1 and θ2\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 — r1r_1 and r2r_2, respectively. Hence

Bob 1

x1b=r1cos⁡θ1x˙1b=r˙1cos⁡θ1−r1θ˙1sin⁡θ1y1b=r1sin⁡θ1y˙1b=r˙1sin⁡θ1+r1θ˙1cos⁡θ1\begin{aligned} x_{1b} &= r_1 \cos{\theta_1} & \dot{x}_{1b} &= \dot{r}_1 \cos{\theta_1} - r_1 \dot{\theta}_1 \sin{\theta_1}\\ y_{1b} &= r_1 \sin{\theta_1} & \dot{y}_{1b} &= \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1 \cos{\theta_1} \end{aligned}

This means that the velocity of the first pendulum bob is

v⃗1b=[x˙1by˙1b]=[r˙1cos⁡θ1−r1θ˙1sin⁡θ1r˙1sin⁡θ1+r1θ˙1cos⁡θ1].\begin{aligned} \vec{v}_{1b} &= \begin{bmatrix} \dot{x}_{1b} \\ \dot{y}_{1b} \end{bmatrix} \\ &= \begin{bmatrix} \dot{r}_1 \cos{\theta_1} - r_1 \dot{\theta}_1 \sin{\theta_1} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1 \cos{\theta_1} \end{bmatrix}. \end{aligned}

Hence

∣v⃗1b∣2=r˙12+r12θ˙12.\begin{aligned} |\vec{v}_{1b}|^2 &= \dot{r}_1^2 + r_1^2 \dot{\theta}_1^2. \end{aligned}

Rod 1

We will calculate the rod positions and velocities two different ways. For calculating the kinetic energy, we will parameterize them in terms of ss — the distance along the pendulum — so that we can integrate over them. For calculating the potential energy and generalized dissipative force, we will calculate them using a centre of mass approach, as we will assume centre of mass approximations will work for these.

Parameterization

x1r=scos⁡θ1x˙1r=s˙cos⁡θ1−sθ˙1sin⁡θ1y1r=ssin⁡θ1y˙1r=s˙sin⁡θ1+sθ˙1cos⁡θ1∣v⃗1r∣2=s˙2+s2θ˙12=s2r˙12r12+s2θ˙12=s2[r˙12r12+θ˙12].\begin{aligned} x_{1r} &= s\cos{\theta_1} & \dot{x}_{1r} = \dot{s}\cos{\theta_1} - s\dot{\theta}_1 \sin{\theta_1} \\ y_{1r} &= s\sin{\theta_1} & \dot{y}_{1r} = \dot{s}\sin{\theta_1} + s\dot{\theta}_1 \cos{\theta_1} \\ |\vec{v}_{1r}|^2 &= \dot{s}^2 + s^2\dot{\theta}_1^2 \\ &= \dfrac{s^2\dot{r}_1^2}{r_1^2} + s^2\dot{\theta}_1^2 \\ &= s^2 \left[\dfrac{\dot{r}_1^2}{r_1^2}+\dot{\theta}_1^2\right]. \end{aligned}

Here we have assumed that the motion of the pendulum rod is uniform.

Centre of mass

x1r=r1cos⁡θ12x˙1r=r1˙cos⁡θ12−r1θ˙1sin⁡θ12y1r=r1sin⁡θ12y˙1r=r1˙sin⁡θ12+r1θ˙1cos⁡θ12v⃗1r=12[r˙1cos⁡θ1−r1θ˙1sin⁡θ1r˙1sin⁡θ1+r1θ˙1cos⁡θ1]∣v⃗1r∣2=14(r˙12+r12θ˙12)=∣v⃗1b∣24.\begin{aligned} x_{1r} &= \dfrac{r_1\cos{\theta_1}}{2} & \dot{x}_{1r} = \dfrac{\dot{r_1}\cos{\theta_1}}{2} - \dfrac{r_1\dot{\theta}_1 \sin{\theta_1}}{2} \\ y_{1r} &= \dfrac{r_1\sin{\theta_1}}{2} & \dot{y}_{1r} = \dfrac{\dot{r_1}\sin{\theta_1}}{2} + \dfrac{r_1\dot{\theta}_1 \cos{\theta_1}}{2} \\ \vec{v}_{1r} &= \dfrac{1}{2} \begin{bmatrix} \dot{r}_1 \cos{\theta_1} - r_1\dot{\theta}_1\sin{\theta_1} \\ \dot{r}_1 \sin{\theta_1} + r_1\dot{\theta}_1\cos{\theta_1} \end{bmatrix} \\ |\vec{v}_{1r}|^2 &= \dfrac{1}{4} \left(\dot{r}_1^2 + r_1^2 \dot{\theta}_1^2\right)\\ &= \dfrac{|\vec{v}_{1b}|^2}{4}. \end{aligned}

Bob 2

x2b=x1b+r2cos⁡θ2x˙2=x˙1b+r˙2cos⁡θ2−r2θ˙2sin⁡θ2y2b=y1b+r2sin⁡θ2y˙2=y˙1b+r˙2sin⁡θ2+r2θ˙2cos⁡θ2.\begin{aligned} x_{2b} &= x_{1b} + r_2\cos{\theta_2} & \dot{x}_2 &= \dot{x}_{1b} + \dot{r}_2\cos{\theta_2} - r_2\dot{\theta}_2 \sin{\theta_2} \\ y_{2b} &= y_{1b} + r_2\sin{\theta_2} & \dot{y}_2 &= \dot{y}_{1b} + \dot{r}_2\sin{\theta_2} + r_2\dot{\theta}_2 \cos{\theta_2}. \end{aligned}

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

v⃗2b=[r˙1cos⁡θ1−r1θ˙1sin⁡θ1+r˙2cos⁡θ2−r2θ˙2sin⁡θ2r˙1sin⁡θ1+r1θ˙1cos⁡θ1+r˙2sin⁡θ2+r2θ˙2cos⁡θ2].\begin{aligned} \vec{v}_{2b} &= \begin{bmatrix} \dot{r}_1 \cos{\theta_1} - r_1 \dot{\theta}_1 \sin{\theta_1} + \dot{r}_2\cos{\theta_2} - r_2\dot{\theta}_2 \sin{\theta_2} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1 \cos{\theta_1} + \dot{r}_2\sin{\theta_2} + r_2\dot{\theta}_2 \cos{\theta_2} \end{bmatrix}. \end{aligned}

Let Δ=θ2−θ1\Delta = \theta_2-\theta_1, then square of the velocity is

∣v⃗2b∣2=r˙12cos⁡2θ1+r12θ˙12sin⁡2θ1+r˙22cos⁡2θ2+r22θ˙22sin⁡2θ2−2r1r˙1θ˙1cos⁡θ1sin⁡θ1+2r˙1r˙2cos⁡θ1cos⁡θ2−2r˙1r2θ˙2cos⁡θ1sin⁡θ2−2r1r˙2θ˙1sin⁡θ1cos⁡θ2+2r1r2θ˙1θ˙2sin⁡θ1sin⁡θ2−2r˙2r2θ˙2cos⁡θ2sin⁡θ2+r˙12sin⁡2θ1+r12θ˙12cos⁡2θ1+r˙22sin⁡2θ2+r22θ˙22cos⁡2θ2+2r1r˙1θ˙1sin⁡θ1cos⁡θ1+2r˙1r˙2sin⁡θ1sin⁡θ2+2r˙1r2θ˙2sin⁡θ1cos⁡θ2+2r1r˙2θ˙1cos⁡θ1sin⁡θ2+2r1r2θ˙1θ˙2cos⁡θ1cos⁡θ2+2r2r˙2θ˙2sin⁡θ2cos⁡θ2=r˙12+r12θ˙12+r˙22+r22θ˙22+2r1r˙1θ˙1(−cos⁡θ1sin⁡θ1+cos⁡θ1sin⁡θ1)+2r˙1r˙2(cos⁡θ1cos⁡θ2+sin⁡θ1sin⁡θ2)+2r˙1r2θ˙2(−cos⁡θ1sin⁡θ2+sin⁡θ1cos⁡θ2)+2r1r˙2θ˙1(−sin⁡θ1cos⁡θ2+cos⁡θ1sin⁡θ2)+2r1r2θ˙1θ˙2(sin⁡θ1sin⁡θ2+cos⁡θ1cos⁡θ2)+2r2r˙2θ˙2(−cos⁡θ2sin⁡θ2+sin⁡θ2cos⁡θ2)=r˙12+r12θ˙12+r˙22+r22θ˙22+2r˙1r˙2cos⁡(θ2−θ1)−2r˙1r2θ2˙sin⁡(θ2−θ1)+2r1r˙2θ˙1sin⁡(θ2−θ1)+2r1r2θ˙1θ˙2cos⁡(θ2−θ1)=r˙12+r12θ˙12+r˙22+r22θ˙22+2r˙1r˙2cos⁡Δ−2r˙1r2θ2˙sin⁡Δ+2r1r˙2θ˙1sin⁡Δ+2r1r2θ˙1θ˙2cos⁡Δ=r˙12+r12θ˙12+r˙22+r22θ˙22+2cos⁡Δ(r˙1r˙2+r1r2θ˙1θ˙2)+2sin⁡Δ(r1r˙2θ˙1−r˙1r2θ˙2).\begin{aligned} |\vec{v}_{2b}|^2 &= \dot{r}_1^2 \cos^2{\theta_1} + r_1^2 \dot{\theta}_1^2\sin^2{\theta_1} + \dot{r}_2^2\cos^2{\theta_2} + r_2^2\dot{\theta}_2^2\sin^2{\theta_2} -2r_1\dot{r}_1\dot{\theta}_1 \cos{\theta}_1\sin{\theta_1} + 2\dot{r}_1\dot{r}_2\cos{\theta_1}\cos{\theta_2} - 2\dot{r}_1r_2\dot{\theta}_2\cos{\theta_1}\sin{\theta_2} - 2r_1\dot{r}_2 \dot{\theta}_1 \sin{\theta_1}\cos{\theta_2} + 2r_1r_2 \dot{\theta}_1\dot{\theta}_2\sin{\theta_1}\sin{\theta_2}\\ &-2\dot{r}_2r_2\dot{\theta}_2\cos{\theta_2}\sin{\theta_2} + \dot{r}_1^2\sin^2{\theta_1} + r_1^2\dot{\theta}_1^2\cos^2{\theta_1} + \dot{r}_2^2\sin^2{\theta_2} + r_2^2\dot{\theta}_2^2\cos^2{\theta_2} + 2r_1\dot{r}_1\dot{\theta}_1\sin{\theta_1}\cos{\theta_1} + 2\dot{r}_1\dot{r}_2\sin{\theta_1}\sin{\theta_2} + 2\dot{r}_1r_2\dot{\theta}_2\sin{\theta_1}\cos{\theta_2} + 2r_1\dot{r}_2\dot{\theta}_1 \cos{\theta_1}\sin{\theta_2}\\ &+2r_1r_2\dot{\theta}_1\dot{\theta}_2\cos{\theta_1}\cos{\theta_2} + 2r_2\dot{r}_2\dot{\theta}_2\sin{\theta_2}\cos{\theta_2} \\ &= \dot{r}_1^2 + r_1^2 \dot{\theta}_1^2 + \dot{r}_2^2 + r_2^2\dot{\theta}_2^2 + 2r_1\dot{r}_1\dot{\theta}_1(-\cos{\theta_1}\sin{\theta_1} + \cos{\theta_1}\sin{\theta_1}) + 2\dot{r}_1\dot{r}_2(\cos{\theta_1}\cos{\theta_2}+\sin{\theta_1}\sin{\theta_2})+2\dot{r}_1r_2\dot{\theta}_2(-\cos{\theta_1}\sin{\theta_2} + \sin{\theta_1}\cos{\theta_2}) \\ &+ 2r_1\dot{r}_2\dot{\theta}_1(-\sin{\theta_1}\cos{\theta_2}+\cos{\theta_1}\sin{\theta_2}) + 2r_1r_2\dot{\theta}_1\dot{\theta}_2(\sin{\theta_1}\sin{\theta_2} + \cos{\theta_1}\cos{\theta_2}) + 2r_2\dot{r}_2\dot{\theta}_2 (-\cos{\theta_2}\sin{\theta_2} + \sin{\theta_2}\cos{\theta_2}) \\ &= \dot{r}_1^2 + r_1^2 \dot{\theta}_1^2 + \dot{r}_2^2 + r_2^2\dot{\theta}_2^2 + 2\dot{r}_1\dot{r}_2 \cos{(\theta_2-\theta_1)} - 2\dot{r}_1r_2\dot{\theta_2}\sin{(\theta_2-\theta_1)} + 2r_1\dot{r}_2 \dot{\theta}_1 \sin{(\theta_2-\theta_1)} + 2r_1r_2\dot{\theta}_1\dot{\theta}_2\cos{(\theta_2-\theta_1)} \\ &= \dot{r}_1^2 + r_1^2 \dot{\theta}_1^2 + \dot{r}_2^2 + r_2^2\dot{\theta}_2^2 + 2\dot{r}_1\dot{r}_2 \cos{\Delta} - 2\dot{r}_1r_2\dot{\theta_2}\sin{\Delta} + 2r_1\dot{r}_2 \dot{\theta}_1 \sin{\Delta} + 2r_1r_2\dot{\theta}_1\dot{\theta}_2\cos{\Delta} \\ &= \dot{r}_1^2 + r_1^2 \dot{\theta}_1^2 + \dot{r}_2^2 + r_2^2\dot{\theta}_2^2 + 2\cos{\Delta}(\dot{r}_1\dot{r}_2 + r_1r_2\dot{\theta}_1\dot{\theta}_2) + 2\sin{\Delta}(r_1\dot{r}_2\dot{\theta}_1-\dot{r}_1r_2\dot{\theta}_2). \end{aligned}

Let us define ∣v⃗2bI∣2=r˙22+r22θ˙22|\vec{v}_{2b}^I|^2 = \dot{r}_2^2+r_2^2\dot{\theta}_2^2, as this will simplify our Lagrangian later. As for the remaining terms in ∣v⃗2b∣2|\vec{v}_{2b}|^2, they are contained with ∣Δv⃗2,1∣2|\Delta \vec{v}_{2,1}|^2.

∣Δv⃗2,1∣2=2cos⁡Δ(r˙1r˙2+r1r2θ˙1θ˙2)+2sin⁡Δ(r1r˙2θ˙1−r˙1r2θ˙2).\begin{aligned} |\Delta \vec{v}_{2,1}|^2 &= 2\cos{\Delta}(\dot{r}_1\dot{r}_2 + r_1r_2\dot{\theta}_1\dot{\theta}_2) + 2\sin{\Delta}(r_1\dot{r}_2\dot{\theta}_1-\dot{r}_1r_2\dot{\theta}_2). \end{aligned}

Rod 2

Parameterization

x2r=r1cos⁡θ1+scos⁡θ2x˙2r=r˙1cos⁡θ1−r1θ˙1sin⁡θ1+s˙cos⁡θ2−sθ˙2sin⁡θ2y2r=r1sin⁡θ1+ssin⁡θ2y˙2r=r˙1sin⁡θ1+r1θ˙1cos⁡θ1+s˙sin⁡θ2+sθ˙2cos⁡θ2\begin{aligned} x_{2r} &= r_1\cos{\theta_1} + s\cos{\theta_2} & \dot{x}_{2r} = \dot{r}_1 \cos{\theta_1} - r_1 \dot{\theta}_1\sin{\theta_1} + \dot{s}\cos{\theta_2} - s\dot{\theta}_2 \sin{\theta_2} \\ y_{2r} &= r_1\sin{\theta_1} + s\sin{\theta_2} & \dot{y}_{2r} = \dot{r}_1 \sin{\theta_1} + r_1 \dot{\theta}_1\cos{\theta_1} + \dot{s}\sin{\theta_2} + s\dot{\theta}_2 \cos{\theta_2} \end{aligned}
∣v⃗2r∣2=r˙12+r12θ˙12+s˙2+s2θ˙22+2cos⁡Δ(r˙1s˙+r1sθ˙1θ˙2)+2sin⁡Δ(r1s˙θ˙1−r˙1sθ˙2)=r˙12+r12θ˙12+r˙22s2r22+s2θ˙22+2cos⁡Δ(r˙1r˙2sr2+r1sθ˙1θ˙2)+2sin⁡Δ(r1sr˙2θ˙1r2−r˙1sθ˙2).\begin{aligned} |\vec{v}_{2r}|^2 &= \dot{r}_1^2 + r_1^2 \dot{\theta}_1^2 + \dot{s}^2 + s^2\dot{\theta}_2^2 + 2\cos{\Delta}\left(\dot{r}_1\dot{s} + r_1s\dot{\theta}_1\dot{\theta}_2\right) + 2\sin{\Delta}\left(r_1\dot{s}\dot{\theta}_1-\dot{r}_1s\dot{\theta}_2\right) \\ &= \dot{r}_1^2 + r_1^2 \dot{\theta}_1^2 + \dfrac{\dot{r}_2^2 s^2}{r_2^2} + s^2\dot{\theta}_2^2 + 2\cos{\Delta}\left(\dfrac{\dot{r}_1\dot{r}_2 s}{r_2} + r_1s\dot{\theta}_1\dot{\theta}_2\right) + 2\sin{\Delta}\left(\dfrac{r_1s\dot{r}_2\dot{\theta}_1}{r_2}-\dot{r}_1s\dot{\theta}_2\right). \end{aligned}

Centre of mass

x2r=r1cos⁡θ1+r2cos⁡θ22x˙2r=r1˙cos⁡θ1−r1θ˙1sin⁡θ1+r2˙cos⁡θ22−r2θ˙2sin⁡θ22y2r=r1sin⁡θ1+r2sin⁡θ22y˙2r=r1˙sin⁡θ1+r1θ˙1cos⁡θ1+r2˙sin⁡θ22+r2θ˙2cos⁡θ22\begin{aligned} x_{2r} &= r_1\cos{\theta_1} + \dfrac{r_2\cos{\theta_2}}{2} & \dot{x}_{2r} = \dot{r_1}\cos{\theta_1} - r_1\dot{\theta}_1 \sin{\theta_1} + \dfrac{\dot{r_2}\cos{\theta_2}}{2} - \dfrac{r_2\dot{\theta}_2 \sin{\theta_2}}{2} \\ y_{2r} &= r_1\sin{\theta_1} + \dfrac{r_2\sin{\theta_2}}{2} & \dot{y}_{2r} = \dot{r_1}\sin{\theta_1} + r_1\dot{\theta}_1 \cos{\theta_1} + \dfrac{\dot{r_2}\sin{\theta_2}}{2} + \dfrac{r_2\dot{\theta}_2 \cos{\theta_2}}{2} \\ \end{aligned}
∣v⃗2r∣2=r˙12cos⁡2θ1−2r˙1r1θ˙1sin⁡θ1cos⁡θ1+r12θ˙12sin⁡2θ1+r˙22cos⁡2θ24+r22θ˙22sin⁡2θ24−r2r˙2θ˙2cos⁡θ2sin⁡θ22+r˙1r˙2cos⁡θ1cos⁡θ2−r˙1r2θ˙2cos⁡θ1sin⁡θ2+r1r˙2θ˙1sin⁡θ1cos⁡θ2+r1r2θ˙1θ˙2sin⁡θ1sin⁡θ2+r˙12sin⁡2θ1+r12θ˙12cos⁡2θ1+r˙22sin⁡2θ24+r22θ˙22cos⁡2θ24+2r˙1r1θ˙1sin⁡θ1cos⁡θ1+r˙1r˙2sin⁡θ1sin⁡θ2+r˙1r2θ˙2sin⁡θ1cos⁡θ2+r1r˙2θ˙1cos⁡θ1sin⁡θ2+r1r2θ1˙θ2˙cos⁡θ1cos⁡θ2+r2r˙2θ˙2sin⁡θ2cos⁡θ22=r˙12+r12θ˙12+r˙22+r22θ˙224+r1r˙2θ˙1sin⁡Δ+r˙1r˙2cos⁡Δ+r1r2θ˙1θ˙2cos⁡Δ−r˙1r1θ˙1sin⁡Δ=r˙12+r12θ˙12+r˙22+r22θ˙224+sin⁡Δ(r1r˙2θ˙1−r˙1r1θ˙1)+cos⁡Δ(r˙1r˙2+r1r2θ˙1θ˙2)=∣v⃗1b∣2+∣v⃗2bI∣24+∣Δv⃗2,1∣22\begin{aligned} |\vec{v}_{2r}|^2 &= \dot{r}_1^2\cos^2{\theta_1} - 2\dot{r}_1r_1\dot{\theta}_1 \sin{\theta_1}\cos{\theta_1} + r_1^2\dot{\theta}_1^2\sin^2{\theta_1} + \dfrac{\dot{r}_2^2\cos^2{\theta_2}}{4} + \dfrac{r_2^2\dot{\theta}_2^2\sin^2{\theta_2}}{4} - \dfrac{r_2\dot{r}_2\dot{\theta}_2 \cos{\theta_2}\sin{\theta_2}}{2} + \\ & \dot{r}_1\dot{r}_2\cos{\theta_1}\cos{\theta_2} - \dot{r}_1r_2\dot{\theta}_2 \cos{\theta_1}\sin{\theta_2} + r_1\dot{r}_2 \dot{\theta}_1 \sin{\theta_1}\cos{\theta_2} + r_1r_2 \dot{\theta}_1 \dot{\theta}_2 \sin{\theta_1}\sin{\theta_2} + \dot{r}_1^2\sin^2{\theta_1} + \\ & r_1^2 \dot{\theta}_1^2\cos^2{\theta_1} + \dfrac{\dot{r}_2^2\sin^2{\theta_2}}{4} + \dfrac{r_2^2\dot{\theta}_2^2\cos^2{\theta_2}}{4} + 2\dot{r}_1r_1 \dot{\theta}_1\sin{\theta_1}\cos{\theta_1} + \dot{r}_1\dot{r}_2\sin{\theta_1}\sin{\theta_2} + \dot{r}_1 r_2\dot{\theta}_2 \sin{\theta_1}\cos{\theta_2} + \\ & r_1\dot{r}_2 \dot{\theta}_1 \cos{\theta_1}\sin{\theta_2} + r_1r_2\dot{\theta_1}\dot{\theta_2}\cos{\theta_1}\cos{\theta_2} + \dfrac{r_2\dot{r}_2\dot{\theta}_2\sin{\theta_2}\cos{\theta_2}}{2}\\ &= \dot{r}_1^2 + r_1^2\dot{\theta}_1^2 + \dfrac{\dot{r}_2^2 + r_2^2\dot{\theta}_2^2}{4} + r_1\dot{r}_2\dot{\theta}_1\sin{\Delta} + \dot{r}_1\dot{r}_2 \cos{\Delta} + r_1r_2 \dot{\theta}_1\dot{\theta}_2 \cos{\Delta} - \dot{r}_1r_1 \dot{\theta}_1\sin{\Delta} \\ &= \dot{r}_1^2 + r_1^2\dot{\theta}_1^2 + \dfrac{\dot{r}_2^2 + r_2^2\dot{\theta}_2^2}{4} + \sin{\Delta}(r_1\dot{r}_2\dot{\theta}_1 - \dot{r}_1r_1 \dot{\theta}_1) + \cos{\Delta}(\dot{r}_1\dot{r}_2 + r_1r_2 \dot{\theta}_1\dot{\theta}_2) \\ &= |\vec{v}_{1b}|^2 + \dfrac{|\vec{v}_{2b}^I|^2}{4} + \dfrac{|\Delta \vec{v}_{2,1}|^2}{2} \end{aligned}

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

F⃗D,j=−(bj+cj∣v⃗j∣)v⃗j.\begin{aligned} \vec{F}_{D,j} &= -(b_j+c_j|\vec{v}_j|)\vec{v}_j. \end{aligned}

Where jj is the pendulum bob of interest, bjb_j and cjc_j are constants.

Kinetic energy

The kinetic energy of the system is given by

T=m1b2∣v⃗1b∣2+∫0r1m1r2r1∣v⃗1r∣2ds+m2b2∣v⃗2b∣2+∫0r2m2r2r2∣v⃗2r∣2ds\begin{aligned} T &= \dfrac{m_{1b}}{2}|\vec{v}_{1b}|^2 + \int_0^{r_1}\dfrac{m_{1r}}{2r_1} |\vec{v}_{1r}|^2 ds + \dfrac{m_{2b}}{2}|\vec{v}_{2b}|^2 + \int_0^{r_2} \dfrac{m_{2r}}{2r_2} |\vec{v}_{2r}|^2 ds \end{aligned}

Rod 1

∫0r1m1r2r1∣v⃗1r∣2ds=∫0r1m1r2r1s2[r˙12r12+θ˙12]ds=m1r2r1[r˙12r12+θ˙12][s33]0r1=m1r2r1[r˙12r12+θ˙12]r133=m1r6[r˙12+r12θ˙12]=m1r6∣v⃗1b∣2.\begin{aligned} \int_0^{r_1}\dfrac{m_{1r}}{2r_1} |\vec{v}_{1r}|^2 ds &= \int_0^{r_1} \dfrac{m_{1r}}{2r_1} s^2 \left[\dfrac{\dot{r}_1^2}{r_1^2}+\dot{\theta}_1^2\right] ds \\ &= \dfrac{m_{1r}}{2r_1} \left[\dfrac{\dot{r}_1^2}{r_1^2}+\dot{\theta}_1^2\right] \left[\dfrac{s^3}{3}\right]_0^{r_1} \\ &= \dfrac{m_{1r}}{2r_1} \left[\dfrac{\dot{r}_1^2}{r_1^2}+\dot{\theta}_1^2\right] \dfrac{r_1^3}{3} \\ &= \dfrac{m_{1r}}{6}\left[\dot{r}_1^2+r_1^2\dot{\theta}_1^2\right] \\ &= \dfrac{m_{1r}}{6} |\vec{v}_{1b}|^2. \end{aligned}

Rod 2

∫0r2m2r2r2∣v⃗2r∣2ds=∫0r2m2r2r2[r˙12+r12θ˙12+r˙22s2r22+s2θ˙22+2cos⁡Δ(r˙1r˙2sr2+r1sθ˙1θ˙2)+2sin⁡Δ(r1sr˙2θ˙1r2−r˙1sθ˙2)]ds=m2r2∣v⃗1b∣2+m2r2r2∫0r2s2(r˙22r22+θ˙22)+2s[cos⁡Δ(r˙1r˙2r2+r1θ˙1θ˙2)+sin⁡Δ(r1r˙2θ˙1r2−r˙1θ˙2)]ds=m2r2∣v⃗1b∣2+m2r2r2[r233(r˙22r22+θ˙22)+r22[cos⁡Δ(r˙1r˙2r2+r1θ˙1θ˙2)+sin⁡Δ(r1r˙2θ˙1r2−r˙1θ˙2)]]=m2r2∣v⃗1b∣2+m2r6∣v⃗2bI∣2+m2r4∣Δv⃗2,1I∣2\begin{aligned} \int_0^{r_2} \dfrac{m_{2r}}{2r_2} |\vec{v}_{2r}|^2 ds &= \int_0^{r_2} \dfrac{m_{2r}}{2r_2} \left[\dot{r}_1^2 + r_1^2 \dot{\theta}_1^2 + \dfrac{\dot{r}_2^2 s^2}{r_2^2} + s^2\dot{\theta}_2^2 + 2\cos{\Delta}\left(\dfrac{\dot{r}_1\dot{r}_2 s}{r_2} + r_1s\dot{\theta}_1\dot{\theta}_2\right) + 2\sin{\Delta}\left(\dfrac{r_1s\dot{r}_2\dot{\theta}_1}{r_2}-\dot{r}_1s\dot{\theta}_2\right)\right] ds \\ &= \dfrac{m_{2r}}{2} |\vec{v}_{1b}|^2 + \dfrac{m_{2r}}{2r_2}\int_0^{r_2} s^2\left(\dfrac{\dot{r}_2^2}{r_2^2} + \dot{\theta}_2^2\right) + 2s\left[\cos{\Delta}\left(\dfrac{\dot{r}_1\dot{r}_2 }{r_2} + r_1\dot{\theta}_1\dot{\theta}_2\right) + \sin{\Delta}\left(\dfrac{r_1\dot{r}_2\dot{\theta}_1}{r_2}-\dot{r}_1\dot{\theta}_2\right)\right]ds \\ &= \dfrac{m_{2r}}{2} |\vec{v}_{1b}|^2 + \dfrac{m_{2r}}{2r_2}\left[\dfrac{r_2^3}{3}\left(\dfrac{\dot{r}_2^2}{r_2^2} + \dot{\theta}_2^2\right) + r_2^2\left[\cos{\Delta}\left(\dfrac{\dot{r}_1\dot{r}_2 }{r_2} + r_1\dot{\theta}_1\dot{\theta}_2\right) + \sin{\Delta}\left(\dfrac{r_1\dot{r}_2\dot{\theta}_1}{r_2}-\dot{r}_1\dot{\theta}_2\right)\right]\right] \\ &= \dfrac{m_{2r}}{2} |\vec{v}_{1b}|^2 + \dfrac{m_{2r}}{6} |\vec{v}_{2b}^{I}|^2 + \dfrac{m_{2r}}{4}|\Delta \vec{v}_{2,1}^I|^2 \end{aligned}

Final result

T=m1b2∣v⃗1b∣2+m1r6∣v⃗1b∣2+m2b2(∣v⃗1b∣2+∣v⃗2bI∣2+∣Δv⃗2,1I∣2)+m2r2∣v⃗1b∣2+m2r6∣v⃗2bI∣2+m2r4∣Δv⃗2,1I∣2=m1b+m1r3+m2b+m2r2∣v⃗1b∣2+m2b+m2r32∣v⃗2bI∣2+m2b+m2r22∣Δv⃗2,1I∣2.\begin{aligned} T &= \dfrac{m_{1b}}{2}|\vec{v}_{1b}|^2 + \dfrac{m_{1r}}{6} |\vec{v}_{1b}|^2 + \dfrac{m_{2b}}{2}(|\vec{v}_{1b}|^2+|\vec{v}_{2b}^I|^2+|\Delta \vec{v}_{2,1}^I|^2) + \dfrac{m_{2r}}{2} |\vec{v}_{1b}|^2 + \dfrac{m_{2r}}{6} |\vec{v}_{2b}^{I}|^2 + \dfrac{m_{2r}}{4}|\Delta \vec{v}_{2,1}^I|^2 \\ &= \dfrac{m_{1b}+\dfrac{m_{1r}}{3}+m_{2b}+m_{2r}}{2} |\vec{v}_{1b}|^2 + \dfrac{m_{2b} + \dfrac{m_{2r}}{3}}{2}|\vec{v}_{2b}^I|^2 + \dfrac{m_{2b}+\dfrac{m_{2r}}{2}}{2}|\Delta \vec{v}_{2,1}^I|^2. \end{aligned}

Let Mj=∑i=jb2rmi(1−2δi,jr3)M_j = \displaystyle \sum_{i=jb}^{2r} m_i \left(1-\dfrac{2\delta_{i,jr}}{3}\right) and μj=∑i=jb2rmi(1−δi,jr2)\mu_j = \displaystyle \sum_{i=jb}^{2r} m_i \left(1-\dfrac{\delta_{i,jr}}{2}\right), where δi,j\delta_{i,j} is the Kronecker delta symbol. For instance, M1=m1b+m1r3+m2b+m2rM_1 = m_{1b} + \dfrac{m_{1r}}{3} + m_{2b} + m_{2r} and μ2=m2b+m2r2\mu_2 = m_{2b} + \dfrac{m_{2r}}{2}. Then

T=M12∣v⃗1b∣2+M22∣v⃗2bI∣2+μ22∣Δv⃗2,1I∣2.\begin{aligned} T &= \dfrac{M_1}{2} |\vec{v}_{1b}|^2 + \dfrac{M_2}{2}|\vec{v}_{2b}^I|^2 + \dfrac{\mu_2}{2}|\Delta \vec{v}_{2,1}^I|^2. \end{aligned}

Potential energy

The potential energy of the system is given by

V=m1bgy1b+m1rgy1r+m2bgy2b+m2rgy2r+k1(r1−l1)2+k2(r2−l2)22=m1bgr1sin⁡θ1+m1rgr1sin⁡θ12+m2rg(r1sin⁡θ1+r2sin⁡θ22)+m2bg(r1sin⁡θ1+r2sin⁡θ2)+k1(r1−l1)2+k2(r2−l2)22=(m1b+m1r2+m2b+m2r)gr1sin⁡θ1+(m2b+m2r2)gr2sin⁡θ2+k1(r1−l1)2+k2(r2−l2)22=μ1gr1sin⁡θ1+μ2gr2sin⁡θ2+k1(r1−l1)2+k2(r2−l2)22\begin{aligned} V &= m_{1b} gy_{1b} + m_{1r}gy_{1r} + m_{2b} gy_{2b} + m_{2r}gy_{2r} + \dfrac{k_1(r_1-l_1)^2+k_2(r_2-l_2)^2}{2}\\ &= m_{1b} gr_1\sin{\theta_1} + \dfrac{m_{1r} gr_1\sin{\theta_1}}{2} + m_{2r}g\left(r_1\sin{\theta_1} + \dfrac{r_2\sin{\theta_2}}{2}\right) + m_{2b}g(r_1\sin{\theta_1} + r_2\sin{\theta_2}) + \dfrac{k_1(r_1-l_1)^2+k_2(r_2-l_2)^2}{2}\\ &= \left(m_{1b} + \dfrac{m_{1r}}{2} + m_{2b} + m_{2r}\right)gr_1\sin{\theta_1} + \left(m_{2b}+\dfrac{m_{2r}}{2}\right)gr_2\sin{\theta_2} + \dfrac{k_1(r_1-l_1)^2+k_2(r_2-l_2)^2}{2} \\ &= \mu_1 gr_1\sin{\theta_1} + \mu_2 gr_2\sin{\theta_2} + \dfrac{k_1(r_1-l_1)^2+k_2(r_2-l_2)^2}{2} \end{aligned}

Lagrangian

Hence the Lagrangian of the system is

L=T−V=M12∣v⃗1b∣2+M22∣v⃗2bI∣2+μ22∣Δv⃗2,1I∣2−(μ1gr1sin⁡θ1+μ2gr2sin⁡θ2+k1(r1−l1)2+k2(r2−l2)22)=M12∣v⃗1b∣2+M22∣v⃗2bI∣2+μ22(∣Δv⃗2,1I∣2−2gr2sin⁡θ2)−μ1gr1sin⁡θ1−k1(r1−l1)2+k2(r2−l2)22.\begin{aligned} \mathcal{L} &= T - V \\ &= \dfrac{M_1}{2} |\vec{v}_{1b}|^2 + \dfrac{M_2}{2}|\vec{v}_{2b}^I|^2 + \dfrac{\mu_2}{2}|\Delta \vec{v}_{2,1}^I|^2 - \left(\mu_1 gr_1\sin{\theta_1} + \mu_2 gr_2\sin{\theta_2} + \dfrac{k_1(r_1-l_1)^2+k_2(r_2-l_2)^2}{2}\right) \\ &= \dfrac{M_1}{2} |\vec{v}_{1b}|^2 + \dfrac{M_2}{2}|\vec{v}_{2b}^I|^2 + \dfrac{\mu_2}{2}(|\Delta \vec{v}_{2,1}^I|^2 - 2gr_2 \sin{\theta_2}) - \mu_1 gr_1 \sin{\theta_1} - \dfrac{k_1(r_1-l_1)^2+k_2(r_2-l_2)^2}{2}. \end{aligned}

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:

∂∣v⃗1b∣2∂r1=2r1θ˙12∂∣v⃗1b∣2∂r2=0∂∣v⃗1b∣2∂θ1=0∂∣v⃗1b∣2∂θ2=0∂∣v⃗1b∣2∂r˙1=2r˙1∂∣v⃗1b∣2∂r˙2=0∂∣v⃗1b∣2∂θ˙1=2r12θ˙1∂∣v⃗1b∣2∂θ˙2=0ddt∂∣v⃗1b∣2∂r˙1=2r¨1ddt∂∣v⃗1b∣2∂r˙2=0ddt∂∣v⃗1b∣2∂θ˙1=2r12θ¨1+4r1r˙1θ˙1ddt∂∣v⃗1b∣2∂θ˙2=0\begin{aligned} \dfrac{\partial |\vec{v}_{1b}|^2}{\partial r_1} &= 2r_1\dot{\theta}_1^2 & \dfrac{\partial |\vec{v}_{1b}|^2}{\partial r_2} &= 0 & \dfrac{\partial |\vec{v}_{1b}|^2}{\partial \theta_1} &= 0 & \dfrac{\partial |\vec{v}_{1b}|^2}{\partial \theta_2} &= 0\\ \dfrac{\partial |\vec{v}_{1b}|^2}{\partial \dot{r}_1} &= 2\dot{r}_1 & \dfrac{\partial |\vec{v}_{1b}|^2}{\partial \dot{r}_2} &= 0 & \dfrac{\partial |\vec{v}_{1b}|^2}{\partial \dot{\theta}_1} &= 2r_1^2\dot{\theta}_1 & \dfrac{\partial |\vec{v}_{1b}|^2}{\partial \dot{\theta}_2} &= 0\\ \dfrac{d}{dt} \dfrac{\partial |\vec{v}_{1b}|^2}{\partial \dot{r}_1} &= 2\ddot{r}_1 & \dfrac{d}{dt}\dfrac{\partial |\vec{v}_{1b}|^2}{\partial \dot{r}_2} &= 0 & \dfrac{d}{dt}\dfrac{\partial |\vec{v}_{1b}|^2}{\partial \dot{\theta}_1} &= 2r_1^2 \ddot{\theta}_1 + 4r_1\dot{r}_1\dot{\theta}_1 & \dfrac{d}{dt}\dfrac{\partial |\vec{v}_{1b}|^2}{\partial \dot{\theta}_2} &= 0 \end{aligned}

Hence the negative functional derivatives are

δ′∣v⃗1b∣2δ′r1=ddt∂∣v⃗1b∣2∂r˙1−∂∣v⃗1b∣2∂r1δ′∣v⃗1b∣2δ′r2=ddt∂∣v⃗1b∣2∂r˙2−∂∣v⃗1b∣2∂r2δ′∣v⃗1b∣2δ′θ1=ddt∂∣v⃗1b∣2∂θ˙1−∂∣v⃗1b∣2∂θ1δ′∣v⃗1b∣2δ′θ2=ddt∂∣v⃗1b∣2∂θ˙2−∂∣v⃗1b∣2∂θ2=2r¨1−2r1θ˙12=0=2r12θ¨1+4r1r˙1θ˙1=0.\begin{aligned} \dfrac{\delta' |\vec{v}_{1b}|^2}{\delta' r_1} &= \dfrac{d}{dt}\dfrac{\partial |\vec{v}_{1b}|^2}{\partial \dot{r}_1} - \dfrac{\partial |\vec{v}_{1b}|^2}{\partial r_1} & \dfrac{\delta' |\vec{v}_{1b}|^2}{\delta' r_2} &= \dfrac{d}{dt}\dfrac{\partial |\vec{v}_{1b}|^2}{\partial \dot{r}_2} - \dfrac{\partial |\vec{v}_{1b}|^2}{\partial r_2} & \dfrac{\delta' |\vec{v}_{1b}|^2}{\delta' \theta_1} &= \dfrac{d}{dt}\dfrac{\partial |\vec{v}_{1b}|^2}{\partial \dot{\theta}_1} - \dfrac{\partial |\vec{v}_{1b}|^2}{\partial \theta_1} & \dfrac{\delta' |\vec{v}_{1b}|^2}{\delta' \theta_2} &= \dfrac{d}{dt}\dfrac{\partial |\vec{v}_{1b}|^2}{\partial \dot{\theta}_2} - \dfrac{\partial |\vec{v}_{1b}|^2}{\partial \theta_2}\\ &= 2\ddot{r}_1 - 2r_1\dot{\theta}_1^2 & &= 0 & &=2r_1^2\ddot{\theta}_1 + 4r_1\dot{r}_1\dot{\theta}_1 & &= 0. \end{aligned}

Second bob's independent velocity squared

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

∂∣v⃗2bI∣2∂r1=0∂∣v⃗2bI∣2∂r2=2r2θ˙22∂∣v⃗2bI∣2∂θ1=0∂∣v⃗2bI∣2∂θ2=0∂∣v⃗2bI∣2∂r˙1=0∂∣v⃗2bI∣2∂r˙2=2r˙2∂∣v⃗2bI∣2∂θ˙1=0∂∣v⃗2bI∣2∂θ˙2=2r22θ˙2\begin{aligned} \dfrac{\partial |\vec{v}_{2b}^I|^2}{\partial r_1} &= 0 & \dfrac{\partial |\vec{v}_{2b}^I|^2}{\partial r_2} &= 2r_2\dot{\theta}_2^2 \\ \dfrac{\partial |\vec{v}_{2b}^I|^2}{\partial \theta_1} &= 0 & \dfrac{\partial |\vec{v}_{2b}^I|^2}{\partial \theta_2} &= 0 \\ \dfrac{\partial |\vec{v}_{2b}^I|^2}{\partial \dot{r}_1} &= 0 & \dfrac{\partial |\vec{v}_{2b}^I|^2}{\partial \dot{r}_2} &= 2\dot{r}_2 \\ \dfrac{\partial |\vec{v}_{2b}^I|^2}{\partial \dot{\theta}_1} &= 0 & \dfrac{\partial |\vec{v}_{2b}^I|^2}{\partial \dot{\theta}_2} &=2r_2^2\dot{\theta}_2 \end{aligned}
ddt∂∣v⃗2bI∣2∂r˙1=0ddt∂∣v⃗2bI∣2∂r˙2=2r¨2ddt∂∣v⃗2bI∣2∂θ˙1=0ddt∂∣v⃗2bI∣2∂θ˙2=2r22θ¨2+4r2r˙2θ˙2.\begin{aligned} \dfrac{d}{dt} \dfrac{\partial |\vec{v}_{2b}^I|^2}{\partial \dot{r}_1} &= 0 & \dfrac{d}{dt} \dfrac{\partial |\vec{v}_{2b}^I|^2}{\partial \dot{r}_2} &= 2\ddot{r}_2 \\ \dfrac{d}{dt} \dfrac{\partial |\vec{v}_{2b}^I|^2}{\partial \dot{\theta}_1} &= 0 & \dfrac{d}{dt} \dfrac{\partial |\vec{v}_{2b}^I|^2}{\partial \dot{\theta}_2} &= 2r_2^2\ddot{\theta}_2 + 4r_2\dot{r}_2\dot{\theta}_2. \end{aligned}

Hence

δr1∣v⃗2bI∣2=0δr2∣v⃗2bI∣2=2(r¨2−r2θ˙22)δθ1∣v⃗2bI∣2=0δθ2∣v⃗2bI∣2=2r22θ¨2+4r2r˙2θ˙2.\begin{aligned} \delta_{r_1} |\vec{v}_{2b}^I|^2 &= 0 & \delta_{r_2} |\vec{v}_{2b}^I|^2 &= 2(\ddot{r}_2 - r_2\dot{\theta}_2^2) \\ \delta_{\theta_1} |\vec{v}_{2b}^I|^2 &=0 & \delta_{\theta_2} |\vec{v}_{2b}^I|^2 &= 2r_2^2 \ddot{\theta}_2 + 4r_2\dot{r}_2\dot{\theta}_2. \end{aligned}

Difference in independent velocities squared

∂∣Δv⃗2,1I∣2∂r1=2r2θ˙1θ˙2cos⁡Δ+2r˙2θ˙1sin⁡Δ∂∣Δv⃗2,1I∣2∂r2=2r1θ˙1θ˙2cos⁡Δ−2r˙1θ˙2sin⁡Δ∂∣Δv⃗2,1I∣2∂θ1=−2sin⁡Δ⋅−1(r˙1r˙2+r1r2θ˙1θ˙2)+2cos⁡Δ⋅−1(r1r˙2θ˙1−r˙1r2θ˙2)∂∣Δv⃗2,1I∣2∂θ2=−2sin⁡Δ(r˙1r˙2+r1r2θ˙1θ˙2)+2cos⁡Δ(r1r˙2θ˙1−r˙1r2θ˙2)=2sin⁡Δ(r˙1r˙2+r1r2θ˙1θ˙2)−2cos⁡Δ(r1r˙2θ˙1−r˙1r2θ˙2)∂∣Δv⃗2,1∣2∂r˙1=2r˙2cos⁡Δ−2r2θ˙2sin⁡Δ∂∣Δv⃗2,1I∣2∂r˙2=2r˙1cos⁡Δ+2r1θ˙1sin⁡Δ∂∣Δv⃗2,1I∣2∂θ˙1=2r1r2θ˙2cos⁡Δ+2r1r˙2sin⁡Δ∂∣Δv⃗2,1I∣2∂θ˙2=2r1r2θ˙1cos⁡Δ−2r˙1r2sin⁡Δ\begin{aligned} \dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial r_1} &= 2r_2\dot{\theta}_1\dot{\theta}_2\cos{\Delta} + 2\dot{r}_2\dot{\theta}_1\sin{\Delta} & \dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial r_2} &= 2r_1\dot{\theta}_1\dot{\theta}_2\cos{\Delta}-2\dot{r}_1\dot{\theta}_2\sin{\Delta}\\ \dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \theta_1} &= -2\sin{\Delta}\cdot -1(\dot{r}_1\dot{r}_2 + r_1r_2\dot{\theta}_1\dot{\theta}_2) + 2\cos{\Delta}\cdot -1(r_1\dot{r}_2\dot{\theta}_1-\dot{r}_1r_2\dot{\theta}_2) & \dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \theta_2} &= -2\sin{\Delta}(\dot{r}_1\dot{r}_2 + r_1r_2\dot{\theta}_1\dot{\theta}_2) + 2\cos{\Delta}(r_1\dot{r}_2\dot{\theta}_1-\dot{r}_1r_2\dot{\theta}_2) \\ &= 2\sin{\Delta}(\dot{r}_1\dot{r}_2 + r_1r_2\dot{\theta}_1\dot{\theta}_2) - 2\cos{\Delta}(r_1\dot{r}_2\dot{\theta}_1-\dot{r}_1r_2\dot{\theta}_2)\\ \dfrac{\partial |\Delta \vec{v}_{2,1}|^2}{\partial \dot{r}_1} &= 2\dot{r}_2\cos{\Delta} - 2r_2\dot{\theta}_2\sin{\Delta} & \dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \dot{r}_2} &= 2\dot{r}_1\cos{\Delta} + 2r_1\dot{\theta}_1\sin{\Delta} \\ \dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \dot{\theta}_1} &= 2r_1r_2\dot{\theta}_2\cos{\Delta} + 2r_1\dot{r}_2\sin{\Delta} & \dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \dot{\theta}_2} &=2r_1r_2\dot{\theta}_1\cos{\Delta} - 2\dot{r}_1r_2\sin{\Delta} \end{aligned}

Let us define Δ˙1=2θ˙1−θ˙2\dot{\Delta}_1 = 2\dot{\theta}_1 - \dot{\theta}_2 and Δ˙2=2θ˙2−θ˙1\dot{\Delta}_2 = 2\dot{\theta}_2 - \dot{\theta}_1.

ddt∂∣Δv⃗2,1I∣2∂r˙1=2r¨2cos⁡Δ−2r˙2Δ˙sin⁡Δ−2r˙2θ˙2sin⁡Δ−2r2θ¨2sin⁡Δ−2r2θ˙2Δ˙cos⁡Δ=2cos⁡Δ(r¨2−r2θ˙2Δ˙)−2sin⁡Δ(r˙2Δ˙+r˙2θ˙2+r2θ¨2)=2cos⁡Δ(r¨2−r2θ˙2Δ˙)−2sin⁡Δ(r˙2(2θ˙2−θ˙1)+r2θ¨2)=2cos⁡Δ(r¨2−r2θ˙2Δ˙)−2sin⁡Δ(r˙2Δ˙2+r2θ¨2)ddt∂∣Δv⃗2,1I∣2∂r˙2=2r¨1cos⁡Δ−2r˙1Δ˙sin⁡Δ+2r˙1θ˙1sin⁡Δ+2r1θ¨1sin⁡Δ+2r1θ˙1Δ˙cos⁡Δ=2cos⁡Δ(r¨1+r1θ˙1Δ˙)+2sin⁡Δ(r1θ¨1+r˙1θ˙1−r˙1Δ˙)=2cos⁡Δ(r¨1+r1θ˙1Δ˙)+2sin⁡Δ(r1θ¨1+r˙1(2θ˙1−θ¨2))=2cos⁡Δ(r¨1+r1θ˙1Δ˙)+2sin⁡Δ(r1θ¨1+r˙1Δ˙1)ddt∂∣Δv⃗2,1I∣2∂θ˙1=2r˙1r2θ˙2cos⁡Δ+2r1r˙2θ˙2cos⁡Δ+2r1r2θ¨2cos⁡Δ−2r1r2θ˙2Δ˙sin⁡Δ+2r˙1r˙2sin⁡Δ+2r1r¨2sin⁡Δ+2r1r˙2Δ˙cos⁡Δ=2cos⁡Δ(r˙1r2θ˙2+r1r˙2(θ˙2+Δ˙)+r1r2θ¨2)+2sin⁡Δ(r˙1r˙2+r1r¨2−r1r2θ˙2Δ˙)=2cos⁡Δ(r˙1r2θ˙2+r1r˙2Δ˙2+r1r2θ¨2)+2sin⁡Δ(r˙1r˙2+r1r¨2−r1r2θ˙2Δ˙)ddt∂∣Δv⃗2,1I∣2∂θ˙2=2r˙1r2θ˙1cos⁡Δ+2r1r˙2θ˙1cos⁡Δ+2r1r2θ¨1cos⁡Δ−2r1r2θ˙1Δ˙sin⁡Δ−2r¨1r2sin⁡Δ−2r˙1r˙2sin⁡Δ−2r˙1r2Δ˙cos⁡Δ=2cos⁡Δ(r˙1r2(θ˙1−Δ˙)+r1r˙2θ˙1+r1r2θ¨1)−2sin⁡Δ(r1r2θ˙1Δ˙+r¨1r2+r˙1r˙2)=2cos⁡Δ(r˙1r2Δ˙1+r1r˙2θ˙1+r1r2θ¨1)−2sin⁡Δ(r1r2θ˙1Δ˙+r¨1r2+r˙1r˙2).\begin{aligned} \dfrac{d}{dt} \dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \dot{r}_1} &= 2\ddot{r}_2\cos{\Delta} - 2\dot{r}_2\dot{\Delta}\sin{\Delta} - 2\dot{r}_2\dot{\theta}_2\sin{\Delta} - 2r_2\ddot{\theta}_2\sin{\Delta} - 2r_2\dot{\theta}_2\dot{\Delta}\cos{\Delta} \\ &= 2\cos{\Delta}(\ddot{r}_2-r_2\dot{\theta}_2\dot{\Delta}) - 2\sin{\Delta}(\dot{r}_2\dot{\Delta} + \dot{r}_2\dot{\theta}_2+r_2\ddot{\theta}_2)\\ &= 2\cos{\Delta}(\ddot{r}_2-r_2\dot{\theta}_2\dot{\Delta}) - 2\sin{\Delta}(\dot{r}_2(2\dot{\theta}_2-\dot{\theta}_1)+r_2\ddot{\theta}_2)\\ &= 2\cos{\Delta}(\ddot{r}_2-r_2\dot{\theta}_2\dot{\Delta}) - 2\sin{\Delta}(\dot{r}_2\dot{\Delta}_2+r_2\ddot{\theta}_2)\\ \dfrac{d}{dt}\dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \dot{r}_2} &= 2\ddot{r}_1 \cos{\Delta} -2\dot{r}_1\dot{\Delta}\sin{\Delta} + 2\dot{r}_1\dot{\theta}_1\sin{\Delta} + 2r_1\ddot{\theta}_1\sin{\Delta} + 2r_1\dot{\theta}_1\dot{\Delta}\cos{\Delta} \\ &= 2\cos{\Delta}(\ddot{r}_1 + r_1\dot{\theta}_1 \dot{\Delta}) + 2\sin{\Delta}(r_1\ddot{\theta}_1 + \dot{r}_1\dot{\theta}_1 - \dot{r}_1\dot{\Delta})\\ &= 2\cos{\Delta}(\ddot{r}_1 + r_1\dot{\theta}_1 \dot{\Delta}) + 2\sin{\Delta}(r_1\ddot{\theta}_1 + \dot{r}_1(2\dot{\theta}_1 - \ddot{\theta}_2))\\ &= 2\cos{\Delta}(\ddot{r}_1 + r_1\dot{\theta}_1 \dot{\Delta}) + 2\sin{\Delta}(r_1\ddot{\theta}_1 + \dot{r}_1\dot{\Delta}_1)\\ \dfrac{d}{dt}\dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \dot{\theta}_1} &= 2\dot{r}_1r_2\dot{\theta}_2\cos{\Delta} + 2r_1\dot{r}_2\dot{\theta}_2\cos{\Delta} + 2r_1r_2\ddot{\theta}_2\cos{\Delta} - 2r_1r_2\dot{\theta}_2\dot{\Delta}\sin{\Delta} + 2\dot{r}_1\dot{r}_2\sin{\Delta} + 2r_1\ddot{r}_2\sin{\Delta} + 2r_1\dot{r}_2\dot{\Delta}\cos{\Delta} \\ &= 2\cos{\Delta}(\dot{r}_1r_2\dot{\theta}_2 + r_1\dot{r}_2 (\dot{\theta}_2+\dot{\Delta})+r_1r_2\ddot{\theta}_2) +2\sin{\Delta}(\dot{r}_1\dot{r}_2+r_1\ddot{r}_2-r_1r_2\dot{\theta}_2\dot{\Delta})\\ &= 2\cos{\Delta}(\dot{r}_1r_2\dot{\theta}_2 + r_1\dot{r}_2 \dot{\Delta}_2+r_1r_2\ddot{\theta}_2) +2\sin{\Delta}(\dot{r}_1\dot{r}_2+r_1\ddot{r}_2-r_1r_2\dot{\theta}_2\dot{\Delta})\\ \dfrac{d}{dt}\dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \dot{\theta}_2} &= 2\dot{r}_1r_2\dot{\theta}_1\cos{\Delta} + 2r_1\dot{r}_2\dot{\theta}_1\cos{\Delta} + 2r_1r_2\ddot{\theta}_1\cos{\Delta} - 2r_1r_2\dot{\theta}_1\dot{\Delta}\sin{\Delta} - 2\ddot{r}_1r_2\sin{\Delta} - 2\dot{r}_1\dot{r}_2\sin{\Delta} - 2\dot{r}_1r_2\dot{\Delta}\cos{\Delta} \\ &=2\cos{\Delta}(\dot{r}_1r_2(\dot{\theta}_1-\dot{\Delta})+r_1\dot{r}_2\dot{\theta}_1+r_1r_2\ddot{\theta}_1)-2\sin{\Delta}(r_1r_2\dot{\theta}_1\dot{\Delta} + \ddot{r}_1r_2 + \dot{r}_1\dot{r}_2) \\ &=2\cos{\Delta}(\dot{r}_1r_2\dot{\Delta}_1+r_1\dot{r}_2\dot{\theta}_1+r_1r_2\ddot{\theta}_1)-2\sin{\Delta}(r_1r_2\dot{\theta}_1\dot{\Delta} + \ddot{r}_1r_2 + \dot{r}_1\dot{r}_2). \end{aligned}

Hence the negative functional derivative for r1r_1 is

δr1′∣Δv⃗2,1I∣2=ddt∂∣Δv⃗2,1I∣2∂r˙1−∂∣Δv⃗2,1I∣2∂r1=2cos⁡Δ(r¨2−r2θ˙2Δ˙)−2sin⁡Δ(r˙2Δ˙2+r2θ¨2)−(2r2θ˙1θ˙2cos⁡Δ+2r˙2θ˙1sin⁡Δ)=2cos⁡Δ(r¨2−r2θ˙2(Δ˙+θ˙1))−2sin⁡Δ(r˙2(Δ˙2+θ˙1)+r2θ¨2).\begin{aligned} \delta'_{r_1} |\Delta \vec{v}_{2,1}^I|^2 &= \dfrac{d}{dt}\dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \dot{r}_1} - \dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial r_1} \\ &= 2\cos{\Delta}(\ddot{r}_2-r_2\dot{\theta}_2\dot{\Delta}) - 2\sin{\Delta}(\dot{r}_2\dot{\Delta}_2+r_2\ddot{\theta}_2) - \left(2r_2\dot{\theta}_1\dot{\theta}_2\cos{\Delta} + 2\dot{r}_2\dot{\theta}_1\sin{\Delta}\right) \\ &= 2\cos{\Delta}(\ddot{r}_2-r_2\dot{\theta}_2(\dot{\Delta}+\dot{\theta}_1)) - 2\sin{\Delta} (\dot{r}_2(\dot{\Delta}_2+\dot{\theta}_1)+r_2\ddot{\theta}_2). \end{aligned}

Where Δ˙+θ˙1=θ˙2−θ˙1+θ˙1=θ˙2\dot{\Delta} + \dot{\theta}_1 = \dot{\theta}_2 - \dot{\theta}_1 + \dot{\theta}_1 = \dot{\theta}_2 and Δ˙2+θ˙1=2θ˙2−θ˙1+θ˙1=2θ˙2\dot{\Delta}_2 + \dot{\theta}_1 = 2\dot{\theta}_2 - \dot{\theta}_1 + \dot{\theta}_1 = 2\dot{\theta}_2. (Confirmed with SymPy)

δr1′∣Δv⃗2,1I∣2=2cos⁡Δ(r¨2−r2θ˙22)−2sin⁡Δ(2r˙2θ˙2+r2θ¨2).\begin{aligned} \delta'_{r_1} |\Delta \vec{v}_{2,1}^I|^2 &= 2\cos{\Delta}(\ddot{r}_2-r_2\dot{\theta}_2^2) - 2\sin{\Delta} (2\dot{r}_2\dot{\theta}_2+r_2\ddot{\theta}_2). \end{aligned}

As for r2r_2 (checked)

δr2′∣Δv⃗2,1I∣2=ddt∂∣Δv⃗2,1I∣2∂r˙2−∂∣Δv⃗2,1I∣2∂r2=2cos⁡Δ(r¨1+r1θ˙1Δ˙)+2sin⁡Δ(r1θ¨1+r˙1Δ˙1)−[2r1θ˙1θ˙2cos⁡Δ−2r˙1θ˙2sin⁡Δ]=2cos⁡Δ(r¨1+r1θ˙1(Δ˙−θ˙2))+2sin⁡Δ(r1θ¨1+r˙1(Δ˙1+θ˙2)).\begin{aligned} \delta'_{r_2} |\Delta \vec{v}_{2,1}^I|^2 &= \dfrac{d}{dt}\dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \dot{r}_2} - \dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial r_2} \\ &= 2\cos{\Delta}(\ddot{r}_1 + r_1\dot{\theta}_1 \dot{\Delta}) + 2\sin{\Delta}(r_1\ddot{\theta}_1 + \dot{r}_1\dot{\Delta}_1) - \left[2r_1\dot{\theta}_1\dot{\theta}_2\cos{\Delta}-2\dot{r}_1\dot{\theta}_2\sin{\Delta}\right]\\ &= 2\cos{\Delta}(\ddot{r}_1 + r_1\dot{\theta}_1( \dot{\Delta}-\dot{\theta}_2))+2\sin{\Delta}(r_1\ddot{\theta}_1 + \dot{r}_1(\dot{\Delta}_1+\dot{\theta}_2)). \end{aligned}

Hence Δ˙−θ˙2=θ˙2−θ˙1−θ˙2=−θ˙1\dot{\Delta}-\dot{\theta}_2 = \dot{\theta}_2-\dot{\theta}_1-\dot{\theta}_2 = -\dot{\theta}_1 and Δ˙1+θ˙2=2θ˙1−θ˙2+θ˙2=2θ˙1\dot{\Delta}_1 + \dot{\theta}_2 = 2\dot{\theta}_1 - \dot{\theta}_2 + \dot{\theta}_2 = 2\dot{\theta}_1.

δr2′∣Δv⃗2,1I∣2=2cos⁡Δ(r¨1−r1θ˙12)+2sin⁡Δ(r1θ¨1+2r˙1θ˙1).\begin{aligned} \delta'_{r_2} |\Delta \vec{v}_{2,1}^I|^2 &= 2\cos{\Delta}(\ddot{r}_1 - r_1\dot{\theta}_1^2)+2\sin{\Delta}(r_1\ddot{\theta}_1 + 2\dot{r}_1\dot{\theta}_1). \end{aligned}

As for θ1\theta_1 (confirmed by SymPy)

δθ1′∣Δv⃗2,1I∣2=ddt∂∣Δv⃗2,1I∣2∂θ˙1−∂∣Δv⃗2,1I∣2∂θ1=2cos⁡Δ(r˙1r2θ˙2+r1r˙2Δ˙2+r1r2θ¨2)+2sin⁡Δ(r˙1r˙2+r1r¨2−r1r2θ˙2Δ˙)−[2sin⁡Δ(r˙1r˙2+r1r2θ˙1θ˙2)−2cos⁡Δ(r1r˙2θ˙1−r˙1r2θ˙2)]=2cos⁡Δ(r˙1r2(θ˙2−θ˙2)+r1r˙2(Δ˙2+θ˙1)+r1r2θ¨2)+2sin⁡Δ(r˙1r˙2−r˙1r˙2+r1r¨2−r1r2θ˙2(Δ˙+θ˙1))=2cos⁡Δ(2r1r˙2θ˙2+r1r2θ¨2)+2sin⁡Δ(r1r¨2−r1r2θ˙22).\begin{aligned} \delta'_{\theta_1} |\Delta \vec{v}_{2,1}^I|^2 &= \dfrac{d}{dt}\dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \dot{\theta}_1} - \dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \theta_1} \\ &= 2\cos{\Delta}(\dot{r}_1r_2\dot{\theta}_2 + r_1\dot{r}_2 \dot{\Delta}_2+r_1r_2\ddot{\theta}_2) +2\sin{\Delta}(\dot{r}_1\dot{r}_2+r_1\ddot{r}_2-r_1r_2\dot{\theta}_2\dot{\Delta}) - \left[2\sin{\Delta}(\dot{r}_1\dot{r}_2 + r_1r_2\dot{\theta}_1\dot{\theta}_2) - 2\cos{\Delta}(r_1\dot{r}_2\dot{\theta}_1-\dot{r}_1r_2\dot{\theta}_2)\right] \\ &= 2\cos{\Delta}(\dot{r}_1r_2(\dot{\theta}_2 - \dot{\theta}_2)+ r_1\dot{r}_2 (\dot{\Delta}_2+\dot{\theta}_1)+r_1r_2\ddot{\theta}_2) +2\sin{\Delta}(\dot{r}_1\dot{r}_2-\dot{r}_1\dot{r}_2+r_1\ddot{r}_2-r_1r_2\dot{\theta}_2(\dot{\Delta}+\dot{\theta}_1)) \\ &= 2\cos{\Delta}(2r_1\dot{r}_2\dot{\theta}_2+r_1r_2\ddot{\theta}_2) +2\sin{\Delta}(r_1\ddot{r}_2-r_1r_2\dot{\theta}_2^2). \end{aligned}

As for θ2\theta_2 (correct, checked with SymPy)

δθ2′∣Δv⃗2,1I∣2=ddt∂∣Δv⃗2,1I∣2∂θ˙2−∂∣Δv⃗2,1I∣2∂θ2=2cos⁡Δ(r˙1r2Δ˙1+r1r˙2θ˙1+r1r2θ¨1)−2sin⁡Δ(r1r2θ˙1Δ˙+r¨1r2+r˙1r˙2)−[−2sin⁡Δ(r˙1r˙2+r1r2θ˙1θ˙2)+2cos⁡Δ(r1r˙2θ˙1−r˙1r2θ˙2)]=2cos⁡Δ(r˙1r2(Δ˙1+θ˙2)+r1r˙2(θ˙1−θ˙1)+r1r2θ¨1)+2sin⁡Δ(r1r2θ˙1(θ˙2−Δ˙)−r¨1r2+r˙1r˙2−r˙1r˙2)=2cos⁡Δ(2r˙1r2θ˙1+r1r2θ¨1)+2sin⁡Δ(r1r2θ˙1(θ˙2−(θ˙2−θ˙1))−r¨1r2)=2cos⁡Δ(2r˙1r2θ˙1+r1r2θ¨1)+2sin⁡Δ(r1r2θ˙12−r¨1r2).\begin{aligned} \delta'_{\theta_2} |\Delta \vec{v}_{2,1}^I|^2 &= \dfrac{d}{dt}\dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \dot{\theta}_2} - \dfrac{\partial |\Delta \vec{v}_{2,1}^I|^2}{\partial \theta_2} \\ &= 2\cos{\Delta}(\dot{r}_1r_2\dot{\Delta}_1+r_1\dot{r}_2\dot{\theta}_1+r_1r_2\ddot{\theta}_1)-2\sin{\Delta}(r_1r_2\dot{\theta}_1\dot{\Delta} + \ddot{r}_1r_2 + \dot{r}_1\dot{r}_2) - \left[-2\sin{\Delta}(\dot{r}_1\dot{r}_2 + r_1r_2\dot{\theta}_1\dot{\theta}_2) + 2\cos{\Delta}(r_1\dot{r}_2\dot{\theta}_1-\dot{r}_1r_2\dot{\theta}_2)\right] \\ &= 2\cos{\Delta}(\dot{r}_1r_2(\dot{\Delta}_1+\dot{\theta}_2)+r_1\dot{r}_2(\dot{\theta}_1-\dot{\theta}_1)+r_1r_2\ddot{\theta}_1)+2\sin{\Delta}(r_1r_2\dot{\theta}_1(\dot{\theta}_2-\dot{\Delta}) - \ddot{r}_1r_2 + \dot{r}_1\dot{r}_2-\dot{r}_1\dot{r}_2) \\ &= 2\cos{\Delta}(2\dot{r}_1r_2\dot{\theta}_1+r_1r_2\ddot{\theta}_1)+2\sin{\Delta}(r_1r_2\dot{\theta}_1(\dot{\theta}_2-(\dot{\theta}_2-\dot{\theta}_1))-\ddot{r}_1r_2) \\ &= 2\cos{\Delta}(2\dot{r}_1r_2\dot{\theta}_1+r_1r_2\ddot{\theta}_1)+2\sin{\Delta}(r_1r_2\dot{\theta}_1^2-\ddot{r}_1r_2). \end{aligned}

Euler-Lagrange equations with dissipation

r1r_1

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

δr1′L=M12δr1′∣v⃗1b∣2+M22δr1′∣v⃗2bI∣2+μ22δr1′∣Δv⃗2,1I∣2+μ1gsin⁡θ1+k1(r1−l1)\begin{aligned} \delta'_{r_1} \mathcal{L} &= \dfrac{M_1}{2} \delta'_{r_1} |\vec{v}_{1b}|^2 + \dfrac{M_2}{2} \delta'_{r_1}|\vec{v}_{2b}^I|^2 + \dfrac{\mu_2}{2}\delta'_{r_1} |\Delta \vec{v}_{2,1}^I|^2 + \mu_1 g\sin{\theta_1} + k_1(r_1-l_1) \end{aligned}

We have deliberately ignored the m2gr2sin⁡θ2m_2gr_2\sin{\theta_2} and −k2(r2−l2)22-\dfrac{k_2(r_2-l_2)^2}{2} as they are independent of r1r_1.

δr1′L=M12(2r¨1−2r1θ˙12)+M22⋅0+μ22(2cos⁡Δ(r¨2−r2θ˙2(Δ˙+θ˙1))−2sin⁡Δ(r˙2(Δ˙2+θ˙1)+r2θ¨2))+μ1gsin⁡θ1+k1(r1−l1)=M1(r¨1−r1θ˙12)+μ2(cos⁡Δ(r¨2−r2θ˙2(Δ˙+θ˙1))−sin⁡Δ(r˙2(Δ˙2+θ˙1)+r2θ¨2))+μ1gsin⁡θ1+k1(r1−l1)=M1(r¨1−r1θ˙12)+μ2(cos⁡Δ(r¨2−r2θ˙22)−sin⁡Δ(2r˙2θ˙2+r2θ¨2))+μ1gsin⁡θ1+k1(r1−l1).\begin{aligned} \delta'_{r_1} \mathcal{L} &= \dfrac{M_1}{2}(2\ddot{r}_1-2r_1\dot{\theta}_1^2) + \dfrac{M_2}{2} \cdot 0 + \dfrac{\mu_2}{2} \left(2\cos{\Delta}(\ddot{r}_2-r_2\dot{\theta}_2(\dot{\Delta}+\dot{\theta}_1)) - 2\sin{\Delta} (\dot{r}_2(\dot{\Delta}_2+\dot{\theta}_1)+r_2\ddot{\theta}_2)\right) + \mu_1g\sin{\theta_1} + k_1(r_1-l_1) \\ &= M_1(\ddot{r}_1-r_1\dot{\theta}_1^2) + \mu_2 \left(\cos{\Delta}(\ddot{r}_2-r_2\dot{\theta}_2(\dot{\Delta}+\dot{\theta}_1)) - \sin{\Delta} (\dot{r}_2(\dot{\Delta}_2+\dot{\theta}_1)+r_2\ddot{\theta}_2)\right) + \mu_1g\sin{\theta_1} + k_1(r_1-l_1) \\ &= M_1(\ddot{r}_1-r_1\dot{\theta}_1^2) + \mu_2 \left(\cos{\Delta}(\ddot{r}_2-r_2\dot{\theta}_2^2) - \sin{\Delta} (2\dot{r}_2\dot{\theta}_2+r_2\ddot{\theta}_2)\right) + \mu_1g\sin{\theta_1} + k_1(r_1-l_1). \end{aligned}

The generalized dissipation force canonical to r1r_1 is hence

Qr1=−(b1b+c1b∣v⃗1b∣)v⃗1b⋅∂r⃗1b∂r1−(b1r+c1r∣v⃗1r∣)v⃗1r⋅∂r⃗1r∂r1−(b2b+c2b∣v⃗2b∣)v⃗2b⋅∂r⃗2b∂r1−(b2r+c2r∣v⃗2r∣)v⃗2r⋅∂r⃗2r∂r1=−(b1b+c1b∣v⃗1b∣)[r˙1cos⁡θ1−r1θ˙1sin⁡θ1r˙1sin⁡θ1+r1θ˙1cos⁡θ1]⋅[cos⁡θ1sin⁡θ1]−(b1r+c1r∣v⃗1r∣)12[r˙1cos⁡θ1−r1θ˙1sin⁡θ1r˙1sin⁡θ1+r1θ˙1cos⁡θ1]⋅12[cos⁡θ1sin⁡θ1]−(b2b+c2b∣v⃗2b∣)[r˙1cos⁡θ1−r1θ˙1sin⁡θ1+r˙2cos⁡θ2−r2θ˙2sin⁡θ2r˙1sin⁡θ1+r1θ˙1cos⁡θ1+r˙2sin⁡θ2+r2θ˙2cos⁡θ2]⋅[cos⁡θ1sin⁡θ1]−(b2r+c2r∣v⃗2r∣)[r˙1cos⁡θ1−r1θ˙1sin⁡θ1+r˙2cos⁡θ2−r2θ˙2sin⁡θ22r˙1sin⁡θ1+r1θ˙1cos⁡θ1+r˙2sin⁡θ2+r2θ˙2cos⁡θ22]⋅[cos⁡θ1sin⁡θ1]=−(b1b+c1b∣v⃗1b∣)[r˙1cos⁡2θ1−r1θ˙1sin⁡θ1cos⁡θ1+r˙1sin⁡2θ1+r1θ˙1cos⁡θ1sin⁡θ1]−14(b1r+c1r∣v⃗1r∣)[r˙1cos⁡2θ1−r1θ˙1sin⁡θ1cos⁡θ1+r˙1sin⁡2θ1+r1θ˙1cos⁡θ1sin⁡θ1]−(b2b+c2b∣v⃗2b∣)[r˙1cos⁡2θ1−r1θ˙1sin⁡θ1cos⁡θ1+r˙2cos⁡θ1cos⁡θ2−r2θ˙2cos⁡θ1sin⁡θ2+r˙1sin⁡2θ1+r1θ˙1cos⁡θ1sin⁡θ1+r˙2sin⁡θ1sin⁡θ2+r2θ˙2sin⁡θ1cos⁡θ2]−(b2r+c2r∣v⃗2r∣)[r˙1cos⁡2θ1−r1θ˙1sin⁡θ1cos⁡θ1+12[r˙2cos⁡θ1cos⁡θ2−r2θ˙2cos⁡θ1sin⁡θ2]+r˙1sin⁡2θ1+r1θ˙1cos⁡θ1sin⁡θ1+r˙2sin⁡θ1sin⁡θ2+r2θ˙2sin⁡θ1cos⁡θ22]=−(b1b+c1b∣v⃗1b∣)[r˙1(cos⁡2θ1+sin⁡2θ1)+r1θ˙1(−sin⁡θ1cos⁡θ1+cos⁡θ1sin⁡θ1)]−14(b1r+c1r∣v⃗1r∣)[r˙1(cos⁡2θ1+sin⁡2θ1)+r1θ˙1(−sin⁡θ1cos⁡θ1+cos⁡θ1sin⁡θ1)]−(b2b+c2b∣v⃗2b∣)[r˙1(cos⁡2θ1+sin⁡2θ1)+r1θ˙1(−sin⁡θ1cos⁡θ1+cos⁡θ1sin⁡θ1)+r˙2(cos⁡θ1cos⁡θ2+sin⁡θ1sin⁡θ2)+r2θ˙2(−cos⁡θ1sin⁡θ2+sin⁡θ1cos⁡θ2)]−(b2r+c2r∣v⃗2r∣)[r˙1(cos⁡2θ1+sin⁡2θ1)+r1θ˙1(−sin⁡θ1cos⁡θ1+cos⁡θ1sin⁡θ1)+12[r˙2(cos⁡θ1cos⁡θ2+sin⁡θ1sin⁡θ2)+r2θ˙2(−cos⁡θ1sin⁡θ2+sin⁡θ1cos⁡θ2)]]=−(b1b+c1b∣v⃗1b∣)r˙1−r˙14(b1r+c1r∣v⃗1r∣)−(b2b+c2b∣v⃗2b∣)[r˙1+r˙2cos⁡Δ−r2θ˙2sin⁡Δ]−(b2r+c2r∣v⃗2r∣)[r˙1+r˙2cos⁡Δ−r2θ˙2sin⁡Δ2].\begin{aligned} Q_{r_1} &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)\vec{v}_{1b}\cdot \dfrac{\partial \vec{r}_{1b}}{\partial r_{1}} -(b_{1r}+c_{1r}|\vec{v}_{1r}|)\vec{v}_{1r}\cdot \dfrac{\partial \vec{r}_{1r}}{\partial r_{1}} - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\vec{v}_{2b}\cdot \dfrac{\partial \vec{r}_{2b}}{\partial r_1} - (b_{2r}+c_{2r}|\vec{v}_{2r}|)\vec{v}_{2r}\cdot \dfrac{\partial \vec{r}_{2r}}{\partial r_1} \\ &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)\begin{bmatrix} \dot{r}_1\cos{\theta_1} - r_1\dot{\theta}_1\sin{\theta_1} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1\cos{\theta_1} \end{bmatrix} \cdot \begin{bmatrix} \cos{\theta_1} \\ \sin{\theta_1} \end{bmatrix} -(b_{1r}+c_{1r}|\vec{v}_{1r}|)\dfrac{1}{2}\begin{bmatrix} \dot{r}_1\cos{\theta_1} - r_1\dot{\theta}_1\sin{\theta_1} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1\cos{\theta_1} \end{bmatrix} \cdot \dfrac{1}{2} \begin{bmatrix} \cos{\theta_1} \\ \sin{\theta_1} \end{bmatrix} - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\begin{bmatrix} \dot{r}_1 \cos{\theta_1} - r_1 \dot{\theta}_1 \sin{\theta_1} + \dot{r}_2\cos{\theta_2} - r_2\dot{\theta}_2 \sin{\theta_2} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1 \cos{\theta_1} + \dot{r}_2\sin{\theta_2} + r_2\dot{\theta}_2 \cos{\theta_2} \end{bmatrix} \cdot \begin{bmatrix} \cos{\theta_1} \\ \sin{\theta_1} \end{bmatrix}- (b_{2r}+c_{2r}|\vec{v}_{2r}|)\begin{bmatrix} \dot{r}_1 \cos{\theta_1} - r_1 \dot{\theta}_1 \sin{\theta_1} + \dfrac{\dot{r}_2\cos{\theta_2} - r_2\dot{\theta}_2 \sin{\theta_2}}{2} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1 \cos{\theta_1} + \dfrac{\dot{r}_2\sin{\theta_2} + r_2\dot{\theta}_2 \cos{\theta_2}}{2} \end{bmatrix} \cdot \begin{bmatrix} \cos{\theta_1} \\ \sin{\theta_1} \end{bmatrix} \\ &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)\left[\dot{r}_1\cos^2{\theta_1} - r_1\dot{\theta}_1\sin{\theta_1}\cos{\theta_1} + \dot{r}_1\sin^2{\theta_1} + r_1\dot{\theta}_1\cos{\theta_1}\sin{\theta_1}\right] \\ &- \dfrac{1}{4}(b_{1r}+c_{1r}|\vec{v}_{1r}|)\left[\dot{r}_1\cos^2{\theta_1} - r_1\dot{\theta}_1\sin{\theta_1}\cos{\theta_1} + \dot{r}_1\sin^2{\theta_1} + r_1\dot{\theta}_1\cos{\theta_1}\sin{\theta_1}\right] \\ &- (b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[\dot{r}_1\cos^2{\theta_1} - r_1\dot{\theta}_1\sin{\theta_1}\cos{\theta_1} + \dot{r}_2\cos{\theta_1}\cos{\theta_2}-r_2\dot{\theta}_2\cos{\theta_1}\sin{\theta_2} \right.\\ &\left.+ \dot{r}_1\sin^2{\theta_1} + r_1\dot{\theta}_1\cos{\theta_1}\sin{\theta_1} + \dot{r}_2\sin{\theta_1}\sin{\theta_2} + r_2\dot{\theta}_2\sin{\theta_1}\cos{\theta_2} \right] \\ &- (b_{2r}+c_{2r}|\vec{v}_{2r}|)\left[\dot{r}_1\cos^2{\theta_1} - r_1\dot{\theta}_1\sin{\theta_1}\cos{\theta_1} + \dfrac{1}{2}\left[\dot{r}_2\cos{\theta_1}\cos{\theta_2}-r_2\dot{\theta}_2\cos{\theta_1}\sin{\theta_2}\right]\right.\\ &\left.+ \dot{r}_1\sin^2{\theta_1} + r_1\dot{\theta}_1\cos{\theta_1}\sin{\theta_1} + \dfrac{\dot{r}_2\sin{\theta_1}\sin{\theta_2} + r_2\dot{\theta}_2\sin{\theta_1}\cos{\theta_2}}{2} \right] \\ &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)\left[\dot{r}_1(\cos^2{\theta_1}+\sin^2{\theta_1}) + r_1\dot{\theta}_1(-\sin{\theta_1}\cos{\theta_1} + \cos{\theta_1}\sin{\theta_1})\right] \\ &- \dfrac{1}{4}(b_{1r}+c_{1r}|\vec{v}_{1r}|)\left[\dot{r}_1(\cos^2{\theta_1}+\sin^2{\theta_1}) + r_1\dot{\theta}_1(-\sin{\theta_1}\cos{\theta_1} + \cos{\theta_1}\sin{\theta_1})\right] \\ &- (b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[\dot{r}_1(\cos^2{\theta_1}+\sin^2{\theta_1}) + r_1\dot{\theta}_1(-\sin{\theta_1}\cos{\theta_1} + \cos{\theta_1}\sin{\theta_1}) \right.\\ &\left.+ \dot{r}_2(\cos{\theta_1}\cos{\theta_2}+\sin{\theta_1}\sin{\theta_2})+r_2\dot{\theta}_2(-\cos{\theta_1}\sin{\theta_2} + \sin{\theta_1}\cos{\theta_2}) \right] \\ &- (b_{2r}+c_{2r}|\vec{v}_{2r}|)\left[\dot{r}_1(\cos^2{\theta_1}+\sin^2{\theta_1}) + r_1\dot{\theta}_1(-\sin{\theta_1}\cos{\theta_1} + \cos{\theta_1}\sin{\theta_1})\right.\\ &\left.+ \dfrac{1}{2}\left[\dot{r}_2(\cos{\theta_1}\cos{\theta_2}+\sin{\theta_1}\sin{\theta_2})+r_2\dot{\theta}_2(-\cos{\theta_1}\sin{\theta_2}+\sin{\theta_1}\cos{\theta_2})\right]\right]\\ &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)\dot{r}_1 - \dfrac{\dot{r}_1}{4}(b_{1r}+c_{1r}|\vec{v}_{1r}|) - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[\dot{r}_1 + \dot{r}_2\cos{\Delta}-r_2\dot{\theta}_2\sin{\Delta} \right] \\ &- (b_{2r}+c_{2r}|\vec{v}_{2r}|)\left[\dot{r}_1+ \dfrac{\dot{r}_2\cos{\Delta}-r_2\dot{\theta}_2\sin{\Delta}}{2}\right]. \end{aligned}

Hence the Euler-Lagrange equation for r1r_1 with dissipative forces is

M1(r¨1−r1θ˙12)+μ2(cos⁡Δ(r¨2−r2θ˙22)−sin⁡Δ(2r˙2θ˙2+r2θ¨2))+μ1gsin⁡θ1+k1(r1−l1)=Qr1.\begin{aligned} M_1(\ddot{r}_1-r_1\dot{\theta}_1^2) + \mu_2 \left(\cos{\Delta}(\ddot{r}_2-r_2\dot{\theta}_2^2) - \sin{\Delta} (2\dot{r}_2\dot{\theta}_2+r_2\ddot{\theta}_2)\right) + \mu_1g\sin{\theta_1} + k_1(r_1-l_1) &= Q_{r_1}. \end{aligned}

So, this line of our matrix equation for q¨\mathbf{\ddot{q}} will be

[M1μ2cos⁡Δ0−μ2r2sin⁡Δ][r¨1r¨2θ¨1θ¨2]=[M1r12θ˙12+μ2r2θ˙22cos⁡Δ+2μ2r˙2θ˙2sin⁡Δ−μ1gsin⁡θ1−k1(r1−l1)+Qr1].\begin{aligned} \begin{bmatrix} M_1 & \mu_2 \cos{\Delta} & 0 & -\mu_2r_2 \sin{\Delta} \end{bmatrix} \begin{bmatrix} \ddot{r}_1 \\ \ddot{r}_2 \\ \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{bmatrix} &= \begin{bmatrix} M_1r_1^2 \dot{\theta}_1^2 + \mu_2 r_2 \dot{\theta}_2^2\cos{\Delta} + 2\mu_2 \dot{r}_2 \dot{\theta}_2\sin{\Delta} - \mu_1g\sin{\theta_1} - k_1(r_1-l_1) + Q_{r_1} \end{bmatrix}. \end{aligned}

r2r_2

As for r2r_2

δr2′L=M12δr2′∣v⃗1b∣2+M22δr2′∣v⃗2bI∣2+μ22(δr2′∣Δv⃗2,1I∣2+2gsin⁡θ2)+k2(r2−l2)=M22(2(r¨2−r2θ˙22))+μ22(2cos⁡Δ(r¨1−r1θ˙12)+2sin⁡Δ(r1θ¨1+2r˙1θ˙1)+2gsin⁡θ2)+k2(r2−l2)=M2(r¨2−r2θ˙22)+μ2(cos⁡Δ(r¨1−r1θ˙12)+sin⁡Δ(r1θ¨1+2r˙1θ˙1)+gsin⁡θ2)+k2(r2−l2)=M2(r¨2−r2θ˙22)+μ2(cos⁡Δ(r¨1−r1θ˙12)+sin⁡Δ(r1θ¨1+2r˙1θ˙1)+gsin⁡θ2)+k2(r2−l2).\begin{aligned} \delta'_{r_2} \mathcal{L} &= \dfrac{M_1}{2} \delta'_{r_2} |\vec{v}_{1b}|^2 + \dfrac{M_2}{2} \delta'_{r_2} |\vec{v}_{2b}^I|^2 + \dfrac{\mu_2}{2}(\delta'_{r_2} |\Delta \vec{v}_{2,1}^I|^2+2g\sin{\theta_2}) + k_2(r_2-l_2) \\ &= \dfrac{M_2}{2} (2(\ddot{r}_2 - r_2\dot{\theta}_2^2)) + \dfrac{\mu_2}{2}(2\cos{\Delta}(\ddot{r}_1 - r_1\dot{\theta}_1^2)+2\sin{\Delta}(r_1\ddot{\theta}_1 + 2\dot{r}_1\dot{\theta}_1)+2g\sin{\theta_2}) + k_2(r_2-l_2) \\ &= M_2 (\ddot{r}_2 - r_2\dot{\theta}_2^2) + \mu_2 (\cos{\Delta}(\ddot{r}_1 - r_1\dot{\theta}_1^2)+\sin{\Delta}(r_1\ddot{\theta}_1 + 2\dot{r}_1\dot{\theta}_1)+g\sin{\theta_2}) + k_2(r_2-l_2) \\ &= M_2 (\ddot{r}_2 - r_2\dot{\theta}_2^2) + \mu_2 (\cos{\Delta}(\ddot{r}_1 - r_1\dot{\theta}_1^2)+\sin{\Delta}(r_1\ddot{\theta}_1 + 2\dot{r}_1\dot{\theta}_1)+g\sin{\theta_2}) + k_2(r_2-l_2). \end{aligned}

Generalized dissipative forces

Qr2=−(b1b+c1b∣v⃗1b∣)v⃗1b⋅∂r⃗1b∂r2−(b1r+c1r∣v⃗1r∣)v⃗1r⋅∂r⃗1r∂r2−(b2b+c2b∣v⃗2b∣)v⃗2b⋅∂r⃗2b∂r2−(b2r+c2r∣v⃗2r∣)v⃗2r⋅∂r⃗2r∂r2=−(b1b+c1b∣v⃗1b∣)[r˙1cos⁡θ1−r1θ˙1sin⁡θ1r˙1sin⁡θ1+r1θ˙1cos⁡θ1]⋅0⃗−(b1r+c1r∣v⃗1r∣)12[r˙1cos⁡θ1−r1θ˙1sin⁡θ1r˙1sin⁡θ1+r1θ˙1cos⁡θ1]⋅0⃗−(b2b+c2b∣v⃗2b∣)[r˙1cos⁡θ1−r1θ˙1sin⁡θ1+r˙2cos⁡θ2−r2θ˙2sin⁡θ2r˙1sin⁡θ1+r1θ˙1cos⁡θ1+r˙2sin⁡θ2+r2θ˙2cos⁡θ2]⋅[cos⁡θ2sin⁡θ2]−(b2r+c2r∣v⃗2r∣)[r˙1cos⁡θ1−r1θ˙1sin⁡θ1+r˙2cos⁡θ2−r2θ˙2sin⁡θ22r˙1sin⁡θ1+r1θ˙1cos⁡θ1+r˙2sin⁡θ2+r2θ˙2cos⁡θ22]⋅12[cos⁡θ2sin⁡θ2]=−(b2b+c2b∣v⃗2b∣)[r˙1cos⁡θ1cos⁡θ2−r1θ˙1sin⁡θ1cos⁡θ2+r˙2cos⁡2θ2−r2θ˙2sin⁡θ2cos⁡θ2+r˙1sin⁡θ1sin⁡θ2+r1θ˙1cos⁡θ1sin⁡θ2+r˙2sin⁡2θ2+r2θ˙2cos⁡θ2sin⁡θ2]−12(b2r+c2r∣v⃗2r∣)[r˙1cos⁡θ1cos⁡θ2−r1θ˙1sin⁡θ1cos⁡θ2+r˙2cos⁡2θ2−r2θ˙2sin⁡θ2cos⁡θ22+r˙1sin⁡θ1sin⁡θ2+r1θ˙1cos⁡θ1sin⁡θ2+r˙2sin⁡2θ2+r2θ˙2cos⁡θ2sin⁡θ22]=−(b2b+c2b∣v⃗2b∣)[r˙1(cos⁡θ1cos⁡θ2+sin⁡θ1sin⁡θ2)+r1θ˙1(−sin⁡θ1cos⁡θ2+cos⁡θ1sin⁡θ2)+r˙2(cos⁡2θ2+sin⁡2θ2)+r2θ˙2(−sin⁡θ2cos⁡θ2+cos⁡θ2sin⁡θ2)]−12(b2r+c2r∣v⃗2r∣)[r˙1(cos⁡θ1cos⁡θ2+sin⁡θ1sin⁡θ2)+r1θ˙1(−sin⁡θ1cos⁡θ2+cos⁡θ1sin⁡θ2)+r˙2(cos⁡2θ2+sin⁡2θ2)+r2θ˙2(−sin⁡θ2cos⁡θ2+cos⁡θ2sin⁡θ2)2]=−(b2b+c2b∣v⃗2b∣)[r˙1cos⁡(θ2−θ1)+r1θ˙1sin⁡(θ2−θ1)+r˙2]−12(b2r+c2r∣v⃗2r∣)[r˙1cos⁡Δ+r1θ˙1sin⁡Δ+r˙22]=−(b2b+c2b∣v⃗2b∣)[r˙1cos⁡Δ+r1θ˙1sin⁡Δ+r˙2]−12(b2r+c2r∣v⃗2r∣)[r˙1cos⁡Δ+r1θ˙1sin⁡Δ+r˙22].\begin{aligned} Q_{r_2} &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)\vec{v}_{1b} \cdot \dfrac{\partial \vec{r}_{1b}}{\partial r_2}-(b_{1r}+c_{1r}|\vec{v}_{1r}|)\vec{v}_{1r} \cdot \dfrac{\partial \vec{r}_{1r}}{\partial r_2} - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\vec{v}_{2b} \cdot \dfrac{\partial \vec{r}_{2b}}{\partial r_{2}} - (b_{2r}+c_{2r}|\vec{v}_{2r}|)\vec{v}_{2r} \cdot \dfrac{\partial \vec{r}_{2r}}{\partial r_{2}} \\ &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)\begin{bmatrix} \dot{r}_1\cos{\theta_1} - r_1\dot{\theta}_1\sin{\theta_1} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1\cos{\theta_1} \end{bmatrix} \cdot \vec{0} -(b_{1r}+c_{1r}|\vec{v}_{1r}|)\dfrac{1}{2}\begin{bmatrix} \dot{r}_1\cos{\theta_1} - r_1\dot{\theta}_1\sin{\theta_1} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1\cos{\theta_1} \end{bmatrix} \cdot \vec{0} - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\begin{bmatrix} \dot{r}_1 \cos{\theta_1} - r_1 \dot{\theta}_1 \sin{\theta_1} + \dot{r}_2\cos{\theta_2} - r_2\dot{\theta}_2 \sin{\theta_2} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1 \cos{\theta_1} + \dot{r}_2\sin{\theta_2} + r_2\dot{\theta}_2 \cos{\theta_2} \end{bmatrix} \cdot \begin{bmatrix} \cos{\theta_2} \\ \sin{\theta_2} \end{bmatrix} - (b_{2r}+c_{2r}|\vec{v}_{2r}|)\begin{bmatrix} \dot{r}_1 \cos{\theta_1} - r_1 \dot{\theta}_1 \sin{\theta_1} + \dfrac{\dot{r}_2\cos{\theta_2} - r_2\dot{\theta}_2 \sin{\theta_2}}{2} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1 \cos{\theta_1} + \dfrac{\dot{r}_2\sin{\theta_2} + r_2\dot{\theta}_2 \cos{\theta_2}}{2} \end{bmatrix} \cdot \dfrac{1}{2}\begin{bmatrix} \cos{\theta_2} \\ \sin{\theta_2} \end{bmatrix} \\ &= - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[\dot{r}_1\cos{\theta_1}\cos{\theta_2} - r_1\dot{\theta}_1\sin{\theta_1}\cos{\theta_2} + \dot{r}_2\cos^2{\theta_2}-r_2\dot{\theta}_2\sin{\theta_2}\cos{\theta_2} + \dot{r}_1\sin{\theta_1}\sin{\theta_2} + r_1\dot{\theta}_1\cos{\theta_1}\sin{\theta_2} + \dot{r}_2\sin^2{\theta_2} + r_2\dot{\theta}_2 \cos{\theta_2}\sin{\theta_2}\right] - \dfrac{1}{2}(b_{2r}+c_{2r}|\vec{v}_{2r}|)\left[\dot{r}_1\cos{\theta_1}\cos{\theta_2} - r_1\dot{\theta}_1\sin{\theta_1}\cos{\theta_2} + \dfrac{\dot{r}_2\cos^2{\theta_2}-r_2\dot{\theta}_2\sin{\theta_2}\cos{\theta_2}}{2} + \dot{r}_1\sin{\theta_1}\sin{\theta_2} + r_1\dot{\theta}_1\cos{\theta_1}\sin{\theta_2} + \dfrac{\dot{r}_2\sin^2{\theta_2} + r_2\dot{\theta}_2 \cos{\theta_2}\sin{\theta_2}}{2}\right] \\ &= - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[\dot{r}_1(\cos{\theta_1}\cos{\theta_2}+\sin{\theta_1}\sin{\theta_2}) + r_1\dot{\theta}_1(-\sin{\theta_1}\cos{\theta_2} + \cos{\theta_1}\sin{\theta_2}) + \dot{r}_2(\cos^2{\theta_2} + \sin^2{\theta_2})+r_2\dot{\theta}_2(-\sin{\theta_2}\cos{\theta_2} + \cos{\theta_2}\sin{\theta_2})\right] - \dfrac{1}{2}(b_{2r}+c_{2r}|\vec{v}_{2r}|)\left[\dot{r}_1(\cos{\theta_1}\cos{\theta_2}+\sin{\theta_1}\sin{\theta_2}) + r_1\dot{\theta}_1(-\sin{\theta_1}\cos{\theta_2} + \cos{\theta_1}\sin{\theta_2}) + \dfrac{\dot{r}_2(\cos^2{\theta_2} + \sin^2{\theta_2})+r_2\dot{\theta}_2(-\sin{\theta_2}\cos{\theta_2} + \cos{\theta_2}\sin{\theta_2})}{2}\right] \\ &= - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[\dot{r}_1\cos{(\theta_2-\theta_1)} + r_1\dot{\theta}_1\sin{(\theta_2-\theta_1)} + \dot{r}_2\right] - \dfrac{1}{2}(b_{2r}+c_{2r}|\vec{v}_{2r}|)\left[\dot{r}_1\cos{\Delta} + r_1\dot{\theta}_1\sin{\Delta} + \dfrac{\dot{r}_2}{2}\right]\\ &= - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[\dot{r}_1\cos{\Delta} + r_1\dot{\theta}_1\sin{\Delta} + \dot{r}_2\right]- \dfrac{1}{2}(b_{2r}+c_{2r}|\vec{v}_{2r}|)\left[\dot{r}_1\cos{\Delta} + r_1\dot{\theta}_1\sin{\Delta} + \dfrac{\dot{r}_2}{2}\right]. \end{aligned}

Hence the Euler-Lagrange equation for r2r_2 with dissipative forces is

M2(r¨2−r2θ˙22)+μ2(cos⁡Δ(r¨1−r1θ˙12)+sin⁡Δ(r1θ¨1+2r˙1θ˙1)+gsin⁡θ2)+k2(r2−l2)=Qr2.\begin{aligned} M_2 (\ddot{r}_2 - r_2\dot{\theta}_2^2) + \mu_2 (\cos{\Delta}(\ddot{r}_1 - r_1\dot{\theta}_1^2)+\sin{\Delta}(r_1\ddot{\theta}_1 + 2\dot{r}_1\dot{\theta}_1)+g\sin{\theta_2}) + k_2(r_2-l_2) &= Q_{r_2}. \end{aligned}

Or, in matrix form

[μ2cos⁡ΔM2μ2r1sin⁡Δ0][r¨1r¨2θ¨1θ¨2]=[M2r2θ˙22+μ2(r1θ˙12cos⁡Δ−2r˙1θ˙1sin⁡Δ−gsin⁡θ2)−k2(r2−l2)+Qr2].\begin{aligned} \begin{bmatrix} \mu_2 \cos{\Delta} & M_2 & \mu_2 r_1\sin{\Delta} & 0 \end{bmatrix} \begin{bmatrix} \ddot{r}_1 \\ \ddot{r}_2 \\ \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{bmatrix} &= \begin{bmatrix} M_2r_2\dot{\theta}_2^2 + \mu_2 (r_1\dot{\theta}_1^2\cos{\Delta}-2\dot{r}_1\dot{\theta}_1\sin{\Delta} - g\sin{\theta_2}) - k_2(r_2-l_2) + Q_{r_2} \end{bmatrix}. \end{aligned}

θ1\theta_1

As for θ1\theta_1

δθ1′L=M12δθ1′∣v⃗1b∣2+M22δθ1′∣v⃗2bI∣2+μ22(δθ1′∣Δv⃗2,1I∣2)+μ1gr1cos⁡θ1=M12(2r12θ¨1+4r1r˙1θ˙1)+μ22(2cos⁡Δ(2r1r˙2θ˙2+r1r2θ¨2)+2sin⁡Δ(r1r¨2−r1r2θ˙22))+μ1gr1cos⁡θ1=M1(r12θ¨1+2r1r˙1θ˙1)+μ2(cos⁡Δ(2r1r˙2θ˙2+r1r2θ¨2)+sin⁡Δ(r1r¨2−r1r2θ˙22))+μ1gr1cos⁡θ1.\begin{aligned} \delta'_{\theta_1} \mathcal{L} &= \dfrac{M_1}{2} \delta'_{\theta_1} |\vec{v}_{1b}|^2 + \dfrac{M_2}{2} \delta'_{\theta_1} |\vec{v}_{2b}^I|^2 + \dfrac{\mu_2}{2}(\delta'_{\theta_1}|\Delta \vec{v}_{2,1}^I|^2) + \mu_1 gr_1 \cos{\theta_1} \\ &= \dfrac{M_1}{2} (2r_1^2 \ddot{\theta}_1 + 4r_1\dot{r}_1\dot{\theta}_1) + \dfrac{\mu_2}{2}(2\cos{\Delta}(2r_1\dot{r}_2\dot{\theta}_2+r_1r_2\ddot{\theta}_2) +2\sin{\Delta}(r_1\ddot{r}_2-r_1r_2\dot{\theta}_2^2)) + \mu_1 gr_1 \cos{\theta_1} \\ &= M_1 (r_1^2 \ddot{\theta}_1 + 2r_1\dot{r}_1\dot{\theta}_1) + \mu_2(\cos{\Delta}(2r_1\dot{r}_2\dot{\theta}_2+r_1r_2\ddot{\theta}_2) +\sin{\Delta}(r_1\ddot{r}_2-r_1r_2\dot{\theta}_2^2)) + \mu_1 gr_1 \cos{\theta_1}. \end{aligned}

Generalized dissipative force

Qθ1=−(b1b+c1b∣v⃗1b∣)v⃗1b⋅∂r⃗1b∂θ1−(b1r+c1r∣v⃗1r∣)v⃗1r⋅∂r⃗1r∂θ1−(b2b+c2b∣v⃗2b∣)v⃗2b⋅∂r⃗2b∂θ1−(b2r+c2r∣v⃗2r∣)v⃗2r⋅∂r⃗2r∂θ1=−(b1b+c1b∣v⃗1b∣)[r˙1cos⁡θ1−r1θ˙1sin⁡θ1r˙1sin⁡θ1+r1θ˙1cos⁡θ1]⋅r1[−sin⁡θ1cos⁡θ1]−14(b1r+c1r∣v⃗1r∣)[r˙1cos⁡θ1−r1θ˙1sin⁡θ1r˙1sin⁡θ1+r1θ˙1cos⁡θ1]⋅r1[−sin⁡θ1cos⁡θ1]−(b2b+c2b∣v⃗2b∣)[r˙1cos⁡θ1−r1θ˙1sin⁡θ1+r˙2cos⁡θ2−r2θ˙2sin⁡θ2r˙1sin⁡θ1+r1θ˙1cos⁡θ1+r˙2sin⁡θ2+r2θ˙2cos⁡θ2]⋅r1[−sin⁡θ1cos⁡θ1]−(b2r+c2r∣v⃗2r∣)[r˙1cos⁡θ1−r1θ˙1sin⁡θ1+r˙2cos⁡θ2−r2θ˙2sin⁡θ22r˙1sin⁡θ1+r1θ˙1cos⁡θ1+r˙2sin⁡θ2+r2θ˙2cos⁡θ22]⋅r1[−sin⁡θ1cos⁡θ1]=−(b1b+c1b∣v⃗1b∣)[−r1r˙1cos⁡θ1sin⁡θ1+r12θ˙1sin⁡2θ1+r1r˙1sin⁡θ1cos⁡θ1+r12θ˙1cos⁡2θ1]−14(b1r+c1r∣v⃗1r∣)[−r1r˙1cos⁡θ1sin⁡θ1+r12θ˙1sin⁡2θ1+r1r˙1sin⁡θ1cos⁡θ1+r12θ˙1cos⁡2θ1]−(b2b+c2b∣v⃗2b∣)[−r1r˙1cos⁡θ1sin⁡θ1+r12θ1˙sin⁡2θ1−r1r˙2sin⁡θ1cos⁡θ2+r1r2θ˙2sin⁡θ1sin⁡θ2+r1r˙1cos⁡θ1sin⁡θ1+r12θ˙1cos⁡2θ1+r1r˙2cos⁡θ1sin⁡θ2+r1r2θ˙2cos⁡θ1cos⁡θ2]−(b2r+c2r∣v⃗2r∣)[−r1r˙1cos⁡θ1sin⁡θ1+r12θ1˙sin⁡2θ1−r1r˙2sin⁡θ1cos⁡θ2+r1r2θ˙2sin⁡θ1sin⁡θ22+r1r˙1cos⁡θ1sin⁡θ1+r12θ˙1cos⁡2θ1+r1r˙2cos⁡θ1sin⁡θ2+r1r2θ˙2cos⁡θ1cos⁡θ22]=−(b1b+c1b∣v⃗1b∣)[r1r˙1(−cos⁡θ1sin⁡θ1+sin⁡θ1cos⁡θ1)+r12θ˙1(sin⁡2θ1+cos⁡2θ1)]−(b2b+c2b∣v⃗2b∣)[r1r˙1(−cos⁡θ1sin⁡θ1+cos⁡θ1sin⁡θ1)+r12θ1˙(sin⁡2θ1+cos⁡2θ1)+r1r˙2(−sin⁡θ1cos⁡θ2+cos⁡θ1sin⁡θ2)+r1r2θ˙2(sin⁡θ1sin⁡θ2+cos⁡θ1cos⁡θ2)]=−(b1b+c1b∣v⃗1b∣)r12θ˙1−14(b1r+c1r∣v⃗1r∣)r12θ˙1−(b2b+c2b∣v⃗2b∣)[r12θ1˙+r1r˙2sin⁡Δ+r1r2θ˙2cos⁡Δ]−(b2r+c2r∣v⃗2r∣)[r12θ1˙+r1r˙2sin⁡Δ+r1r2θ˙2cos⁡Δ2].\begin{aligned} Q_{\theta_1} &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)\vec{v}_{1b} \cdot \dfrac{\partial \vec{r}_{1b}}{\partial \theta_1} -(b_{1r}+c_{1r}|\vec{v}_{1r}|)\vec{v}_{1r} \cdot \dfrac{\partial \vec{r}_{1r}}{\partial \theta_1} - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\vec{v}_{2b} \cdot \dfrac{\partial \vec{r}_{2b}}{\partial \theta_1}- (b_{2r}+c_{2r}|\vec{v}_{2r}|)\vec{v}_{2r} \cdot \dfrac{\partial \vec{r}_{2r}}{\partial \theta_1} \\ &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)\begin{bmatrix} \dot{r}_1\cos{\theta_1} - r_1\dot{\theta}_1\sin{\theta_1} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1\cos{\theta}_1 \end{bmatrix} \cdot r_1\begin{bmatrix} -\sin{\theta_1} \\ \cos{\theta_1} \end{bmatrix} -\dfrac{1}{4}(b_{1r}+c_{1r}|\vec{v}_{1r}|)\begin{bmatrix} \dot{r}_1\cos{\theta_1} - r_1\dot{\theta}_1\sin{\theta_1} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1\cos{\theta}_1 \end{bmatrix} \cdot r_1\begin{bmatrix} -\sin{\theta_1} \\ \cos{\theta_1} \end{bmatrix} \\ &- (b_{2b}+c_{2b}|\vec{v}_{2b}|) \begin{bmatrix} \dot{r}_1 \cos{\theta_1} - r_1 \dot{\theta}_1 \sin{\theta_1} + \dot{r}_2\cos{\theta_2} - r_2\dot{\theta}_2 \sin{\theta_2} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1 \cos{\theta_1} + \dot{r}_2\sin{\theta_2} + r_2\dot{\theta}_2 \cos{\theta_2} \end{bmatrix} \cdot r_1\begin{bmatrix} -\sin{\theta_1} \\ \cos{\theta_1} \end{bmatrix} \\ &- (b_{2r}+c_{2r}|\vec{v}_{2r}|) \begin{bmatrix} \dot{r}_1 \cos{\theta_1} - r_1 \dot{\theta}_1 \sin{\theta_1} + \dfrac{\dot{r}_2\cos{\theta_2} - r_2\dot{\theta}_2 \sin{\theta_2}}{2} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1 \cos{\theta_1} + \dfrac{\dot{r}_2\sin{\theta_2} + r_2\dot{\theta}_2 \cos{\theta_2}}{2} \end{bmatrix} \cdot r_1\begin{bmatrix} -\sin{\theta_1} \\ \cos{\theta_1} \end{bmatrix} &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)\left[-r_1\dot{r}_1\cos{\theta_1}\sin{\theta_1}+r_1^2\dot{\theta}_1\sin^2{\theta_1} + r_1\dot{r}_1\sin{\theta_1}\cos{\theta_1}+r_1^2\dot{\theta}_1\cos^2{\theta_1}\right] \\ &-\dfrac{1}{4}(b_{1r}+c_{1r}|\vec{v}_{1r}|)\left[-r_1\dot{r}_1\cos{\theta_1}\sin{\theta_1}+r_1^2\dot{\theta}_1\sin^2{\theta_1} + r_1\dot{r}_1\sin{\theta_1}\cos{\theta_1}+r_1^2\dot{\theta}_1\cos^2{\theta_1}\right] \\ &- (b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[-r_1\dot{r}_1\cos{\theta_1}\sin{\theta_1} + r_1^2\dot{\theta_1}\sin^2{\theta_1} -r_1\dot{r}_2\sin{\theta_1}\cos{\theta_2}+r_1r_2\dot{\theta}_2\sin{\theta_1}\sin{\theta_2} \right.\\ &\left.+r_1\dot{r}_1\cos{\theta_1}\sin{\theta_1}+r_1^2\dot{\theta}_1\cos^2{\theta_1} + r_1\dot{r}_2\cos{\theta_1}\sin{\theta_2}+r_1r_2\dot{\theta}_2\cos{\theta_1}\cos{\theta_2}\right]\\ &- (b_{2r}+c_{2r}|\vec{v}_{2r}|)\left[-r_1\dot{r}_1\cos{\theta_1}\sin{\theta_1} + r_1^2\dot{\theta_1}\sin^2{\theta_1} -\dfrac{r_1\dot{r}_2\sin{\theta_1}\cos{\theta_2}+r_1r_2\dot{\theta}_2\sin{\theta_1}\sin{\theta_2}}{2} \right.\\ &\left.+r_1\dot{r}_1\cos{\theta_1}\sin{\theta_1}+r_1^2\dot{\theta}_1\cos^2{\theta_1} + \dfrac{r_1\dot{r}_2\cos{\theta_1}\sin{\theta_2}+r_1r_2\dot{\theta}_2\cos{\theta_1}\cos{\theta_2}}{2}\right] &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)\left[r_1\dot{r}_1(-\cos{\theta_1}\sin{\theta_1}+\sin{\theta_1}\cos{\theta_1}) + r_1^2\dot{\theta}_1(\sin^2{\theta_1} +\cos^2{\theta_1})\right] - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[r_1\dot{r}_1(-\cos{\theta_1}\sin{\theta_1} + \cos{\theta_1}\sin{\theta_1}) +r_1^2\dot{\theta_1}(\sin^2{\theta_1}+\cos^2{\theta_1})\right.\\ &\left.+r_1\dot{r}_2(-\sin{\theta_1}\cos{\theta_2}+\cos{\theta_1}\sin{\theta_2})+r_1r_2\dot{\theta}_2(\sin{\theta_1}\sin{\theta_2}+\cos{\theta_1}\cos{\theta_2})\right] \\ &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)r_1^2\dot{\theta}_1 -\dfrac{1}{4}(b_{1r}+c_{1r}|\vec{v}_{1r}|)r_1^2\dot{\theta}_1 - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[ r_1^2\dot{\theta_1}+r_1\dot{r}_2\sin{\Delta}+r_1r_2\dot{\theta}_2\cos{\Delta}\right] \\ &- (b_{2r}+c_{2r}|\vec{v}_{2r}|)\left[ r_1^2\dot{\theta_1}+\dfrac{r_1\dot{r}_2\sin{\Delta}+r_1r_2\dot{\theta}_2\cos{\Delta}}{2}\right]. \end{aligned}

Hence Equation (1) is

M1(r12θ¨1+2r1r˙1θ˙1)+μ2(cos⁡Δ(2r1r˙2θ˙2+r1r2θ¨2)+sin⁡Δ(r1r¨2−r1r2θ˙22))+μ1gr1cos⁡θ1=Qθ1.\begin{aligned} M_1 (r_1^2 \ddot{\theta}_1 + 2r_1\dot{r}_1\dot{\theta}_1) + \mu_2(\cos{\Delta}(2r_1\dot{r}_2\dot{\theta}_2+r_1r_2\ddot{\theta}_2) +\sin{\Delta}(r_1\ddot{r}_2-r_1r_2\dot{\theta}_2^2)) + \mu_1 gr_1 \cos{\theta_1} &= Q_{\theta_1}. \end{aligned}

Hence the corresponding line of our matrix equation is

[0μ2r1sin⁡ΔM1r12μ2r1r2cos⁡Δ][r¨1r¨2θ¨1θ¨2]=[−2M1r1r˙1θ˙1+μ2(r1r2θ˙22sin⁡Δ−2r1r˙2θ˙2cos⁡Δ)−μ1gr1cos⁡θ1+Qθ1].\begin{aligned} \begin{bmatrix} 0 & \mu_2 r_1\sin{\Delta} & M_1r_1^2 & \mu_2 r_1r_2\cos{\Delta} \end{bmatrix} \begin{bmatrix} \ddot{r}_1 \\ \ddot{r}_2 \\ \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{bmatrix} &= \begin{bmatrix} -2M_1r_1\dot{r}_1 \dot{\theta}_1 + \mu_2(r_1r_2\dot{\theta}_2^2 \sin{\Delta} -2r_1\dot{r}_2 \dot{\theta}_2\cos{\Delta}) - \mu_1 gr_1\cos{\theta_1} + Q_{\theta_1} \end{bmatrix}. \end{aligned}

θ2\theta_2

As for θ2\theta_2

δθ2′L=M12δθ2′∣v⃗1b∣2+M22δθ2′∣v⃗2bI∣2+μ22(δθ2′∣Δv⃗2,1I∣2+2gr2cos⁡θ2)=M22(2r22θ¨2+4r2r˙2θ˙2)+μ22(2cos⁡Δ(2r˙1r2θ˙1+r1r2θ¨1)+2sin⁡Δ(r1r2θ˙12−r¨1r2)+2gr2cos⁡θ2)=M2(r22θ¨2+2r2r˙2θ˙2)+μ2(cos⁡Δ(2r˙1r2θ˙1+r1r2θ¨1)+sin⁡Δ(r1r2θ˙12−r¨1r2)+gr2cos⁡θ2).\begin{aligned} \delta'_{\theta_2} \mathcal{L} &= \dfrac{M_1}{2} \delta'_{\theta_2} |\vec{v}_{1b}|^2 + \dfrac{M_2}{2} \delta'_{\theta_2} |\vec{v}_{2b}^I|^2 + \dfrac{\mu_2}{2}(\delta'_{\theta_2} |\Delta \vec{v}_{2,1}^I|^2+2gr_2\cos{\theta_2}) \\ &= \dfrac{M_2}{2}(2r_2^2 \ddot{\theta}_2 + 4r_2 \dot{r}_2 \dot{\theta}_2) + \dfrac{\mu_2}{2}\left(2\cos{\Delta}(2\dot{r}_1r_2\dot{\theta}_1+r_1r_2\ddot{\theta}_1)+2\sin{\Delta}(r_1r_2\dot{\theta}_1^2-\ddot{r}_1r_2)+2gr_2\cos{\theta_2}\right) \\ &= M_2(r_2^2 \ddot{\theta}_2 + 2r_2 \dot{r}_2 \dot{\theta}_2) + \mu_2 (\cos{\Delta}(2\dot{r}_1r_2\dot{\theta}_1+r_1r_2\ddot{\theta}_1)+\sin{\Delta}(r_1r_2\dot{\theta}_1^2-\ddot{r}_1r_2)+gr_2\cos{\theta_2}). \end{aligned}
Qθ2=−(b1b+c1b∣v⃗1b∣)v⃗1b⋅∂r⃗1b∂θ2−(b1r+c1r∣v⃗1r∣)v⃗1r⋅∂r⃗1r∂θ2−(b2b+c2b∣v⃗2b∣)v⃗2b⋅∂r⃗2b∂θ2−(b2r+c2r∣v⃗2r∣)v⃗2r⋅∂r⃗2r∂θ2=−(b2b+c2b∣v⃗2b∣)[r˙1cos⁡θ1−r1θ˙1sin⁡θ1+r˙2cos⁡θ2−r2θ˙2sin⁡θ2r˙1sin⁡θ1+r1θ˙1cos⁡θ1+r˙2sin⁡θ2+r2θ˙2cos⁡θ2]⋅r2[−sin⁡θ2cos⁡θ2]−(b2r+c2r∣v⃗2r∣)[r˙1cos⁡θ1−r1θ˙1sin⁡θ1+r˙2cos⁡θ2−r2θ˙2sin⁡θ22r˙1sin⁡θ1+r1θ˙1cos⁡θ1+r˙2sin⁡θ2+r2θ˙2cos⁡θ22]⋅r22[−sin⁡θ2cos⁡θ2]=−(b2b+c2b∣v⃗2b∣)[−r˙1r2cos⁡θ1sin⁡θ2+r1r2θ˙1sin⁡θ1sin⁡θ2−r2r˙2cos⁡θ2sin⁡θ2+r22θ˙2sin⁡2θ2+r˙1r2sin⁡θ1cos⁡θ2+r1r2θ˙1cos⁡θ1cos⁡θ2+r2r˙2sin⁡θ2cos⁡θ2+r22θ˙2cos⁡2θ2]−12(b2r+c2r∣v⃗2r∣)[−r˙1r2cos⁡θ1sin⁡θ2+r1r2θ˙1sin⁡θ1sin⁡θ2−r2r˙2cos⁡θ2sin⁡θ2+r22θ˙2sin⁡2θ22+r˙1r2sin⁡θ1cos⁡θ2+r1r2θ˙1cos⁡θ1cos⁡θ2+r2r˙2sin⁡θ2cos⁡θ2+r22θ˙2cos⁡2θ22]=−(b2b+c2b∣v⃗2b∣)[r˙1r2(−cos⁡θ1sin⁡θ2+sin⁡θ1cos⁡θ2)+r1r2θ˙1(sin⁡θ1sin⁡θ2+cos⁡θ1cos⁡θ2)+r2r˙2(−cos⁡θ2sin⁡θ2+sin⁡θ2cos⁡θ2)+r22θ˙2(sin⁡2θ2+cos⁡2θ2)]−12(b2r+c2r∣v⃗2r∣)[r˙1r2(−cos⁡θ1sin⁡θ2+sin⁡θ1cos⁡θ2)+r1r2θ˙1(sin⁡θ1sin⁡θ2+cos⁡θ1cos⁡θ2)+r2r˙2(−cos⁡θ2sin⁡θ2+sin⁡θ2cos⁡θ2)+r22θ˙2(sin⁡2θ2+cos⁡2θ2)2]=−(b2b+c2b∣v⃗2b∣)[r22θ˙2−r˙1r2sin⁡Δ+r1r2θ˙1cos⁡Δ]−12(b2r+c2r∣v⃗2r∣)[r22θ˙22−r˙1r2sin⁡Δ+r1r2θ˙1cos⁡Δ].\begin{aligned} Q_{\theta_2} &= -(b_{1b}+c_{1b}|\vec{v}_{1b}|)\vec{v}_{1b} \cdot \dfrac{\partial \vec{r}_{1b}}{\partial \theta_2}-(b_{1r}+c_{1r}|\vec{v}_{1r}|)\vec{v}_{1r} \cdot \dfrac{\partial \vec{r}_{1r}}{\partial \theta_2} - (b_{2b}+c_{2b}|\vec{v}_{2b}|)\vec{v}_{2b}\cdot \dfrac{\partial \vec{r}_{2b}}{\partial \theta_2} -(b_{2r}+c_{2r}|\vec{v}_{2r}|)\vec{v}_{2r}\cdot \dfrac{\partial \vec{r}_{2r}}{\partial \theta_2} \\ &= -(b_{2b}+c_{2b}|\vec{v}_{2b}|)\begin{bmatrix} \dot{r}_1 \cos{\theta_1} - r_1 \dot{\theta}_1 \sin{\theta_1} + \dot{r}_2\cos{\theta_2} - r_2\dot{\theta}_2 \sin{\theta_2} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1 \cos{\theta_1} + \dot{r}_2\sin{\theta_2} + r_2\dot{\theta}_2 \cos{\theta_2} \end{bmatrix} \cdot r_2\begin{bmatrix} -\sin{\theta_2}\\ \cos{\theta_2} \end{bmatrix} -(b_{2r}+c_{2r}|\vec{v}_{2r}|)\begin{bmatrix} \dot{r}_1 \cos{\theta_1} - r_1 \dot{\theta}_1 \sin{\theta_1} + \dfrac{\dot{r}_2\cos{\theta_2} - r_2\dot{\theta}_2 \sin{\theta_2}}{2} \\ \dot{r}_1\sin{\theta_1} + r_1\dot{\theta}_1 \cos{\theta_1} + \dfrac{\dot{r}_2\sin{\theta_2} + r_2\dot{\theta}_2 \cos{\theta_2}}{2} \end{bmatrix} \cdot \dfrac{r_2}{2}\begin{bmatrix} -\sin{\theta_2}\\ \cos{\theta_2} \end{bmatrix} \\ &= -(b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[-\dot{r}_1r_2\cos{\theta}_1\sin{\theta_2} + r_1r_2\dot{\theta}_1\sin{\theta_1}\sin{\theta_2} - r_2\dot{r}_2\cos{\theta_2}\sin{\theta_2} + r_2^2\dot{\theta}_2\sin^2{\theta_2} + \dot{r}_1r_2\sin{\theta_1}\cos{\theta_2} + r_1r_2\dot{\theta}_1\cos{\theta_1}\cos{\theta_2} + r_2\dot{r}_2\sin{\theta_2}\cos{\theta_2} + r_2^2\dot{\theta}_2\cos^2{\theta_2}\right]\\ &-\dfrac{1}{2}(b_{2r}+c_{2r}|\vec{v}_{2r}|)\left[-\dot{r}_1r_2\cos{\theta}_1\sin{\theta_2} + r_1r_2\dot{\theta}_1\sin{\theta_1}\sin{\theta_2} - \dfrac{r_2\dot{r}_2\cos{\theta_2}\sin{\theta_2} + r_2^2\dot{\theta}_2\sin^2{\theta_2}}{2} + \dot{r}_1r_2\sin{\theta_1}\cos{\theta_2} + r_1r_2\dot{\theta}_1\cos{\theta_1}\cos{\theta_2} + \dfrac{r_2\dot{r}_2\sin{\theta_2}\cos{\theta_2} + r_2^2\dot{\theta}_2\cos^2{\theta_2}}{2}\right]\\ &= -(b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[\dot{r}_1r_2(-\cos{\theta}_1\sin{\theta_2} + \sin{\theta_1}\cos{\theta_2}) + r_1r_2\dot{\theta}_1(\sin{\theta_1}\sin{\theta_2}+\cos{\theta_1}\cos{\theta_2}) + r_2\dot{r}_2(-\cos{\theta_2}\sin{\theta_2} + \sin{\theta_2}\cos{\theta_2}) + r_2^2\dot{\theta}_2(\sin^2{\theta_2} +\cos^2{\theta_2})\right]\\ &-\dfrac{1}{2}(b_{2r}+c_{2r}|\vec{v}_{2r}|)\left[\dot{r}_1r_2(-\cos{\theta}_1\sin{\theta_2} + \sin{\theta_1}\cos{\theta_2}) + r_1r_2\dot{\theta}_1(\sin{\theta_1}\sin{\theta_2}+\cos{\theta_1}\cos{\theta_2}) + \dfrac{r_2\dot{r}_2(-\cos{\theta_2}\sin{\theta_2} + \sin{\theta_2}\cos{\theta_2}) + r_2^2\dot{\theta}_2(\sin^2{\theta_2} +\cos^2{\theta_2})}{2}\right]\\ &= -(b_{2b}+c_{2b}|\vec{v}_{2b}|)\left[ r_2^2\dot{\theta}_2-\dot{r}_1r_2\sin{\Delta} + r_1r_2\dot{\theta}_1\cos{\Delta}\right]-\dfrac{1}{2}(b_{2r}+c_{2r}|\vec{v}_{2r}|)\left[\dfrac{r_2^2\dot{\theta}_2}{2}-\dot{r}_1r_2\sin{\Delta} + r_1r_2\dot{\theta}_1\cos{\Delta}\right].\\ \end{aligned}

Hence Equation (1) for θ2\theta_2 is

M2(r22θ¨2+2r2r˙2θ˙2)+μ2(cos⁡Δ(2r˙1r2θ˙1+r1r2θ¨1)+sin⁡Δ(r1r2θ˙12−r¨1r2)+gr2cos⁡θ2)=Qθ2.\begin{aligned} M_2(r_2^2 \ddot{\theta}_2 + 2r_2 \dot{r}_2 \dot{\theta}_2) + \mu_2 (\cos{\Delta}(2\dot{r}_1r_2\dot{\theta}_1+r_1r_2\ddot{\theta}_1)+\sin{\Delta}(r_1r_2\dot{\theta}_1^2-\ddot{r}_1r_2)+gr_2\cos{\theta_2}) &= Q_{\theta_2}. \end{aligned}

Or, in matrix form

[−μ2r2sin⁡Δ0μ2r1r2cos⁡ΔM2r22][r¨1r¨2θ¨1θ¨2]=[−2M2r2r˙2θ˙2−μ2(2r˙1r2θ˙1cos⁡Δ+r1r2θ˙12sin⁡Δ+gr2cos⁡θ2)+Qθ2].\begin{aligned} \begin{bmatrix} -\mu_2r_2 \sin{\Delta} & 0 & \mu_2r_1 r_2\cos{\Delta} & M_2r_2^2 \end{bmatrix} \begin{bmatrix} \ddot{r}_1\\ \ddot{r}_2\\ \ddot{\theta}_1\\ \ddot{\theta}_2\\ \end{bmatrix} &= \begin{bmatrix} -2M_2 r_2\dot{r}_2 \dot{\theta}_2 - \mu_2(2\dot{r}_1r_2\dot{\theta}_1\cos{\Delta} + r_1r_2\dot{\theta}_1^2\sin{\Delta} + gr_2\cos{\theta_2}) + Q_{\theta_2} \end{bmatrix}. \end{aligned}

Matrix form

[M1μ2cos⁡Δ0−μ2r2sin⁡Δμ2cos⁡ΔM2μ2r1sin⁡Δ00μ2r1sin⁡ΔM1r12μ2r1r2cos⁡Δ−μ2r2sin⁡Δ0μ2r1r2cos⁡ΔM2r22][r¨1r¨2θ¨1θ¨2]=[M1r1θ˙12+μ2(r2θ˙22cos⁡Δ+2r˙2θ˙2sin⁡Δ)−μ1gsin⁡θ1−k1(r1−l1)+Qr1M2r2θ˙22+μ2(r1θ˙12cos⁡Δ−2r˙1θ˙1sin⁡Δ−gsin⁡θ2)−k2(r2−l2)+Qr2−2M1r1r˙1θ˙1+μ2(r1r2θ˙22sin⁡Δ−2r1r˙2θ˙2cos⁡Δ)−μ1gr1cos⁡θ1+Qθ1−2M2r2r˙2θ˙2−μ2(2r˙1r2θ˙1cos⁡Δ+r1r2θ˙12sin⁡Δ+gr2cos⁡θ2)+Qθ2]\begin{aligned} \begin{bmatrix} M_1 & \mu_2 \cos{\Delta} & 0 & -\mu_2r_2 \sin{\Delta}\\ \mu_2 \cos{\Delta} & M_2 & \mu_2 r_1\sin{\Delta} & 0 \\ 0 & \mu_2 r_1\sin{\Delta} & M_1r_1^2 & \mu_2 r_1r_2\cos{\Delta}\\ -\mu_2r_2 \sin{\Delta} & 0 & \mu_2r_1 r_2\cos{\Delta} & M_2r_2^2 \end{bmatrix} \begin{bmatrix} \ddot{r}_1\\ \ddot{r}_2\\ \ddot{\theta}_1\\ \ddot{\theta}_2\\ \end{bmatrix} &= \begin{bmatrix} M_1r_1\dot{\theta}_1^2 + \mu_2 (r_2 \dot{\theta}_2^2\cos{\Delta} + 2 \dot{r}_2 \dot{\theta}_2\sin{\Delta}) - \mu_1g\sin{\theta_1} - k_1(r_1-l_1) + Q_{r_1} \\ M_2r_2\dot{\theta}_2^2 + \mu_2 (r_1\dot{\theta}_1^2\cos{\Delta}-2\dot{r}_1\dot{\theta}_1\sin{\Delta} - g\sin{\theta_2}) - k_2(r_2-l_2) + Q_{r_2}\\ -2M_1r_1\dot{r}_1 \dot{\theta}_1 + \mu_2(r_1r_2\dot{\theta}_2^2 \sin{\Delta} -2r_1\dot{r}_2 \dot{\theta}_2\cos{\Delta}) - \mu_1 gr_1\cos{\theta_1} + Q_{\theta_1} \\ -2M_2 r_2\dot{r}_2 \dot{\theta}_2 - \mu_2(2\dot{r}_1r_2\dot{\theta}_1\cos{\Delta} + r_1r_2\dot{\theta}_1^2\sin{\Delta} + gr_2\cos{\theta_2}) + Q_{\theta_2} \end{bmatrix} \end{aligned}
© Brenton Horne. I do not claim to be infallible, so feel free to open issues or pull requests on GitHub if you notice flaws or things that can be improved.
I favour Oxford spelling, so if you notice I generally use different spelling to Australian spelling, that is probably why.
Homepage. Source code: https://github.com/fusion809/fusion809.github.io.
Website built with Franklin.jl.

Mathematical solver and related pages

Partial differential equation solver (non-interactive)

Heat equation solver

Ordinary differential equation solvers (interactive)

Complex analysis pages

Numerical mathematics

Lagrangian mechanics pages

Other mathematics

Education, my personal opinion largely

Linux


Last modified: August 13, 2025.