Equations of motion for double spherical pendulum simply?

Question

I am attempting to simulate a double spherical pendulum, i.e. a combination of the spherical pendulum and the double pendulum.

I understand that the equations of motion can be derived via the Lagrangian and the Euler-Lagrange equations. However this method rapidly becomes very messy.

Is there an alternate method that could be used to simplify the calculations required?

-UPDATE-

This is an update to specify the Lagrangian of the system as requested in the comments.

We have two spherical pendulums, denoted $1$ and $2$. The lengths of each pendulum $l_1=l_2=1$ whilst the masses $m_1=m_2=1$.

There are two parameters that describe the location of the pendulum mass $\theta$, the angle from the vertical, and $\phi$ the azimuthal angle about the vertical axis. For a diagram see here.

The locations of the masses $r_1 , r_2$ are therefore given by (for a vertical $z$ axis): $$x_1 = sin(\theta_1)cos(\phi_1)$$ $$y_1 = sin(\theta_1)sin(\phi_1)$$ $$z_1 = -cos(\theta_1)$$

and

$$x_2 = x_1+sin(\theta_2)cos(\phi_2)$$ $$y_2 = y_1+sin(\theta_2)sin(\phi_2)$$ $$z_2 = z_1-cos(\theta_2)$$

The Lagrangian is given by the difference between the kinetic and potential energies $L = E_k -E_p$

$$E_k = \frac{1}{2} \left( \dot{r}_1^2 +\dot{r}_2^2 \right)$$ $$E_p = g (z_1 +z_2)$$

The Lagrangian then follows simply from taking the time derivatives of $r_1$ and $r_2$. (Note that this leads me to a very long and complicated form the Lagrangian)

Answer 1

Not that I know of. However if you're fine with considering only small oscillations, then you can replace $\sin \theta$ by $\theta$ and $\cos \theta$ by $1-\frac{\theta^2}{2}$. This might make things simpler although the solution you get will be acceptable for small angles.