27 Example 5
27.1 A Constrained Van der Pol System and a Dummy-Variable Reformulation
The Van der Pol equation x'' - \mu(1-x^2)x' + x = 0 represents a classic example of a nonlinear oscillator with nonlinear damping, widely used in the study of dynamical systems and electrical circuits.
27.2 Formulation as a Stiff DAE
To construct a constrained DAE system [1], we introduce:
- The variable y = x' to obtain a first-order system
- An additional state z
The resulting DAE is: \left\{ \begin{aligned} x' &= y \\ y' &= \mu(1-x^2)y - x \\ 0 &= x^2 - z^2 + 5 \end{aligned} \right.
In this original formulation, x and y are differential variables and z is algebraic. Since \partial(x^2-z^2+5)/\partial z=-2z\neq0 on the constraint manifold z^2=x^2+5, this is a regular semi-explicit index-1 DAE. The dummy-variable construction below is used to illustrate how one obtains an ODE representation together with its hidden invariant; it does not change the index of the original DAE.
Dummy-Variable Reformulation
Introduce dummy variables y_{\text{dot}}, z_{\text{dot}}. In this augmented representation, the system decomposes into differential (blue) and algebraic (red) components:
\left\{ \color{blue} \begin{aligned} x' & = y\\ y' & = y_{\text{dot}} \\ z' & = z_{\text{dot}} \end{aligned} \right. \quad \bigoplus \quad \left\{ \color{DarkRed} \begin{aligned} 0 & = {\color{blue}y_{\text{dot}}} - \mu(1-x^2)y + x & \text{(Constraint 1)}\\ 0 & = x^2 - z^2 + 5& \text{(Constraint 2)} \end{aligned} \right.
For the augmented system we have:
- Differential variables: x, y, z.
- Algebraic variables: y_{\text{dot}}, z_{\text{dot}}.
Index Reduction Procedure
The equations determine y_{\text{dot}} but not z_{\text{dot}}. To obtain the ODE representation, we differentiate the geometric constraint once.
Step 1: Differentiate algebraic constraints
We differentiate only the geometric constraint: this introduces z_{\text{dot}} without differentiating an algebraic variable again.
We construct a dependency table to identify which constraints to differentiate:
\begin{array}{c|cc} & y_{\text{dot}} & z_{\text{dot}} \\[1em] \hline \\[0.1em] \text{(Constraint 1)} & 1 & 0 \\[1em] \text{(Constraint 2)} & 0 & 0 \end{array}
A 0 means that the function does not depend on the variable in the column. A 1 means dependence. The last constraint is the candidate for differentiation because it does not introduce a new algebraic variable. Differentiate it and substitute known derivatives from the differential part:
\text{(Constraint 2)}' = 2xx' - 2 z z' = 2 x y - 2 z z_{\text{dot}}
The augmented DAE becomes:
\left\{ \color{blue} \begin{aligned} x' & = y\\ y' & = y_{\text{dot}} \\ z' & = z_{\text{dot}} \end{aligned} \right. \quad \bigoplus \quad \left\{ \color{DarkRed} \begin{aligned} 0 & = y_{\text{dot}} - \mu(1-x^2)y + x & \text{(Constraint 1)}\\ 0 & = 2 x y - 2 z z_{\text{dot}} & \text{(Constraint 2)}' \\[1em] \color{DarkGreen} 0 & \color{DarkGreen} = x^2 - z^2 + 5 & \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} & y_{\text{dot}} & z_{\text{dot}} \\[1em] \hline \\[0.1em] \text{(Constraint 1)} & 1 & 0 \\[1em] \text{(Constraint 2)}' & 0 & 1 \end{array}
The table shows that the nonlinear system is structurally non singular and algebraic constraints can be solved. From \text{(Constraint 2)}' solve for z_{\text{dot}} and from \text{(Constraint 1)} solve for y_{\text{dot}}
z_{\text{dot}} = \dfrac{xy}{z}, \qquad y_{\text{dot}} = \mu(1-x^2)y - x
The ODE
Thus DAE now can be reduced to an ODE and can be written as:
\left\{ \color{blue} \begin{aligned} x' & = y\\ y' & = y_{\text{dot}} = \mu(1-x^2)y - x \\ z' & = z_{\text{dot}} = \dfrac{xy}{z} \end{aligned} \right. \quad \bigoplus \quad \color{DarkGreen} x^2 - z^2 + 5=0,\quad\text{(Constraint 1)}
27.3 Baumgarte Stabilization
Baumgarte stabilization is a classical technique used to prevent this drift. Instead of enforcing a constraint like:
z(t) = 0
we instead modify it to a differential equation:
\frac{dz}{dt} + \frac{1}{\tau} z = 0
This form ensures that any initial deviation from the constraint decays exponentially with a rate controlled by the parameter \tau. The smaller \tau is, the stronger the correction applied.
The exact solution is:
z(t) = z_0 \, e^{-t/\tau}
In the following example, we simulate the stabilized equation for different values of \tau and visualize how quickly the constraint error decays.
The plot clearly shows the impact of the stabilization parameter \tau:
- Small \tau (e.g., 0.01): The solution quickly decays to zero, meaning constraint violations are rapidly corrected. However, this can introduce numerical stiffness.
- Large \tau (e.g., 5): The decay is slow, allowing more deviation from the constraint over time.
This trade-off is at the heart of Baumgarte stabilization: small \tau improves accuracy but increases stiffness and numerical cost, while large \tau is easier to integrate but less effective at controlling constraint drift.
In real DAE systems, choosing \tau often involves tuning based on system stiffness and desired constraint enforcement.
For this DAE the hidden constraints can be used to stabilize the algebraic constraint \text{(Constraint 1)}
\color{Olive}\text{(Constraint 1)}_{\text{stab}} = \text{(Constraint 1)}' + {\color{DarkOrange}\frac{1}{\tau}} \text{(Constraint 1)}= 0
and substituting (diving also by 2)
\color{Olive}\text{(Constraint 1)}_{\text{stab}} = x\,y-z\,z_{\text{dot}}+ {\color{DarkOrange}\frac{1}{\tau}}\dfrac{x^2-z^2+5}{2}
and using the stabilized constraints
\left\{ \color{blue} \begin{aligned} x' & = y\\ y' & = \mu(1-x^2)y - x \\ z' & = \dfrac{xy}{z}+{\color{DarkOrange}\frac{1}{\tau}\dfrac{x^2-z^2+5}{2z}} \end{aligned} \right. \quad \bigoplus \quad \color{DarkGreen} x^2 - z^2 + 5=0,\quad\text{(Constraint 1)}