I followed the lead of "Theoretische Physik", 1e, 2015 by Bartelmann et al. (pp. 171 - 174) to form the set of constituting Lagrange equations of the 1st kind for the double pendulum: eight 1st order ODEs of the form
$$ \dot{x} = f(x, \lambda(x)), $$
where $x$ are states, and $\lambda$ are lagrange multipliers.
The defining equations for the lagrange multipliers are
$$ \sum_{i=1}^{2N} \frac{1}{m_i} \frac{\partial f_a}{\partial x_i} \sum_{b=1}^r \lambda_b \frac{\partial f_b}{\partial x_i} = g_a - \sum_{i=1}^{2N}\frac{1}{m_i} \frac{\partial f_a}{\partial x_i} F_i, \; \; a = 1..n $$
where $n$ is the number of constraints, $N$ is the number of mass points, $f_a$ are the (holonomic) constraints, and $g_a$ are terms obtained from taking the $2^{\text{nd}}$ total time derivative of the constraints:
$$ 0 = \ddot{f_a} = \sum_{i=1}^{2N} \frac{\partial f_a}{\partial x_i} \ddot{x_i} - g_a(t, x, \dot{x}). $$
Dynamics are then computed using
$$ m_i \ddot{x}_i = \sum_{a=1}^{r} \lambda_a \frac{\partial f_a}{\partial x_i} + F_i, \; \; i = 1..2N. $$
I observed the following:
While the integration of the mass points' ODEs yielded plausible results, the two masses were slightly off the constraints of circular motion about the pivot points which I stated as constraints.
Given my equations are right, could approximation errors in the ODE solver, or the linear equations involving the lagrange multipliers, be the cause of the non-circular motion? Could the constraints being considered only as second deriatives over time be the reason?
The constraints are as follows:
$$ \begin{eqnarray} f_1 & = & x_1^2 + y_1^2 - (l_1/2)^2 = 0 \\ f_2 & = & (x_2 - 2x_1)^2 + (y_2 - 2y_1)^2 - (l_2/2)^2 = 0 \end{eqnarray} $$
They state that mass $m_1$ is on a circle about $(0, 0)$ with radius $l_1/2$, and that mass $m_2$ is on a circle about the joint at position $(2x_1, 2x_2)$ with radius $l_2/2$.
The intial conditions are
$$ \begin{eqnarray} x_0 & = & (x_{10}, \; \dot{x}_{10}, \; y_{10}, \; \dot{y}_{10}, \; x_{20}, \; \dot{x}_{20}, \; y_{20}, \; \dot{y}_{20}) \\ & = & (0, \; 0, \; -l_1/2, \; 0, \; l_2/2, \; 0, \; -l_1, \; 0), \end{eqnarray} $$
which are consistent with the constraints.
The below given trajectory of $m_1$ in the $(x, y)$ plane shows that $m_1$ is off the circle about $(0, 0)$. For completeness, a time series plot of $m_1$'s $x$- and $y$-coordinates are shown.
Parameters for integration are $l_1 = l_2 = m_1 = m_2 = g = 1$.
Strangely, a similar drifting effect occurs when the movement of the two masses is constrained by the following set of equations. They represent a crossbar movement of point masses where mass $1$ moves horizontally, mass $2$ moves vertically, and both masses are connected by a rigid rod of length $L$:
$$ \begin{eqnarray} f_1 & = & x_1^2 + y_2^2 - L^2 = 0 \\ f_2 & = & x_2 = 0 \\ f_3 & = & y_1 = 0 \end{eqnarray} $$
Again, the initial conditions are consistent with the constraints:
$$ \begin{eqnarray} x_0 & = & (x_{10}, \; \dot{x}_{10}, \; y_{10}, \; \dot{y}_{10}, \; x_{20}, \; \dot{x}_{20}, \; y_{20}, \; \dot{y}_{20}) \\ & = & (-L, \; 0, \; 0, \; 0, \; 0, \; 0, \; 0, \; 0). \end{eqnarray} $$
Below is a plot showing the change of the length of the rod over time which clearly violates the first constraint $f_1$:
