The cartesian coordinates $x_1, y_1, x_2, y_2$ are given by:

$x_1 = L_1\sin\theta_1$		$y_1 = -L_1\cos\theta_1$
$x_2 = x_1 + L_2\sin\theta_2$		$y_2 = y_1 - L_2\cos\theta_2$

Their derivatives (needed later) are:

$\dot x_1 = L_1\cos\theta_1$		$\dot y_1 = L_1\sin\theta_1$
$\dot x_2 = \dot x_1 + L_2\cos\theta_2$		$\dot y_2 = \dot y_1 + L_2\sin\theta_2$

The total kinetic energy is given by:

$KE = \frac{1}{2} m_1 (\dot x_1^2 + \dot y_1^2) + \frac{1}{2} m_2 (\dot x_2^2 + \dot y_2^2)$

which after some algebra reduces to: $$KE = \frac{1}{2} M L_1^2 \dot\theta_1^2 + \frac{1}{2} m_2 L_2^2 \dot\theta_2^2 + m_2 L_1 L_2 \dot\theta_1\dot\theta_2\cos\Delta\theta\label{KE}$$ where $M\equiv m_1 + m_2$ and $\Delta\theta\equiv \theta_1 - \theta_2$.

The potential energy due to gravity is the usual:

$PE = m_1 g y_1 + m_2 g y_2$

which reduces to

$$PE = -MgL_1\cos\theta_1 - m_2gL_2\cos\theta_2\label{PE}$$

The full Lagrangian is therefore: $$\mathcal{L} = \frac{1}{2} M L_1^2 \dot\theta_1^2 + \frac{1}{2} m_2 L_2^2 \dot\theta_2^2 + m_2 L_1 L_2 \dot\theta_1\dot\theta_2\cos\Delta\theta + MgL_1\cos\theta_1 + m_2gL_2\cos\theta_2\label{Lagrangian}$$ After applying the Euler-Lagrange equations to $(\dot\theta_1,\theta_1)$ and $(\dot\theta_2,\theta_2)$ and simplying a bit we get the following two equations:

$$ML_1\ddot\theta_1 + m_2L_2\ddot\theta_2\cos\Delta\theta + m_2L_2\dot\theta_2^2\sin\Delta\theta + Mg\sin\theta_1 = 0\label{eqth1}$$ $$L_1\ddot\theta_1\cos\Delta\theta + L_2\ddot\theta_2 - L_1\dot\theta_1^2\sin\Delta\theta + g\sin\theta_2 = 0\label{eqth2}$$

These are of course very non-linear coupled systems, not solvable analytically.

It is useful first to take a few limits. For instance, setting $m_2 = 0$, $L_2 = 0$, and $\theta_1 = \theta_2$, both equations become the simple pendulum equation $\ddot\theta_1 + \omega^2\sin\theta_1 = 0$ where $\omega^2 \equiv g/L$.

To solve these equations numerically in a simulation, we first have to rearrange into two equations, each of which have only 1 second derivative in the time. After a lot of algebra you should get the following two equations: $$\ddot\theta_1 = [-\sin\Delta\theta(m_2L_1\dot\theta_1^2\cos\Delta\theta + m_2L_2\dot\theta_2^2) - g(M\sin\theta_1 - m_2\sin\theta_2\cos\Delta\theta)]/L_1\alpha\label{eq1}$$ $$\ddot\theta_2 = [+\sin\Delta\theta(ML_1\dot\theta_1^2 + m_2L_2\dot\theta_2^2\cos\Delta\theta) + g(M\sin\theta_1\cos\Delta\theta - M\sin\theta_2)]/L_2\alpha\label{eq2}$$ and where $\alpha \equiv m_1 + m_2\sin^2\Delta\theta$.

In the simulation below, we use 3 common methods for the numerical integration: Euler's method; the modified Euler-Cromer; and Runge-Kutta (order 2, RK2). All 3 start with the 2 basic equations, ($\ref{eq1}$) and ($\ref{eq2}$), written in the following way:

$\ddot\theta_1 = f_1(\theta_1,\theta_2,\dot\theta_1,\dot\theta_2)$
$\ddot\theta_2 = f_2(\theta_1,\theta_2,\dot\theta_1,\dot\theta_2)$

where $f_1$ is given by equation ($\ref{eq1}$) and $f_2$ by equation ($\ref{eq2}$).

What we want is a way to use our knowledge of $\theta_1(t)$ to get $\theta_1(t+\dt)$ (and the same for $\theta_2, \dot\theta_1, and \dot\theta_2$). To do this, we first consider $\theta_1$ and make the substitution:

$\alpha\equiv\theta_1$
$\beta\equiv\dot\theta_1$

This allows us to write equation ($\ref{eq1}$), a second order differential equation, in terms of 2 first order differential equations:

$\dot\beta = f_1(\theta_1,\theta_2,\dot\theta_1,\dot\theta_2)$
$\dot\alpha = f_1'(\beta)=\beta$

The Euler technique uses the above 2 equations to numerically find $\alpha(t+\dt)$ and $\beta(t+\dt)$ based on knowledge of $\alpha(t)$, $\beta(t)$, is to use the approximation (based on calculus definitions of derivatives): $$\alpha(t+\dt)=\alpha(t)+\dt\cdot f_1'(\beta)\label{eqa1}$$ $$\beta(t+\dt)=\beta(t)+\dt\cdot f_1\label{eqa2}$$ This is basically just using the definition of derivatives before taking the limit:

$\frac{d\alpha}{dt}=f_1'=\lim_{\dt\to 0}\frac{\alpha(t+\dt)-\alpha(t)}{\dt}$

and solving for $\alpha(t+\dt)$ to get equation ($\ref{eqa1}$). And similarly for $\beta(t+\dt)$.

In words, what we are doing here is using the knowledge of $\alpha(t)$, the value at the beginning of the interval, and "swimming" to $\alpha(t+\dt)$ using the value of the derivative at the beginning of the interval. So to be explicit: $$\alpha(t+\dt)=\alpha(t)+\dt\cdot f_1'(\beta(t))$$ $$\beta(t+\dt)=\beta(t)+\dt\cdot f_1(\alpha(t),\dot\alpha(t),\beta(t),\dot\beta(t))$$ The Euler-Cromer technique is just a variation where we use the slope (derivative) not at the beginning, but at the end of the interval: $$\alpha(t+\dt)=\alpha(t)+\dt\cdot f_1'(\beta(t+\dt))$$ $$\beta(t+\dt)=\beta(t)+\dt\cdot f_1(\alpha(t+\dt),\dot\alpha(t+\dt),\beta(t+\dt),\dot\beta(t+\dt))$$ And finally, the Runge-Kutta (2nd order, RK2) uses the value of the slope in the middle of the interval: $$\alpha(t+\dt)=\alpha(t)+\dt\cdot f_1'(\beta_h)$$ $$\beta(t+\dt)=\beta(t)+\dt\cdot f_1(\alpha_h,\dot\alpha_h,\beta_h,\dot\beta_h)$$ where $$\alpha_h\equiv\alpha(t)+\half\dt\cdot \dot\alpha(t)$$ $$\dot\alpha_h\equiv\dot\alpha(t)+\half\dt\cdot f_1'(\beta(t))$$ and etc for $\beta$ and $\dot\beta$. The RK2 really amounts to starting at the beginning of the interval and using Euler to calculate the slope at the midpoint of the interval, $f_1(\alpha_h,\dot\alpha_h,\beta_h,\dot\beta_h)$, to find to get us to the end of the interval. For more, surf to

In the simulation below, you can choose which 3 numerical integration methods you want to use, and change the time interval $\dt$, the two masses $m_1$ and $m_2$, and the two pendulum lengths $L_1$ and $L_2$ by clicking on the arrows (hold the arrows down and the masses and pendlums and it will change continuously). You can even change the masses and lengths while the thing is in motion, just for fun!

To get a feel for how good each method is as an approximation, the bottom plot labeled "Energy (arbitrary scale..." is of the total energy E = KE + PE. Any deviation from a constant value shows how good (or bad) an approximation the method is. As you can see quite readily, Euler's approximation is pretty poor for such a non-linear system!