# Do these Lagrange equations of 1st kind exhibit numerical instabilities?

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$$: • The initial conditions must fulfilled the constraint equations, did you check it ?
– Eli
Aug 10 at 18:27
• Yes, I did. Thanks for the reminder! Aug 10 at 19:36
• Can you plot it over the simulation time
– Eli
Aug 10 at 19:40
• Sure. Will provide them tomorrow. Aug 10 at 20:27
• Suggestion: Consider to include the definition of 1st kind as it varies from author to author. Aug 16 at 5:32

Yes, integration of the crossbar system with Runge-Kutta scheme of order 4 and 5 does result in numerical instabilities. One can tell by removing the constant states $$x_2$$ and $$y_1$$ (and its timely derivatives) from the set of ordinary differential equations. Integrating the reduced system

$$\begin{eqnarray} m_1 \ddot{x}_1 & = & - \frac{\dot{x}_1^2 + \dot{y}_2^2 - y_2g}{x_1^2/m_1 + y_2^2/m_2}x_1, \\ m_2 \ddot{y}_2 & = & - \frac{\dot{x}_1^2 + \dot{y}_2^2 - y_2g}{x_1^2/m_1 + y_2^2/m_2}y_2 - m_2 g, \end{eqnarray}$$

with $$\textbf{x}_0 = (-L, \; 0, \; 0, \; 0)$$, and $$m_1 = m_2 = g = L = 1$$ makes a better approximation, but still violates constraint $$f_1$$ (blue: full order system, orange: reduced order system): Important to note that the choice of the integration method Radau reduced the error an order of magnitude of 3: This is surprising, since the two masses oscillate with the same frequency. Hence, no stiff behavior and no need for a stiff solver.

As for the double pendulum system, it is clear that the full set of eight ordinary differential equations somehow overdetermines the dynamics, since only four are required (two angles and angular velocities).

But it is still an open question to me what exactly goes wrong here.