10 Semi-Explicit Form and Its Classification
Any DAE \bm F(\bm x,\dot{\bm x},t)=\bm0 can be brought to the semi-explicit form
\dot{\bm x}=\bm y, \qquad \bm 0=\bm F(\bm x,\bm y,t),
by the single substitution \bm y=\dot{\bm x}. The transformation is free, and it is worth making because the semi-explicit form exposes immediately which equations are differential and which are algebraic — and, as this chapter shows, allows the index to be determined by inspecting which variables each algebraic equation depends on, without performing any reduction.
10.1 From an Implicit DAE to Semi-Explicit Form
An implicit differential-algebraic equation (DAE) is written as
\bm F\bigl(\bm x(t),\bm x'(t),t\bigr)=\bm 0, \qquad \bm x(0)=\bm x_0,
where \bm x(t)\in\mathbb R^n. If the Jacobian \bm F_{\bm x'} is nonsingular at a solution point, the implicit function theorem gives, locally,
\bm x'(t)=\bm f\bigl(\bm x(t),t\bigr).
The system is then locally an ODE. When \bm F_{\bm x'} is singular, the equation cannot in general be solved for the full derivative; the remaining relations constrain the state and lead to a DAE. Singular Jacobians alone do not determine an index: the derivative constraints and their rank conditions also matter.
Semi-explicit form
A common structure separates differential variables \bm x(t)\in\mathbb R^n from algebraic variables \bm z(t)\in\mathbb R^m:
\begin{aligned} \bm x'(t) &= \bm f\bigl(\bm x(t),\bm z(t),t\bigr),\\ \bm 0 &= \bm g\bigl(\bm x(t),\bm z(t),t\bigr). \end{aligned}
The first equation supplies the differential dynamics; the second imposes an algebraic constraint. This is called semi-explicit form. The variables \bm z need not be constant: their values are determined at each time by the constraint and the evolving differential state.
From an implicit DAE to semi-explicit form
If the equations can be partitioned so that some equations explicitly define the derivatives of selected variables and the remaining equations contain no derivatives, they can be arranged in semi-explicit form by reordering variables and equations. For a general equation
\bm F\bigl(\bm x(t),\bm x'(t),t\bigr)=\bm 0,
one can also introduce an auxiliary variable \bm v=\bm x' and write
\begin{aligned} \bm x' &= \bm v,\\ \bm 0 &= \bm F(\bm x,\bm v,t). \end{aligned}
This is a semi-explicit representation in the broad sense, but it may not be the usual differential/algebraic split with an invertible algebraic Jacobian. Under the regularity assumptions for this order-lifting construction, it typically raises the differentiation index by one. It is therefore not an index-preserving simplification in general.
The index-one Jacobian condition
For the semi-explicit system, suppose \bm g is continuously differentiable and the square Jacobian
\bm g_{\bm z}(\bm x,\bm z,t) =\frac{\partial\bm g}{\partial\bm z}
is nonsingular at a point satisfying the constraint. The implicit function theorem then gives a local function \bm z=\bm\phi(\bm x,t) such that \bm g(\bm x,\bm\phi(\bm x,t),t)=\bm 0. Substitution reduces the differential equation locally to
\bm x'(t)=\bm f\bigl(\bm x(t),\bm\phi(\bm x(t),t),t\bigr).
This regular semi-explicit case is called index one under the differentiation-index convention used in the next chapter. If \bm g_{\bm z} is singular, this argument does not apply; one may need to differentiate constraints to reveal additional relations. The definition and worked calculations are given in Section 9.1.
Theorem 10.1 (Implicit Function Theorem) Let \bm g:U\subset\mathbb R^{n+m+1}\to\mathbb R^m be continuously differentiable, and let (\bm x_0,\bm z_0,t_0)\in U satisfy \bm g(\bm x_0,\bm z_0,t_0)=\bm 0. If \bm g_{\bm z}(\bm x_0,\bm z_0,t_0) is nonsingular, then there are neighborhoods of (\bm x_0,t_0) and \bm z_0 and a unique continuously differentiable function \bm z=\bm\phi(\bm x,t) such that, locally,
\bm g(\bm x,\bm z,t)=\bm 0 \quad\Longleftrightarrow\quad \bm z=\bm\phi(\bm x,t).
Its derivatives satisfy
\bm\phi_{\bm x} =-\bm g_{\bm z}^{-1}\bm g_{\bm x}, \qquad \bm\phi_t=-\bm g_{\bm z}^{-1}\bm g_t,
with Jacobians evaluated along \bm z=\bm\phi(\bm x,t).
Example: an index-one system
Consider
\begin{aligned} x'(t)&=-y(t),\\ 0&=x(t)+y(t)-\sin t. \end{aligned}
Here g_y=1, so the constraint gives y=\sin t-x. The system reduces to
x'(t)=x(t)-\sin t, \qquad y(t)=\sin t-x(t).
Equivalently, differentiating the constraint once gives x'+y'-\cos t=0, and using x'=-y yields y'=y+\cos t. Thus one differentiation supplies an explicit ODE for (x,y), in agreement with the index-one Jacobian test.
10.2 Classification of Semi-Explicit DAEs
Once a DAE has been written in the semi-explicit form of Section 10.1, its index can be read off from which variables the algebraic equations depend on, without performing the reduction. The mechanism is always the same: differentiate the constraint, look for the algebraic variable in the result, and if it is absent differentiate again.
Throughout this section \bm x denotes the differential variables and \bm y the algebraic ones, and all Jacobians are assumed to exist and to be continuous.
Index one
The constraint contains the algebraic variable:
\begin{aligned} \dot{\bm x} &= \bm f(\bm x,\bm y,t), \\ \bm 0 &= \bm g(\bm x,\bm y,t). \end{aligned} \tag{10.1}
Differentiating the constraint once and substituting \dot{\bm x}=\bm f,
\bm 0 =\frac{\partial\bm g}{\partial\bm x}\dot{\bm x} +\frac{\partial\bm g}{\partial\bm y}\dot{\bm y} +\frac{\partial\bm g}{\partial t} =\frac{\partial\bm g}{\partial\bm x}\bm f +\frac{\partial\bm g}{\partial\bm y}\dot{\bm y} +\frac{\partial\bm g}{\partial t}.
Definition 10.1 (Index-one condition) Equation 10.1 is a DAE of index one if
\boxed{\; \frac{\partial\bm g}{\partial\bm y}(\bm x,\bm y,t) \quad\text{is nonsingular.} \;} \tag{10.2}
Under Equation 10.2 the derivative of the constraint can be solved for \dot{\bm y},
\dot{\bm y} =\bm h(\bm x,\bm y,t) :=-\left(\frac{\partial\bm g}{\partial\bm y}\right)^{-1} \left(\frac{\partial\bm g}{\partial t} +\frac{\partial\bm g}{\partial\bm x}\bm f\right),
and one differentiation suffices:
\begin{aligned} \dot{\bm x} &= \bm f(\bm x,\bm y,t), \\ \dot{\bm y} &= \bm h(\bm x,\bm y,t), \end{aligned} \qquad\oplus\qquad \bm g(\bm x,\bm y,t)=\bm0 \ \ \text{(invariant)}.
By the implicit function theorem (Theorem 10.1), condition Equation 10.2 also guarantees the existence of \bm y=\bm\varphi(\bm x,t) with \bm g(\bm x,\bm\varphi(\bm x,t),t)\equiv\bm0, so that Equation 10.1 becomes the genuine ODE \dot{\bm x}=\bm f(\bm x,\bm\varphi(\bm x,t),t) on the constraint manifold. This is the trivial case of the coordinate partitioning technique of Coordinate Partitioning. It is of limited practical use, because \bm\varphi is rarely available in closed form; the derivative formula above is what an integrator actually uses.
Index two
The algebraic variable disappears from the constraint:
\begin{aligned} \dot{\bm x} &= \bm f(\bm x,\bm y,t), \\ \bm 0 &= \bm g(\bm x,t). \end{aligned} \tag{10.3}
Now \partial\bm g/\partial\bm y=\bm0 identically, so Equation 10.2 fails and one differentiation cannot produce \dot{\bm y}. It produces instead a new algebraic equation,
\bm 0 =\frac{\partial\bm g}{\partial\bm x}(\bm x,t)\,\bm f(\bm x,\bm y,t) +\frac{\partial\bm g}{\partial t}(\bm x,t) =:\bm h(\bm x,\bm y,t), \tag{10.4}
in which \bm y does appear, through \bm f. The pair \{\dot{\bm x}=\bm f,\ \bm0=\bm h\} has the shape of Equation 10.1, and is therefore of index one exactly when \partial\bm h/\partial\bm y is nonsingular. Differentiating Equation 10.4 with respect to \bm y gives
\frac{\partial\bm h}{\partial\bm y} =\frac{\partial\bm g}{\partial\bm x} \frac{\partial\bm f}{\partial\bm y},
the first instance of the pattern that governs the whole classification.
Definition 10.2 (Index-two condition) Equation 10.3 is a DAE of index two if
\boxed{\; \frac{\partial\bm g}{\partial\bm x}(\bm x,t)\, \frac{\partial\bm f}{\partial\bm y}(\bm x,\bm y,t) \quad\text{is nonsingular.} \;} \tag{10.5}
Two differentiations are then needed, and the reduced ODE carries two invariants: the original constraint \bm g=\bm0 and its first derivative \bm h=\bm0. Note that each is weaker than the previous one: a numerical trajectory can satisfy \bm h=\bm0 to machine precision and still drift away from \bm g=\bm0.
Index three
The constraint reaches the algebraic variable only through two intermediate differential equations:
\begin{aligned} \dot{\bm x} &= \bm f(\bm x,\bm y,\bm z,t), \\ \dot{\bm y} &= \bm g(\bm x,\bm y,t), \\ \bm 0 &= \bm h(\bm y,t), \end{aligned} \tag{10.6}
with \bm z algebraic. Differentiating the constraint,
\bm 0 =\frac{\partial\bm h}{\partial\bm y}\dot{\bm y} +\frac{\partial\bm h}{\partial t} =\frac{\partial\bm h}{\partial\bm y}\bm g(\bm x,\bm y,t) +\frac{\partial\bm h}{\partial t} =:\bm a(\bm x,\bm y,t), \tag{10.7}
which still does not contain \bm z. Collecting \bm w=(\bm x,\bm y) and \bm F(\bm w,\bm z,t)=\bigl(\bm f(\bm x,\bm y,\bm z,t),\, \bm g(\bm x,\bm y,t)\bigr), the system
\dot{\bm w}=\bm F(\bm w,\bm z,t), \qquad \bm0=\bm a(\bm w,t)
is precisely Equation 10.3, hence of index two when (\partial\bm a/\partial\bm w)(\partial\bm F/\partial\bm z) is nonsingular. Since
\frac{\partial\bm a}{\partial\bm w} =\left(\frac{\partial\bm a}{\partial\bm x}, \frac{\partial\bm a}{\partial\bm y}\right), \qquad \frac{\partial\bm F}{\partial\bm z} =\begin{pmatrix}\dfrac{\partial\bm f}{\partial\bm z}\\[0.6em]\bm0\end{pmatrix}, \qquad \frac{\partial\bm a}{\partial\bm x} =\frac{\partial\bm h}{\partial\bm y} \frac{\partial\bm g}{\partial\bm x},
the product collapses to a chain of three Jacobians.
Definition 10.3 (Index-three condition) Equation 10.6 is a DAE of index three if
\boxed{\; \frac{\partial\bm h}{\partial\bm y}(\bm y,t)\, \frac{\partial\bm g}{\partial\bm x}(\bm x,\bm y,t)\, \frac{\partial\bm f}{\partial\bm z}(\bm x,\bm y,\bm z,t) \quad\text{is nonsingular.} \;} \tag{10.8}
Hessenberg form
The three conditions Equation 10.2, Equation 10.5 and Equation 10.8 are the first terms of one chain. The systems for which this is exactly true have a name.
Definition 10.4 (Hessenberg form of size r) A semi-explicit DAE is in Hessenberg form of size r if it can be written
\begin{aligned} \dot{\bm x}_i &= \bm f_i(\bm x_1,\ldots,\bm x_{i+1},t), && i=1,\ldots,r-1, \\ \bm 0 &= \bm g(\bm x_1,t), \end{aligned} \tag{10.9}
that is, each differential equation may depend only on the variables up to the next block, the constraint depends only on \bm x_1, and the chain of Jacobians
\boxed{\; \frac{\partial\bm g}{\partial\bm x_1} \frac{\partial\bm f_1}{\partial\bm x_2} \frac{\partial\bm f_2}{\partial\bm x_3} \cdots \frac{\partial\bm f_{r-1}}{\partial\bm x_r} \quad\text{is nonsingular.} \;} \tag{10.10}
Theorem 10.2 (Index of a Hessenberg DAE) A DAE in Hessenberg form of size r has differentiation index exactly r.
The name refers to the staircase sparsity of the Jacobian of the whole system with respect to the derivatives: each block couples only to its neighbour, exactly as in a Hessenberg matrix. Setting r=2 in Equation 10.10 reproduces Equation 10.5 and r=3 reproduces Equation 10.8, with the identification \bm x_1=\bm y, \bm x_2=\bm x, \bm x_3=\bm z in the latter. Index one is the degenerate case in which the constraint already contains the algebraic variable, \bm0=\bm g(\bm x_1,\bm x_2,t) with \partial\bm g/\partial\bm x_2 nonsingular.
Hessenberg form is not a curiosity: it is the structure of constrained mechanics, and therefore of most of the DAEs that arise in this course. For a mechanical system with generalized coordinates \bm q, mass matrix \bm M(\bm q), applied forces \bm F, holonomic constraints \bm\Phi(\bm q)=\bm0 and constraint Jacobian \bm G=\partial\bm\Phi/\partial\bm q, the equations of motion are
\dot{\bm q}=\bm v, \qquad \bm M(\bm q)\dot{\bm v}=\bm F(\bm q,\bm v,t)-\bm G(\bm q)^T\bm\lambda, \qquad \bm 0=\bm\Phi(\bm q),
which is Equation 10.9 with r=3, \bm x_1=\bm q, \bm x_2=\bm v, \bm x_3=\bm\lambda. The chain Equation 10.10 becomes
\bm G\,\bm I\,\bigl(-\bm M^{-1}\bm G^T\bigr) =-\bm G\bm M^{-1}\bm G^T,
nonsingular precisely when \bm M is positive definite and \bm G has full row rank, that is when the mass matrix is physical and the constraints are independent. Every holonomically constrained mechanical system is an index-three DAE, and the two conditions that make it so are the two conditions a modeller checks anyway.
Knowing r in advance is worth the effort: it tells how many differentiations the reduction of Chapter 11 will need, how many invariants the reduced ODE will carry, and — as Order Reduction in DAE Integration shows — how much accuracy a Runge–Kutta method loses on the algebraic components.
10.3 Can a Runge–Kutta Method Be Applied Directly?
Before developing index reduction it is worth asking whether it is needed at all. A Runge–Kutta method can be written for a fully implicit ODE, and the fully implicit DAE \bm F(\bm x,\dot{\bm x},t)=\bm0 looks like one. The answer is instructive: the formulas exist, they can be implemented, and they do not solve the problem.
Implicit form of a Runge–Kutta method
For the explicit ODE \dot{\bm x}=\bm g(\bm x,t), an s-stage Runge–Kutta method advances as
\bm x_{k+1}=\bm x_k+\sum_{j=1}^{s}b_j\bm K_j, \qquad \bm K_j=\Delta t\, \bm g\Bigl(\bm x_k+\sum_{\ell=1}^{s}a_{j\ell}\bm K_\ell,\; t_k+c_j\Delta t\Bigr),
so that each stage satisfies \bm K_j/\Delta t=\bm g(\,\cdot\,). The quantity \bm K_j/\Delta t is therefore an approximation of the derivative at the stage point. If the equation is given implicitly as \bm F(\bm x,\dot{\bm x},t)=\bm0, the same stage conditions read
\bm F\Bigl( \bm x_k+\sum_{\ell=1}^{s}a_{j\ell}\bm K_\ell,\; \frac{\bm K_j}{\Delta t},\; t_k+c_j\Delta t \Bigr)=\bm0, \qquad j=1,\ldots,s, \tag{10.11}
a nonlinear system of s\,n equations in the s\,n unknowns \bm K_1,\ldots,\bm K_s, to be solved with the Newton iteration of Solving the Stage Equations. Nothing in Equation 10.11 requires \partial\bm F/\partial\dot{\bm x} to be invertible.
Why this is not enough
Equation 10.11 can be written down for any \bm F, which makes it tempting to conclude that any DAE can be integrated by an implicit Runge–Kutta method. The conclusion is false, for three separate reasons.
- The problem may have no solution at all. A DAE can be inconsistent, as the singular-pencil example of Section 8.4 shows: there, no differentiable \bm x satisfies the equations for the given forcing. A convergent iteration on Equation 10.11 cannot produce one.
- The initial condition may be inconsistent. Even for a solvable DAE, an arbitrary \bm x_k need not satisfy the hidden constraints, and then Equation 10.11 has no solution near \bm K_j=\bm0. This is why consistent initialization (Consistent Initialization) is a separate step and not a detail.
- The accuracy is not the one advertised. Even when everything is consistent, the classical order of a Runge–Kutta method does not carry over to the algebraic components of a high-index DAE. The loss depends on the index and on the method (is the subject of Order Reduction in DAE Integration).
The practical conclusion is the one that organizes the rest of this part:
Rather than applying a general-purpose solver to \bm F(\bm x,\dot{\bm x},t)=\bm0 and hoping, one first determines the index, reduces it to one or zero, records the hidden constraints as invariants, and only then integrates — with a method chosen in view of the remaining structure.
10.4 Examples
An index-three nonlinear system
Example from [1] equation (3.1)–(3.3)
\left\{ \begin{aligned} y' &= t-x \\ y'+z' &= 1 - (\alpha x+\eta) y \\ 0 &=y(2-y)+2z \end{aligned} \right.
Semi-Explicit Form
For the previous DAE system introduce dummy variable x_{\text{dot}} and the system decomposes into differential (blue) and algebraic (red) components:
\left\{ \color{blue} \begin{aligned} x' &= x_{\text{dot}}\\ y' &= t-x \\ z' &= 1-t+x - (\alpha x+\eta)y \end{aligned} \right. \quad \bigoplus \quad \color{DarkRed} 0 = y(1-y/2)+z \qquad \text{(Constraint 1)}
For this DAE we have:
- Differential variables: x, y, z.
- Algebraic variables: x_{\text{dot}}.
Reduction Procedure
We attempt to solve for x_{\text{dot}}, but the system is not yet solvable explicitly. We proceed by differentiating.
Step 1: Differentiate Constraints
Differentiate the constraint
\begin{aligned} \text{(Constraint 1)}' &= 1 - yy'+z' \\ &= 1 - y(t-x)+1-t+x - (\alpha x+\eta)y\\ &= 2 + (x-t)(y+1) - (\alpha x+\eta)y \end{aligned}
so that the DAE become:
\left\{ \color{blue} \begin{aligned} x' &= x_{\text{dot}}\\ y' &= t-x \\ z' &= 1-t+x - (\alpha x+\eta)y \end{aligned} \right. \quad \bigoplus \quad \left\{ \color{DarkRed} \begin{aligned} 0 & = 2 + (x-t)(y+1) - (\alpha x+\eta)y & \text{(Constraint 1)}' \\[0.5em] \color{DarkGreen} 0 & \color{DarkGreen} = y(1-y/2)+z & \color{DarkGreen}\text{(Constraint 1)}\\ \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 x_{\text{dot}} the algebraic variable must be differentiated again.
Step 2: Differentiate Again
\begin{aligned} \text{(Constraint 1)}'' &= (x'-1)(y+1)+(x-t)y'-\alpha(x'y+xy')+\eta y' \\ &= (x_{\text{dot}}-1)(y+1)+(x-t)(t-x)-\alpha(x_{\text{dot}}y+x(t-x))+\eta(t-x) \\ &= x_{\text{dot}}(y(1-\alpha)+1) -y-1-(x-t)^2+(\alpha x+\eta)(t-x) = 0 \end{aligned}
The dummy variable x_{\text{dot}} can be solved:
x_{\text{dot}} = \frac{y+1+(x-t+\alpha x+\eta)(x-t)}{y(1-\alpha)+1}
and the DAE is reduced to an index-1 DAE. Having done two derivation of the constraint the original DAE was of index 3. Removing the dummy variable we can obtain an ODE with invariants:
\left\{ \color{blue} \begin{aligned} x' &= \frac{y+1+(x-t+\alpha x+\eta)(x-t)}{y(1-\alpha)+1}\\ y' &= t-x \\ z' &= 1-t+x - (\alpha x+\eta)y \end{aligned} \right. \quad \bigoplus \quad \left\{ \color{DarkGreen} \begin{aligned} 0 & = 2 + (x-t)(y+1) - (\alpha x+\eta)y & \text{(Constraint 1)}' \\[0.5em] 0 & = y(1-y/2)+z & \text{(Constraint 1)} \end{aligned} \right.
Baumgarte Stabilization
As usual we want to stabilize the constraints to avoid the drift on the hidden constraints:
\text{(Constraints 1)}''+{\color{DarkOrange}2\eta\omega}\text{(Constraints 1)}' +{\color{DarkOrange}\omega^2}\text{(Constraints 1)}=0
and substituting the expression
\begin{aligned} & x_{\text{dot}}(y(1-\alpha)+1) -y-1-(x-t)^2+(\alpha x+\eta)(t-x) \\ & +{\color{DarkOrange}2\eta\omega}(2 + (x-t)(y+1) - (\alpha x+\eta)y) \\ & +{\color{DarkOrange}\omega^2}(y(1-y/2)+z)=0 \end{aligned}
so that
\begin{aligned} x_{\text{dot}} & = f(x,y,z,t) \\ f(x,y,z,t) & = \frac{a(x,y,z,t) -{\color{DarkOrange}2\eta\omega}b(x,y,z,t) -{\color{DarkOrange}\omega^2} c(x,y,z,t) }{y(1-\alpha)+1} \\ a(x,y,z,t) &= y+1+(x-t+\alpha x+\eta)(x-t) \\ b(x,y,z,t) &= 2 + (x-t)(y+1) - (\alpha x+\eta)y \\ c(x,y,z,t) &= y(1-y/2)+z \end{aligned}
and the stabilized ODE become
\left\{ \color{blue} \begin{aligned} x' &= f(x,y,z,t)\\ y' &= t-x \\ z' &= 1-t+x - (\alpha x+\eta)y \end{aligned} \right. \quad \bigoplus \quad \left\{ \color{DarkGreen} \begin{aligned} 0 & = 2 + (x-t)(y+1) - (\alpha x+\eta)y & \text{(Constraint 1)}' \\[0.5em] 0 & = y(1-y/2)+z & \text{(Constraint 1)} \end{aligned} \right.