From the Lagrangian I've got the following equations of motion for the double pendulum in 2D. (The masses are different but the lengths of the two pendula are equal.) Let $m_2$ be the lowest-hanging mass.
$$(m_1+m_2)\ddot{\theta_1}+2m_2\ddot\theta_2\cos(\theta_2-\theta_1)=\\ -2m_2\dot\theta_1\dot\theta_2\sin(\theta_1-\theta_2)-(m_1+m_2)g/l\sin(\theta_1)$$
and
$$m_2\ddot{\theta_1}+2m_2\ddot\theta_2\cos(\theta_2-\theta_1)=\\ 2m_2\dot\theta_1\dot\theta_2\sin(\theta_1-\theta_2)-m_2g/l\sin(\theta_1)$$
In the small angle approximation these become, respectively
$$(m_1+m_2)\ddot{\theta_1}+2m_2\ddot\theta_2= -2m_2\dot\theta_1\dot\theta_2(\theta_1-\theta_2)-\theta_1(m_1+m_2)g/l$$
and
$$m_2\ddot{\theta_1}+2m_2\ddot\theta_2= 2m_2\dot\theta_1\dot\theta_2(\theta_1-\theta_2)-\theta_1m_2g/l$$.
Most sources don't have the terms of order $\dot\theta$. This is because they apply the small angle approximation to the Lagrangian before taking the derivatives, thereby ignoring terms of order $\theta.$ What justification do we have for getting rid of these terms?