Consider the continuous-time LTI state-space model (from Chapter 8):

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

Stability of the internal (state) dynamics is determined by the homogeneous system obtained by setting \( \mathbf{u}(t)=\mathbf{0} \):

\[ \dot{\mathbf{x}}(t) = \mathbf{A}\mathbf{x}(t), \quad \mathbf{x}(0)=\mathbf{x}_0. \]

Using the stability concepts from Lesson 1, we will characterize:

• Asymptotic (and exponential) stability: whether \( \mathbf{x}(t) → \mathbf{0} \) as \( t → \infty \).
• Lyapunov (internal) stability: whether small \( \|\mathbf{x}_0\| \) implies uniformly small \( \|\mathbf{x}(t)\| \) for all \( t \ge 0 \).
• Instability: whether some initial conditions yield growth in time.

The central object will be the matrix exponential \( e^{\mathbf{A}t} \), whose behavior is governed by the spectrum \( \sigma(\mathbf{A})=\{\lambda_1,\dots,\lambda_n\} \).

The unique solution of \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \) is

\[ \mathbf{x}(t) = e^{\mathbf{A}t}\mathbf{x}_0, \quad e^{\mathbf{A}t} \triangleq \sum_{k=0}^{\infty}\frac{(\mathbf{A}t)^k}{k!}. \]

The series definition is always convergent (for any finite matrix \( \mathbf{A} \)) and provides analytic dependence on \( t \). The following property justifies the solution formula.

Proposition 2.1 (Derivative property). For all real \( t \),

\[ \frac{d}{dt}e^{\mathbf{A}t} = \mathbf{A}e^{\mathbf{A}t} = e^{\mathbf{A}t}\mathbf{A}. \]

Proof. Differentiate the defining series term-by-term (justified by uniform convergence on bounded intervals):

\[ \frac{d}{dt}\left(\sum_{k=0}^{\infty}\frac{(\mathbf{A}t)^k}{k!}\right) = \sum_{k=1}^{\infty}\frac{k(\mathbf{A}t)^{k-1}\mathbf{A}}{k!} = \sum_{k=1}^{\infty}\frac{(\mathbf{A}t)^{k-1}}{(k-1)!}\mathbf{A} = \left(\sum_{m=0}^{\infty}\frac{(\mathbf{A}t)^m}{m!}\right)\mathbf{A} = e^{\mathbf{A}t}\mathbf{A}. \]

Since \( \mathbf{A} \) commutes with its powers, \( e^{\mathbf{A}t}\mathbf{A}=\mathbf{A}e^{\mathbf{A}t} \). Finally, with \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0 \), we obtain \( \dot{\mathbf{x}}(t)=\mathbf{A}e^{\mathbf{A}t}\mathbf{x}_0=\mathbf{A}\mathbf{x}(t) \) and \( \mathbf{x}(0)=e^{\mathbf{A}0}\mathbf{x}_0=\mathbf{I}\mathbf{x}_0=\mathbf{x}_0 \). □

Therefore, stability of the origin is equivalent to boundedness/decay properties of \( e^{\mathbf{A}t} \) as \( t → \infty \).

Suppose \( \mathbf{A} \) is diagonalizable over \( \mathbb{C} \), i.e., there exists an invertible \( \mathbf{V} \) and diagonal \( \boldsymbol{\Lambda}=\mathrm{diag}(\lambda_1,\dots,\lambda_n) \) such that \( \mathbf{A}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1} \). Then

\[ e^{\mathbf{A}t} = e^{\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1}t} = \mathbf{V}e^{\boldsymbol{\Lambda}t}\mathbf{V}^{-1}, \quad e^{\boldsymbol{\Lambda}t}=\mathrm{diag}\left(e^{\lambda_1 t},\dots,e^{\lambda_n t}\right). \]

Hence the solution decomposes into exponentials (modes). If we write \( \mathbf{c}=\mathbf{V}^{-1}\mathbf{x}_0 \), then

\[ \mathbf{x}(t) = \mathbf{V}\,\mathrm{diag}\left(e^{\lambda_1 t},\dots,e^{\lambda_n t}\right)\mathbf{c} = \sum_{i=1}^{n} c_i e^{\lambda_i t}\mathbf{v}_i, \]

where \( \mathbf{v}_i \) is the \( i \)-th column of \( \mathbf{V} \). This shows that each eigenvalue generates an exponential component.

flowchart TD
  X0["Initial state x0"] --> DECOMP["Compute c = V^{-1} x0"]
  DECOMP --> MODES["Form modes: c_i * exp(lambda_i * t) * v_i"]
  MODES --> SUM["Sum all modes -> x(t)"]
  SUM --> OUT["State trajectory x(t)"]
        

The decay/growth of \( e^{\lambda_i t} \) is determined by \( \Re(\lambda_i) \). In particular:

\[ \Re(\lambda_i) < 0 \Rightarrow e^{\lambda_i t} → 0 \text{ as } t → \infty, \qquad \Re(\lambda_i) > 0 \Rightarrow |e^{\lambda_i t}| → \infty. \]

We now make this into a rigorous stability statement (including matrix conditioning).

Theorem 3.1 (Diagonalizable sufficient condition). If \( \mathbf{A} \) is diagonalizable and \( \max_i \Re(\lambda_i) < 0 \), then the origin is exponentially stable: there exist constants \( M > 0 \) and \( \alpha > 0 \) such that

\[ \|\mathbf{x}(t)\| \le M e^{-\alpha t}\|\mathbf{x}_0\|, \quad \forall t \ge 0. \]

Proof. Let \( \mathbf{A}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1} \). Then \( \|e^{\mathbf{A}t}\|\le \|\mathbf{V}\|\,\|e^{\boldsymbol{\Lambda}t}\|\,\|\mathbf{V}^{-1}\| \). Using the induced matrix norm compatible with vector norms,

\[ \|e^{\boldsymbol{\Lambda}t}\| = \max_i |e^{\lambda_i t}| = \exp\!\left(t\max_i \Re(\lambda_i)\right). \]

Choose \( \alpha \) such that \( 0 < \alpha < -\max_i \Re(\lambda_i) \). Then \( \exp(t\max_i \Re(\lambda_i)) \le e^{-\alpha t} \) for all \( t\ge 0 \). Let \( M=\|\mathbf{V}\|\,\|\mathbf{V}^{-1}\| \). We obtain \( \|e^{\mathbf{A}t}\| \le M e^{-\alpha t} \), and hence the claim. □

Diagonalizability is not guaranteed. To obtain a complete criterion, use the Jordan decomposition (from linear algebra): there exists an invertible \( \mathbf{P} \) such that

\[ \mathbf{A} = \mathbf{P}\mathbf{J}\mathbf{P}^{-1}, \quad \mathbf{J}=\mathrm{blkdiag}(\mathbf{J}_1,\dots,\mathbf{J}_r), \]

where each Jordan block associated to eigenvalue \( \lambda \) has the form

\[ \mathbf{J}_\lambda = \lambda\mathbf{I} + \mathbf{N}, \quad \mathbf{N} = \begin{bmatrix} 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & & \ddots & \ddots & \vdots\\ 0 & \cdots & 0 & 0 & 1\\ 0 & \cdots & \cdots & 0 & 0 \end{bmatrix}, \quad \mathbf{N}^m=\mathbf{0} \text{ for an } m\times m \text{ block}. \]

The exponential of a Jordan block can be computed explicitly using nilpotency of \( \mathbf{N} \).

Lemma 4.1 (Jordan-block exponential). For an \( m\times m \) Jordan block, \( e^{\mathbf{J}_\lambda t} = e^{\lambda t}\sum_{k=0}^{m-1}\frac{(\mathbf{N}t)^k}{k!} \).

Proof. Since \( \lambda\mathbf{I} \) commutes with \( \mathbf{N} \), we have \( e^{(\lambda\mathbf{I}+\mathbf{N})t}=e^{\lambda t}e^{\mathbf{N}t} \). But \( \mathbf{N}^m=\mathbf{0} \), so the series truncates:

\[ e^{\mathbf{N}t} = \sum_{k=0}^{\infty}\frac{(\mathbf{N}t)^k}{k!} = \sum_{k=0}^{m-1}\frac{(\mathbf{N}t)^k}{k!}. \]

Multiplying yields the claimed expression. □

This reveals an essential phenomenon: Jordan blocks introduce polynomial factors \( t^k \) multiplying exponentials. Consequently:

• If \( \Re(\lambda) < 0 \), then \( t^k e^{\lambda t} → 0 \) for any finite \( k \).
• If \( \Re(\lambda)=0 \) and the block size \( m > 1 \), then \( \|e^{\mathbf{J}_\lambda t}\| grows at least like a polynomial (terms such as \( t^{m-1} \)), causing instability.

The matrix \( \mathbf{A} \) is called Hurwitz if all eigenvalues satisfy \( \Re(\lambda_i) < 0 \).

Theorem 5.1 (Hurwitz ⇔ exponential stability). For the continuous-time LTI system \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \), the following are equivalent:

1. \( \mathbf{A} \) is Hurwitz (i.e., \( \max_i \Re(\lambda_i) < 0 \)).
2. The origin is exponentially stable: there exist \( M > 0 \) and \( \alpha > 0 \) such that \( \|e^{\mathbf{A}t}\| \le M e^{-\alpha t} \) for all \( t\ge 0 \).

Proof (1 ⇒ 2). Write the Jordan decomposition \( \mathbf{A}=\mathbf{P}\mathbf{J}\mathbf{P}^{-1} \). Then \( e^{\mathbf{A}t}=\mathbf{P}e^{\mathbf{J}t}\mathbf{P}^{-1} \) and \( \|e^{\mathbf{A}t}\|\le \|\mathbf{P}\|\,\|e^{\mathbf{J}t}\|\,\|\mathbf{P}^{-1}\| \). Since \( \mathbf{J} \) is block diagonal, it suffices to bound each Jordan block. For a block of size \( m \) with eigenvalue \( \lambda \), Lemma 4.1 gives

\[ \|e^{\mathbf{J}_\lambda t}\| \le |e^{\lambda t}|\sum_{k=0}^{m-1}\frac{\|\mathbf{N}\|^k t^k}{k!} \le e^{\Re(\lambda)t}\,p_{m-1}(t), \]

where \( p_{m-1}(t) \) is a polynomial of degree \( m-1 \). Let \( \beta = -\max_i \Re(\lambda_i) \); Hurwitz means \( \beta > 0 \). Pick any \( 0 < \alpha < \beta \). Then for each block, \( e^{\Re(\lambda)t}p_{m-1}(t) \le K e^{-\alpha t} \) for some finite constant \( K \) (because exponential decay dominates any fixed-degree polynomial). Taking the maximum over finitely many blocks yields \( \|e^{\mathbf{J}t}\|\le K_J e^{-\alpha t} \). Hence \( \|e^{\mathbf{A}t}\|\le (\|\mathbf{P}\|\,\|\mathbf{P}^{-1}\|K_J)e^{-\alpha t} \). Set \( M=\|\mathbf{P}\|\,\|\mathbf{P}^{-1}\|K_J \). □

Proof (2 ⇒ 1). Suppose exponential stability holds but \( \mathbf{A} \) has an eigenvalue \( \lambda \) with \( \Re(\lambda) \ge 0 \). Let \( \mathbf{v}\ne \mathbf{0} \) be an eigenvector: \( \mathbf{A}\mathbf{v}=\lambda\mathbf{v} \). Then \( e^{\mathbf{A}t}\mathbf{v}=e^{\lambda t}\mathbf{v} \) (because the exponential series preserves eigenvectors). Taking norms:

\[ \|e^{\mathbf{A}t}\mathbf{v}\| = |e^{\lambda t}|\|\mathbf{v}\| = e^{\Re(\lambda)t}\|\mathbf{v}\|. \]

If \( \Re(\lambda) > 0 \), the right side grows unbounded, contradicting \( \|e^{\mathbf{A}t}\mathbf{v}\| \le \|e^{\mathbf{A}t}\|\,\|\mathbf{v}\| \le M e^{-\alpha t}\|\mathbf{v}\| \). If \( \Re(\lambda)=0 \), the right side is constant and cannot be upper-bounded by \( M e^{-\alpha t}\|\mathbf{v}\| \) for all large \( t \) since the latter decays to zero. Hence no eigenvalue can satisfy \( \Re(\lambda)\ge 0 \). Therefore \( \mathbf{A} \) is Hurwitz. □

Exponential stability is the most important case in control engineering because it implies robustness to small perturbations and yields useful bounds. However, linear systems also admit a clean classification of non-decaying cases.

Theorem 6.1 (Stability classification for \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \)). Let \( \lambda_i \) be eigenvalues of \( \mathbf{A} \).

• Asymptotically (exponentially) stable iff \( \Re(\lambda_i) < 0 \) for all \( i \).
• Lyapunov stable but not asymptotically stable iff \( \Re(\lambda_i) \le 0 \) for all \( i \) and every eigenvalue with \( \Re(\lambda_i)=0 \) is semisimple (its Jordan blocks are all size 1).
• Unstable otherwise (i.e., some \( \Re(\lambda_i) > 0 \), or some \( \Re(\lambda_i)=0 \) has a Jordan block of size \( > 1 \)).

Proof sketch (based on Jordan blocks).

• If any \( \Re(\lambda) > 0 \), then \( e^{\lambda t} \) grows exponentially along its eigenvector/generalized eigenvector direction, implying instability.
• If \( \Re(\lambda)=0 \) with a Jordan block size \( m > 1 \), Lemma 4.1 yields polynomial growth (terms like \( t^{m-1} \)), implying instability.
• If all \( \Re(\lambda) \le 0 \) and all \( \Re(\lambda)=0 \) are semisimple, then \( e^{\mathbf{A}t} \) remains bounded (oscillatory terms from purely imaginary eigenvalues can persist), so trajectories remain bounded and the origin is Lyapunov stable, but not asymptotically stable unless all real parts are strictly negative.

This theorem explains an important engineering warning: purely imaginary eigenvalues are not automatically “safe”; the Jordan structure decides whether bounded oscillations or unbounded drift occurs.

From Chapter 8 (state-space) and Chapter 6 (transfer functions), the transfer matrix of an LTI realization is

\[ \mathbf{G}(s) = \mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B} + \mathbf{D}. \]

The resolvent \( (s\mathbf{I}-\mathbf{A})^{-1} \) exists exactly when \( \det(s\mathbf{I}-\mathbf{A})\ne 0 \). Thus the eigenvalues of \( \mathbf{A} \) are the roots of the characteristic polynomial:

\[ \det(s\mathbf{I}-\mathbf{A}) = 0 \quad \Longleftrightarrow \quad s \in \sigma(\mathbf{A}). \]

Moreover, using the adjugate identity, \( (s\mathbf{I}-\mathbf{A})^{-1} = \frac{\mathrm{adj}(s\mathbf{I}-\mathbf{A})}{\det(s\mathbf{I}-\mathbf{A})} \), each entry of \( \mathbf{G}(s) \) is a rational function with denominator dividing \( \det(s\mathbf{I}-\mathbf{A}) \). Therefore:

• Poles of \( \mathbf{G}(s) \) are a subset of \( \sigma(\mathbf{A}) \).
• If the realization is minimal (no pole-zero cancellations), then the poles of \( \mathbf{G}(s) \) coincide with the eigenvalues of \( \mathbf{A} \).

This explains consistency with Lesson 2 (pole locations in the \( s \)-plane): eigenvalues provide an internal-state perspective on the same stability mechanism.

In computation, we typically use eigendecomposition (or Schur decomposition) to assess stability. The robust classification uses:

• Spectral test: compute \( \lambda_i \) and check \( \Re(\lambda_i) \).
• If any \( \Re(\lambda_i) \) is near zero (within tolerance), test semisimplicity (Jordan defects), e.g., via rank/nullity checks on \( \mathbf{A}-\lambda\mathbf{I} \) and its powers.
• Beware conditioning: nearly defective matrices can make eigenvectors ill-conditioned, inflating transient response even when \( \Re(\lambda_i) < 0 \).

flowchart TD
  A["Given A matrix"] --> EIG["Compute eigenvalues lambda_i"]
  EIG --> TEST1["Any Re(lambda_i) > 0 ?"]
  TEST1 -->|yes| UNST["Unstable"]
  TEST1 -->|no| TEST2["Any Re(lambda_i) < 0 for all i ?"]
  TEST2 -->|yes| EXPST["Exponentially stable"]
  TEST2 -->|no| TEST3["Eigenvalues with Re(lambda)=0 exist"]
  TEST3 --> SEMI["Check Jordan blocks (semisimple?)"]
  SEMI -->|yes| LYAP["Lyapunov stable (bounded)"]
  SEMI -->|no| UNST
        

A useful quantitative measure is the spectral abscissa \( \alpha(\mathbf{A}) \triangleq \max_i \Re(\lambda_i) \). Then:

\[ \alpha(\mathbf{A}) < 0 \Rightarrow \text{exponential stability}, \quad \alpha(\mathbf{A}) > 0 \Rightarrow \text{instability}. \]

The boundary case \( \alpha(\mathbf{A})=0 \) requires the semisimplicity test described above.

Recommended libraries for this lesson: numpy, scipy.linalg (matrix exponential, Schur forms), control (state-space utilities), and sympy (symbolic eigen-analysis when needed).

Example system: a standard second-order mass–spring–damper written in state variables \( \mathbf{x}=[x,\dot{x}]^\top \) (modeled earlier in Chapters 3 and 8). For mass \( m \), damping \( c \), spring \( k \):

\[ \dot{\mathbf{x}} = \underbrace{\begin{bmatrix} 0 & 1\\ -k/m & -c/m \end{bmatrix}}_{\mathbf{A}} \mathbf{x} + \underbrace{\begin{bmatrix} 0\\ 1/m\end{bmatrix}}_{\mathbf{B}}u, \quad y = \underbrace{\begin{bmatrix}1 & 0\end{bmatrix}}_{\mathbf{C}}\mathbf{x}. \]

import numpy as np
from numpy.linalg import eig
from scipy.linalg import expm
import math

def classify_stability(A, tol=1e-9):
    # Returns: "exp_stable", "lyap_stable", or "unstable"
    lam = eig(A)[0]
    re = np.real(lam)

    if np.any(re > tol):
        return "unstable", lam

    if np.all(re < -tol):
        return "exp_stable", lam

    # Boundary case: some Re(lambda) approximately zero.
    # A conservative semisimplicity test: compare geometric vs algebraic multiplicity
    # for eigenvalues near the imaginary axis (Re ~ 0).
    # (Numerical caution: nearly defective matrices can be hard to classify.)
    near = [l for l in lam if abs(np.real(l)) <= 10*tol]
    for l in near:
        # algebraic multiplicity (count close eigenvalues)
        alg = sum(abs(l - lj) <= 1e-6 for lj in lam)

        # geometric multiplicity = nullity of (A - l I)
        M = A - l*np.eye(A.shape[0])
        s = np.linalg.svd(M, compute_uv=False)
        rank = np.sum(s > 1e-7)
        geo = A.shape[0] - rank

        if geo < alg:
            return "unstable", lam

    return "lyap_stable", lam

# Mass-spring-damper parameters
m = 1.0
c = 0.6
k = 4.0

A = np.array([[0.0, 1.0],
              [-k/m, -c/m]])

cls, lam = classify_stability(A)
print("Eigenvalues:", lam)
print("Classification:", cls)

# Simulate homogeneous response x(t) = expm(A t) x0
x0 = np.array([1.0, 0.0])
ts = np.linspace(0.0, 10.0, 400)
X = np.array([expm(A*t) @ x0 for t in ts])  # states over time

# Report a stability-consistent metric
alpha = np.max(np.real(lam))
print("Spectral abscissa alpha(A):", alpha)

# Optional: compare with python-control (state-space utilities)
try:
    import control
    sys = control.ss(A, np.array([[0.0],[1.0/m]]), np.array([[1.0, 0.0]]), np.array([[0.0]]))
    print("control.ss poles:", control.poles(sys))
except Exception as e:
    print("python-control not available or failed:", e)

Notes:

• The classification is numerically conservative near \( \Re(\lambda)=0 \), because defect detection is ill-conditioned. In high-stakes settings, use a Schur-based analysis and verify behavior by bounding \( \|e^{\mathbf{A}t}\| \) or simulating multiple initial conditions.
• Transient growth can occur even when exponentially stable if eigenvectors are ill-conditioned; exponential stability is asymptotic, not necessarily monotonic decay.

Recommended libraries: Eigen (linear algebra; includes eigenvalues and an optional matrix exponential in unsupported modules), and optionally Boost for ODE integration (not required here).

#include <iostream>
#include <complex>
#include <Eigen/Dense>

// If available in your Eigen installation:
//   #include <unsupported/Eigen/MatrixFunctions>

int main() {
    using Eigen::Matrix2d;
    using Eigen::Vector2d;

    double m = 1.0, c = 0.6, k = 4.0;

    Matrix2d A;
    A << 0.0, 1.0,
        -k/m, -c/m;

    Eigen::EigenSolver<Matrix2d> es(A);
    auto evals = es.eigenvalues();

    std::cout << "Eigenvalues:\n";
    for (int i = 0; i < evals.size(); ++i) {
        std::cout << evals(i) << "\n";
    }

    // Classification by spectral abscissa
    double alpha = -1e300;
    for (int i = 0; i < evals.size(); ++i) {
        alpha = std::max(alpha, evals(i).real());
    }
    std::cout << "alpha(A) = " << alpha << "\n";
    if (alpha < 0.0) std::cout << "Exponentially stable (Hurwitz)\n";
    else if (alpha > 0.0) std::cout << "Unstable\n";
    else std::cout << "Boundary case: check Jordan/semisimplicity\n";

    // Exact propagation (if MatrixFunctions is available):
    // double t = 2.0;
    // Eigen::Matrix2d Phi = (A * t).exp();
    // Vector2d x0(1.0, 0.0);
    // Vector2d xt = Phi * x0;
    // std::cout << "x(t)=" << xt.transpose() << "\n";

    return 0;
}

If your toolchain lacks unsupported/Eigen/MatrixFunctions, you can implement a basic scaling-and-squaring matrix exponential (later refined in numerical simulation lessons) or propagate using a small-step ODE integrator. The stability criterion itself depends only on eigenvalues/Jordan structure, not on the chosen numerical propagator.

Recommended libraries: EJML (efficient linear algebra, eigenvalue decomposition) or Apache Commons Math (linear algebra utilities). Below is an EJML-based example for eigenvalues (propagation can be done using a matrix exponential implementation or numerical integration; we focus on eigenvalue classification here).

import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.factory.DecompositionFactory_DDRM;
import org.ejml.interfaces.decomposition.EigenDecomposition_F64;

public class StabilityCheck {
    public static void main(String[] args) {
        double m = 1.0, c = 0.6, k = 4.0;

        DMatrixRMaj A = new DMatrixRMaj(new double[][]{
                {0.0,  1.0},
                {-k/m, -c/m}
        });

        EigenDecomposition_F64<DMatrixRMaj> eig =
                DecompositionFactory_DDRM.eig(A.numRows, false);

        if (!eig.decompose(A)) {
            throw new RuntimeException("Eigen decomposition failed");
        }

        double alpha = -1e300;
        System.out.println("Eigenvalues:");
        for (int i = 0; i < A.numRows; i++) {
            var val = eig.getEigenvalue(i);
            System.out.println(val);
            alpha = Math.max(alpha, val.getReal());
        }

        System.out.println("alpha(A) = " + alpha);
        if (alpha < 0.0) System.out.println("Exponentially stable (Hurwitz)");
        else if (alpha > 0.0) System.out.println("Unstable");
        else System.out.println("Boundary case: check semisimplicity (Jordan defects)");
    }
}

For exact state propagation in Java, a standard approach is to implement scaling-and-squaring with Padé approximation for the matrix exponential (the same mathematical object as in Section 2). Many numerical libraries also provide matrix exponential utilities; however, the stability decision should always be based on eigen-structure, not on a single numerical trajectory.

MATLAB provides direct access to eigenvalues and matrix exponentials, and the Control System Toolbox supports state-space objects. Simulink provides a State-Space block for simulation.

% Mass-spring-damper example
m = 1.0; c = 0.6; k = 4.0;

A = [0 1; -k/m -c/m];
B = [0; 1/m];
C = [1 0];
D = 0;

lambda = eig(A);
disp('Eigenvalues:'); disp(lambda);

alpha = max(real(lambda));
disp('alpha(A):'); disp(alpha);

if alpha < 0
    disp('Exponentially stable (Hurwitz).');
elseif alpha > 0
    disp('Unstable.');
else
    disp('Boundary case: check Jordan defects (semisimplicity).');
end

% Exact homogeneous solution x(t) = expm(A t) x0
x0 = [1; 0];
t = linspace(0,10,400);
X = zeros(length(t),2);
for i=1:length(t)
    X(i,:) = (expm(A*t(i))*x0).';
end

% State-space object and pole check
sys = ss(A,B,C,D);
disp('Poles from ss object:'); disp(pole(sys));

Simulink guidance (no images):

• Open Simulink and create a new model.
• Add a State-Space block and set A, B, C, D as above.
• Set initial conditions to \( \mathbf{x}_0 \) (block parameter).
• Connect the output to a Scope block; run and observe boundedness/decay.

You may also create the model programmatically (useful for reproducibility):

% Programmatic Simulink model creation (basic)
model = 'EigenStabilityDemo';
new_system(model); open_system(model);

add_block('simulink/Continuous/State-Space', [model '/StateSpace']);
set_param([model '/StateSpace'], 'A', mat2str(A), 'B', mat2str(B), 'C', mat2str(C), 'D', mat2str(D));
set_param([model '/StateSpace'], 'X0', mat2str(x0));

add_block('simulink/Sinks/Scope', [model '/Scope']);

add_line(model, 'StateSpace/1', 'Scope/1');

set_param(model, 'StopTime', '10');
save_system(model);

Mathematica directly supports eigenvalues and matrix exponentials. Use Eigenvalues and MatrixExp.

m = 1.0; c = 0.6; k = 4.0;

A = { {0, 1}, {-k/m, -c/m} };
lambda = Eigenvalues[A]
alpha = Max[Re /@ lambda]

If[alpha < 0, "Exponentially stable (Hurwitz)",
 If[alpha > 0, "Unstable",
  "Boundary case: check Jordan defects (semisimplicity)"]]

x0 = {1, 0};
x[t_] := MatrixExp[A t].x0

(* Sample trajectory values *)
Table[{tt, x[tt]}, {tt, 0, 5, 0.5}]

Problem 1 (Hurwitz test by eigenvalues): Consider \( \mathbf{A}=\begin{bmatrix}-2 & 1\\ 0 & -3\end{bmatrix} \). Classify the stability of \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \).

Solution: The eigenvalues of an upper triangular matrix are its diagonal entries: \( \lambda_1=-2 \), \( \lambda_2=-3 \). Hence \( \Re(\lambda_i) < 0 \) for all \( i \), so \( \mathbf{A} \) is Hurwitz and the origin is exponentially stable by Theorem 5.1.

Problem 2 (Jordan block on the imaginary axis): Let \( \mathbf{A}=\begin{bmatrix}0 & 1\\ 0 & 0\end{bmatrix} \). Determine whether the origin is Lyapunov stable, asymptotically stable, or unstable.

Solution: The characteristic polynomial is \( \lambda^2=0 \), so \( \lambda=0 \) has algebraic multiplicity 2. However, \( \mathbf{A} \) is already a Jordan block of size 2 (defective), so the eigenvalue on the imaginary axis is not semisimple. Using Lemma 4.1 with \( \lambda=0 \) and \( \mathbf{N}=\mathbf{A} \),

\[ e^{\mathbf{A}t} = \mathbf{I} + \mathbf{A}t = \begin{bmatrix}1 & t\\ 0 & 1\end{bmatrix}. \]

Then for \( \mathbf{x}_0=[0,1]^\top \), we get \( \mathbf{x}(t)=[t,1]^\top \), which grows unbounded. Therefore the origin is unstable.

Problem 3 (Marginally stable oscillation): Consider \( \mathbf{A}=\begin{bmatrix}0 & -\omega\\ \omega & 0\end{bmatrix} \) with \( \omega \ne 0 \). Classify stability and describe trajectories.

Solution: The eigenvalues are \( \lambda_{1,2}=\pm j\omega \), so \( \Re(\lambda_{1,2})=0 \). The matrix is diagonalizable over \( \mathbb{C} \) (no defective Jordan blocks), hence the imaginary-axis eigenvalues are semisimple. By Theorem 6.1, the origin is Lyapunov stable but not asymptotically stable. Trajectories are bounded rotations (constant energy) in the plane.

Problem 4 (Exponential decay bound): Suppose \( \mathbf{A} \) is diagonalizable and \( \max_i \Re(\lambda_i) = -\beta \) with \( \beta > 0 \). Show that there exist \( M > 0 \) and \( \alpha \in (0,\beta) \) such that \( \|e^{\mathbf{A}t}\| \le M e^{-\alpha t} \).

Solution: Using Section 3, write \( \mathbf{A}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1} \). Then \( \|e^{\mathbf{A}t}\|\le \|\mathbf{V}\|\,\|\mathbf{V}^{-1}\|\,\|e^{\boldsymbol{\Lambda}t}\| \) and \( \|e^{\boldsymbol{\Lambda}t}\|=\exp(t\max_i \Re(\lambda_i))=\exp(-\beta t) \). Choose any \( \alpha \in (0,\beta) \). Then \( e^{-\beta t} \le e^{-\alpha t} \) for all \( t\ge 0 \). Set \( M=\|\mathbf{V}\|\,\|\mathbf{V}^{-1}\| \). This yields the bound.

Problem 5 (State-space vs transfer-function poles): Let \( \mathbf{G}(s)=\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B} \) (assume \( \mathbf{D}=\mathbf{0} \)). Prove that any pole of any entry of \( \mathbf{G}(s) \) must satisfy \( \det(s\mathbf{I}-\mathbf{A})=0 \).

Solution: Using the adjugate identity,

\[ (s\mathbf{I}-\mathbf{A})^{-1} = \frac{\mathrm{adj}(s\mathbf{I}-\mathbf{A})}{\det(s\mathbf{I}-\mathbf{A})}. \]

Therefore, \( \mathbf{G}(s)=\frac{\mathbf{C}\,\mathrm{adj}(s\mathbf{I}-\mathbf{A})\,\mathbf{B}}{\det(s\mathbf{I}-\mathbf{A})} \). Each entry is a rational function whose denominator divides \( \det(s\mathbf{I}-\mathbf{A}) \). Hence any pole must be a root of the denominator, and thus satisfies \( \det(s\mathbf{I}-\mathbf{A})=0 \), i.e., it is an eigenvalue of \( \mathbf{A} \).

We derived eigenvalue-based stability criteria for LTI state-space models by analyzing the matrix exponential \( e^{\mathbf{A}t} \). The key results are:

• \( \Re(\lambda_i) < 0 \) for all eigenvalues ⇔ \( \mathbf{A} \) is Hurwitz ⇔ exponential stability.
• Boundary cases with \( \Re(\lambda)=0 \) require a Jordan-block (semisimplicity) test to decide Lyapunov stability vs instability.
• Transfer-function poles are tied to eigenvalues through \( \det(s\mathbf{I}-\mathbf{A}) \), matching earlier pole-location stability insights.

In the next lesson, we will examine how eigenvalues/poles move under parameter variations, connecting stability margins to robustness insight.