26 Example 4
26.1 Double Pendulum
We consider a constrained formulation of a planar double pendulum as a differential-algebraic equation (DAE) system:
\left\{ \color{blue} \begin{aligned} x'_1 & = u_1\\ y'_1 & = v_1\\ x'_2 & = u_2\\ y'_2 & = v_2\\ u'_1 & = \dfrac{x_1\lambda_1 + (x_1-x_2)\lambda_2}{m_1} \\ v'_1 & = \dfrac{y_1\lambda_1 + (y_1-y_2)\lambda_2}{m_1} - g \\ u'_2 & = \dfrac{(x_2-x_1)\lambda_2}{m_2} \\ v'_2 & = \dfrac{(y_2-y_1)\lambda_2}{m_2} - g \\ \color{DarkRed} \ell_1^2 & \color{DarkRed} = x_1^2 + y_1^2\\ \color{DarkRed} \ell_2^2 & \color{DarkRed} = (x_1-x_2)^2 + (y_1-y_2)^2 \end{aligned} \right.
| Parameter | Value | Parameter | Value | Parameter | Value |
|---|---|---|---|---|---|
| m_1 | 5 | m_2 | 1 | \ell_1 | 1 |
| \ell_2 | 1 | g | 9.81 | \theta_1 | \pi/2 |
| \theta_2 | \pi/4 | \omega_1 | 0 | \omega_2 | -0.1 |
We define the variable vector and matrices:
\bm{E}(\bm{Z},t) = \begin{pmatrix} 1 & . & . & . & . & . & . & . & . & . \\[0.5em] . & 1 & . & . & . & . & . & . & . & . \\[0.5em] . & . & 1 & . & . & . & . & . & . & . \\[0.5em] . & . & . & 1 & . & . & . & . & . & . \\[0.5em] . & . & . & . & 1 & . & . & . & . & . \\[0.5em] . & . & . & . & . & 1 & . & . & . & . \\[0.5em] . & . & . & . & . & . & 1 & . & . & . \\[0.5em] . & . & . & . & . & . & . & 1 & . & . \\[0.5em] . & . & . & . & . & . & . & . & . & . \\[0.5em] . & . & . & . & . & . & . & . & . & . \end{pmatrix}, \quad \bm{Z} = \begin{pmatrix} x_1 \\[0.5em] y_1 \\[0.5em] x_2 \\[0.5em] y_2 \\[0.5em] u_1 \\[0.5em] v_1 \\[0.5em] u_2 \\[0.5em] v_2 \\[0.5em] \lambda_1 \\[0.5em] \lambda_2 \end{pmatrix},
\bm{G}(\bm{Z},t) = \begin{pmatrix} u_1 \\[0.5em] v_2 \\[0.5em] u_2 \\[0.5em] v_2 \\[0.5em] \dfrac{x_1\lambda_1 + (x_1-x_2)\lambda_2}{m_1} \\[0.5em] \dfrac{y_1\lambda_1 + (y_1-y_2)\lambda_2}{m_1} - g \\[0.5em] \dfrac{(x_2-x_1)\lambda_2}{m_2} \\[0.5em] \dfrac{(y_2-y_1)\lambda_2}{m_2} - g \\[0.5em] x_1^2 + y_1^2 - \ell_1^2\\[0.5em] (x_1-x_2)^2 + (y_1-y_2)^2 - \ell_2^2 \end{pmatrix}
So that the DAE reads:
\bm{E}(\bm{Z},t)\,\bm{Z}^\prime(t) = \bm{G}(\bm{Z},t)
26.2 Index reduction
The table
\bm{T} = \left[\begin{array}{c|c} \bm{E}(\bm{Z},t) & \bm{G}(\bm{Z},t) \end{array} \right] = \left[ \begin{array}{rrrrrrrrrr|c} 1 & . & . & . & . & . & . & . & . & . & u_1 \\[0.5em] . & 1 & . & . & . & . & . & . & . & . & v_2 \\[0.5em] . & . & 1 & . & . & . & . & . & . & . & u_2 \\[0.5em] . & . & . & 1 & . & . & . & . & . & . & v_2 \\[0.5em] . & . & . & . & 1 & . & . & . & . & . & \dfrac{x_1\lambda_1 + (x_1-x_2)\lambda_2}{m_1} \\[0.5em] . & . & . & . & . & 1 & . & . & . & . & \dfrac{y_1\lambda_1 + (y_1-y_2)\lambda_2}{m_1} - g \\[0.5em] . & . & . & . & . & . & 1 & . & . & . & \dfrac{(x_2-x_1)\lambda_2}{m_2} \\[0.5em] . & . & . & . & . & . & . & 1 & . & . & \dfrac{(y_2-y_1)\lambda_2}{m_2} - g \\[0.5em] . & . & . & . & . & . & . & . & . & . & x_1^2 + y_1^2 - \ell_1^2\\[0.5em] . & . & . & . & . & . & . & . & . & . & (x_1-x_2)^2 + (y_1-y_2)^2 - \ell_2^2 \end{array} \right]
is used to memorize the operation done on the DAE.
Remark 26.1 (semi-explicit). Notice that
\bm{E}(\bm{Z},t) = \begin{pmatrix} \bm{I} & \bm{0} \\ \bm{0} & \bm{0} \end{pmatrix}
this permits to simplify index reduction working only on the last two rows.
First reduction of the index by one
Derive the algebraic constraint
\begin{aligned} \color{blue}\text{[alg1]} &= \begin{pmatrix} x_1^2 + y_1^2 - \ell_1^2 \\[0.5em] (x_1-x_2)^2 + (y_1-y_2)^2 - \ell_2^2 \end{pmatrix} \\ \dfrac{\mathrm{d}\text{[alg1]}}{\mathrm{d}t} &= \begin{pmatrix} 2x_1x'_1+ 2y_1y'_1 \\[0.5em] 2(x_1-x_2)(x'_1-x'_2) + 2(y_1-y_2)(y'_1-y'_2) \end{pmatrix} \end{aligned}
Update reduced table \widetilde{\bm{T}} of the last two rows
\widetilde{\bm{T}} \leftarrow \left[ \begin{array}{rrrrrrrrrr|c} 2 x_1 & 2 y_1 & . & . & . & . & . & . & . & . & 0 \\[0.5em] 2 (x_1-x_2) & 2 (y_1-y_2) & -2 (x_1-x_2) & -2 (y_1-y_2) & . & . & . & . & . & . & 0 \end{array} \right]
apply elimination by subtract a combination of the fist four row to eliminate 2x_1, 2y_1 and so on. However this is equivalent to substitute x'_1 with u_1, y'_1 with v_1 and so on., So that
\text{[alg2]} = \dfrac{\mathrm{d}\text{[alg1]}}{\mathrm{d}t} = \begin{pmatrix} 2x_1u_1+ 2y_1v_1 \\[0.5em] 2(x_1-x_2)(u_1-u_2) + 2(y_1-y_2)(v_1-v_2) \end{pmatrix}
Second reduction of the index by one
Derive the constraint and apply substitution (using \delta x = x_1-x_2, \delta y = y_1-y_2, \delta u = u_1-u_2, \delta v = v_1-v_2)
\begin{aligned} \text{[alg3]}= \dfrac{\mathrm{d}\text{[alg2]}}{\mathrm{d}t} &= \dfrac{\mathrm{d}}{\mathrm{d}t} \begin{pmatrix} 2x_1u_1+ 2y_1v_1 \\[0.5em] 2(x_1-x_2)(u_1-u_2) + 2(y_1-y_2)(v_1-v_2) \end{pmatrix} \\ &= \dfrac{\mathrm{d}}{\mathrm{d}t} \begin{pmatrix} 2x_1u_1+ 2y_1v_1 \\[0.5em] 2\,\delta x\, \delta u_1 + 2\,\delta y\,\delta v \end{pmatrix} \\ & = \begin{pmatrix} 2x'_1u_1+ 2y'_1v_1 + 2x_1u'_1+ 2y_1v'_1 \\[0.5em] 2\,\delta x'\,\delta u + 2\,\delta y'\,\delta v + 2\,\delta x \,\delta u' + 2\,\delta y\,\delta v' \end{pmatrix} \\ & = \begin{pmatrix} 2u_1^2+ 2v_1^2 + 2x_1u'_1+ 2y_1v'_1 \\[0.5em] 2\,\delta u^2 + 2\,\delta v^2 + 2\,\delta x\,\delta u' + 2\,\delta y\,\delta v' \end{pmatrix} \end{aligned}
defining 1/\widetilde{m}=1/m_1+1/m_2
\begin{aligned} \delta u' = u'_1-u'_2 & = \dfrac{x_1\lambda_1 + (x_1-x_2)\lambda_2}{m_1} - \dfrac{(x_2-x_1)\lambda_2}{m_2} \\ & = \dfrac{x_1\lambda_1}{m_1} + \delta x\lambda_2\left(\dfrac{1}{m_1}+\dfrac{1}{m_2}\right) \\ & = \dfrac{x_1\lambda_1}{m_1} + \dfrac{\delta x\lambda_2}{\widetilde{m}} \\ \delta v' = v'_1-v'_2 & = \dfrac{y_1\lambda_1 + (y_1-y_2)\lambda_2}{m_1} - g - \dfrac{(y_2-y_1)\lambda_2}{m_2} + g \\ & = \dfrac{y_1\lambda_1}{m_1} + \delta y\lambda_2\left(\dfrac{1}{m_1}+\dfrac{1}{m_2}\right) \\ & = \dfrac{y_1\lambda_1}{m_1} + \dfrac{\delta y\lambda_2}{\widetilde{m}} \end{aligned}
The new constraint become
\begin{aligned} \text{[alg3]} & = \begin{pmatrix} \text{[alg3]}_1\\[0.5em] \text{[alg3]}_2 \end{pmatrix} \\ \text{[alg3]}_1 & = 2(u_1^2+v_1^2) + \dfrac{2}{m_1}\lambda_1(x_1^2+y_1^2) + \dfrac{2}{m_1}\lambda_2(x_1\delta x+y_1\delta y) - 2gy_1 \\ \text{[alg3]}_2 & = 2(\delta u^2+\delta v^2) + \dfrac{2\lambda_1}{m_1}(x_1\delta x+y_1\delta y) + \dfrac{2\lambda_2}{\widetilde{m}} \left(\delta x^2 + \delta y^2\right) \end{aligned}
DAE reduced to index 1 DAE
Remark 26.2. At this point \lambda_1 and \lambda_2 appears. Thus after another derivation the DAE is reduced to an ODE and, thus, the index is 3.
The derivation produce a complex ODE and is not reccomeneded. Instead, \lambda_1 and \lambda_2 can be solved, or better written as a function of position and velocity.
Elimination of the multiplier
Rewrite \text{[alg3]} as a linear system
\dfrac{1}{m_1} \begin{pmatrix} x_1^2+y_1^2 & x_1\delta x+y_1\delta y \\[0.5em] x_1\delta x+y_1\delta y & \dfrac{m_1}{\widetilde{m}} (\delta x^2 + \delta y^2) \end{pmatrix} \begin{pmatrix} \lambda_1 \\[0.5em] \lambda_2 \end{pmatrix} = \begin{pmatrix} gy_1-(u_1^2+v_1^2) \\[0.5em] -(\delta u^2+\delta v^2) \end{pmatrix}
the solution of the linear system are the function
\begin{pmatrix} \lambda_1 \\[0.5em] \lambda_2 \end{pmatrix} = m_1 \begin{pmatrix} x_1^2+y_1^2 & x_1\delta x+y_1\delta y \\[0.5em] x_1\delta x+y_1\delta y & \dfrac{m_1}{\widetilde{m}} (\delta x^2 + \delta y^2) \end{pmatrix}^{-1} \begin{pmatrix} gy_1-(u_1^2+v_1^2) \\[0.5em] -(\delta u^2+\delta v^2) \end{pmatrix}
The final ODE is
\left\{ \begin{aligned} x'_1 & = u_1\\ y'_1 & = v_1\\ x'_2 & = u_2\\ y'_2 & = v_2\\ u'_1 & = \dfrac{x_1}{m_1}\lambda_1(x_1,y_1,x_2,y_2,u_1,v_1,u_2,v_2)+ \dfrac{(x_1-x_2)}{m_1}\lambda_2(x_1,y_1,x_2,y_2,u_1,v_1,u_2,v_2) \\ v'_1 & = \dfrac{y_1}{m_1} \lambda_1(x_1,y_1,x_2,y_2,u_1,v_1,u_2,v_2)+ \dfrac{(y_1-y_2)}{m_1}\lambda_2(x_1,y_1,x_2,y_2,u_1,v_1,u_2,v_2) - g \\ u'_2 & = \dfrac{(x_2-x_1)}{m_2}\lambda_2(x_1,y_1,x_2,y_2,u_1,v_1,u_2,v_2) \\ v'_2 & = \dfrac{(y_2-y_1)}{m_2}\lambda_2(x_1,y_1,x_2,y_2,u_1,v_1,u_2,v_2) - g \end{aligned} \right.
| Parameter | Parameter | Parameter |
|---|---|---|
| m_1=5 | m_2=1 | \ell_1=1 |
| \ell_2=1 | g=9.81 | \theta_1=\pi/2 |
| \theta_2=\pi/4 | \omega_1=0 | \omega_2=-0.1 |
| x_1(0)=\ell_1 \sin(\theta_1) | y_1(0)=-\ell_1 \cos(\theta_1) | x_2(0)=x_1(0) + \ell_2 \sin(\theta_2) |
| y_2(0)=x_2(0) - \ell_2 \cos(\theta_2) | u_1(0)=0 | v_1(0)=0 |
| u_2(0)=0 | v_2(0)=0 |
Baumgarte Stabilization
To stabilize the index-3 DAE system of the double pendulum and avoid constraint drift during numerical integration, we introduce Baumgarte stabilization, a method that replaces the original constraints with second-order damped versions of them.
We define the stabilized constraints by enforcing a second-order differential equation for each component:
\dfrac{\mathrm{d}^2\text{[alg1]}}{\mathrm{d}t^2}+ 2\eta\omega \dfrac{\mathrm{d}\text{[alg1]}}{\mathrm{d}t} + \omega^2\text{[alg1]} = \bm{0}
or equivalently
\text{[alg3]}+ 2\eta\omega\text{[alg2]} + \omega^2\text{[alg1]} = \bm{0}
that become
\begin{aligned} 2(u_1^2+v_1^2) + \dfrac{2}{m_1}\lambda_1(x_1^2+y_1^2) + \dfrac{2}{m_1}\lambda_2(x_1\delta x+y_1\delta y) - 2gy_1 & + \\[0.5em] {\color{blue}2\eta\omega}(2x_1u_1+ 2y_1v_1) + {\color{blue}\omega^2}(x_1^2 + y_1^2 - \ell_1^2) & = 0 \\[2em] 2(\delta u^2+\delta v^2) + \dfrac{2\lambda_1}{m_1}(x_1\delta x+y_1\delta y) + \dfrac{2\lambda_2}{\widetilde{m}} \left(\delta x^2 + \delta y^2\right) & +\\[0.5em] {\color{blue}2\eta\omega}(2\delta x\,\delta u + 2\delta y\,\delta v) + {\color{blue}\omega^2}(\delta x^2 + \delta y^2 - \ell_2^2) & = 0 \end{aligned}
that become the linear system
\dfrac{1}{m_1} \begin{pmatrix} x_1^2+y_1^2 & x_1\delta x+y_1\delta y \\[0.5em] x_1\delta x+y_1\delta y & \dfrac{m_1}{\widetilde{m}} (\delta x^2 + \delta y^2) \end{pmatrix} \begin{pmatrix} \lambda_1 \\[0.5em] \lambda_2 \end{pmatrix} = \begin{pmatrix} r_1 \\[0.5em] r_2 \end{pmatrix}
where
\begin{aligned} r_1 &= gy_1-(u_1^2+v_1^2) -{2\color{blue}\eta\omega}(x_1u_1+y_1v_1) -\frac{1}{2}{\color{blue}\omega^2}(x_1^2 + y_1^2 - \ell_1^2)\\ r_2 &= -(\delta u^2+\delta v^2) -{2\color{blue}\eta\omega}(\delta x\,\delta u + \delta y\,\delta v) - \frac{1}{2}{\color{blue}\omega^2}(\delta x^2 + \delta y^2 - \ell_2^2) \end{aligned}
Parameters for the double pendulum solution of the reduced DAE with stabilization
| Parameter | Parameter |
|---|---|
| m_1=5 | m_2=1 |
| \ell_1=1 | \ell_2=1 |
| g=9.81 | \theta_1=\pi/2 |
| \theta_2=\pi/4 | \omega_1=0 |
| \omega_2=-0.1 | x_1(0)=\ell_1 \sin(\theta_1) |
| y_1(0)=-\ell_1 \cos(\theta_1) | x_2(0)=x_1(0) + \ell_2 \sin(\theta_2) |
| y_2(0)=x_2(0) - \ell_2 \cos(\theta_2) | u_1(0)=0 |
| v_1(0)=0 | u_2(0)=0 |
| v_2(0)=0 | \eta=0.7 |
| \Delta t=0.05 | \omega=2.4/\Delta t |
26.3 Index reduction using dummy variable
Variable classification
- Differential variables: x_1, y_1, x_2, y_2, u_1, v_1, u_2, v_2
- Algebraic variables: \lambda_1, \lambda_2.
We attempt to solve for \lambda_1 and \lambda_2, but the system is not yet solvable explicitly. We proceed by differentiating.
Step 1: Differentiate algebraic constraints
Now we must derivate the constraints in order to reduce the index. Here we want to differentiate only the algebraic equation that can reduce the index and do not introduce new differential equation. Here we want to avoid to introduce derivatove of the algebraic variable unless is strictly necessary.
We construct a dependency table to identify which constraints to differentiate:
\begin{array}{c|cc|} & \lambda_1 & \lambda_2 \\[0.5em] \hline \\ \text{(Constraint 1)} & 0 & 0 \\[1em] \text{(Constraint 2)} & 0 & 0 \end{array}
a 0 means the the function do not depend on the variable on the column. A 1 means dependence. Both constraints are candidate for differentiation as it doesn’t introduce new algebraic variables. Differentiate thus the constraints and substitute known variable in differential form using the ODE part
\begin{aligned} \text{(Constraint 1)}' & = 2x_1x'_1 + 2y_1y'_1 = 2(x_1u_1+y_1v_1) \\ \text{(Constraint 2)}' & = 2(x_1-x_2)(x'_1-x'_2) + 2(y_1-y_2)(y'_1-y'_2) \\ & = 2((x_1-x_2)(u_1-u_2)+(y_1-y_2)(v_1-v_2)) \end{aligned}
Substitute known differential variables and use shortcuts \delta x = x_1-x_2, \delta y = y_1-y_2, \delta u = u_1-u_2, \delta v = v_1-v_2 so that the DAE become: \left\{ \color{blue} \begin{aligned} x'_1 & = u_1\\ y'_1 & = v_1\\ x'_2 & = u_2\\ y'_2 & = v_2\\ u'_1 & = \dfrac{x_1\lambda_1 + (x_1-x_2)\lambda_2}{m_1} \\ v'_1 & = \dfrac{y_1\lambda_1 + (y_1-y_2)\lambda_2}{m_1} - g \\ u'_2 & = \dfrac{(x_2-x_1)\lambda_2}{m_2} \\ v'_2 & = \dfrac{(y_2-y_1)\lambda_2}{m_2} - g \end{aligned} \right. \quad \bigoplus \quad \left\{ \color{DarkRed} \begin{aligned} 0 & = 2(x_1u_1+y_1v_1) && (\text{Constraint 1})'\\ 0 & = 2(\delta x \delta u+\delta y \delta v) && (\text{Constraint 2})'\\[1em] \color{DarkGreen}\ell_1^2 &\color{DarkGreen} = x_1^2 + y_1^2 &&\color{DarkGreen} (\text{Constraint 1})\\ \color{DarkGreen}\ell_2^2 &\color{DarkGreen} = \delta x^2 + \delta y^2 &&\color{DarkGreen} (\text{Constraint 2}) \end{aligned} \right.
The red algebraic equations are the new algebraic equation after the derivation. The green algebraic equation is the hidden constraints that is removed from the DAE system but that is satisfied by choosing the right initial conditions.
The red part cannot determine the algebraic variable and must be differentiated again.
Step 2: Differentiate (again) the algebraic constraints
To find the right candidate(s) here we write a table of the dependencies \begin{array}{c|cc} & \lambda_1 & \lambda_2 \\[1em] \hline \\[0.1em] \text{(Constraint 1)}' & 0 & 0 \\[1em] \text{(Constraint 2)}' & 0 & 0 \end{array}
Both constraints are candidate for differentiation as it doesn’t introduce new algebraic variables. Differentiate thus the constraints and substitute known variable in differential form using the ODE part
Step 3: Differentiate algebraic constraints
\color{DarkGreen} \begin{aligned} \text{(Constraint 1)}'' & = 2x'_1 u_1 + 2 y'_1 v_1+ 2x_1 + 2 y_1 v'_1 =0\\ \text{(Constraint 2)}'' &= 2\delta x' \delta u + 2 \delta y' \delta v + 2\delta x \delta u' + 2 \delta y \delta v' =0 \end{aligned}
Substitute known differential variables
\begin{aligned} 2(u_1^2+v_1^2)+ 2\dfrac{x_1^2\lambda_1 + x_1\delta x\lambda_2}{m_1} + 2\dfrac{y_1^2\lambda_1 + y_1\delta y\lambda_2}{m_1} -2 y_1 g &=0\\ 2(\delta u^2 + \delta v^2) + 2\delta x \delta u' + 2 \delta y \delta v' &=0 \end{aligned}
using 1/\widetilde{m}=1/m_1+1/m_2
\begin{aligned} \delta u' & = \dfrac{x_1\lambda_1 + \delta x\lambda_2}{m_1} +\dfrac{\delta x\lambda_2}{m_2} = \dfrac{x_1}{m_1}\lambda_1+ \dfrac{\delta x}{\widetilde{m}}\lambda_2 \\ \delta v' & = \dfrac{y_1\lambda_1 + \delta y\lambda_2}{m_1} + \dfrac{\delta y\lambda_2}{m_2} = \dfrac{y_1}{m_1}\lambda_1+ \dfrac{\delta y}{\widetilde{m}}\lambda_2 \end{aligned}
so that constraints after substitution become
\color{DarkGreen} \begin{aligned} \text{(Constraint 1)}'' &= 2(u_1^2+v_1^2)+ \dfrac{2}{m_1} ( (x_1^2+y_1^2)\lambda_1 + (x_1\delta x+y_1\delta y)\lambda_2) -2 y_1 g =0\\ \text{(Constraint 2)}'' &= 2(\delta u^2 + \delta v^2) + 2\delta x\left( \dfrac{x_1}{m_1}\lambda_1+ \dfrac{\delta x}{\widetilde{m}}\lambda_2 \right) + 2 \delta y\left(\dfrac{y_1}{m_1}\lambda_1+ \dfrac{\delta y}{\widetilde{m}}\lambda_2\right) =0 \end{aligned}
and now \lambda_1 and \lambda_2 appears linearly:
\begin{pmatrix} x_1^2+y_1^2 & x_1\delta x+y_1\delta y \\[1em] x_1\delta x+y_1\delta y & \dfrac{m_1}{\widetilde{m}}(\delta x^2+\delta y^2) \end{pmatrix} \begin{pmatrix} \lambda_1 \\[1em] \lambda_2 \end{pmatrix} = m_1 \begin{pmatrix} y_1 g - (u_1^2+v_1^2) \\[1em] - (\delta u^2 + \delta v^2) \end{pmatrix}
The DAE of index 1
Thus DAE now is of index 1 and can be written as: \left\{ \color{blue} \begin{aligned} x'_1 & = u_1\\ y'_1 & = v_1\\ x'_2 & = u_2\\ y'_2 & = v_2\\ u'_1 & = \dfrac{x_1\lambda_1 + (x_1-x_2)\lambda_2}{m_1} \\ v'_1 & = \dfrac{y_1\lambda_1 + (y_1-y_2)\lambda_2}{m_1} - g \\ u'_2 & = \dfrac{(x_2-x_1)\lambda_2}{m_2} \\ v'_2 & = \dfrac{(y_2-y_1)\lambda_2}{m_2} - g \\ \color{DarkRed} \begin{pmatrix} \lambda_1 \\[1em] \lambda_2 \end{pmatrix} &\color{DarkRed}=\bm{A}(t,\bm{Z})^{-1} \bm{b}(t,\bm{Z}) \end{aligned} \right.
where
\bm{Z} = \begin{pmatrix} x_1 \\ y_1 \\ x_2 \\ y_2 \\ u_1 \\ v_1 \\ u_2 \\ v_2 \end{pmatrix}, \quad \begin{aligned} \bm{A}(t,\bm{Z}) &= \begin{pmatrix} x_1^2+y_1^2 & x_1\delta x+y_1\delta y \\[1em] x_1\delta x+y_1\delta y & \dfrac{m_1}{\widetilde{m}}(\delta x^2+\delta y^2) \end{pmatrix}, \\[4em] \bm{b}(t,\bm{Z})&= m_1 \begin{pmatrix} y_1 g - (u_1^2+v_1^2) \\[1em] - (\delta u^2 + \delta v^2) \end{pmatrix} \end{aligned}
Moreover using the right initial condition satisfy the hidden constraints.
\left\{ \color{DarkGreen} \begin{aligned} 0&\color{DarkGreen} = x_1^2 + y_1^2 -\ell_1^2& (\text{Constraint 1})\\ 0&\color{DarkGreen} = \delta x^2 + \delta y^2-\ell_2^2 & (\text{Constraint 2})\\ 0 & = 2(x_1u_1+y_1v_1) & (\text{Constraint 1})'\\ 0 & = 2(\delta x \delta u+\delta y \delta v) & (\text{Constraint 2})' \end{aligned} \right.
Baumgarte Stabilization
For this DAE the hidden constraints can be used to stabilize the algebric constraints: \text{(Constraint 1)} and \text{(Constraint 2)}
\begin{aligned} \color{Olive}\text{(Constraint 1)}_{\text{stab}} &= \text{(Constraint 1)}'' + {\color{DarkOrange}2\eta_1\omega_1} \text{(Constraint 1)}' + {\color{DarkOrange}\omega_1^2} \text{(Constraint 1)} = 0 \\ \color{Olive}\text{(Constraint 2)}_{\text{stab}} &= \text{(Constraint 2)}'' + {\color{DarkOrange}2\eta_2\omega_2} \text{(Constraint 2)}' + {\color{DarkOrange}\omega_2^2} \text{(Constraint 2)} = 0 \end{aligned}
and substituting (diving also by 2)
\begin{aligned} \color{Olive}\text{(Constraint 1)}_{\text{stab}} &= u_1^2+v_1^2+ \dfrac{(x_1^2+y_1^2)\lambda_1 + (x_1\delta x+y_1\delta y)\lambda_2}{m_1} -y_1 g +{\color{DarkOrange}(\star)_1} \\ \color{Olive}\text{(Constraint 2)}_{\text{stab}} &= \delta u^2 + \delta v^2 + \dfrac{(x_1\delta x+y_1\delta y)\lambda_1}{m_1}+ \dfrac{(\delta x^2+\delta y^2)\lambda_2}{\widetilde{m}} + {\color{DarkOrange}(\star)_2} \\ \color{DarkOrange}(\star)_1 &= {\color{DarkOrange}2\eta_1\omega_1}(x_1u_1+y_1v_1)+ {\color{DarkOrange}\omega_1^2}\dfrac{x_1^2 + y_1^2 -\ell_1^2}{2} \\ \color{DarkOrange}(\star)_2 &= {\color{DarkOrange}2\eta_2\omega_2}(\delta x \delta u+\delta y \delta v)+ {\color{DarkOrange}\omega_2^2}\dfrac{\delta x^2 + \delta y^2-\ell_2^2}{2} \end{aligned}
and using the stabilized constraints
\bm{b}_{\text{stab}}(t,\bm{Z})= m_1 \begin{pmatrix} y_1 g - (u_1^2+v_1^2) - {\color{DarkOrange}(\star)_1}\\[1em] - (\delta u^2 + \delta v^2) - {\color{DarkOrange}(\star)_2} \end{pmatrix}