In earlier chapters, a continuous-time LTI system was written as \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \). Descriptor systems generalize this by allowing a matrix multiplier on the derivative term:

\[ \mathbf{E}\dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t), \qquad \mathbf{y}(t)=\mathbf{C}\mathbf{x}(t)+\mathbf{D}\mathbf{u}(t). \]

Here \( \mathbf{E},\mathbf{A}\in\mathbb{R}^{n\times n} \), \( \mathbf{B}\in\mathbb{R}^{n\times m} \), \( \mathbf{C}\in\mathbb{R}^{p\times n} \), and \( \mathbf{D}\in\mathbb{R}^{p\times m} \). The model is called descriptor because the vector \( \mathbf{x} \) may contain true dynamic states and auxiliary variables used to describe constraints. It is called singular when

\[ \operatorname{rank}(\mathbf{E}) < n. \]

A singular \( \mathbf{E} \) means that not every row of the model contains a derivative. Some rows impose instantaneous algebraic relations among \( \mathbf{x} \) and \( \mathbf{u} \). Therefore, not every initial vector \( \mathbf{x}(0) \) is physically or mathematically admissible.

flowchart TD
  A["Model equations"] --> B["Check rank of E"]
  B --> C["rank E equals n?"]
  C -->|"yes"| D["Standard state-space ODE"]
  C -->|"no"| E["Descriptor or singular model"]
  E --> F["Check regularity of pencil"]
  F --> G["Separate differential and algebraic parts"]
  G --> H["Find consistent initial conditions"]
  H --> I["Analyze finite modes and constraint modes"]
        

Descriptor dynamics are governed not by the ordinary eigenvalue problem of \( \mathbf{A} \), but by the matrix pencil \( \lambda\mathbf{E}-\mathbf{A} \). The pencil is called regular if its determinant is not the zero polynomial:

\[ \det(\lambda\mathbf{E}-\mathbf{A}) \not\equiv 0. \]

Regularity is the descriptor analogue of having a well-defined square state model. If the pencil is regular, the finite generalized eigenvalues are the roots of

\[ \det(\lambda\mathbf{E}-\mathbf{A})=0. \]

If \( \mathbf{E} \) is nonsingular, then \( \lambda\mathbf{E}-\mathbf{A} \) is regular and the generalized eigenvalue problem reduces to the ordinary one:

\[ \lambda\mathbf{E}\mathbf{v}=\mathbf{A}\mathbf{v} \quad \Longleftrightarrow \quad \lambda\mathbf{v}=\mathbf{E}^{-1}\mathbf{A}\mathbf{v}. \]

When \( \mathbf{E} \) is singular, the degree of \( \det(\lambda\mathbf{E}-\mathbf{A}) \) may be smaller than \( n \). The missing eigenvalues are interpreted as infinite eigenvalues. These correspond to algebraic constraint structure rather than ordinary exponential modes.

\[ \deg\det(\lambda\mathbf{E}-\mathbf{A}) = r \quad \Longrightarrow \quad n-r \text{ infinite eigenvalues, counted with algebraic multiplicity}. \]

A useful first case occurs when coordinates can be ordered so that \( \mathbf{x}=\begin{bmatrix}\mathbf{x}_d'&\mathbf{x}_a'\end{bmatrix}' \), where \( \mathbf{x}_d \) contains differential variables and \( \mathbf{x}_a \) contains algebraic variables. In an index-1 form, the descriptor matrix can be written as

\[ \mathbf{E}=\begin{bmatrix}\mathbf{I}_r & \mathbf{0}\\ \mathbf{0} & \mathbf{0}\end{bmatrix},\qquad \mathbf{x}=\begin{bmatrix}\mathbf{x}_d\\\mathbf{x}_a\end{bmatrix}. \]

Partitioning \( \mathbf{A} \) and \( \mathbf{B} \) conformably gives

\[ \begin{aligned} \dot{\mathbf{x} }_d &= \mathbf{A}_{11}\mathbf{x}_d+ \mathbf{A}_{12}\mathbf{x}_a+\mathbf{B}_1\mathbf{u},\\ \mathbf{0} &= \mathbf{A}_{21}\mathbf{x}_d+ \mathbf{A}_{22}\mathbf{x}_a+\mathbf{B}_2\mathbf{u}. \end{aligned} \]

If \( \mathbf{A}_{22} \) is nonsingular, the algebraic variable is uniquely determined by \( \mathbf{x}_d \) and \( \mathbf{u} \):

\[ \mathbf{x}_a=-\mathbf{A}_{22}^{-1}\mathbf{A}_{21}\mathbf{x}_d -\mathbf{A}_{22}^{-1}\mathbf{B}_2\mathbf{u}. \]

Substitution gives a reduced ordinary state equation:

\[ \dot{\mathbf{x} }_d=\left(\mathbf{A}_{11}- \mathbf{A}_{12}\mathbf{A}_{22}^{-1}\mathbf{A}_{21}\right)\mathbf{x}_d +\left(\mathbf{B}_{1}-\mathbf{A}_{12}\mathbf{A}_{22}^{-1}\mathbf{B}_{2}\right)\mathbf{u}. \]

flowchart TD
  A["Descriptor model"] --> B["Reorder variables: dynamic and algebraic"]
  B --> C["Partition E, A, and B"]
  C --> D["Is A22 invertible?"]
  D -->|"yes"| E["Solve algebraic equation for xa"]
  E --> F["Substitute into dynamic equations"]
  F --> G["Obtain reduced standard ODE"]
  D -->|"no"| H["Higher index or \nsingular constraint case"]
        

In a standard ODE, any initial condition \( \mathbf{x}(0)=\mathbf{x}_0 \) is admissible. In a descriptor model, the algebraic constraints must already hold at the initial time. For the partitioned index-1 case, consistency requires

\[ \mathbf{0}=\mathbf{A}_{21}\mathbf{x}_{d0}+ \mathbf{A}_{22}\mathbf{x}_{a0}+\mathbf{B}_{2}\mathbf{u}(0). \]

When \( \mathbf{A}_{22} \) is nonsingular, one may freely choose \( \mathbf{x}_{d0} \), but the algebraic component is fixed by

\[ \mathbf{x}_{a0}=-\mathbf{A}_{22}^{-1}\mathbf{A}_{21}\mathbf{x}_{d0} -\mathbf{A}_{22}^{-1}\mathbf{B}_2\mathbf{u}(0). \]

If this condition is violated, a numerical DAE solver may generate an impulsive correction, fail to initialize, or silently project the initial condition onto a consistent manifold. For control design, consistency is not a numerical detail; it is part of the model definition.

For the homogeneous descriptor system \( \mathbf{E}\dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \), finite modes are determined by the generalized eigenvalues of \( \lambda\mathbf{E}-\mathbf{A} \). For a regular, impulse-free continuous-time descriptor system, asymptotic stability is characterized by the location of finite eigenvalues:

\[ \operatorname{Re}(\lambda_i)<0 \quad \text{for every finite generalized eigenvalue } \lambda_i. \]

This condition mirrors the standard state-space stability test, but it must be read through the generalized eigenvalue problem. Infinite modes describe constraint dynamics; if they are not structurally benign, they may produce impulsive behavior. A first lesson should therefore retain the following separation:

\[ \begin{aligned} \text{finite eigenvalues} \;&\Rightarrow\; \text{ordinary exponential modes}, \\ \text{infinite eigenvalues} \;&\Rightarrow\; \text{algebraic constraint structure}. \end{aligned} \]

Consider the descriptor system

\[ \mathbf{E}=\begin{bmatrix}1&0&0\\0&1&0\\0&0&0\end{bmatrix},\quad \mathbf{A}=\begin{bmatrix}0&1&0\\-2&-3&1\\1&0&1\end{bmatrix},\quad \mathbf{B}=\begin{bmatrix}0\\1\\-1\end{bmatrix}. \]

The descriptor equations are

\[ \begin{aligned} \dot{x}_1 &= x_2,\\ \dot{x}_2 &= -2x_1-3x_2+x_3+u,\\ 0 &= x_1+x_3-u. \end{aligned} \]

The algebraic equation gives \( x_3=-x_1+u \). Hence the reduced ODE is

\[ \begin{bmatrix}\dot{x}_1\\\dot{x}_2\end{bmatrix} =\begin{bmatrix}0&1\\-3&-3\end{bmatrix} \begin{bmatrix}x_1\\x_2\end{bmatrix} +\begin{bmatrix}0\\2\end{bmatrix}u. \]

The characteristic polynomial of the reduced dynamics is

\[ \det\left(s\mathbf{I}-\begin{bmatrix}0&1\\-3&-3\end{bmatrix}\right) =s^2+3s+3, \]

so the finite modes are stable complex conjugates:

\[ s_{1,2}=\frac{-3\pm j\sqrt{3} }{2}. \]

Theorem 1 (Reduction of an index-1 descriptor system):

Suppose the partitioned descriptor system in Section 3 has nonsingular \( \mathbf{A}_{22} \). Then every solution satisfying the algebraic constraint is equivalent to a solution of the reduced ODE for \( \mathbf{x}_d \), with \( \mathbf{x}_a \) reconstructed algebraically.

Proof:

The algebraic row is \( \mathbf{0}=\mathbf{A}_{21}\mathbf{x}_d+\mathbf{A}_{22}\mathbf{x}_a+\mathbf{B}_2\mathbf{u} \). Since \( \mathbf{A}_{22} \) is nonsingular, solving this linear equation gives a unique \( \mathbf{x}_a \). Substituting that expression into the differential row yields the reduced ODE in Section 3. Conversely, any solution of the reduced ODE, together with the reconstructed \( \mathbf{x}_a \), satisfies both the differential and algebraic rows. Therefore, the two descriptions are equivalent on the consistent constraint manifold.

\[ \boxed{\mathbf{A}_{22}\text{ nonsingular } \Longrightarrow \text{index-1 descriptor model reduces to an ordinary state equation}.} \]

Theorem 2 (Standard systems as a special case):

If \( \mathbf{E} \) is nonsingular, then the descriptor model is exactly a standard state-space model after premultiplication by \( \mathbf{E}^{-1} \).

Proof:

\[ \mathbf{E}\dot{\mathbf{x} }=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \quad \Longrightarrow \quad \dot{\mathbf{x} }=\mathbf{E}^{-1}\mathbf{A}\mathbf{x}+ \mathbf{E}^{-1}\mathbf{B}\mathbf{u}. \]

Therefore, all ordinary state-space results apply with \( \mathbf{A}_s=\mathbf{E}^{-1}\mathbf{A} \) and \( \mathbf{B}_s=\mathbf{E}^{-1}\mathbf{B} \).

Python support for descriptor systems is usually built from numerical linear algebra and DAE tools. Useful libraries include numpy, scipy.linalg.eig for generalized eigenvalues, scipy.integrate.solve_ivp for the reduced ODE, matplotlib for visualization, and specialized packages such as Slycot or Python-Control for broader state-space workflows.

"""
Chapter30_Lesson1.py
Descriptor (singular) system example: E xdot = A x + B u.

This script demonstrates:
1. generalized eigenvalue analysis of the pencil lambda E - A,
2. index-1 algebraic reduction when E = diag(I, 0),
3. simulation of the reduced ODE and reconstruction of the algebraic state.

Required packages: numpy, scipy, matplotlib
"""

import numpy as np
from scipy.linalg import eig
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Descriptor model with two differential states and one algebraic state.
# E xdot = A x + B u
E = np.diag([1.0, 1.0, 0.0])
A = np.array([
    [0.0,  1.0, 0.0],
    [-2.0, -3.0, 1.0],
    [1.0,  0.0, 1.0],
])
B = np.array([[0.0], [1.0], [-1.0]])

# Generalized eigenvalues solve det(lambda E - A) = 0.
finite_and_infinite_eigs = eig(A, E, right=False)
print("Generalized eigenvalues of (A, E):")
print(finite_and_infinite_eigs)

# Partition x = [xd; xa], where xd = [x1, x2]^T and xa = x3.
A11 = A[:2, :2]
A12 = A[:2, 2:3]
A21 = A[2:3, :2]
A22 = A[2:3, 2:3]
B1 = B[:2, :]
B2 = B[2:3, :]

# Algebraic equation: 0 = A21 xd + A22 xa + B2 u.
# If A22 is nonsingular, xa = -A22^{-1} A21 xd - A22^{-1} B2 u.
A22_inv = np.linalg.inv(A22)
Ar = A11 - A12 @ A22_inv @ A21
Br = B1 - A12 @ A22_inv @ B2

print("Reduced ODE matrix Ar:")
print(Ar)
print("Reduced input matrix Br:")
print(Br)
print("Eigenvalues of the reduced dynamics:")
print(np.linalg.eigvals(Ar))


def u_of_t(t: float) -> float:
    """Piecewise-smooth input."""
    return 1.0 if t >= 1.0 else 0.0


def rhs(t: float, xd: np.ndarray) -> np.ndarray:
    u = np.array([[u_of_t(t)]])
    return (Ar @ xd.reshape(2, 1) + Br @ u).ravel()


def reconstruct_xa(xd: np.ndarray, u: float) -> float:
    xd_col = xd.reshape(2, 1)
    xa = -A22_inv @ A21 @ xd_col - A22_inv @ B2 * u
    return float(xa[0, 0])

# Consistent initial condition: choose xd(0), then compute xa(0) from constraint.
xd0 = np.array([0.2, 0.0])
xa0 = reconstruct_xa(xd0, u_of_t(0.0))
print("Consistent initial x0 =", np.r_[xd0, xa0])

sol = solve_ivp(rhs, (0.0, 8.0), xd0, t_eval=np.linspace(0.0, 8.0, 401))
x1 = sol.y[0, :]
x2 = sol.y[1, :]
x3 = np.array([reconstruct_xa(sol.y[:, k], u_of_t(sol.t[k])) for k in range(sol.t.size)])
constraint_residual = A21 @ sol.y + A22 @ x3.reshape(1, -1) + B2 @ np.array([[u_of_t(t) for t in sol.t]])
print("Maximum algebraic constraint residual:", np.max(np.abs(constraint_residual)))

plt.figure(figsize=(8, 4.8))
plt.plot(sol.t, x1, label="x1 differential")
plt.plot(sol.t, x2, label="x2 differential")
plt.plot(sol.t, x3, label="x3 algebraic")
plt.xlabel("time")
plt.ylabel("state")
plt.title("Descriptor system reduced to an index-1 ODE")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

The following implementations avoid external control libraries. They use the reduced 2-state ODE, reconstruct the algebraic state, and verify the algebraic constraint residual during simulation.

Chapter30_Lesson1.cpp

/*
Chapter30_Lesson1.cpp
Index-1 descriptor system reduction for E xdot = A x + B u.

This standalone C++ example avoids external libraries. It uses the reduced
ODE obtained from the algebraic constraint and simulates it with RK4.
*/

#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>

using Vec2 = std::array<double, 2>;
using Mat2 = std::array<std::array<double, 2>, 2>;

constexpr Mat2 Ar = { { {0.0, 1.0}, {-3.0, -3.0} } };
constexpr Vec2 Br = {0.0, 2.0};

double u_of_t(double t) {
    return (t >= 1.0) ? 1.0 : 0.0;
}

Vec2 f(double t, const Vec2& x) {
    const double u = u_of_t(t);
    return {
        Ar[0][0] * x[0] + Ar[0][1] * x[1] + Br[0] * u,
        Ar[1][0] * x[0] + Ar[1][1] * x[1] + Br[1] * u
    };
}

Vec2 add_scaled(const Vec2& x, const Vec2& k, double scale) {
    return {x[0] + scale * k[0], x[1] + scale * k[1]};
}

Vec2 rk4_step(double t, const Vec2& x, double h) {
    Vec2 k1 = f(t, x);
    Vec2 k2 = f(t + 0.5 * h, add_scaled(x, k1, 0.5 * h));
    Vec2 k3 = f(t + 0.5 * h, add_scaled(x, k2, 0.5 * h));
    Vec2 k4 = f(t + h, add_scaled(x, k3, h));
    return {
        x[0] + (h / 6.0) * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]),
        x[1] + (h / 6.0) * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1])
    };
}

double algebraic_x3(const Vec2& xd, double u) {
    // Constraint: 0 = x1 + x3 - u, hence x3 = -x1 + u.
    return -xd[0] + u;
}

void print_reduced_eigenvalues() {
    // Eigenvalues of a 2 by 2 matrix using the quadratic formula.
    const double trace = Ar[0][0] + Ar[1][1];
    const double determinant = Ar[0][0] * Ar[1][1] - Ar[0][1] * Ar[1][0];
    const double discriminant = trace * trace - 4.0 * determinant;
    std::cout << "Reduced characteristic polynomial: lambda^2 - (trace) lambda + det" << std::endl;
    if (discriminant >= 0.0) {
        const double r1 = 0.5 * (trace + std::sqrt(discriminant));
        const double r2 = 0.5 * (trace - std::sqrt(discriminant));
        std::cout << "Reduced eigenvalues: " << r1 << ", " << r2 << std::endl;
    } else {
        const double real = 0.5 * trace;
        const double imag = 0.5 * std::sqrt(-discriminant);
        std::cout << "Reduced eigenvalues: " << real << " + " << imag << "i, "
                  << real << " - " << imag << "i" << std::endl;
    }
}

int main() {
    print_reduced_eigenvalues();

    Vec2 xd = {0.2, 0.0};
    const double h = 0.01;
    const int steps = 800;

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "time,x1,x2,x3,constraint_residual" << std::endl;

    for (int k = 0; k <= steps; ++k) {
        const double t = k * h;
        const double u = u_of_t(t);
        const double x3 = algebraic_x3(xd, u);
        const double residual = xd[0] + x3 - u;

        if (k % 40 == 0) {
            std::cout << t << "," << xd[0] << "," << xd[1] << "," << x3 << "," << residual << std::endl;
        }
        xd = rk4_step(t, xd, h);
    }
    return 0;
}

Chapter30_Lesson1.java

/*
Chapter30_Lesson1.java
Index-1 descriptor system reduction for E xdot = A x + B u.
Compile: javac Chapter30_Lesson1.java
Run:     java Chapter30_Lesson1
*/

public class Chapter30_Lesson1 {
    static final double[][] Ar = {
        {0.0, 1.0},
        {-3.0, -3.0}
    };
    static final double[] Br = {0.0, 2.0};

    static double uOfT(double t) {
        return t >= 1.0 ? 1.0 : 0.0;
    }

    static double[] f(double t, double[] x) {
        double u = uOfT(t);
        return new double[] {
            Ar[0][0] * x[0] + Ar[0][1] * x[1] + Br[0] * u,
            Ar[1][0] * x[0] + Ar[1][1] * x[1] + Br[1] * u
        };
    }

    static double[] addScaled(double[] x, double[] k, double scale) {
        return new double[] {x[0] + scale * k[0], x[1] + scale * k[1]};
    }

    static double[] rk4Step(double t, double[] x, double h) {
        double[] k1 = f(t, x);
        double[] k2 = f(t + 0.5 * h, addScaled(x, k1, 0.5 * h));
        double[] k3 = f(t + 0.5 * h, addScaled(x, k2, 0.5 * h));
        double[] k4 = f(t + h, addScaled(x, k3, h));
        return new double[] {
            x[0] + (h / 6.0) * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]),
            x[1] + (h / 6.0) * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1])
        };
    }

    static double algebraicX3(double[] xd, double u) {
        // Constraint: 0 = x1 + x3 - u, hence x3 = -x1 + u.
        return -xd[0] + u;
    }

    static void printReducedEigenvalues() {
        double trace = Ar[0][0] + Ar[1][1];
        double determinant = Ar[0][0] * Ar[1][1] - Ar[0][1] * Ar[1][0];
        double discriminant = trace * trace - 4.0 * determinant;
        if (discriminant >= 0.0) {
            double r1 = 0.5 * (trace + Math.sqrt(discriminant));
            double r2 = 0.5 * (trace - Math.sqrt(discriminant));
            System.out.println("Reduced eigenvalues: " + r1 + ", " + r2);
        } else {
            double real = 0.5 * trace;
            double imag = 0.5 * Math.sqrt(-discriminant);
            System.out.println("Reduced eigenvalues: " + real + " + " + imag + "i, "
                    + real + " - " + imag + "i");
        }
    }

    public static void main(String[] args) {
        printReducedEigenvalues();
        double[] xd = {0.2, 0.0};
        double h = 0.01;
        int steps = 800;

        System.out.println("time,x1,x2,x3,constraint_residual");
        for (int k = 0; k <= steps; k++) {
            double t = k * h;
            double u = uOfT(t);
            double x3 = algebraicX3(xd, u);
            double residual = xd[0] + x3 - u;
            if (k % 40 == 0) {
                System.out.printf("%.6f,%.6f,%.6f,%.6f,%.6e%n", t, xd[0], xd[1], x3, residual);
            }
            xd = rk4Step(t, xd, h);
        }
    }
}

MATLAB directly supports generalized eigenvalues through eig(A,E). For simulation, an index-1 descriptor model can be reduced explicitly as shown below, or represented in Simulink with Integrator blocks for differential states and an Algebraic Constraint block for the constraint equation.

% Chapter30_Lesson1.m
% Descriptor (singular) system example: E xdot = A x + B u.
% This script computes generalized eigenvalues, performs an index-1 reduction,
% simulates the reduced ODE, and reconstructs the algebraic state.

clear; clc; close all;

E = diag([1 1 0]);
A = [ 0  1  0;
     -2 -3  1;
      1  0  1];
B = [0; 1; -1];

lambda = eig(A, E);
disp('Generalized eigenvalues of (A,E):');
disp(lambda);

A11 = A(1:2, 1:2);
A12 = A(1:2, 3);
A21 = A(3, 1:2);
A22 = A(3, 3);
B1 = B(1:2);
B2 = B(3);

% Algebraic equation: 0 = A21*xd + A22*xa + B2*u.
% Since A22 is nonsingular, xa = -A22^(-1)*A21*xd - A22^(-1)*B2*u.
Ar = A11 - A12 * (A22 \ A21);
Br = B1 - A12 * (A22 \ B2);

disp('Reduced matrix Ar:');
disp(Ar);
disp('Reduced input matrix Br:');
disp(Br);
disp('Reduced eigenvalues:');
disp(eig(Ar));

u = @(t) double(t >= 1.0);
rhs = @(t, xd) Ar*xd + Br*u(t);
reconstruct_xa = @(xd, input_value) -A22 \ (A21*xd + B2*input_value);

xd0 = [0.2; 0.0];
xa0 = reconstruct_xa(xd0, u(0));
disp('Consistent initial condition [xd; xa]:');
disp([xd0; xa0]);

[t, xd] = ode45(rhs, [0 8], xd0);
x3 = zeros(length(t), 1);
residual = zeros(length(t), 1);
for k = 1:length(t)
    x3(k) = reconstruct_xa(xd(k, :).', u(t(k)));
    residual(k) = A21*xd(k, :).'+ A22*x3(k) + B2*u(t(k));
end
fprintf('Maximum algebraic constraint residual: %.3e\n', max(abs(residual)));

figure;
plot(t, xd(:, 1), 'DisplayName', 'x1 differential'); hold on;
plot(t, xd(:, 2), 'DisplayName', 'x2 differential');
plot(t, x3, 'DisplayName', 'x3 algebraic');
grid on;
xlabel('time');
ylabel('state');
title('Descriptor system reduced to an index-1 ODE');
legend('Location', 'best');

% Simulink note:
% For a block diagram implementation, represent the first two equations by
% Integrator blocks and implement 0 = x1 + x3 - u with an Algebraic Constraint
% block or by explicitly computing x3 = -x1 + u for this index-1 example.

Mathematica can compute generalized eigenvalues symbolically or numerically and can simulate the reduced ODE with NDSolveValue.

Notebook[{
 Cell["Chapter30_Lesson1.nb", "Title"],
 Cell["Descriptor (Singular) Systems: E xdot = A x + B u", "Subtitle"],
 Cell["This notebook computes generalized eigenvalues, performs index-1 reduction, and simulates the reduced ODE.", "Text"],
 Cell["E = DiagonalMatrix[{1, 1, 0}];\nA = { {0, 1, 0}, {-2, -3, 1}, {1, 0, 1} };\nB = { {0}, {1}, {-1} };", "Input"],
 Cell["generalizedEigenvalues = Eigenvalues[{A, E}]", "Input"],
 Cell["A11 = A[[1 ;; 2, 1 ;; 2]];\nA12 = A[[1 ;; 2, 3 ;; 3]];\nA21 = A[[3 ;; 3, 1 ;; 2]];\nA22 = A[[3 ;; 3, 3 ;; 3]];\nB1 = B[[1 ;; 2]];\nB2 = B[[3 ;; 3]];\nAr = A11 - A12 . Inverse[A22] . A21;\nBr = B1 - A12 . Inverse[A22] . B2;\n{Ar, Br, Eigenvalues[Ar]}", "Input"],
 Cell["u[t_] := Piecewise[{ {0, t < 1}, {1, t >= 1} }];\nsol = NDSolveValue[{x1'[t] == x2[t], x2'[t] == -3 x1[t] - 3 x2[t] + 2 u[t], x1[0] == 0.2, x2[0] == 0}, {x1, x2}, {t, 0, 8}];", "Input"],
 Cell["x3[t_] := -sol[[1]][t] + u[t];\nPlot[Evaluate[{sol[[1]][t], sol[[2]][t], x3[t]}], {t, 0, 8}, PlotLegends -> {\"x1 differential\", \"x2 differential\", \"x3 algebraic\"}, GridLines -> Automatic]", "Input"],
 Cell["constraintResidual[t_] := sol[[1]][t] + x3[t] - u[t];\nMax[Abs[constraintResidual /@ Subdivide[0, 8, 100]]]", "Input"]
}]

Problem 1 (Regularity): Let \( \mathbf{E}=\begin{bmatrix}1&0\\0&0\end{bmatrix} \) and \( \mathbf{A}=\begin{bmatrix}0&1\\1&0\end{bmatrix} \). Determine whether the pencil \( \lambda\mathbf{E}-\mathbf{A} \) is regular.

Solution:

\[ \lambda\mathbf{E}-\mathbf{A}=\begin{bmatrix}\lambda&-1\\-1&0\end{bmatrix}, \qquad \det(\lambda\mathbf{E}-\mathbf{A})=-1. \]

Since the determinant is not the zero polynomial, the pencil is regular. It has no finite eigenvalue because the determinant is constant. Both generalized eigenvalues are therefore infinite.

Problem 2 (Consistent initial condition): For the worked example in Section 6, suppose \( u(0)=0 \) and \( x_1(0)=2 \). Find the required value of \( x_3(0) \).

Solution:

\[ 0=x_1(0)+x_3(0)-u(0). \]

Substituting \( x_1(0)=2 \) and \( u(0)=0 \) gives \( x_3(0)=-2 \). Thus \( [2,\;x_2(0),\;-2]' \) is consistent for any chosen \( x_2(0) \).

Problem 3 (Index-1 reduction): Consider

\[ \dot{x}_1=x_2+z,\qquad 0=x_1+2z+u. \]

Reduce the system to an ODE in \( x_1 \) and \( x_2 \) assuming \( \dot{x}_2=-x_1-x_2 \).

Solution:

\[ z=-\frac{1}{2}x_1-\frac{1}{2}u. \]

Therefore,

\[ \dot{x}_1=x_2-\frac{1}{2}x_1-\frac{1}{2}u, \qquad \dot{x}_2=-x_1-x_2. \]

In matrix form,

\[ \begin{bmatrix}\dot{x}_1\\\dot{x}_2\end{bmatrix} =\begin{bmatrix}-1/2&1\\-1&-1\end{bmatrix} \begin{bmatrix}x_1\\x_2\end{bmatrix} +\begin{bmatrix}-1/2\\0\end{bmatrix}u. \]

Problem 4 (Stability of the reduced example): Show that the reduced matrix \( \mathbf{A}_r=\begin{bmatrix}0&1\\-3&-3\end{bmatrix} \) is asymptotically stable.

Solution:

\[ \det(s\mathbf{I}-\mathbf{A}_r)=s^2+3s+3. \]

The roots are \( s_{1,2}=(-3\pm j\sqrt{3})/2 \). Both have negative real part, so the finite dynamics are asymptotically stable.

Problem 5 (Why not invert E?): Explain why the formula \( \dot{\mathbf{x} }=\mathbf{E}^{-1}\mathbf{A}\mathbf{x}+\mathbf{E}^{-1}\mathbf{B}\mathbf{u} \) is invalid for a singular descriptor system.

Solution:

If \( \operatorname{rank}(\mathbf{E})<n \), then \( \mathbf{E}^{-1} \) does not exist. More importantly, the zero rows of \( \mathbf{E} \) encode algebraic constraints, not missing dynamics. Treating the model as an ordinary ODE would discard the consistency manifold and can create spurious modes or invalid simulations.

Descriptor systems extend the standard state-space form by allowing a singular derivative matrix \( \mathbf{E} \). This introduces algebraic constraints, generalized eigenvalues, finite and infinite modes, and consistency requirements for initial conditions. In the important index-1 case, the algebraic variables can be eliminated when \( \mathbf{A}_{22} \) is nonsingular, producing a reduced standard ODE while preserving the original constraint relation. Later advanced courses can build on these ideas for descriptor controllability, observability, feedback regularization, robust control, and numerical DAE simulation.