In Lesson 1, stability was defined by the behavior of the free response near the equilibrium. For a continuous-time autonomous LTI system, the equilibrium at the origin is governed by \( \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t) \), and the solution is

\[ \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}(0), \qquad t\ge 0. \]

Thus, stability is determined by whether the matrix exponential \( e^{\mathbf{A}t} \) remains bounded or decays to zero. For a discrete-time LTI system \( \mathbf{x}[k+1]=\mathbf{A}\mathbf{x}[k] \), the free response is

\[ \mathbf{x}[k]=\mathbf{A}^k\mathbf{x}[0], \qquad k=0,1,2,\dots. \]

Therefore, the central question is spectral: how do the eigenvalues and Jordan structure of \( \mathbf{A} \) determine the long-time behavior of \( e^{\mathbf{A}t} \) or \( \mathbf{A}^k \)?

flowchart TD
  A["State matrix A"] --> B["Compute eigenvalues lambda_i"]
  B --> C["Continuous time: \ninspect \nRe(lambda_i)"]
  B --> D["Discrete time: \ninspect \nabs(lambda_i)"]
  C --> E["All \nRe(lambda_i) \nnegative"]
  C --> F["Some \nRe(lambda_i) \npositive"]
  C --> G["Boundary \nRe(lambda_i) \nzero"]
  D --> H["All \nabs(lambda_i) \nless than 1"]
  D --> I["Some \nabs(lambda_i) \ngreater than 1"]
  D --> J["Boundary \nabs(lambda_i) \nequal to 1"]
  E --> K["Asymptotically stable"]
  H --> K
  F --> L["Unstable"]
  I --> L
  G --> M["Check Jordan block sizes"]
  J --> M
  M --> N["Semisimple boundary \nmodes: stable, \nnot asymptotically stable"]
  M --> O["Defective boundary \nmodes: unstable"]
        

Suppose \( \lambda \) is an eigenvalue of \( \mathbf{A} \) with eigenvector \( \mathbf{v}\neq \mathbf{0} \), so that \( \mathbf{A}\mathbf{v}=\lambda\mathbf{v} \). If the initial condition is aligned with this eigenvector, the continuous-time solution becomes a single modal response:

\[ \mathbf{x}(0)=c\mathbf{v} \quad \Longrightarrow \quad \mathbf{x}(t)=c e^{\lambda t}\mathbf{v}. \]

The magnitude of this mode is controlled by \( e^{\operatorname{Re}(\lambda)t} \). Hence modes with \( \operatorname{Re}(\lambda)<0 \) decay, modes with \( \operatorname{Re}(\lambda)>0 \) grow, and modes on the imaginary axis neither decay nor grow exponentially.

For discrete-time systems, the corresponding modal response is

\[ \mathbf{x}[0]=c\mathbf{v} \quad \Longrightarrow \quad \mathbf{x}[k]=c\lambda^k\mathbf{v}. \]

The magnitude is controlled by \( |\lambda|^k \). Thus, the unit circle replaces the imaginary axis: modes inside the unit circle decay, modes outside the unit circle grow, and modes on the unit circle require additional Jordan-structure analysis.

Consider the autonomous continuous-time system \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \), where \( \mathbf{A}\in\mathbb{R}^{n\times n} \). Let \( \lambda_i(\mathbf{A}) \) denote all eigenvalues of \( \mathbf{A} \), counted with algebraic multiplicity.

\[ \boxed{ \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \text{ is asymptotically stable} \iff \operatorname{Re}(\lambda_i(\mathbf{A}))<0 \text{ for every } i. } \]

This condition says that every eigenvalue must lie strictly in the open left half-plane. The phrase strictly is essential: eigenvalues on the imaginary axis produce persistent oscillations or polynomial growth, depending on the size of their Jordan blocks.

Proof idea using Jordan form. Let

\[ \mathbf{A}=\mathbf{T}\mathbf{J}\mathbf{T}^{-1}, \qquad e^{\mathbf{A}t}=\mathbf{T}e^{\mathbf{J}t}\mathbf{T}^{-1}. \]

A Jordan block associated with eigenvalue \( \lambda \) has the form \( \mathbf{J}_\lambda=\lambda\mathbf{I}+\mathbf{N} \), where \( \mathbf{N} \) is nilpotent. Its exponential is

\[ e^{\mathbf{J}_\lambda t} = e^{\lambda t} \left( \mathbf{I}+t\mathbf{N}+\frac{t^2}{2!}\mathbf{N}^2+ \cdots+ \frac{t^{m-1} }{(m-1)!}\mathbf{N}^{m-1} \right), \]

where \( m \) is the Jordan block size. If \( \operatorname{Re}(\lambda)<0 \), then the exponential decay \( e^{\operatorname{Re}(\lambda)t} \) dominates every polynomial factor in \( t \), so \( e^{\mathbf{A}t}\to\mathbf{0} \). If at least one eigenvalue has positive real part, the corresponding mode grows exponentially, so the origin is unstable.

For stability without asymptotic stability, all eigenvalues must satisfy \( \operatorname{Re}(\lambda_i)\le 0 \), and every eigenvalue satisfying \( \operatorname{Re}(\lambda_i)=0 \) must be semisimple, meaning its algebraic and geometric multiplicities agree.

\[ \boxed{ \begin{aligned} \text{Stable, not asymptotically stable} &\iff \operatorname{Re}(\lambda_i)\le 0 \text{ for all } i,\\ &\quad \text{and imaginary-axis eigenvalues are semisimple.} \end{aligned} } \]

For the discrete-time autonomous system \( \mathbf{x}[k+1]=\mathbf{A}\mathbf{x}[k] \), the stability boundary is the unit circle. The asymptotic stability condition is

\[ \boxed{ \mathbf{x}[k+1]=\mathbf{A}\mathbf{x}[k] \text{ is asymptotically stable} \iff |\lambda_i(\mathbf{A})|<1 \text{ for every } i. } \]

The proof follows from powers of a Jordan block:

\[ \mathbf{J}_\lambda^k = (\lambda\mathbf{I}+\mathbf{N})^k = \sum_{r=0}^{m-1} \binom{k}{r}\lambda^{k-r}\mathbf{N}^r. \]

If \( |\lambda|<1 \), the exponential decay of \( |\lambda|^k \) dominates the polynomial growth of \( \binom{k}{r} \). If \( |\lambda|>1 \), at least one mode grows. On the boundary \( |\lambda|=1 \), any nontrivial Jordan block causes polynomial growth, so boundary eigenvalues must be semisimple for Lyapunov stability.

\[ \boxed{ \begin{gathered} \text{Discrete-time stability without asymptotic stability} \\ \iff |\lambda_i|\le 1 \text{ for all } i \text{ and unit-circle eigenvalues are semisimple.} \end{gathered} } \]

Boundary eigenvalues are the most common source of confusion. Consider two continuous-time matrices with the same eigenvalue \( \lambda=0 \):

\[ \mathbf{A}_s= \begin{bmatrix} 0&0\\ 0&-1 \end{bmatrix}, \qquad \mathbf{A}_u= \begin{bmatrix} 0&1\\ 0&0 \end{bmatrix}. \]

The first matrix has a semisimple zero eigenvalue and a stable eigenvalue \( -1 \); its response is bounded but does not converge to the origin for every initial condition. The second matrix has a defective zero eigenvalue:

\[ e^{\mathbf{A}_u t} = \begin{bmatrix} 1&t\\ 0&1 \end{bmatrix}, \]

which grows linearly with \( t \). Therefore, the second system is unstable even though no eigenvalue has positive real part. The same phenomenon occurs in discrete time for defective eigenvalues on the unit circle.

flowchart TD
  A["Boundary eigenvalue"] --> B["Algebraic multiplicity equals geometric multiplicity?"]
  B -->|yes| C["Semisimple mode"]
  B -->|no| D["Defective Jordan block"]
  C --> E["Bounded oscillation or constant mode"]
  D --> F["Polynomial factor appears"]
  E --> G["Stable but not asymptotically stable"]
  F --> H["Unstable"]
        

If a continuous-time system \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \) is sampled with period \( T_s>0 \) and zero input, the corresponding discrete transition matrix is

\[ \mathbf{A}_d=e^{\mathbf{A}T_s}. \]

If \( \lambda_i \) is an eigenvalue of \( \mathbf{A} \), then \( e^{\lambda_i T_s} \) is an eigenvalue of \( \mathbf{A}_d \). Consequently,

\[ \operatorname{Re}(\lambda_i)<0 \quad \Longleftrightarrow \quad |e^{\lambda_i T_s}|<1. \]

This mapping explains why the open left half-plane in continuous time corresponds to the open unit disk in discrete time.

In exact mathematics, the criteria are strict. In numerical computation, small floating-point errors can move eigenvalues slightly across the boundary. Therefore, a practical routine uses a tolerance \( \varepsilon \). For continuous time, a robust engineering decision rule is

\[ \alpha(\mathbf{A})=\max_i \operatorname{Re}(\lambda_i(\mathbf{A})), \qquad \begin{cases} \alpha(\mathbf{A})<-\varepsilon & \text{asymptotically stable},\\ \alpha(\mathbf{A})>\varepsilon & \text{unstable},\\ |\alpha(\mathbf{A})|\le \varepsilon & \text{borderline; inspect Jordan structure}. \end{cases} \]

For discrete time, define the spectral radius \( \rho(\mathbf{A})=\max_i |\lambda_i(\mathbf{A})| \). Then:

\[ \begin{cases} \rho(\mathbf{A})<1-\varepsilon & \text{asymptotically stable},\\ \rho(\mathbf{A})>1+\varepsilon & \text{unstable},\\ |\rho(\mathbf{A})-1|\le \varepsilon & \text{borderline; inspect Jordan structure}. \end{cases} \]

In control software, the borderline case should not be ignored. A repeated eigenvalue at zero real part or on the unit circle can be harmless if semisimple, but unstable if defective.

Python workflows typically use numpy and scipy for eigenvalues, matrix exponentials, powers, and norms. The python-control package is also useful for state-space modeling, but the eigenvalue criterion itself only needs linear algebra.

Chapter9_Lesson2.py

# Chapter9_Lesson2.py
"""
Eigenvalue-based stability criteria for continuous-time and discrete-time LTI systems.

Required packages:
    numpy
    scipy

Optional package:
    python-control (not required here, but useful for state-space workflows)
"""

import numpy as np
from scipy.linalg import eigvals, expm, norm


def classify_continuous_lti(A: np.ndarray, tol: float = 1e-10) -> str:
    """
    Classify x_dot = A x using eigenvalue criteria.

    For a complete rigorous marginal-stability test, one must also verify
    semisimplicity of eigenvalues on the imaginary axis. This script estimates
    that by comparing algebraic multiplicity and geometric multiplicity.
    """
    lam = eigvals(A)
    max_real = np.max(np.real(lam))

    if max_real < -tol:
        return "asymptotically stable"

    if max_real > tol:
        return "unstable"

    if imaginary_axis_modes_are_semisimple(A, tol):
        return "stable but not asymptotically stable"

    return "unstable"


def classify_discrete_lti(A: np.ndarray, tol: float = 1e-10) -> str:
    """
    Classify x[k+1] = A x[k] using eigenvalue criteria.
    """
    lam = eigvals(A)
    rho = np.max(np.abs(lam))

    if rho < 1.0 - tol:
        return "asymptotically stable"

    if rho > 1.0 + tol:
        return "unstable"

    if unit_circle_modes_are_semisimple(A, tol):
        return "stable but not asymptotically stable"

    return "unstable"


def algebraic_multiplicity(values: np.ndarray, target: complex, tol: float) -> int:
    return int(np.sum(np.abs(values - target) < 50.0 * tol))


def geometric_multiplicity(A: np.ndarray, target: complex, tol: float) -> int:
    M = A.astype(complex) - target * np.eye(A.shape[0], dtype=complex)
    singular_values = np.linalg.svd(M, compute_uv=False)
    rank = int(np.sum(singular_values > tol))
    return A.shape[0] - rank


def unique_eigenvalues(values: np.ndarray, tol: float) -> list[complex]:
    unique = []
    for value in values:
        if not any(abs(value - old) < 50.0 * tol for old in unique):
            unique.append(value)
    return unique


def imaginary_axis_modes_are_semisimple(A: np.ndarray, tol: float) -> bool:
    lam = eigvals(A)
    for value in unique_eigenvalues(lam, tol):
        if abs(np.real(value)) <= tol:
            am = algebraic_multiplicity(lam, value, tol)
            gm = geometric_multiplicity(A, value, tol)
            if gm < am:
                return False
    return True


def unit_circle_modes_are_semisimple(A: np.ndarray, tol: float) -> bool:
    lam = eigvals(A)
    for value in unique_eigenvalues(lam, tol):
        if abs(abs(value) - 1.0) <= tol:
            am = algebraic_multiplicity(lam, value, tol)
            gm = geometric_multiplicity(A, value, tol)
            if gm < am:
                return False
    return True


def print_report(name: str, A: np.ndarray, system_type: str) -> None:
    lam = eigvals(A)
    print("=" * 72)
    print(name)
    print("A =")
    print(A)
    print("eigenvalues =", np.array2string(lam, precision=6))
    if system_type == "continuous":
        print("classification:", classify_continuous_lti(A))
        sample_times = [0, 1, 2, 5, 10]
        print("||exp(A t)||_2 samples:")
        for t in sample_times:
            print(f"  t={t:>2}: {norm(expm(A * t), 2):.6f}")
    elif system_type == "discrete":
        print("classification:", classify_discrete_lti(A))
        sample_steps = [0, 1, 2, 5, 10]
        print("||A^k||_2 samples:")
        for k in sample_steps:
            print(f"  k={k:>2}: {norm(np.linalg.matrix_power(A, k), 2):.6f}")
    else:
        raise ValueError("system_type must be 'continuous' or 'discrete'")


if __name__ == "__main__":
    # Continuous-time asymptotically stable: eigenvalues -1 and -2
    A1 = np.array([[-1.0, 0.0],
                   [ 0.0,-2.0]])

    # Continuous-time stable but not asymptotically stable:
    # pure rotation, eigenvalues +j and -j, no growing Jordan block.
    A2 = np.array([[0.0, -1.0],
                   [1.0,  0.0]])

    # Continuous-time unstable because the Jordan block at lambda=0 grows linearly.
    A3 = np.array([[0.0, 1.0],
                   [0.0, 0.0]])

    # Discrete-time asymptotically stable: spectral radius 0.8
    Ad1 = np.array([[0.8, 0.0],
                    [0.0, 0.5]])

    # Discrete-time unstable: eigenvalue magnitude larger than 1
    Ad2 = np.array([[1.02, 0.0],
                    [0.0, 0.7]])

    print_report("Continuous example 1", A1, "continuous")
    print_report("Continuous example 2", A2, "continuous")
    print_report("Continuous example 3", A3, "continuous")
    print_report("Discrete example 1", Ad1, "discrete")
    print_report("Discrete example 2", Ad2, "discrete")

In C++, the header-only Eigen library provides dense matrices and eigenvalue solvers suitable for state-space analysis. The example below classifies systems using eigenvalues and prints the modal spectrum.

Chapter9_Lesson2.cpp

// Chapter9_Lesson2.cpp
// Eigenvalue-based stability criteria using the Eigen C++ library.
//
// Build example:
//   g++ -std=c++17 Chapter9_Lesson2.cpp -I /path/to/eigen -O2 -o Chapter9_Lesson2
//
// Eigen is header-only and is widely used for state-space computations.

#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <complex>
#include <iostream>
#include <string>
#include <vector>
#include <cmath>

using Matrix = Eigen::MatrixXd;
using Complex = std::complex<double>;

std::string classifyContinuous(const Matrix& A, double tol = 1e-10) {
    Eigen::EigenSolver<Matrix> solver(A);
    auto eigenvalues = solver.eigenvalues();

    double maxReal = -1.0e300;
    for (int i = 0; i < eigenvalues.size(); ++i) {
        maxReal = std::max(maxReal, eigenvalues[i].real());
    }

    if (maxReal < -tol) return "asymptotically stable";
    if (maxReal > tol) return "unstable";

    return "borderline: check semisimplicity of imaginary-axis eigenvalues";
}

std::string classifyDiscrete(const Matrix& A, double tol = 1e-10) {
    Eigen::EigenSolver<Matrix> solver(A);
    auto eigenvalues = solver.eigenvalues();

    double spectralRadius = 0.0;
    for (int i = 0; i < eigenvalues.size(); ++i) {
        spectralRadius = std::max(spectralRadius, std::abs(eigenvalues[i]));
    }

    if (spectralRadius < 1.0 - tol) return "asymptotically stable";
    if (spectralRadius > 1.0 + tol) return "unstable";

    return "borderline: check semisimplicity of unit-circle eigenvalues";
}

void printEigenReport(const std::string& name, const Matrix& A, bool continuous) {
    Eigen::EigenSolver<Matrix> solver(A);

    std::cout << "============================================================\n";
    std::cout << name << "\nA =\n" << A << "\n";
    std::cout << "Eigenvalues:\n";
    for (int i = 0; i < solver.eigenvalues().size(); ++i) {
        std::cout << "  " << solver.eigenvalues()[i] << "\n";
    }

    if (continuous) {
        std::cout << "Continuous-time classification: "
                  << classifyContinuous(A) << "\n";
    } else {
        std::cout << "Discrete-time classification: "
                  << classifyDiscrete(A) << "\n";
    }
}

int main() {
    Matrix A1(2, 2);
    A1 << -1.0, 0.0,
           0.0, -2.0;

    Matrix A2(2, 2);
    A2 << 0.0, -1.0,
           1.0,  0.0;

    Matrix A3(2, 2);
    A3 << 0.0, 1.0,
           0.0, 0.0;

    Matrix Ad1(2, 2);
    Ad1 << 0.8, 0.0,
            0.0, 0.5;

    printEigenReport("Continuous asymptotically stable example", A1, true);
    printEigenReport("Continuous oscillatory borderline example", A2, true);
    printEigenReport("Continuous defective zero-eigenvalue example", A3, true);
    printEigenReport("Discrete asymptotically stable example", Ad1, false);

    return 0;
}

In Java, Apache Commons Math supplies matrix and eigenvalue routines. For large-scale control problems, more specialized numerical libraries may be preferable, but the following implementation is enough for small- and medium-size educational examples.

Chapter9_Lesson2.java

// Chapter9_Lesson2.java
// Eigenvalue-based stability criteria using Apache Commons Math.
//
// Compile/run example:
//   javac -cp commons-math3-3.6.1.jar Chapter9_Lesson2.java
//   java  -cp .;commons-math3-3.6.1.jar Chapter9_Lesson2
//
// On Linux/macOS use ":" instead of ";" in the classpath.

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.EigenDecomposition;
import org.apache.commons.math3.linear.RealMatrix;

public class Chapter9_Lesson2 {
    static final double TOL = 1e-10;

    public static String classifyContinuous(double[][] data) {
        RealMatrix A = new Array2DRowRealMatrix(data);
        EigenDecomposition eig = new EigenDecomposition(A);

        double maxReal = -Double.MAX_VALUE;
        for (int i = 0; i < A.getRowDimension(); i++) {
            maxReal = Math.max(maxReal, eig.getRealEigenvalue(i));
        }

        if (maxReal < -TOL) return "asymptotically stable";
        if (maxReal > TOL) return "unstable";
        return "borderline: check semisimplicity of imaginary-axis eigenvalues";
    }

    public static String classifyDiscrete(double[][] data) {
        RealMatrix A = new Array2DRowRealMatrix(data);
        EigenDecomposition eig = new EigenDecomposition(A);

        double spectralRadius = 0.0;
        for (int i = 0; i < A.getRowDimension(); i++) {
            double real = eig.getRealEigenvalue(i);
            double imag = eig.getImagEigenvalue(i);
            double magnitude = Math.hypot(real, imag);
            spectralRadius = Math.max(spectralRadius, magnitude);
        }

        if (spectralRadius < 1.0 - TOL) return "asymptotically stable";
        if (spectralRadius > 1.0 + TOL) return "unstable";
        return "borderline: check semisimplicity of unit-circle eigenvalues";
    }

    public static void printReport(String name, double[][] data, boolean continuous) {
        RealMatrix A = new Array2DRowRealMatrix(data);
        EigenDecomposition eig = new EigenDecomposition(A);

        System.out.println("====================================================");
        System.out.println(name);
        System.out.println("Eigenvalues:");
        for (int i = 0; i < A.getRowDimension(); i++) {
            System.out.printf("  %.8f %+ .8fi%n",
                    eig.getRealEigenvalue(i), eig.getImagEigenvalue(i));
        }

        if (continuous) {
            System.out.println("Classification: " + classifyContinuous(data));
        } else {
            System.out.println("Classification: " + classifyDiscrete(data));
        }
    }

    public static void main(String[] args) {
        double[][] A1 = {
            {-1.0, 0.0},
            { 0.0,-2.0}
        };

        double[][] A2 = {
            {0.0, -1.0},
            {1.0,  0.0}
        };

        double[][] A3 = {
            {0.0, 1.0},
            {0.0, 0.0}
        };

        double[][] Ad1 = {
            {0.8, 0.0},
            {0.0, 0.5}
        };

        printReport("Continuous asymptotically stable example", A1, true);
        printReport("Continuous oscillatory borderline example", A2, true);
        printReport("Continuous defective zero-eigenvalue example", A3, true);
        printReport("Discrete asymptotically stable example", Ad1, false);
    }
}

MATLAB is standard in control engineering for eigenvalue analysis. The commands eig, expm, ss, and isstable form a minimal state-space stability workflow. In Simulink, the same matrices can be placed in a State-Space block, and the response can be inspected by applying an initial condition.

Chapter9_Lesson2.m

% Chapter9_Lesson2.m
% Eigenvalue-based stability criteria for LTI systems in MATLAB/Octave.
% Related MATLAB functions: eig, jordan (Symbolic Math Toolbox), expm, ss, isstable.

clear; clc;

A1 = [-1 0; 0 -2];
A2 = [0 -1; 1 0];
A3 = [0 1; 0 0];

Ad1 = [0.8 0; 0 0.5];
Ad2 = [1.02 0; 0 0.7];

reportContinuous("Continuous asymptotically stable example", A1);
reportContinuous("Continuous oscillatory borderline example", A2);
reportContinuous("Continuous defective zero-eigenvalue example", A3);
reportDiscrete("Discrete asymptotically stable example", Ad1);
reportDiscrete("Discrete unstable example", Ad2);

% If Control System Toolbox is available:
% sys = ss(A1, eye(2), eye(2), zeros(2));
% isstable(sys)

function reportContinuous(name, A)
    fprintf("\n============================================================\n");
    fprintf("%s\n", name);
    disp("A ="); disp(A);
    lambda = eig(A);
    disp("eigenvalues ="); disp(lambda);
    fprintf("classification: %s\n", classifyContinuous(A));

    sampleTimes = [0 1 2 5 10];
    fprintf("norm(expm(A*t),2) samples:\n");
    for t = sampleTimes
        fprintf("  t = %2g: %.8f\n", t, norm(expm(A*t), 2));
    end
end

function reportDiscrete(name, A)
    fprintf("\n============================================================\n");
    fprintf("%s\n", name);
    disp("A ="); disp(A);
    lambda = eig(A);
    disp("eigenvalues ="); disp(lambda);
    fprintf("classification: %s\n", classifyDiscrete(A));

    sampleSteps = [0 1 2 5 10];
    fprintf("norm(A^k,2) samples:\n");
    for k = sampleSteps
        fprintf("  k = %2g: %.8f\n", k, norm(A^k, 2));
    end
end

function label = classifyContinuous(A)
    tol = 1e-10;
    lambda = eig(A);
    maxReal = max(real(lambda));

    if maxReal < -tol
        label = "asymptotically stable";
    elseif maxReal > tol
        label = "unstable";
    else
        label = "borderline: check semisimplicity of imaginary-axis eigenvalues";
    end
end

function label = classifyDiscrete(A)
    tol = 1e-10;
    lambda = eig(A);
    rho = max(abs(lambda));

    if rho < 1 - tol
        label = "asymptotically stable";
    elseif rho > 1 + tol
        label = "unstable";
    else
        label = "borderline: check semisimplicity of unit-circle eigenvalues";
    end
end

Simulink workflow: add a State-Space block, set \( \mathbf{A} \), \( \mathbf{B} \), \( \mathbf{C} \), and \( \mathbf{D} \), then use an initial condition block or model configuration to excite the free response. If all eigenvalues of \( \mathbf{A} \) are in the open left half-plane, the displayed state/output response decays to zero.

Wolfram Mathematica is useful for symbolic eigenvalue analysis, exact Jordan forms, and visualizing spectra in the complex plane.

Chapter9_Lesson2.nb

(* Chapter9_Lesson2.nb *)
(* Wolfram Mathematica code for eigenvalue-based stability criteria. *)

ClearAll[classifyContinuous, classifyDiscrete, reportContinuous, reportDiscrete];

classifyContinuous[A_] := Module[
  {lambda = Eigenvalues[A], maxReal},
  maxReal = Max[Re[lambda]];
  Which[
    maxReal < 0, "asymptotically stable",
    maxReal > 0, "unstable",
    True, "borderline: check semisimplicity of imaginary-axis eigenvalues"
  ]
];

classifyDiscrete[A_] := Module[
  {lambda = Eigenvalues[A], rho},
  rho = Max[Abs[lambda]];
  Which[
    rho < 1, "asymptotically stable",
    rho > 1, "unstable",
    True, "borderline: check semisimplicity of unit-circle eigenvalues"
  ]
];

reportContinuous[name_, A_] := Module[
  {lambda = Eigenvalues[A], times = {0, 1, 2, 5, 10} },
  Print["============================================================"];
  Print[name];
  Print["A = ", MatrixForm[A]];
  Print["eigenvalues = ", lambda];
  Print["classification = ", classifyContinuous[A]];
  Print["Norm[MatrixExp[A t], 2] samples:"];
  Table[Print["t = ", t, ": ", N[Norm[MatrixExp[A t], 2]]], {t, times}]
];

reportDiscrete[name_, A_] := Module[
  {lambda = Eigenvalues[A], steps = {0, 1, 2, 5, 10} },
  Print["============================================================"];
  Print[name];
  Print["A = ", MatrixForm[A]];
  Print["eigenvalues = ", lambda];
  Print["classification = ", classifyDiscrete[A]];
  Print["Norm[MatrixPower[A, k], 2] samples:"];
  Table[Print["k = ", k, ": ", N[Norm[MatrixPower[A, k], 2]]], {k, steps}]
];

A1 = { {-1, 0}, {0, -2} };
A2 = { {0, -1}, {1, 0} };
A3 = { {0, 1}, {0, 0} };
Ad1 = { {0.8, 0}, {0, 0.5} };
Ad2 = { {1.02, 0}, {0, 0.7} };

reportContinuous["Continuous asymptotically stable example", A1];
reportContinuous["Continuous oscillatory borderline example", A2];
reportContinuous["Continuous defective zero-eigenvalue example", A3];
reportDiscrete["Discrete asymptotically stable example", Ad1];
reportDiscrete["Discrete unstable example", Ad2];

(* Optional visualization of continuous-time eigenvalues in the complex plane. *)
ListPlot[
  ReIm /@ Eigenvalues[A2],
  AxesLabel -> {"Re(lambda)", "Im(lambda)"},
  PlotLabel -> "Eigenvalues of the rotation matrix"
]

Problem 1 (Continuous-Time Classification): Classify the stability of \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \) for

\[ \mathbf{A}= \begin{bmatrix} -2&0&0\\ 0&-1&4\\ 0&-4&-1 \end{bmatrix}. \]

Solution: The eigenvalues are \( -2 \) and \( -1\pm 4j \). Their real parts are \( -2 \) and \( -1 \), so every eigenvalue lies in the open left half-plane. Therefore, the equilibrium is asymptotically stable.

Problem 2 (Imaginary-Axis Eigenvalues): Determine whether the system with

\[ \mathbf{A}= \begin{bmatrix} 0&-3\\ 3&0 \end{bmatrix} \]

is stable, asymptotically stable, or unstable.

Solution: The eigenvalues are \( \lambda_1=3j \) and \( \lambda_2=-3j \). They are on the imaginary axis and are distinct, so they are semisimple. The state transition matrix is a rotation, hence bounded for all \( t\ge0 \). Therefore, the origin is stable but not asymptotically stable.

Problem 3 (Defective Boundary Eigenvalue): Let

\[ \mathbf{A}= \begin{bmatrix} 0&1\\ 0&0 \end{bmatrix}. \]

The only eigenvalue is \( \lambda=0 \). Is the system stable?

Solution: The zero eigenvalue has algebraic multiplicity two but geometric multiplicity one. Therefore, it is not semisimple. Since

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

a bounded initial condition can produce a state component growing linearly with \( t \). The system is unstable.

Problem 4 (Discrete-Time Classification): Classify \( \mathbf{x}[k+1]=\mathbf{A}\mathbf{x}[k] \) for

\[ \mathbf{A}= \begin{bmatrix} 0.6&0\\ 0&-0.8 \end{bmatrix}. \]

Solution: The eigenvalues are \( 0.6 \) and \( -0.8 \). Their magnitudes are strictly less than one: \( |0.6|<1 \) and \( |-0.8|<1 \). Hence \( \mathbf{A}^k\to\mathbf{0} \), and the discrete-time system is asymptotically stable.

Problem 5 (Sampling and Stability): Suppose a continuous-time system has eigenvalues \( -2 \) and \( -1\pm 5j \). If sampled with period \( T_s=0.1 \), show that the discrete-time eigenvalues are inside the unit circle.

Solution: Discrete eigenvalues are \( z_i=e^{\lambda_iT_s} \). Their magnitudes are

\[ |e^{\lambda_iT_s}| = e^{\operatorname{Re}(\lambda_i)T_s}. \]

For \( \lambda=-2 \), \( |z|=e^{-0.2}<1 \). For \( \lambda=-1\pm5j \), \( |z|=e^{-0.1}<1 \). Therefore, the sampled system is asymptotically stable.

Eigenvalues provide a complete stability test for autonomous LTI systems. In continuous time, asymptotic stability requires all eigenvalues of \( \mathbf{A} \) to lie in the open left half-plane. In discrete time, asymptotic stability requires all eigenvalues to lie inside the unit circle. Boundary eigenvalues require Jordan-structure inspection: semisimple boundary modes are stable but not asymptotically stable, while defective boundary modes are unstable.