Differential-Algebraic Equations, differential equations lesson

A bob on a rigid rod of length $1$ swings on a circle around the pivot. Newton’s second law gives rates of change for its position and velocity, but the rod adds a condition on the position alone, with no derivative in it, that must hold at every instant. Can the bob start anywhere on the circle with any velocity, and what decides how hard the rod pulls? Throughout, $t$ is the independent variable and a prime means $d/dt$ . This lesson shows how such models join rate equations with algebraic constraints, which initial values they accept, and what differentiating a constraint reveals.

Key terms

Differential-algebraic equation

Definition Let $X(t)$ be an unknown vector function with $n$ entries. A differential-algebraic equation, or DAE, is a system of $n$ equations

F(t, X, X') = 0

whose Jacobian matrix $\partial F/\partial X'$ , taken with respect to the entries of $X'$ , is singular at every point of the region considered. Near a point where $F = 0$ and this matrix is invertible, with $F$ continuously differentiable, the implicit function theorem would let one solve for $X'$ and obtain an ordinary first-order system $X' = f(t, X)$ , so the singular Jacobian is what sets a DAE apart. A common case is the linear-implicit form

M X' = \varphi(t, X)

with a constant singular $n \times n$ matrix $M$ , often called the mass matrix. Because $M$ is singular, some nonzero row vector $u^{T}$ satisfies $u^{T}M = 0$ , and multiplying the system by it gives $u^{T}\varphi(t, X) = 0$ , an equation without derivatives that has to hold at every time.

In plain words A DAE mixes rates of change with conditions that must hold at every instant. It cannot be solved for every derivative at once, so the conditions are part of the problem from the start, and they limit both the starting values and the motion.

Example For the unknown $X = (x, z)^{T}$ , the system

\begin{aligned} x' &= -x + z, \\ 0 &= x + z - 1 \end{aligned}

has the linear-implicit form with

\begin{gathered} M = \begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix}, \\[0.4em] \varphi(t, X) = \begin{bmatrix} -x + z \\ x + z - 1 \end{bmatrix}. \end{gathered}

The matrix $M$ has determinant $0$ , and the row vector $u^{T} = (0, 1)$ picks out the second row, the relation $x + z - 1 = 0$ with no derivative in it. Written as $F(t, X, X') = 0$ , the system has $F = (x' + x - z,\ x + z - 1)^{T}$ , and its Jacobian with respect to $(x', z')$ is this same matrix $M$ . No derivative of $z$ appears anywhere, so the system as written gives no formula for $z'$ .

Differential and algebraic variables

Definition A semi-explicit DAE has the form

\begin{aligned} x' &= f(t, x, z), \\ 0 &= g(t, x, z), \end{aligned}

where $x(t) \in \mathbb{R}^m$ and $z(t) \in \mathbb{R}^k$ with $k \ge 1$ are columns that together make up the unknown $X$ . The entries of $x$ are the differential variables, because each has its own rate equation. The entries of $z$ are the algebraic variables, because no derivative of $z$ appears in the system. Each of the $k$ equations $0 = g_i(t, x, z)$ is an algebraic constraint, a relation without derivatives that must hold at every time in the interval of the solution. The Jacobian of this system with respect to $X'$ has an identity block for $x'$ and zero rows for the constraints, so it is singular. A solution on an open interval $I$ is a pair with $x$ differentiable and $z$ continuous on $I$ such that both lines hold for every $t$ in $I$ . The labels describe the system as written, so rewriting the model can change which variables count as algebraic.

In plain words The model hands each differential variable a rate. An algebraic variable gets no rate of its own, and in a well-posed model its values are whatever the constraints require at each instant, sometimes only after the constraints are differentiated. Those values can still change from moment to moment.

Example In the system $x' = -x + z$ , $0 = x + z - 1$ , the variable $x$ is differential and $z$ is algebraic, with $m = k = 1$ . The constraint gives $z = 1 - x$ , and substituting it leaves the linear equation $x' + 2x = 1$ . With the integrating factor $e^{2t}$ it reads $(e^{2t}x)' = e^{2t}$ , so $e^{2t}x = \tfrac{1}{2}e^{2t} + C$ and $x = \tfrac{1}{2} + Ce^{-2t}$ . The start $x(0) = 0$ gives $C = -\tfrac{1}{2}$ , and the constraint returns $z$ , so for all real $t$

\begin{gathered} x = \tfrac{1}{2}\bigl(1 - e^{-2t}\bigr), \\[0.4em] z = \tfrac{1}{2}\bigl(1 + e^{-2t}\bigr). \end{gathered}

A check of every equation confirms this. The derivative $x' = e^{-2t}$ equals $-x + z$ , the sum $x + z$ is $1$ , and the values at $t = 0$ are $x = 0$ and $z = 1$ . The algebraic variable falls from $1$ toward $\tfrac{1}{2}$ , so algebraic does not mean constant.

Consistent initial values

Definition Initial values $X(t_0) = X_0$ of a DAE are consistent if the DAE has a solution on some open interval containing $t_0$ that takes these values at $t_0$ . Consistent values must satisfy every algebraic constraint at $t_0$ , together with every hidden constraint. For a semi-explicit system with $f$ and $g$ continuously differentiable, suppose that

g(t_0, x_0, z_0) = 0

and that the $k \times k$ Jacobian matrix $\partial g/\partial z$ is invertible at $(t_0, x_0, z_0)$ . Then the values $x(t_0) = x_0$ and $z(t_0) = z_0$ are consistent, and any two solutions through them agree on some open interval containing $t_0$ . In this case the differential values $x_0$ are the data one prescribes, and $z_0$ must then solve the constraint at $t_0$ .

In plain words A starting point is allowed only when the model can actually move from it. Values that break a constraint, written or hidden, start no solution at all, and a computation started from them is working on a different problem.

Example For $x' = -x + z$ , $0 = x + z - 1$ , the derivative $\partial g/\partial z$ equals $1$ , so any $x_0$ may be prescribed and the constraint forces $z_0 = 1 - x_0$ . With $x(0) = 0$ the only consistent value is $z(0) = 1$ , and the pair $x(0) = z(0) = 0$ is inconsistent, since $0 + 0 - 1 \ne 0$ . From a general consistent start the solution is

\begin{gathered} x = \tfrac{1}{2} + \bigl(x_0 - \tfrac{1}{2}\bigr)e^{-2t}, \\[0.4em] z = \tfrac{1}{2} - \bigl(x_0 - \tfrac{1}{2}\bigr)e^{-2t}, \end{gathered}

which keeps $x + z = 1$ for every $t$ and sends both variables to $\tfrac{1}{2}$ . One free number, the value $x_0$ , selects the solution.

Hidden constraint

Definition Along a solution on an interval $I$ , a constraint $0 = g(t, x, z)$ holds at every $t$ in $I$ , so its derivative with respect to $t$ along the solution is also zero wherever that derivative exists. Substituting the rate equations into this derivative gives a derived relation that every such solution obeys although it is not written in the system. When a derived relation contains no derivatives after this substitution and does not follow from the written constraints, it is a hidden constraint. It restricts the states a solution can pass through, so consistent initial values must satisfy it too. Differentiating again can reveal further hidden constraints. A derived relation adds to the original constraint. It can stand in for it only when the initial values satisfy the original constraint, because a zero derivative only says that $g$ stays constant along the solution, not that the constant is $0$ .

In plain words A constraint has to keep holding as time runs, and that requirement can impose conditions nobody wrote down. A rod of fixed length, for instance, forbids any motion toward or away from the pivot, so a condition on the position hides a condition on the velocity.

Example Consider $x' = z$ , $0 = x - \sin t$ , with $x$ differential and $z$ algebraic. The constraint says nothing about $z$ directly. Differentiating it gives $x' - \cos t = 0$ , and substituting $x' = z$ gives the hidden constraint

z = \cos t.

So the only consistent initial values at $t = 0$ are $x(0) = 0$ and $z(0) = 1$ , and the solution is $x = \sin t$ , $z = \cos t$ . Here not even the differential variable has a free starting value. Compare the first system, $x' = -x + z$ with $0 = x + z - 1$ . Differentiating its constraint gives $x' + z' = 0$ , so the two rates are equal and opposite, and substituting $x' = -x + z$ gives $z' = x - z$ . That derived relation fixes the rate of $z$ but places no new condition on the starting values.

Differentiation index

Definition Differentiate the equations of a DAE with respect to $t$ along solutions, all of them or only some, as many times as needed. The differentiation index is the smallest number $\nu$ such that the original equations together with their derivatives up to order $\nu$ determine $X'$ as a continuous function of $t$ and $X$ by algebraic steps. The resulting ordinary system $X' = \psi(t, X)$ is called the underlying ODE, and an ordinary system $X' = f(t, X)$ has index $0$ . For a semi-explicit system with continuously differentiable $f$ and $g$ , write $g_t$ , $g_x$ , $g_z$ and $f_z$ for the partial derivatives and Jacobian matrices. If $g_z$ is invertible along the solution, one differentiation of the constraint gives

g_t + g_x\,f + g_z\,z' = 0,

which can be solved for $z'$ , so the index is $1$ . If $g$ does not depend on $z$ , that first derivative contains no $z'$ , so the index is at least $2$ , and it equals $2$ when the matrix $g_x f_z$ is invertible. The index belongs to a formulation and to this convention. Rewriting a model can change it, and other notions, such as the perturbation index, are defined differently and need not give the same number for every system.

In plain words The index counts how many rounds of differentiating the constraints it takes before every unknown, the algebraic ones included, has a rate equation. A higher index means more differentiations between the model as written and an ordinary system, and numerical methods generally find such problems harder.

Example In $x' = -x + z$ , $0 = x + z - 1$ the matrix $g_z$ is the number $1$ , and one differentiation gives $z' = x - z$ . The underlying ODE is

\begin{aligned} x' &= -x + z, \\ z' &= x - z, \end{aligned}

so the index is $1$ . In $x' = z$ , $0 = x - \sin t$ the constraint does not contain $z$ . The first differentiation gives the hidden constraint $z = \cos t$ and the second gives $z' = -\sin t$ , so the index is $2$ , in agreement with $g_x f_z = 1$ . The underlying ODE can have more solutions than the DAE. Adding the two rows of the first underlying ODE gives $(x + z)' = 0$ , so each of its solutions keeps $x + z$ constant, but only those with $x + z = 1$ solve the DAE.

Constrained pendulum

Definition Scale the units so that the mass, the rod length and the gravitational acceleration all equal $1$ . Put the pivot at the origin, write the position of the bob as $(x, y)$ with $y$ positive upward and its velocity as $(v, w)$ , and let the rod act on the bob with the force $-\lambda(x, y)$ . The constrained pendulum is the DAE

\begin{aligned} x' &= v, \\ y' &= w, \\ v' &= -\lambda x, \\ w' &= -1 - \lambda y, \\ 0 &= x^2 + y^2 - 1. \end{aligned}

The differential variables are $x$ , $y$ , $v$ and $w$ , and the multiplier $\lambda$ is the algebraic variable. The rod force has magnitude $\lvert \lambda \rvert$ , and a positive $\lambda$ means the rod pulls the bob toward the pivot. Differentiating the constraint and using $x' = v$ , $y' = w$ gives $2(xv + yw) = 0$ , the hidden constraint

xv + yw = 0.

Differentiating that and inserting the accelerations gives $v^2 + w^2 - y - \lambda(x^2 + y^2) = 0$ , which on the circle reads

\lambda = v^2 + w^2 - y.

A third differentiation gives

\begin{aligned} \lambda' &= -2\lambda(xv + yw) - 3w \\[0.4em] &= -3w, \end{aligned}

which supplies the last missing rate, so this formulation has index $3$ . Consistent initial values put the bob on the circle, make its velocity tangent to the circle, and fix $\lambda$ by the formula above. Of the five unknowns only two numbers are free, the starting point on the circle and the signed speed along it. Other formulations of the same motion have other indices. With the velocity constraint $xv + yw = 0$ in place of the position constraint the index is $2$ , although the bob must still be started on the circle, and in the angle $\theta$ with $x = \sin\theta$ and $y = -\cos\theta$ the motion is the ordinary equation $\theta'' + \sin\theta = 0$ , of index $0$ .

In plain words The rod pulls exactly as hard as needed to keep the bob on the circle. That pull is not an extra starting value to choose, because the position and velocity already decide it, and the fixed length hides a condition on the velocity.

Example At rest at the bottom, $(x, y, v, w) = (0, -1, 0, 0)$ . The bob is on the circle, the tangency condition $xv + yw = 0$ holds, and the formula gives $\lambda = 0 + 0 - (-1) = 1$ . The accelerations are then $v' = -1 \cdot 0 = 0$ and $w' = -1 - 1 \cdot (-1) = 0$ , so the rod force balances gravity and the bob stays at rest, as it should. This small check also catches a sign error in the multiplier convention. Released from rest with the rod horizontal, at $(1, 0, 0, 0)$ , the bob has $\lambda = 0$ , so at that instant the rod exerts no force and the accelerations are $v' = 0$ and $w' = -1$ , the start of free fall.

Check yourself

In the system $x' = z$ , $0 = z + 3x - 6$ , which variable is differential and which is algebraic? With $x(0) = 0$ , find the consistent value of $z(0)$ and solve the problem.

Show the answer

The derivative of $x$ appears, so $x$ is differential, and $z$ is algebraic. The constraint has $\partial g/\partial z = 1$ , so the index is $1$ and $x(0)$ may be prescribed freely. At $t = 0$ the constraint gives $z(0) = 6 - 3 \cdot 0 = 6$ , the only consistent value. Eliminating $z = 6 - 3x$ leaves $x' + 3x = 6$ . With the integrating factor $e^{3t}$ this is $(e^{3t}x)' = 6e^{3t}$ , so $x = 2 + Ce^{-3t}$ , and $x(0) = 0$ gives $C = -2$ . Hence

\begin{gathered} x = 2 - 2e^{-3t}, \\[0.4em] z = 6e^{-3t}. \end{gathered}

As a check, $x' = 6e^{-3t} = z$ , the expression $z + 3x - 6 = 6e^{-3t} + 6 - 6e^{-3t} - 6$ is $0$ , and the values at $t = 0$ are $0$ and $6$ .

A classmate replaces the constraint of $x' = -x + z$ , $0 = x + z - 1$ by its derivative, which turns the problem into the ordinary system $x' = -x + z$ , $z' = x - z$ , and starts it from $x(0) = 0$ and $z(0) = \tfrac{1}{2}$ . Solve the classmate’s problem and explain why its solution does not solve the DAE.

Show the answer

Adding the two equations gives $(x + z)' = 0$ , so $x + z$ keeps its starting value $\tfrac{1}{2}$ . Then $x' = -x + (\tfrac{1}{2} - x) = \tfrac{1}{2} - 2x$ , so $x = \tfrac{1}{4} + Ce^{-2t}$ with $C = -\tfrac{1}{4}$ , and

\begin{gathered} x = \tfrac{1}{4}\bigl(1 - e^{-2t}\bigr), \\[0.4em] z = \tfrac{1}{4}\bigl(1 + e^{-2t}\bigr). \end{gathered}

Both rate equations hold, since $x' = \tfrac{1}{2}e^{-2t} = -x + z$ and $z' = -\tfrac{1}{2}e^{-2t} = x - z$ . But $x + z - 1 = -\tfrac{1}{2}$ at every time, so the original constraint fails throughout, and both variables approach $\tfrac{1}{4}$ instead of $\tfrac{1}{2}$ . The starting values were inconsistent, and the differentiated constraint preserved the wrong sum instead of correcting it. From the consistent start $z(0) = 1$ the same ordinary system gives the true solution, which is why the original constraint is kept beside its derivative and checked during a computation.

For the constrained pendulum in scaled units, decide which of the starting states $(x, y, v, w) = (0.6, -0.8, 0.8, 0.6)$ and $(x, y, v, w) = (0.6, -0.8, 0.6, 0.8)$ is consistent. For the consistent one, find $\lambda(0)$ and the initial accelerations $v'(0)$ and $w'(0)$ .

Show the answer

Both states put the bob on the circle, since $0.36 + 0.64 = 1$ . For the first state the hidden constraint holds, because $xv + yw = 0.48 - 0.48 = 0$ . The formula then gives $\lambda(0) = 0.64 + 0.36 + 0.8 = 1.8$ , and

\begin{aligned} v'(0) &= -1.8 \cdot 0.6 \\ &= -1.08, \\[0.4em] w'(0) &= -1 + 1.8 \cdot 0.8 \\ &= 0.44. \end{aligned}

As a check, the component of this acceleration along the unit vector $(x, y)$ is $(-1.08)(0.6) + (0.44)(-0.8) = -1$ , which is the inward acceleration $-(v^2 + w^2) = -1$ needed to move on the unit circle at speed $1$ . For the second state, $xv + yw = 0.36 - 0.64 = -0.28 \ne 0$ , so its velocity has a component along the rod. It violates the hidden constraint, and no choice of $\lambda(0)$ makes it consistent.

Differential-Algebraic Equations

Key terms

Differential-algebraic equation

Differential and algebraic variables

Consistent initial values

Hidden constraint

Differentiation index

Constrained pendulum

Check yourself

How this lesson connects

Builds on