In classical control, stability is often introduced through closed-loop poles or through the denominator of a transfer function. In modern control, the more primitive object is the internal state trajectory \( \mathbf{x}(t) \in \mathbb{R}^n \). For the autonomous continuous-time linear time-invariant system

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

the equilibrium at the origin is the reference object. The system is stable if initial states that start sufficiently close to the origin remain close for all future time; it is asymptotically stable if they also converge to the origin; it is unstable if this closeness property fails.

The solution learned in earlier lessons is expressed through the state transition matrix:

\[ \mathbf{x}(t)=\boldsymbol{\Phi}(t)\mathbf{x}_0 =e^{\mathbf{A}t}\mathbf{x}_0. \]

Therefore, stability is fundamentally a statement about the boundedness and limiting behavior of the linear operator \( e^{\mathbf{A}t} \) as \( t\ge 0 \).

flowchart TD
  X0["Initial condition x0 near origin"] --> FLOW["State flow x(t)=Phi(t)x0"]
  FLOW --> Q1["Does x(t) remain near origin for all t >= 0?"]
  Q1 -->|no| U["Unstable"]
  Q1 -->|yes| Q2["Does x(t) approach zero as t grows?"]
  Q2 -->|yes| AS["Asymptotically stable"]
  Q2 -->|no| S["Stable but not asymptotically stable"]
        

Let \( \|\cdot\| \) be any norm on \( \mathbb{R}^n \). In finite-dimensional vector spaces, all norms are equivalent, so the qualitative stability classification does not depend on the particular norm used.

Stability in the sense of Lyapunov. The origin is stable if

\[ \forall \varepsilon > 0,\; \exists \delta > 0 \text{ such that } \|\mathbf{x}_0\| < \delta \Rightarrow \|\mathbf{x}(t;\mathbf{x}_0)\| < \varepsilon \quad \forall t\ge 0. \]

Asymptotic stability. The origin is asymptotically stable if it is stable and, additionally,

\[ \exists r>0 \text{ such that } \|\mathbf{x}_0\|<r \Rightarrow \lim_{t\to\infty}\mathbf{x}(t;\mathbf{x}_0)=\mathbf{0}. \]

For a linear time-invariant system, local asymptotic stability and global asymptotic stability are equivalent. The reason is homogeneity: if \( \mathbf{x}(t;\mathbf{x}_0)=e^{\mathbf{A}t}\mathbf{x}_0 \) converges for all sufficiently small initial states, scaling \( \mathbf{x}_0 \) scales the entire trajectory linearly.

Instability. The origin is unstable if it is not stable:

\[ \exists \varepsilon_0>0 \text{ such that for every } \delta>0, \exists \mathbf{x}_0 \text{ with } \|\mathbf{x}_0\|<\delta \text{ and } \exists t_1\ge0 \text{ such that } \|\mathbf{x}(t_1;\mathbf{x}_0)\|\ge \varepsilon_0. \]

For \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \), the entire stability question can be expressed through \( e^{\mathbf{A}t} \).

Theorem 1 (bounded transition matrix and stability). The origin is stable if and only if there exists a finite constant \( M \) such that

\[ \|e^{\mathbf{A}t}\|\le M,\qquad \forall t\ge0. \]

Proof. If the above bound holds, then

\[ \|\mathbf{x}(t)\| = \|e^{\mathbf{A}t}\mathbf{x}_0\| \le \|e^{\mathbf{A}t}\|\,\|\mathbf{x}_0\| \le M\|\mathbf{x}_0\|. \]

Given \( \varepsilon>0 \), choose \( \delta=\varepsilon/M \). Then \( \|\mathbf{x}_0\|<\delta \) implies \( \|\mathbf{x}(t)\|<\varepsilon \) for all \( t\ge0 \). Conversely, if the origin is stable, the family of linear maps \( e^{\mathbf{A}t} \) is pointwise bounded on every sufficiently small ball; finite-dimensional linearity then implies uniform boundedness over directions, so the induced operator norm remains bounded for all future time.

Theorem 2 (vanishing transition matrix and asymptotic stability). The origin is asymptotically stable if and only if

\[ \lim_{t\to\infty} e^{\mathbf{A}t} = \mathbf{0}. \]

Proof. If \( e^{\mathbf{A}t}\to\mathbf{0} \), then \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0\to\mathbf{0} \) for every initial condition, and boundedness of the transition matrix gives stability. Conversely, if the origin is asymptotically stable, then every basis vector \( \mathbf{e}_i \) satisfies \( e^{\mathbf{A}t}\mathbf{e}_i\to\mathbf{0} \). Since the columns of \( e^{\mathbf{A}t} \) are exactly these images, the full matrix converges to the zero matrix.

Stability is not the same as convergence. A frictionless oscillator has bounded trajectories but does not decay to the origin. For example,

\[ \mathbf{A}= \begin{bmatrix} 0 & 1\\ -1 & 0 \end{bmatrix}, \qquad e^{\mathbf{A}t}= \begin{bmatrix} \cos t & \sin t\\ -\sin t & \cos t \end{bmatrix}. \]

This transition matrix is a rotation. Its norm remains bounded, but it does not approach the zero matrix. Hence the origin is stable but not asymptotically stable.

By contrast, if a component grows like \( e^{0.2t} \), then arbitrarily small initial conditions along that direction eventually leave any fixed neighborhood of the origin. That is instability. If all components decay, the origin is asymptotically stable.

\[ \begin{array}{c|c|c} \text{State behavior} & \text{Bounded?} & \text{Converges to zero?} \\ \hline \|\mathbf{x}(t)\|\le C\|\mathbf{x}_0\| & \text{yes} & \text{not necessarily} \\ \|\mathbf{x}(t)\|\le C e^{-\alpha t}\|\mathbf{x}_0\|,\; \alpha>0 & \text{yes} & \text{yes} \\ \|\mathbf{x}(t_k)\|\to\infty \text{ for some sequence } t_k & \text{no} & \text{no} \end{array} \]

In many engineering analyses, asymptotic stability is strengthened to exponential stability. The origin is exponentially stable if there exist constants \( M\ge1 \) and \( \alpha>0 \) such that

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

For finite-dimensional continuous-time LTI systems, asymptotic stability and exponential stability are equivalent. This fact is one reason LTI stability theory is so powerful: convergence is not merely eventual, but can be bounded by a decaying exponential envelope.

The value of \( \alpha \) is related to the slowest decaying mode, while \( M \) captures coordinate scaling, non-normal transient growth, and the geometry of eigenvectors. Even if a system is asymptotically stable, a non-normal matrix can temporarily amplify \( \|\mathbf{x}(t)\| \) before eventual decay.

\[ \sup_{t\ge0} \|e^{\mathbf{A}t}\| \text{ can be large even when } \lim_{t\to\infty}e^{\mathbf{A}t}=\mathbf{0}. \]

A useful state-space viewpoint is to measure the size of the state by a positive definite quadratic function \( V(\mathbf{x})=\mathbf{x}^T\mathbf{P}\mathbf{x} \), where \( \mathbf{P}=\mathbf{P}^T>0 \). Along trajectories,

\[ \dot{V} = \dot{\mathbf{x} }^T\mathbf{P}\mathbf{x} + \mathbf{x}^T\mathbf{P}\dot{\mathbf{x} } = \mathbf{x}^T(\mathbf{A}^T\mathbf{P}+\mathbf{P}\mathbf{A})\mathbf{x}. \]

If one can find \( \mathbf{P}=\mathbf{P}^T>0 \) and \( \mathbf{Q}=\mathbf{Q}^T>0 \) satisfying

\[ \mathbf{A}^T\mathbf{P}+\mathbf{P}\mathbf{A}=-\mathbf{Q}, \]

then \( \dot{V}=-\mathbf{x}^T\mathbf{Q}\mathbf{x}<0 \) for all nonzero states, proving asymptotic stability. The full Lyapunov theorem will be used more heavily in later design chapters; here it gives an energy interpretation of what it means for every state direction to decay.

Computationally, one usually combines state trajectories, eigenvalue inspection, transition matrix behavior, and Lyapunov equation checks. Lesson 2 develops eigenvalue criteria rigorously; this lesson uses them as diagnostic evidence while emphasizing the state-space definitions.

flowchart TD
  A["Given state matrix A"] --> B["Compute sample \ntrajectories"]
  A --> C["Compute eigenvalue \nreal parts"]
  A --> D["Inspect exp(A t) or \nsimulate Phi(t)"]
  C --> E["All real parts negative?"]
  E -->|yes| F["Expect decay to origin"]
  E -->|no| G["Positive real part or borderline case"]
  G --> H["Check boundedness and Jordan structure"]
  F --> I["Optional: solve Lyapunov equation"]
  H --> J["Classify as stable, marginal, or unstable"]
  I --> J
        

Python implementations typically use NumPy for matrices, SciPy for expm, eigenvalues, and Lyapunov equations, and python-control for higher-level state-space workflows.

Chapter9_Lesson1.py

# Chapter9_Lesson1.py
"""
Modern Control — Chapter 9, Lesson 1
Stability, Asymptotic Stability, and Instability (State-Space View)

This script studies continuous-time LTI systems:
    x_dot = A x
through the state transition matrix exp(A t), trajectories, eigenvalue
abscissa, and a quadratic Lyapunov equation check.
"""

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


def classify_ct_lti(A: np.ndarray, tol: float = 1e-9) -> str:
    """
    Basic continuous-time LTI classification from eigenvalues.

    For a complete rigorous marginal-stability decision, imaginary-axis
    eigenvalues must also be semisimple. This function flags the borderline
    case as "marginal-candidate" because Lesson 2 develops the full test.
    """
    lam = eigvals(A)
    alpha = np.max(np.real(lam))

    if alpha < -tol:
        return "asymptotically stable"
    if alpha > tol:
        return "unstable"

    return "marginal-candidate: inspect Jordan structure for imaginary-axis modes"


def exact_trajectory(A: np.ndarray, x0: np.ndarray, t_grid: np.ndarray) -> np.ndarray:
    """Compute x(t)=exp(A t)x0 on a grid."""
    return np.vstack([expm(A * t) @ x0 for t in t_grid])


def lyapunov_certificate(A: np.ndarray):
    """
    Solve A^T P + P A = -I.
    If P is symmetric positive definite, then V=x^T P x proves asymptotic stability.
    """
    Q = np.eye(A.shape[0])
    P = solve_continuous_lyapunov(A.T, -Q)
    eig_P = np.linalg.eigvalsh((P + P.T) / 2.0)
    return P, eig_P


def report_system(name: str, A: np.ndarray, x0: np.ndarray) -> None:
    print(f"\n=== {name} ===")
    print("A =")
    print(A)

    lam = eigvals(A)
    print("eigenvalues =", lam)
    print("classification =", classify_ct_lti(A))

    t_grid = np.linspace(0, 10, 101)
    X = exact_trajectory(A, x0, t_grid)
    norms = np.linalg.norm(X, axis=1)
    print("||x(0)|| =", norms[0])
    print("||x(10)|| =", norms[-1])
    print("max_t ||x(t)|| on [0,10] =", np.max(norms))

    if classify_ct_lti(A).startswith("asymptotically"):
        P, eig_P = lyapunov_certificate(A)
        print("Lyapunov P solving A^T P + P A = -I:")
        print(P)
        print("eigenvalues(P) =", eig_P)


if __name__ == "__main__":
    x0 = np.array([1.0, -0.5])

    A_stable = np.array([[-1.0, 2.0],
                         [-3.0, -2.0]])

    A_center = np.array([[0.0, 1.0],
                         [-1.0, 0.0]])

    A_unstable = np.array([[0.2, 1.0],
                           [0.0, -1.0]])

    report_system("Stable spiral", A_stable, x0)
    report_system("Marginal center", A_center, x0)
    report_system("Unstable saddle/source component", A_unstable, x0)

For large-scale C++ control software, Eigen, Armadillo, or LAPACK wrappers are common choices. The implementation below is self-contained for two-dimensional systems and uses trace-determinant classification plus a fourth-order Runge-Kutta trajectory simulation.

Chapter9_Lesson1.cpp

// Chapter9_Lesson1.cpp
// Modern Control — Chapter 9, Lesson 1
// Stability, Asymptotic Stability, and Instability (State-Space View)
//
// Self-contained 2-by-2 implementation.
// For larger systems, production C++ projects commonly use Eigen, Armadillo,
// or LAPACK wrappers for eigenvalues and matrix exponentials.

#include <cmath>
#include <iostream>
#include <string>
#include <algorithm>

struct Matrix2 {
    double a11, a12, a21, a22;

    double trace() const {
        return a11 + a22;
    }

    double det() const {
        return a11 * a22 - a12 * a21;
    }

    void multiply(const double x[2], double y[2]) const {
        y[0] = a11 * x[0] + a12 * x[1];
        y[1] = a21 * x[0] + a22 * x[1];
    }
};

double norm2(const double x[2]) {
    return std::sqrt(x[0] * x[0] + x[1] * x[1]);
}

std::string classify2x2(const Matrix2& A, double tol = 1e-10) {
    double tr = A.trace();
    double determinant = A.det();
    double discriminant = tr * tr - 4.0 * determinant;

    if (discriminant >= 0.0) {
        double root = std::sqrt(discriminant);
        double lambda1 = 0.5 * (tr + root);
        double lambda2 = 0.5 * (tr - root);

        if (lambda1 < -tol && lambda2 < -tol) return "asymptotically stable";
        if (lambda1 > tol || lambda2 > tol) return "unstable";
        return "borderline: inspect semisimplicity";
    }

    double realPart = 0.5 * tr;
    if (realPart < -tol) return "asymptotically stable";
    if (realPart > tol) return "unstable";
    return "marginal center candidate";
}

void addScaled(const double x[2], const double k[2], double scale, double y[2]) {
    y[0] = x[0] + scale * k[0];
    y[1] = x[1] + scale * k[1];
}

void rk4Step(const Matrix2& A, const double x[2], double h, double xNext[2]) {
    double k1[2], k2[2], k3[2], k4[2], temp[2];

    A.multiply(x, k1);
    addScaled(x, k1, 0.5 * h, temp);
    A.multiply(temp, k2);
    addScaled(x, k2, 0.5 * h, temp);
    A.multiply(temp, k3);
    addScaled(x, k3, h, temp);
    A.multiply(temp, k4);

    xNext[0] = x[0] + h * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0;
    xNext[1] = x[1] + h * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0;
}

void reportSystem(const std::string& name, const Matrix2& A, const double x0[2]) {
    std::cout << "
=== " << name << " ===
";
    std::cout << "A = [[" << A.a11 << ", " << A.a12 << "], ["
              << A.a21 << ", " << A.a22 << "]]
";
    std::cout << "trace = " << A.trace() << ", determinant = " << A.det() << "
";
    std::cout << "classification = " << classify2x2(A) << "
";

    double h = 0.01;
    int steps = 1000;
    double x[2] = {x0[0], x0[1]};
    double maxNorm = norm2(x);

    for (int k = 0; k < steps; ++k) {
        double xNext[2];
        rk4Step(A, x, h, xNext);
        x[0] = xNext[0];
        x[1] = xNext[1];
        maxNorm = std::max(maxNorm, norm2(x));
    }

    std::cout << "||x(0)|| = " << norm2(x0) << "
";
    std::cout << "||x(10)|| approx = " << norm2(x) << "
";
    std::cout << "max ||x(t)|| approx on [0,10] = " << maxNorm << "
";
}

int main() {
    double x0[2] = {1.0, -0.5};

    Matrix2 stable{-1.0, 2.0, -3.0, -2.0};
    Matrix2 center{0.0, 1.0, -1.0, 0.0};
    Matrix2 unstable{0.2, 1.0, 0.0, -1.0};

    reportSystem("Stable spiral", stable, x0);
    reportSystem("Marginal center", center, x0);
    reportSystem("Unstable saddle/source component", unstable, x0);

    return 0;
}

Java does not include advanced numerical linear algebra in the standard library. Production projects often use EJML, Apache Commons Math, or ojAlgo. The self-contained code below handles two-dimensional systems and demonstrates trace-determinant classification with RK4 simulation.

Chapter9_Lesson1.java

// Chapter9_Lesson1.java
// Modern Control — Chapter 9, Lesson 1
// Stability, Asymptotic Stability, and Instability (State-Space View)
//
// This self-contained example handles 2-by-2 continuous-time systems.
// It classifies by trace/determinant eigenvalue information and simulates
// x_dot = A x using RK4.

public class Chapter9_Lesson1 {

    static class Matrix2 {
        double a11, a12, a21, a22;

        Matrix2(double a11, double a12, double a21, double a22) {
            this.a11 = a11;
            this.a12 = a12;
            this.a21 = a21;
            this.a22 = a22;
        }

        double trace() {
            return a11 + a22;
        }

        double det() {
            return a11 * a22 - a12 * a21;
        }

        double[] multiply(double[] x) {
            return new double[] {
                a11 * x[0] + a12 * x[1],
                a21 * x[0] + a22 * x[1]
            };
        }
    }

    static double norm(double[] x) {
        return Math.sqrt(x[0] * x[0] + x[1] * x[1]);
    }

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

    static double[] rk4Step(Matrix2 A, double[] x, double h) {
        double[] k1 = A.multiply(x);
        double[] k2 = A.multiply(add(x, k1, 0.5 * h));
        double[] k3 = A.multiply(add(x, k2, 0.5 * h));
        double[] k4 = A.multiply(add(x, k3, h));

        return new double[] {
            x[0] + h * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
            x[1] + h * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0
        };
    }

    static String classify2x2(Matrix2 A, double tol) {
        double tr = A.trace();
        double det = A.det();
        double disc = tr * tr - 4.0 * det;

        if (disc >= 0.0) {
            double sqrtDisc = Math.sqrt(disc);
            double lambda1 = 0.5 * (tr + sqrtDisc);
            double lambda2 = 0.5 * (tr - sqrtDisc);

            if (lambda1 < -tol && lambda2 < -tol) return "asymptotically stable";
            if (lambda1 > tol || lambda2 > tol) return "unstable";
            return "borderline: inspect semisimplicity";
        } else {
            double realPart = 0.5 * tr;
            if (realPart < -tol) return "asymptotically stable";
            if (realPart > tol) return "unstable";
            return "marginal center candidate";
        }
    }

    static void reportSystem(String name, Matrix2 A, double[] x0) {
        System.out.println("\n=== " + name + " ===");
        System.out.printf("A = [[%.3f, %.3f], [%.3f, %.3f]]%n", A.a11, A.a12, A.a21, A.a22);
        System.out.printf("trace = %.6f, determinant = %.6f%n", A.trace(), A.det());
        System.out.println("classification = " + classify2x2(A, 1e-10));

        double h = 0.01;
        int steps = 1000;
        double[] x = new double[] {x0[0], x0[1]};
        double maxNorm = norm(x);

        for (int k = 0; k < steps; k++) {
            x = rk4Step(A, x, h);
            maxNorm = Math.max(maxNorm, norm(x));
        }

        System.out.printf("||x(0)|| = %.6f%n", norm(x0));
        System.out.printf("||x(10)|| approx = %.6f%n", norm(x));
        System.out.printf("max ||x(t)|| approx on [0,10] = %.6f%n", maxNorm);
    }

    public static void main(String[] args) {
        double[] x0 = {1.0, -0.5};

        Matrix2 stable = new Matrix2(-1.0, 2.0, -3.0, -2.0);
        Matrix2 center = new Matrix2(0.0, 1.0, -1.0, 0.0);
        Matrix2 unstable = new Matrix2(0.2, 1.0, 0.0, -1.0);

        reportSystem("Stable spiral", stable, x0);
        reportSystem("Marginal center", center, x0);
        reportSystem("Unstable saddle/source component", unstable, x0);
    }
}

MATLAB is especially direct for this lesson because eig, expm, lyap, and state-space/Simulink tools are native to the Control System Toolbox and Simulink ecosystem.

Chapter9_Lesson1.m

% Chapter9_Lesson1.m
% Modern Control — Chapter 9, Lesson 1
% Stability, Asymptotic Stability, and Instability (State-Space View)
%
% Demonstrates continuous-time LTI stability using eig(A), expm(A*t),
% Lyapunov equation checks, and a programmatic Simulink state-space model.

clear; clc; close all;

A_stable = [-1  2;
            -3 -2];

A_center = [ 0  1;
            -1  0];

A_unstable = [0.2  1;
              0   -1];

x0 = [1; -0.5];

reportSystem("Stable spiral", A_stable, x0);
reportSystem("Marginal center", A_center, x0);
reportSystem("Unstable saddle/source component", A_unstable, x0);

% Plot exact trajectories using the matrix exponential.
t = linspace(0, 10, 500);
Xstable = zeros(2, numel(t));
Xcenter = zeros(2, numel(t));
Xunstable = zeros(2, numel(t));

for k = 1:numel(t)
    Xstable(:, k) = expm(A_stable * t(k)) * x0;
    Xcenter(:, k) = expm(A_center * t(k)) * x0;
    Xunstable(:, k) = expm(A_unstable * t(k)) * x0;
end

figure;
plot(t, vecnorm(Xstable), "LineWidth", 1.5); hold on;
plot(t, vecnorm(Xcenter), "LineWidth", 1.5);
plot(t, vecnorm(Xunstable), "LineWidth", 1.5);
grid on;
xlabel("time t");
ylabel("state norm ||x(t)||");
legend("stable", "center", "unstable", "Location", "best");
title("State norm behavior for stability classes");

% Simulink-oriented construction of x_dot = A x with output y = x.
% This section requires Simulink. It creates a model with a State-Space block.
model = "Chapter9_Lesson1_StateSpace_Model";
if license("test", "Simulink")
    if bdIsLoaded(model)
        close_system(model, 0);
    end
    new_system(model);
    open_system(model);

    add_block("simulink/Continuous/State-Space", model + "/State-Space");
    set_param(model + "/State-Space", ...
        "A", mat2str(A_stable), ...
        "B", mat2str(zeros(2, 1)), ...
        "C", mat2str(eye(2)), ...
        "D", mat2str(zeros(2, 1)), ...
        "X0", mat2str(x0));

    add_block("simulink/Sources/Constant", model + "/ZeroInput");
    set_param(model + "/ZeroInput", "Value", "0");

    add_block("simulink/Sinks/Scope", model + "/Scope");

    add_line(model, "ZeroInput/1", "State-Space/1");
    add_line(model, "State-Space/1", "Scope/1");

    set_param(model, "StopTime", "10");
    save_system(model);
    disp("Simulink model created and saved: " + model + ".slx");
else
    disp("Simulink license not available; skipped model creation.");
end

function reportSystem(name, A, x0)
    fprintf("\n=== %s ===\n", name);
    disp("A ="); disp(A);

    lambda = eig(A);
    disp("eigenvalues ="); disp(lambda);

    alpha = max(real(lambda));
    tol = 1e-10;

    if alpha < -tol
        classification = "asymptotically stable";
    elseif alpha > tol
        classification = "unstable";
    else
        classification = "marginal-candidate: inspect Jordan structure";
    end

    disp("classification = " + classification);

    x10 = expm(A * 10) * x0;
    fprintf("||x(0)|| = %.6f\n", norm(x0));
    fprintf("||x(10)|| = %.6f\n", norm(x10));

    if classification == "asymptotically stable"
        P = lyap(A', eye(size(A)));
        disp("P solving A'P + P A = -I:");
        disp(P);
        disp("eig(P) =");
        disp(eig(P));
    end
end

Mathematica is useful for symbolic and exact matrix exponentials, as well as high-level visualization of state norm behavior.

Chapter9_Lesson1.nb

Notebook[{
    Cell["Chapter9_Lesson1.nb", "Title"],
    Cell["Modern Control — Chapter 9, Lesson 1: Stability, Asymptotic Stability, and Instability", "Subtitle"],
    Cell[BoxData[
        RowBox[{
            RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}], ";"}]], "Input"],
    Cell[BoxData[
        RowBox[{
            RowBox[{"AStable", "=", RowBox[{"{ {", RowBox[{
                RowBox[{"-", "1"}], ",", "2", "},", "{", RowBox[{
                RowBox[{"-", "3"}], ",", RowBox[{"-", "2"}]}], "} }"}]}], ";"}]], "Input"],
    Cell[BoxData[
        RowBox[{
            RowBox[{"ACenter", "=", RowBox[{"{ {", RowBox[{
                "0", ",", "1", "},", "{", RowBox[{
                RowBox[{"-", "1"}], ",", "0"}], "} }"}]}], ";"}]], "Input"],
    Cell[BoxData[
        RowBox[{
            RowBox[{"AUnstable", "=", RowBox[{"{ {", RowBox[{
                "0.2", ",", "1", "},", "{", RowBox[{
                "0", ",", RowBox[{"-", "1"}]}], "} }"}]}], ";"}]], "Input"],
    Cell[BoxData[
        RowBox[{
            RowBox[{"x0", "=", RowBox[{"{", RowBox[{"1", ",", RowBox[{"-", "0.5"}]}], "}"}]}], ";"}]], "Input"],
    Cell[BoxData[
        RowBox[{
            RowBox[{"classify", "[", "A_", "]"}], ":=", RowBox[{"Module", "[",
                RowBox[{
                    RowBox[{"{", RowBox[{"lam", ",", "alpha", ",", "tol"}], "}"}], ",",
                    RowBox[{
                        RowBox[{"tol", "=", RowBox[{"10", "^", RowBox[{"-", "10"}]}]}], ";",
                        RowBox[{"lam", "=", RowBox[{"Eigenvalues", "[", "A", "]"}]}], ";",
                        RowBox[{"alpha", "=", RowBox[{"Max", "[", RowBox[{"Re", "[", "lam", "]"}], "]"}]}], ";",
                        RowBox[{"Which", "[",
                            RowBox[{
                                RowBox[{"alpha", "<", RowBox[{"-", "tol"}]}], ",", "\"asymptotically stable\"", ",",
                                RowBox[{"alpha", ">", "tol"}], ",", "\"unstable\"", ",",
                                "True", ",", "\"marginal-candidate: inspect Jordan structure\""}], "]"}]}]}], "]"}]}]], "Input"],
    Cell[BoxData[
        RowBox[{
            RowBox[{"trajectory", "[", RowBox[{"A_", ",", "t_"}], "]"}], ":=",
            RowBox[{RowBox[{"MatrixExp", "[", RowBox[{"A", " ", "t"}], "]"}], ".", "x0"}]}]], "Input"],
    Cell[BoxData[
        RowBox[{"Table", "[",
            RowBox[{
                RowBox[{"{", RowBox[{"name", ",", RowBox[{"Eigenvalues", "[", "A", "]"}], ",", RowBox[{"classify", "[", "A", "]"}], ",", RowBox[{"Norm", "[", RowBox[{"trajectory", "[", RowBox[{"A", ",", "10"}], "]"}], "]"}]}], "}"}], ",",
                RowBox[{"{", RowBox[{"{", RowBox[{"name", ",", "A"}], "}", ",",
                    RowBox[{"{", RowBox[{
                        RowBox[{"{", RowBox[{"\"Stable spiral\"", ",", "AStable"}], "}"}], ",",
                        RowBox[{"{", RowBox[{"\"Marginal center\"", ",", "ACenter"}], "}"}], ",",
                        RowBox[{"{", RowBox[{"\"Unstable\"", ",", "AUnstable"}], "}"}]}], "}"}]}], "}"}]}], "]"}]], "Input"],
    Cell[BoxData[
        RowBox[{"Plot", "[",
            RowBox[{
                RowBox[{"Evaluate", "[", RowBox[{"{",
                    RowBox[{
                        RowBox[{"Norm", "[", RowBox[{"trajectory", "[", RowBox[{"AStable", ",", "t"}], "]"}], "]"}], ",",
                        RowBox[{"Norm", "[", RowBox[{"trajectory", "[", RowBox[{"ACenter", ",", "t"}], "]"}], "]"}], ",",
                        RowBox[{"Norm", "[", RowBox[{"trajectory", "[", RowBox[{"AUnstable", ",", "t"}], "]"}], "]"}]}],
                    "}"}], "]"}], ",",
                RowBox[{"{", RowBox[{"t", ",", "0", ",", "10"}], "}"}], ",",
                RowBox[{"PlotLegends", "->", RowBox[{"{", RowBox[{"\"stable\"", ",", "\"center\"", ",", "\"unstable\""}], "}"}]}], ",",
                RowBox[{"AxesLabel", "->", RowBox[{"{", RowBox[{"\"t\"", ",", "\"||x(t)||\""}], "}"}]}]}], "]"}]], "Input"]
}]

Problem 1 (Definition Check): Consider an autonomous system whose trajectories satisfy \( \|\mathbf{x}(t)\|\le 4\|\mathbf{x}_0\| \) for all \( t\ge0 \). Prove that the origin is stable.

Solution: Given \( \varepsilon>0 \), choose \( \delta=\varepsilon/4 \). If \( \|\mathbf{x}_0\|<\delta \), then \( \|\mathbf{x}(t)\|\le4\|\mathbf{x}_0\|<\varepsilon \) for all \( t\ge0 \). This is exactly Lyapunov stability.

Problem 2 (Rotation System): For \( \mathbf{A}=\begin{bmatrix}0&1\\-1&0\end{bmatrix} \), show that the origin is stable but not asymptotically stable.

Solution: The matrix exponential is

\[ e^{\mathbf{A}t}= \begin{bmatrix} \cos t & \sin t\\ -\sin t & \cos t \end{bmatrix}. \]

This matrix is orthogonal, so \( \|e^{\mathbf{A}t}\mathbf{x}_0\|_2=\|\mathbf{x}_0\|_2 \) for all \( t \). Hence the origin is stable. However, unless \( \mathbf{x}_0=\mathbf{0} \), the state rotates forever and does not converge to zero. Therefore, it is not asymptotically stable.

Problem 3 (Scalar State System): Classify the origin of \( \dot{x}=ax \) for \( a<0 \), \( a=0 \), and \( a>0 \).

Solution: The solution is \( x(t)=e^{at}x_0 \). If \( a<0 \), then \( e^{at} o0 \), so the origin is asymptotically stable. If \( a=0 \), then \( x(t)=x_0 \), so the origin is stable but not asymptotically stable. If \( a>0 \), then any nonzero initial condition grows exponentially, so the origin is unstable.

Problem 4 (Transition Matrix Bound): Suppose \( \|e^{\mathbf{A}t}\|\le Me^{-\alpha t} \) for constants \( M\ge1 \) and \( \alpha>0 \). Prove asymptotic stability.

Solution: The solution satisfies \( \|\mathbf{x}(t)\|\le Me^{-\alpha t}\|\mathbf{x}_0\| \). Stability follows by choosing \( \delta=\varepsilon/M \). Since \( e^{-\alpha t} o0 \), the right-hand side converges to zero, so \( \mathbf{x}(t) o\mathbf{0} \).

Problem 5 (Quadratic Energy Derivative): Let \( V(\mathbf{x})=\mathbf{x}^T\mathbf{P}\mathbf{x} \) with \( \mathbf{P}=\mathbf{P}^T>0 \). Derive \( \dot{V} \) for \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \).

Solution: Applying the product rule,

\[ \dot{V} = \dot{\mathbf{x} }^T\mathbf{P}\mathbf{x} + \mathbf{x}^T\mathbf{P}\dot{\mathbf{x} } = (\mathbf{A}\mathbf{x})^T\mathbf{P}\mathbf{x} + \mathbf{x}^T\mathbf{P}\mathbf{A}\mathbf{x} = \mathbf{x}^T(\mathbf{A}^T\mathbf{P}+\mathbf{P}\mathbf{A})\mathbf{x}. \]

Stability in state space is a statement about the internal trajectory, not merely a transfer-function denominator. For LTI systems, the state transition matrix \( e^{\mathbf{A}t} \) determines whether the origin is stable, asymptotically stable, or unstable. Stable systems keep nearby states nearby; asymptotically stable systems also force states to the origin; unstable systems allow arbitrarily small perturbations to leave a fixed neighborhood. The next lesson turns these definitions into practical eigenvalue-based criteria.