The kinetic & potential energies of the system are given by $$T=\frac{m}{2}\bigg(v_1^2+v_2^2\bigg)+\frac{I}{2} \bigg(\dot{\theta}_1^2+\dot{\theta}_2^2\bigg),$$ $$\text{and}$$ $$V=mg\cdot(y_1+y_2),$$ respectively. From this, one can construct the Lagrangian $$\mathcal{L}=T-V.$$ The Euler-Lagrange equation reads $$\frac{\partial\mathcal L}{\partial q}-\frac{d}{dt} \frac{\partial\mathcal L}{\partial\dot q}=0$$ and leads to the following equations of motion: $$\dot{\theta}_1=\frac{6}{ml^2}\cdot \frac{2p_{\theta_1}-3\cos(\theta_1-\theta_2) \cdot p_{\theta_2} }{16-9\cos^2(\theta_1-\theta_2}$$ $$\dot{\theta}_2=\frac{6}{ml^2}\cdot \frac{8p_{\theta_2}-3\cos(\theta_1-\theta_2) \cdot p_{\theta_1} }{16-9\cos^2(\theta_1-\theta_2}$$ $$\dot{p}_{\theta_1}= -\frac{ml^2}{2} \bigg(\dot{\theta}_1\dot{\theta}_2 \sin(\theta_1-\theta_2) +3\frac{g}{l}\sin(\theta_1)\bigg)$$ $$\dot{p}_{\theta_1}= -\frac{ml^2}{2} \bigg(-\dot{\theta}_1\dot{\theta}_2 \sin(\theta_1-\theta_2) +\frac{g}{l}\sin(\theta_2)\bigg)$$ For the above animation, these equations were integrated utilizing a 4th order Runge-Kutta scheme.