24 Example 2
24.1 3D Pendulum
Given the following definitions:
\begin{aligned} T &= \frac{m}{2}(v_x^2 + v_y^2 + v_z^2) \quad && \text{(Kinetic energy)} \\[0.5em] V &= -mgz \quad && \text{(Potential energy)} \\[0.5em] \Phi(x, y, z) &= x^2 + y^2 + z^2 - \ell^2 = 0 \quad && \text{(Holonomic constraint)} \end{aligned}
The constraint enforces that the pendulum mass lies on the surface of a sphere of radius \ell, centered at the origin.
The corresponding system of Differential-Algebraic Equations (DAE) is:
\left\{ \begin{aligned} x'(t) &= v_x(t)\\ y'(t) &= v_y(t)\\ z'(t) &= v_z(t)\\ m\, v_x'(t) &= 2\, x(t)\,\lambda(t)\\ m\, v_y'(t) &= 2\, y(t)\,\lambda(t)\\ m\, v_z'(t) &= 2\, z(t)\,\lambda(t)-mg\\ x^2+y^2+z^2 &= \ell^2 \end{aligned} \right.
This is a DAE of index 3 due to the algebraic constraint \Phi(x, y, z) = 0, and the fact that the constraint must be differentiated three times to obtain an explicit expression for the highest derivative (\lambda').
We define
\bm{Z} = \begin{pmatrix} x \\[0.5em] y \\[0.5em] z \\[0.5em] v_x \\[0.5em] v_y \\[0.5em] v_z \\[0.5em] \lambda \end{pmatrix}, \quad \bm{E}(\bm{Z},t) = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\[0.5em] 0 & 1 & 0 & 0 & 0 & 0 & 0 \\[0.5em] 0 & 0 & 1 & 0 & 0 & 0 & 0 \\[0.5em] 0 & 0 & 0 & 1 & 0 & 0 & 0 \\[0.5em] 0 & 0 & 0 & 0 & 1 & 0 & 0 \\[0.5em] 0 & 0 & 0 & 0 & 0 & 1 & 0 \\[0.5em] 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \quad \bm{G}(\bm{Z},t) = \begin{pmatrix} v_x(t)\\[0.75em] v_y(t)\\[0.75em] v_z(t)\\[0.75em] \dfrac{2}{m} x(t)\,\lambda(t)\\[0.75em] \dfrac{2}{m} y(t)\,\lambda(t)\\[0.75em] \dfrac{2}{m} z(t)\,\lambda(t)-g\\[0.75em] x^2+y^2+z^2-\ell^2 \end{pmatrix}
so that the DAE can be written as
\bm{E}(\bm{Z},t)\,\bm{Z}^\prime(t) = \bm{G}(\bm{Z},t)
24.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}{rrrrrrr|c} 1 & 0 & 0 & 0 & 0 & 0 & 0 & v_x(t)\\[0.25em] 0 & 1 & 0 & 0 & 0 & 0 & 0 & v_y(t)\\[0.25em] 0 & 0 & 1 & 0 & 0 & 0 & 0 & v_z(t)\\[0.25em] 0 & 0 & 0 & 1 & 0 & 0 & 0 & \frac{2}{m} x(t)\,\lambda(t)\\[0.25em] 0 & 0 & 0 & 0 & 1 & 0 & 0 & \frac{2}{m} y(t)\,\lambda(t)\\[0.25em] 0 & 0 & 0 & 0 & 0 & 1 & 0 & \frac{2}{m} z(t)\,\lambda(t)-g\\[0.25em] 0 & 0 & 0 & 0 & 0 & 0 & 0 & x^2+y^2+z^2-\ell^2 \end{array} \right]
is used to memorize the operation done on the DAE. Let \textrm{index}\leftarrow 0.
The algebraic constraint is
\color{blue} \textrm{[alg1]} = x^2+y^2+z^2-\ell^2 = 0
Reduce by one the index by differentiating \textrm{[alg1]} (\textrm{index}\leftarrow =1)
\color{blue} \textrm{[alg2]} = \dfrac{\mathrm{d}\textrm{[alg1]}}{\mathrm{d}t} = 2\left(xx'+yy'+zz'\right) = 2\left(xv_x+yv_y+zv_z\right)
Reduce by one the index by differentiating \textrm{[alg2]} (\textrm{index}\leftarrow =2)
\begin{aligned} {\color{blue} \textrm{[alg3]}} = \dfrac{\mathrm{d}\textrm{[alg2]}}{\mathrm{d}t} &= 2\left(x'v_x+y'v_y+z'v_z\right)+ 2\left(xv'_x+yv'_y+zv'_z\right)\\ &= \color{blue}2\left(v_x^2+v_y^2+v_z^2\right)+ \frac{4\lambda}{m}\left(x^2+y^2+z^2\right)-2gz \end{aligned}
Now \lambda can be solved from \textrm{[alg3]}=0
\color{DarkRed} \lambda = \frac{m}{2}\dfrac{gz-\left(v_x^2+v_y^2+v_z^2\right)}{x^2+y^2+z^2}
DAE reduced to index 1 DAE
using this solution the index 1 DAE (after the reduction) can be reduced to an ODE by substituting \lambda as follows:
\left\{ \begin{aligned} x' &= v_x\\ y' &= v_y\\ z' &= v_z\\ v_x' &= x\left[\dfrac{gz-\left(v_x^2+v_y^2+v_z^2\right)}{x^2+y^2+z^2}\right]\\ v_y' &= y\left[\dfrac{gz-\left(v_x^2+v_y^2+v_z^2\right)}{x^2+y^2+z^2}\right]\\ v_z' &= z\left[\dfrac{gz-\left(v_x^2+v_y^2+v_z^2\right)}{x^2+y^2+z^2}\right]-mg \end{aligned} \right.
DAE reduced to ODE
If you want to include also \lambda in the reduced ODE derive \textrm{[alg3]}=0 so that \textrm{index}\leftarrow =3
\begin{aligned} \dfrac{\mathrm{d}\textrm{[alg3]}}{\mathrm{d}t} &= \dfrac{\mathrm{d}}{\mathrm{d}t} \left[ 2\left(v_x^2+v_y^2+v_z^2\right)+ \frac{4\lambda}{m}\left(x^2+y^2+z^2\right)-2gz \right]\\ &= 4\left(v_xv'_x+v_yv'_y+v_zv'_z\right)+ \frac{4\lambda'}{m}\left(x^2+y^2+z^2\right)\\ &+\frac{8\lambda}{m}\left(xx'+yy'+zz'\right) -2gz'\\ &= \dfrac{8\lambda}{m}\left(xv_x+yv_y+zv_z\right)-4v_zg+ \frac{4\lambda'}{m}\left(x^2+y^2+z^2\right)\\ &+\frac{8\lambda}{m}\left(xv_x+yv_y+zv_z\right) -2v_zg\\ &= \dfrac{16\lambda}{m}\left(xv_x+yv_y+zv_z\right)-6v_zg+ \frac{4\lambda'}{m}\left(x^2+y^2+z^2\right) = 0 \end{aligned}
Solving for \lambda'
\lambda' = \frac{1}{2} \dfrac{3\,m\,v_z\,g-8\lambda\left(xv_x+yv_y+zv_z\right)}{x^2+y^2+z^2}
The DAE reduced to the ODE become
\left\{ \begin{aligned} x' &= v_x\\[0.5em] y' &= v_y\\[0.5em] z' &= v_z\\[0.5em] v_x' &= \dfrac{2}{m} x\lambda\\[0.5em] v_y' &= \dfrac{2}{m} y\lambda\\[0.5em] v_z' &= \dfrac{2}{m} z\lambda-g \\[0.5em] \lambda' &= \frac{1}{2} \dfrac{3\,m\,v_z\,g-8\lambda\left(xv_x+yv_y+zv_z\right)}{x^2+y^ 2+z^2} \end{aligned} \right.
with the hidden constraints
\color{blue} \begin{aligned} \textrm{[alg1]}\qquad && x^2+y^2+z^2-\ell^2 &= 0 \\[0.5em] \textrm{[alg2]}\qquad && 2\left(xv_x+yv_y+zv_z\right) &= 0 \\[0.5em] \textrm{[alg3]}\qquad && 2\left(v_x^2+v_y^2+v_z^2\right)+ \frac{4\lambda}{m}\left(x^2+y^2+z^2\right)-2gz &= 0 \end{aligned}
24.3 Index reduction using dummy variable
Variable classification
- Differential variables: x, y, z, v_x, v_y, v_z
- Algebraic variables: \lambda
We attempt to solve for \lambda, the system is not yet solvable explicitly. We proceed by differentiating.
Step 1: Differentiate algebraic constraints
Differentiate \text{(Constraint 1)}=x^2+y^2+z^2-\ell^2:
\dfrac{\mathrm{d}\text{(Constraint 1)}}{\mathrm{d}t}= 2xx'+2yy'+2zz' = 0 \tag{24.1}
Substitute known differential variables:
\color{DarkRed}\text{(Constraint 2)} = 2xv_x+2yv_y+2zv_z = 0
Still no \lambda. We differentiate once more.
Step 2: Differentiate algebraic constraint
Differentiate Equation 24.1:
\dfrac{\mathrm{d}\text{(Constraint 2)}}{\mathrm{d}t} = 2x'v_x+2y'v_y+2z'v_z+ 2xv'_x+2yv'_y+2zv'_z = 0
Substitute known derivatives:
\color{DarkRed} \text{(Constraint 3)}= 2(v_x^2+v_y^2+v_z^2)+ \dfrac{4\lambda}{m} (x^2+y^2+z^2) - 2 z g = 0
Now \lambda appears in the DAE. We can solve for \lambda
\lambda = \dfrac{m}{2} \cdot \dfrac{z g -(v_x^2+v_y^2+v_z^2)}{x^2+y^2+z^2}
now we can substitute \lambda in the DAE, but this means reduce the DAE by eliminating a variable. Better define
\widetilde{\lambda}(x,y,z,v_x,v_y,v_z) = \dfrac{m}{2} \cdot \dfrac{z g -(v_x^2+v_y^2+v_z^2)}{x^2+y^2+z^2}
and obtain the index 1 DAE
\left\{ \begin{aligned} \color{blue}x' &\color{blue}= v_x\\ \color{blue}y' &\color{blue}= v_y\\ \color{blue}z' &\color{blue}= v_z\\ \color{blue}v_x' &\color{blue}= \dfrac{2}{m} x\,\lambda\\ \color{blue}v_y' &\color{blue}= \dfrac{2}{m} y\,\lambda\\ \color{blue}v_z' &\color{blue}= \dfrac{2}{m} z\,\lambda-g\\ \color{DarkRed}\lambda & \color{DarkRed}=\widetilde{\lambda}(x,y,z,v_x,v_y,v_z) \end{aligned} \right.
having done two differetiation to obtain an index 1 DAE he original DAE was of index 3.