After completing this lesson, a student should be able to:

• state Lyapunov, asymptotic, global, and exponential stability precisely;
• verify positive definiteness, radial unboundedness, and sign definiteness of a derivative;
• prove the basic direct-method stability theorem using invariant sublevel sets;
• apply LaSalle's invariance principle when the derivative is only negative semidefinite;
• construct and interpret quadratic certificates for stable linear systems;
• distinguish pointwise stability from uniform stability under time-varying uncertainty; and
• recognize why a state-only candidate often fails for adaptive error dynamics.

In fixed-parameter control, the controller is designed first and stability is checked afterward. In Lyapunov-based adaptive control, the analysis and the update mechanism are usually coupled: one proposes a Lyapunov candidate, differentiates it along the uncertain closed-loop dynamics, and later chooses an adaptation law so that troublesome mixed terms are removed or dominated. This lesson reviews the certificate side of that procedure. The construction of a combined state–parameter candidate is intentionally reserved for Lesson 2.

flowchart TD
  A["Closed-loop state or tracking-error dynamics"] --> B["Choose a scalar candidate V"]
  B --> C["Check V is positive definite"]
  C --> D["Differentiate V along trajectories"]
  D --> E{"Sign of V_dot?"}
  E -->|"negative definite"| F["Conclude asymptotic or \nexponential stability"]
  E -->|"negative semidefinite"| G["Find the largest invariant \nzero-derivative set"]
  E -->|"indefinite"| H["Candidate or design \nis insufficient"]
  G --> I["Apply invariance reasoning"]
  H --> J["Redesign candidate or \nlater choose an \nadaptive update"]
        

Begin with the autonomous nonlinear system \( \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}) \), where \( \mathbf{x}\in\mathbb{R}^{n} \), \( \mathbf{f}(\mathbf{0})=\mathbf{0} \), and \( \mathbf{f} \) is locally Lipschitz. The local Lipschitz condition guarantees a unique local solution through each initial condition. The origin is therefore an equilibrium.

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

\[ \forall\,\varepsilon > 0\;\exists\,\delta(\varepsilon)>0:\quad \|\mathbf{x}(0)\|<\delta \;\Longrightarrow\; \|\mathbf{x}(t)\|<\varepsilon,\quad \forall t\ge 0. \]

Stability is a bounded-deviation statement. It does not by itself require the state to converge to the origin.

Attractivity. The origin is locally attractive if there is a radius \( r>0 \) such that

\[ \|\mathbf{x}(0)\|<r \;\Longrightarrow\; \mathbf{x}(t)→\mathbf{0}\quad\text{as}\quad t→\infty. \]

The equilibrium is locally asymptotically stable when it is both stable and locally attractive. It is globally asymptotically stable when stability holds and every solution in \( \mathbb{R}^{n} \) converges to the origin.

Exponential stability. The origin is locally exponentially stable if

\[ \|\mathbf{x}(t)\|\le M e^{-\lambda t}\|\mathbf{x}(0)\|, \qquad M\ge 1,\quad \lambda>0, \]

for all sufficiently small initial conditions. Exponential stability gives an explicit decay rate and is stronger than asymptotic stability.

A continuously differentiable scalar function \( V:\mathcal{D}\subset\mathbb{R}^{n}\rightarrow\mathbb{R} \) is positive definite about the origin if

\[ V(\mathbf{0})=0,\qquad V(\mathbf{x})>0 \quad\forall\mathbf{x}\ne\mathbf{0}. \]

It is negative definite when \( -V \) is positive definite. A function is positive semidefinite if \( V(\mathbf{x})\ge0 \) but may vanish away from the origin.

A useful comparison representation is

\[ \alpha_{1}(\|\mathbf{x}\|)\le V(\mathbf{x}) \le \alpha_{2}(\|\mathbf{x}\|), \]

where \( \alpha_{1} \) and \( \alpha_{2} \) are continuous, strictly increasing, and zero at the origin. Such bounds convert statements about the scalar \( V \) into statements about the state norm.

The sublevel set associated with \( c>0 \) is

\[ \Omega_{c}=\{\mathbf{x}\in\mathcal{D}:V(\mathbf{x})\le c\}. \]

If \( \dot V\le0 \) along trajectories, every sublevel set lying inside the model domain is positively invariant: once a trajectory enters it, the trajectory cannot leave through a boundary on which \( V=c \).

For global results, one commonly also requires radial unboundedness:

\[ \|\mathbf{x}\|→\infty \quad\Longrightarrow\quad V(\mathbf{x})→\infty. \]

Radial unboundedness makes all finite sublevel sets bounded and prevents the certificate from losing information far from the origin.

Theorem 1 (local stability). Suppose there is a continuously differentiable function \( V \) on a neighborhood \( \mathcal{D} \) of the origin such that

\[ V(\mathbf{0})=0,\qquad V(\mathbf{x})>0\;\text{for}\;\mathbf{x}\ne\mathbf{0}, \qquad \dot V(\mathbf{x})=\nabla V(\mathbf{x})^{T}\mathbf{f}(\mathbf{x})\le0. \]

Then the origin is stable.

Proof. Fix an arbitrary \( \varepsilon>0 \) such that the closed ball \( \overline B_{\varepsilon} \) is contained in \( \mathcal{D} \). On the compact sphere \( \|\mathbf{x}\|=\varepsilon \), continuity and positive definiteness imply that the minimum

\[ c_{\varepsilon}= \min_{\|\mathbf{x}\|=\varepsilon}V(\mathbf{x}) \]

exists and satisfies \( c_{\varepsilon}>0 \). Since \( V(\mathbf{0})=0 \) and \( V \) is continuous, there is a \( \delta>0 \) such that

\[ \|\mathbf{x}(0)\|<\delta \quad\Longrightarrow\quad V(\mathbf{x}(0))<c_{\varepsilon}. \]

Because \( \dot V\le0 \), the value of \( V \) cannot increase:

\[ V(\mathbf{x}(t))\le V(\mathbf{x}(0))<c_{\varepsilon}, \qquad t\ge0. \]

A trajectory cannot cross the sphere \( \|\mathbf{x}\|=\varepsilon \), because every point on that sphere has \( V\ge c_{\varepsilon} \). Therefore \( \|\mathbf{x}(t)\|<\varepsilon \) for all \( t\ge0 \). This is Lyapunov stability. □

The proof is important for adaptive control because it reveals the real mechanism: a nonincreasing scalar certificate traps the complete augmented trajectory inside a sublevel set. No closed-form solution is required.

Theorem 2 (asymptotic stability). Under the assumptions of Theorem 1, if

\[ \dot V(\mathbf{x})<0 \quad\forall\mathbf{x}\ne\mathbf{0}, \]

then the origin is locally asymptotically stable.

A proof can be obtained by considering a compact annulus \( \eta\le\|\mathbf{x}\|\le\varepsilon \). On that annulus, the continuous positive-definite function \( W(\mathbf{x})=-\dot V(\mathbf{x}) \) has a strictly positive minimum. Hence a trajectory cannot remain forever in the annulus, because \( V \) would decrease below zero in finite time. Since the inner radius \( \eta \) is arbitrary, the state approaches the origin.

A particularly useful exponential certificate is the existence of constants \( c_{1},c_{2},c_{3}>0 \) satisfying

\[ c_{1}\|\mathbf{x}\|^{2}\le V(\mathbf{x}) \le c_{2}\|\mathbf{x}\|^{2}, \qquad \dot V(\mathbf{x})\le-c_{3}\|\mathbf{x}\|^{2}. \]

Since \( \|\mathbf{x}\|^{2}\ge V/c_{2} \),

\[ \dot V\le-\frac{c_{3}}{c_{2}}V. \]

Multiplication by the integrating factor gives

\[ V(t)\le V(0)\exp\!\left(-\frac{c_{3}}{c_{2}}t\right). \]

Using the quadratic bounds once more yields

\[ \|\mathbf{x}(t)\| \le \sqrt{\frac{c_{2}}{c_{1}}} \exp\!\left(-\frac{c_{3}}{2c_{2}}t\right) \|\mathbf{x}(0)\|. \]

In adaptive systems, negative definiteness in the full augmented state is often unavailable. One may instead prove boundedness and convergence of a tracking error by combining semidefinite Lyapunov decay with an invariance or uniform-continuity argument. The latter argument is developed later with Barbalat's lemma.

Negative semidefiniteness does not imply that every state component converges to zero. LaSalle's invariance principle resolves this by examining the dynamics restricted to the zero-derivative set.

Theorem 3 (LaSalle). Let \( \Omega \) be compact and positively invariant. If \( V \) is continuously differentiable and

\[ \dot V(\mathbf{x})\le0\quad\text{on}\quad\Omega, \]

define

\[ E=\{\mathbf{x}\in\Omega:\dot V(\mathbf{x})=0\}. \]

Every trajectory in \( \Omega \) approaches the largest invariant subset \( M \) contained in \( E \). Therefore, if \( M=\{\mathbf{0}\} \), the origin is asymptotically stable within \( \Omega \).

Worked example: damped nonlinear oscillator.

\[ \dot q=v,\qquad \dot v=-kq-cv-\beta q^{3}, \qquad k>0,\;c>0,\;\beta>0. \]

Choose the total-energy candidate

\[ V(q,v)=\frac{1}{2}v^{2}+\frac{k}{2}q^{2}+\frac{\beta}{4}q^{4}. \]

Its derivative is

\[ \begin{aligned} \dot V &=v\dot v+kq\dot q+\beta q^{3}\dot q\\ &=v(-kq-cv-\beta q^{3})+kqv+\beta q^{3}v\\ &=-cv^{2}\le0. \end{aligned} \]

The zero-derivative set is \( E=\{(q,v):v=0\} \). A trajectory can remain in \( E \) only if \( \dot v=-kq-\beta q^{3}=0 \), which implies \( q=0 \). Thus the largest invariant subset is the origin. Because \( V \) is radially unbounded, the equilibrium is globally asymptotically stable.

For \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \), choose \( V=\mathbf{x}^{T}\mathbf{P}\mathbf{x} \), where \( \mathbf{P}=\mathbf{P}^{T}>0 \). Then

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

Theorem 4. The matrix \( \mathbf{A} \) is Hurwitz if and only if, for every \( \mathbf{Q}=\mathbf{Q}^{T}>0 \), there is a unique \( \mathbf{P}=\mathbf{P}^{T}>0 \) satisfying

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

Constructive proof of the forward direction. If \( \mathbf{A} \) is Hurwitz, define

\[ \mathbf{P}=\int_{0}^{\infty} e^{\mathbf{A}^{T}\tau}\mathbf{Q}e^{\mathbf{A}\tau}\,d\tau. \]

Exponential decay of \( e^{\mathbf{A}\tau} \) guarantees convergence. For any nonzero \( \mathbf{x} \),

\[ \mathbf{x}^{T}\mathbf{P}\mathbf{x} =\int_{0}^{\infty} (e^{\mathbf{A}\tau}\mathbf{x})^{T} \mathbf{Q} (e^{\mathbf{A}\tau}\mathbf{x})\,d\tau>0, \]

so \( \mathbf{P}>0 \). Furthermore,

\[ \begin{aligned} \mathbf{A}^{T}\mathbf{P}+\mathbf{P}\mathbf{A} &=\int_{0}^{\infty} \frac{d}{d\tau} \left(e^{\mathbf{A}^{T}\tau}\mathbf{Q}e^{\mathbf{A}\tau}\right)d\tau\\ &=\left[e^{\mathbf{A}^{T}\tau}\mathbf{Q}e^{\mathbf{A}\tau}\right]_{0}^{\infty} =-\mathbf{Q}. \end{aligned} \]

For the converse direction, suppose a positive-definite solution exists. Then \( V=\mathbf{x}^{T}\mathbf{P}\mathbf{x} \) satisfies \( \dot V=-\mathbf{x}^{T}\mathbf{Q}\mathbf{x} \). Quadratic eigenvalue bounds give exponential stability of the origin, so every eigenvalue of \( \mathbf{A} \) must have a negative real part; hence \( \mathbf{A} \) is Hurwitz.

This equation appears repeatedly in adaptive control because a desired or nominal error matrix is deliberately selected to be Hurwitz, and its Lyapunov equation supplies the state-error weighting matrix used in a later augmented candidate.

Adaptive closed loops are generally time-varying even when the physical plant is time invariant, because estimated parameters evolve. For

\[ \dot{\mathbf{x}}=\mathbf{f}(t,\mathbf{x}), \]

a time-dependent candidate satisfies

\[ \dot V(t,\mathbf{x}) =\frac{\partial V}{\partial t} +\frac{\partial V}{\partial\mathbf{x}}\mathbf{f}(t,\mathbf{x}). \]

Uniform stability requires the admissible initial radius to be independent of the initial time. A standard sufficient condition is the existence of time-independent comparison functions such that

\[ \alpha_{1}(\|\mathbf{x}\|)\le V(t,\mathbf{x}) \le\alpha_{2}(\|\mathbf{x}\|), \qquad \dot V(t,\mathbf{x})\le-W(\mathbf{x}), \]

uniformly for all admissible \( t \). Uniformity matters because a certificate that deteriorates with starting time gives no reusable guarantee for an online controller.

A related robust viewpoint considers a parameterized family

\[ \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\boldsymbol{\theta}), \qquad \boldsymbol{\theta}\in\Theta. \]

If one function \( V(\mathbf{x}) \) satisfies the same positive-definite bounds and

\[ \nabla V(\mathbf{x})^{T} \mathbf{f}(\mathbf{x},\boldsymbol{\theta}) \le-W(\mathbf{x}) \quad\forall\boldsymbol{\theta}\in\Theta, \]

it is a common Lyapunov function for the family. The resulting stability statement is uniform with respect to every parameter value in \( \Theta \).

9.1 Why a state-only candidate may be inconclusive

Consider the scalar tracking-error model

\[ \dot e=-\lambda e+ ilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}(t), \qquad \lambda>0, \]

where \( \tilde{\boldsymbol{\theta}} \) is a parameter mismatch and \( \boldsymbol{\phi}(t) \) is a known regressor. With the state-only candidate \( V_{e}=\tfrac12 e^{2} \),

\[ \dot V_{e} =-\lambda e^{2} +e\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}(t). \]

The mixed term has no fixed sign. This does not prove that the system is unstable; it proves only that this candidate, with the information used so far, is insufficient. In the next lesson, the parameter mismatch is included in an augmented Lyapunov function. Later lessons choose an update law that shapes the mixed terms.

flowchart TD
  A["Nominal stable error term"] --> C["Differentiate the state-only candidate"]
  B["Unknown parameter mismatch times regressor"] --> C
  C --> D{"Does an uncertain \nmixed term remain?"}
  D -->|"no"| E["Use the current certificate"]
  D -->|"yes"| F["The state-only certificate is inconclusive"]
  F --> G["Next lesson: add parameter-error energy"]
  G --> H["Later: select update dynamics to shape cross terms"]
        

9.2 Uniformly stable uncertain scalar family

Now consider

\[ \dot e=-a(t)e, \qquad a(t)\ge a_{\min}>0. \]

For \( V=\tfrac12e^{2} \),

\[ \dot V=-a(t)e^{2}\le-a_{\min}e^{2}=-2a_{\min}V. \]

Therefore

\[ V(t)\le V(0)e^{-2a_{\min}t}, \qquad |e(t)|\le e^{-a_{\min}t}|e(0)|. \]

The result is uniform over every measurable coefficient satisfying the same lower bound. The numerical implementations below use

\[ a(t)=1+0.35\sin(0.8t)+0.15\cos(0.21t), \]

for which the conservative analytical bound \( a(t)\ge0.5 \) follows directly from the triangle inequality.

Each implementation performs the same fourth-order Runge–Kutta simulation, computes \( V \), \( \dot V \), and the comparison bound, and checks three numerical conditions:

\[ V_{k+1}\le V_{k},\qquad \dot V_{k}\le-2a_{\min}V_{k},\qquad V_{k}\le V(0)e^{-2a_{\min}t_{k}}. \]

• Python: NumPy for arrays, Matplotlib for plots, and the standard CSV library.
• C++: C++17 standard library only; results are exported to CSV.
• Java: JDK standard numerical and file APIs; results are exported to CSV.
• MATLAB/Simulink: base MATLAB for RK4 and plotting, plus an optional programmatically generated Simulink model.
• Wolfram Mathematica: NDSolveValue, symbolic function definitions, plotting, and CSV export.

Numerical checks support an analytical proof but do not replace it. A finite simulation samples only selected times and initial conditions, whereas a Lyapunov inequality quantifies an entire continuum of trajectories.

Chapter3_Lesson1.py

"""Chapter 3, Lesson 1: Lyapunov verification for a time-varying uncertain system.

System:
    e_dot = -a(t) e
    a(t) = 1 + 0.35 sin(0.8 t) + 0.15 cos(0.21 t)

Lyapunov candidate:
    V(e) = 0.5 e^2
    V_dot = -a(t) e^2 <= -2 a_min V,  a_min = 0.5

The script integrates the system with RK4, checks the Lyapunov inequalities,
exports a CSV file, and plots the state, V, its exponential bound, and V_dot.
"""

from __future__ import annotations

import csv
import math
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


def a_of_t(t: float) -> float:
    """Time-varying coefficient with the guaranteed bound a(t) >= 0.5."""
    return 1.0 + 0.35 * math.sin(0.8 * t) + 0.15 * math.cos(0.21 * t)


def rhs(t: float, e: float) -> float:
    """Scalar error dynamics."""
    return -a_of_t(t) * e


def rk4_step(t: float, e: float, dt: float) -> float:
    """One classical fourth-order Runge-Kutta step."""
    k1 = rhs(t, e)
    k2 = rhs(t + 0.5 * dt, e + 0.5 * dt * k1)
    k3 = rhs(t + 0.5 * dt, e + 0.5 * dt * k2)
    k4 = rhs(t + dt, e + dt * k3)
    return e + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def main() -> None:
    dt = 0.002
    final_time = 20.0
    initial_error = 2.5
    a_min = 0.5

    time = np.arange(0.0, final_time + dt, dt)
    error = np.zeros_like(time)
    error[0] = initial_error

    for k in range(time.size - 1):
        error[k + 1] = rk4_step(float(time[k]), float(error[k]), dt)

    coefficient = np.array([a_of_t(float(t)) for t in time])
    lyapunov = 0.5 * error**2
    lyapunov_dot = -coefficient * error**2
    exponential_bound = lyapunov[0] * np.exp(-2.0 * a_min * time)

    tolerance = 1.0e-10
    monotonic_violations = int(np.count_nonzero(np.diff(lyapunov) > tolerance))
    derivative_violations = int(
        np.count_nonzero(lyapunov_dot > -2.0 * a_min * lyapunov + tolerance)
    )
    bound_violations = int(
        np.count_nonzero(lyapunov > exponential_bound + 5.0e-8)
    )

    print(f"minimum sampled a(t): {coefficient.min():.8f}")
    print(f"maximum sampled V_dot: {lyapunov_dot.max():.8e}")
    print(f"monotonicity violations: {monotonic_violations}")
    print(f"differential-inequality violations: {derivative_violations}")
    print(f"exponential-bound violations: {bound_violations}")
    print(f"final error: {error[-1]:.8e}")

    output_path = Path("Chapter3_Lesson1.csv")
    with output_path.open("w", newline="", encoding="utf-8") as csv_file:
        writer = csv.writer(csv_file)
        writer.writerow(["time", "a", "error", "V", "V_dot", "V_bound"])
        writer.writerows(
            zip(time, coefficient, error, lyapunov, lyapunov_dot, exponential_bound)
        )
    print(f"wrote {output_path.resolve()}")

    plt.figure()
    plt.plot(time, error)
    plt.xlabel("time")
    plt.ylabel("e(t)")
    plt.title("State error")
    plt.grid(True)

    plt.figure()
    plt.semilogy(time, np.maximum(lyapunov, np.finfo(float).tiny), label="V(t)")
    plt.semilogy(time, exponential_bound, "--", label="guaranteed bound")
    plt.xlabel("time")
    plt.ylabel("Lyapunov value")
    plt.title("Lyapunov decay and comparison bound")
    plt.legend()
    plt.grid(True)

    plt.figure()
    plt.plot(time, lyapunov_dot)
    plt.xlabel("time")
    plt.ylabel("V_dot(t)")
    plt.title("Derivative of the Lyapunov candidate")
    plt.grid(True)

    plt.show()


if __name__ == "__main__":
    main()

Chapter3_Lesson1.cpp

// Chapter 3, Lesson 1: Lyapunov verification for a time-varying uncertain system.
// Build: g++ -std=c++17 -O2 Chapter3_Lesson1.cpp -o Chapter3_Lesson1

#include <algorithm>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>

namespace {

double a_of_t(double t) {
    return 1.0 + 0.35 * std::sin(0.8 * t) + 0.15 * std::cos(0.21 * t);
}

double rhs(double t, double e) {
    return -a_of_t(t) * e;
}

double rk4_step(double t, double e, double dt) {
    const double k1 = rhs(t, e);
    const double k2 = rhs(t + 0.5 * dt, e + 0.5 * dt * k1);
    const double k3 = rhs(t + 0.5 * dt, e + 0.5 * dt * k2);
    const double k4 = rhs(t + dt, e + dt * k3);
    return e + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
}

}  // namespace

int main() {
    try {
        constexpr double dt = 0.002;
        constexpr double final_time = 20.0;
        constexpr double initial_error = 2.5;
        constexpr double a_min = 0.5;
        constexpr double tolerance = 1.0e-10;

        const std::size_t sample_count =
            static_cast<std::size_t>(std::llround(final_time / dt)) + 1U;

        std::vector<double> time(sample_count);
        std::vector<double> coefficient(sample_count);
        std::vector<double> error(sample_count);
        std::vector<double> lyapunov(sample_count);
        std::vector<double> lyapunov_dot(sample_count);
        std::vector<double> exponential_bound(sample_count);

        error.front() = initial_error;
        for (std::size_t k = 0; k < sample_count; ++k) {
            time[k] = static_cast<double>(k) * dt;
            if (k + 1U < sample_count) {
                error[k + 1U] = rk4_step(time[k], error[k], dt);
            }
        }

        for (std::size_t k = 0; k < sample_count; ++k) {
            coefficient[k] = a_of_t(time[k]);
            lyapunov[k] = 0.5 * error[k] * error[k];
            lyapunov_dot[k] = -coefficient[k] * error[k] * error[k];
            exponential_bound[k] =
                lyapunov.front() * std::exp(-2.0 * a_min * time[k]);
        }

        std::size_t monotonic_violations = 0U;
        std::size_t derivative_violations = 0U;
        std::size_t bound_violations = 0U;

        for (std::size_t k = 0; k + 1U < sample_count; ++k) {
            if (lyapunov[k + 1U] > lyapunov[k] + tolerance) {
                ++monotonic_violations;
            }
        }
        for (std::size_t k = 0; k < sample_count; ++k) {
            if (lyapunov_dot[k] > -2.0 * a_min * lyapunov[k] + tolerance) {
                ++derivative_violations;
            }
            if (lyapunov[k] > exponential_bound[k] + 5.0e-8) {
                ++bound_violations;
            }
        }

        const auto min_a = *std::min_element(coefficient.begin(), coefficient.end());
        const auto max_v_dot =
            *std::max_element(lyapunov_dot.begin(), lyapunov_dot.end());

        std::cout << std::setprecision(10)
                  << "minimum sampled a(t): " << min_a << '\n'
                  << "maximum sampled V_dot: " << max_v_dot << '\n'
                  << "monotonicity violations: " << monotonic_violations << '\n'
                  << "differential-inequality violations: "
                  << derivative_violations << '\n'
                  << "exponential-bound violations: " << bound_violations << '\n'
                  << "final error: " << error.back() << '\n';

        std::ofstream csv("Chapter3_Lesson1.csv");
        if (!csv) {
            throw std::runtime_error("Could not create Chapter3_Lesson1.csv");
        }
        csv << "time,a,error,V,V_dot,V_bound\n";
        csv << std::setprecision(16);
        for (std::size_t k = 0; k < sample_count; ++k) {
            csv << time[k] << ',' << coefficient[k] << ',' << error[k] << ','
                << lyapunov[k] << ',' << lyapunov_dot[k] << ','
                << exponential_bound[k] << '\n';
        }
        std::cout << "wrote Chapter3_Lesson1.csv\n";
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "error: " << error.what() << '\n';
        return 1;
    }
}

Chapter3_Lesson1.java

// Chapter 3, Lesson 1: Lyapunov verification for a time-varying uncertain system.
// Build and run:
//   javac Chapter3_Lesson1.java
//   java Chapter3_Lesson1

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.charset.StandardCharsets;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Locale;

public final class Chapter3_Lesson1 {
    private Chapter3_Lesson1() {
        // Utility class.
    }

    private static double aOfT(double t) {
        return 1.0 + 0.35 * Math.sin(0.8 * t) + 0.15 * Math.cos(0.21 * t);
    }

    private static double rhs(double t, double e) {
        return -aOfT(t) * e;
    }

    private static double rk4Step(double t, double e, double dt) {
        double k1 = rhs(t, e);
        double k2 = rhs(t + 0.5 * dt, e + 0.5 * dt * k1);
        double k3 = rhs(t + 0.5 * dt, e + 0.5 * dt * k2);
        double k4 = rhs(t + dt, e + dt * k3);
        return e + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);

        final double dt = 0.002;
        final double finalTime = 20.0;
        final double initialError = 2.5;
        final double aMin = 0.5;
        final double tolerance = 1.0e-10;
        final int sampleCount = (int) Math.round(finalTime / dt) + 1;

        double[] time = new double[sampleCount];
        double[] coefficient = new double[sampleCount];
        double[] error = new double[sampleCount];
        double[] lyapunov = new double[sampleCount];
        double[] lyapunovDot = new double[sampleCount];
        double[] exponentialBound = new double[sampleCount];

        error[0] = initialError;
        for (int k = 0; k < sampleCount; ++k) {
            time[k] = k * dt;
            if (k + 1 < sampleCount) {
                error[k + 1] = rk4Step(time[k], error[k], dt);
            }
        }

        for (int k = 0; k < sampleCount; ++k) {
            coefficient[k] = aOfT(time[k]);
            lyapunov[k] = 0.5 * error[k] * error[k];
            lyapunovDot[k] = -coefficient[k] * error[k] * error[k];
            exponentialBound[k] = lyapunov[0] * Math.exp(-2.0 * aMin * time[k]);
        }

        int monotonicViolations = 0;
        int derivativeViolations = 0;
        int boundViolations = 0;
        double minimumA = Double.POSITIVE_INFINITY;
        double maximumVDot = Double.NEGATIVE_INFINITY;

        for (int k = 0; k < sampleCount; ++k) {
            minimumA = Math.min(minimumA, coefficient[k]);
            maximumVDot = Math.max(maximumVDot, lyapunovDot[k]);
            if (lyapunovDot[k] > -2.0 * aMin * lyapunov[k] + tolerance) {
                derivativeViolations++;
            }
            if (lyapunov[k] > exponentialBound[k] + 5.0e-8) {
                boundViolations++;
            }
            if (k + 1 < sampleCount && lyapunov[k + 1] > lyapunov[k] + tolerance) {
                monotonicViolations++;
            }
        }

        System.out.printf("minimum sampled a(t): %.8f%n", minimumA);
        System.out.printf("maximum sampled V_dot: %.8e%n", maximumVDot);
        System.out.printf("monotonicity violations: %d%n", monotonicViolations);
        System.out.printf("differential-inequality violations: %d%n", derivativeViolations);
        System.out.printf("exponential-bound violations: %d%n", boundViolations);
        System.out.printf("final error: %.8e%n", error[sampleCount - 1]);

        Path output = Path.of("Chapter3_Lesson1.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(
                output, StandardCharsets.UTF_8)) {
            writer.write("time,a,error,V,V_dot,V_bound");
            writer.newLine();
            for (int k = 0; k < sampleCount; ++k) {
                writer.write(String.format(
                        Locale.ROOT,
                        "%.16g,%.16g,%.16g,%.16g,%.16g,%.16g%n",
                        time[k], coefficient[k], error[k], lyapunov[k],
                        lyapunovDot[k], exponentialBound[k]));
            }
            System.out.println("wrote " + output.toAbsolutePath());
        } catch (IOException exception) {
            System.err.println("Could not write CSV: " + exception.getMessage());
            System.exit(1);
        }
    }
}

Chapter3_Lesson1.m

%% Chapter3_Lesson1.m
% Lyapunov verification for the time-varying uncertain system
%   e_dot = -a(t)e,
%   a(t) = 1 + 0.35 sin(0.8t) + 0.15 cos(0.21t).
%
% Candidate: V = 0.5 e^2.  Since a(t) >= 0.5,
%   V_dot = -a(t)e^2 <= -2(0.5)V.
%
% The first part is a MATLAB RK4 simulation.  The optional second part
% programmatically builds an equivalent Simulink model.

clear; clc; close all;

dt = 0.002;
finalTime = 20.0;
initialError = 2.5;
aMin = 0.5;
tolerance = 1.0e-10;

time = (0:dt:finalTime).';
errorSignal = zeros(size(time));
errorSignal(1) = initialError;

aOfT = @(tau) 1.0 + 0.35 .* sin(0.8 .* tau) + 0.15 .* cos(0.21 .* tau);
rhs = @(tau, value) -aOfT(tau) .* value;

for k = 1:numel(time)-1
    t = time(k);
    e = errorSignal(k);
    k1 = rhs(t, e);
    k2 = rhs(t + 0.5*dt, e + 0.5*dt*k1);
    k3 = rhs(t + 0.5*dt, e + 0.5*dt*k2);
    k4 = rhs(t + dt, e + dt*k3);
    errorSignal(k+1) = e + (dt/6.0) * (k1 + 2*k2 + 2*k3 + k4);
end

coefficient = aOfT(time);
V = 0.5 .* errorSignal.^2;
Vdot = -coefficient .* errorSignal.^2;
Vbound = V(1) .* exp(-2.0 .* aMin .* time);

monotonicViolations = nnz(diff(V) > tolerance);
derivativeViolations = nnz(Vdot > -2.0 .* aMin .* V + tolerance);
boundViolations = nnz(V > Vbound + 5.0e-8);

fprintf('minimum sampled a(t): %.8f\n', min(coefficient));
fprintf('maximum sampled V_dot: %.8e\n', max(Vdot));
fprintf('monotonicity violations: %d\n', monotonicViolations);
fprintf('differential-inequality violations: %d\n', derivativeViolations);
fprintf('exponential-bound violations: %d\n', boundViolations);
fprintf('final error: %.8e\n', errorSignal(end));

results = table(time, coefficient, errorSignal, V, Vdot, Vbound, ...
    'VariableNames', {'time','a','error','V','V_dot','V_bound'});
writetable(results, 'Chapter3_Lesson1.csv');

figure;
plot(time, errorSignal, 'LineWidth', 1.2);
xlabel('time'); ylabel('e(t)'); title('State error'); grid on;

figure;
semilogy(time, max(V, realmin), 'LineWidth', 1.2); hold on;
semilogy(time, Vbound, '--', 'LineWidth', 1.2);
xlabel('time'); ylabel('Lyapunov value');
title('Lyapunov decay and comparison bound');
legend('V(t)', 'guaranteed bound', 'Location', 'best'); grid on;

figure;
plot(time, Vdot, 'LineWidth', 1.2);
xlabel('time'); ylabel('Vdot(t)');
title('Derivative of the Lyapunov candidate'); grid on;

%% Optional programmatic Simulink implementation
% Set this flag to true when Simulink is installed.
buildSimulink = false;

if buildSimulink
    if ~license('test', 'Simulink')
        error('Simulink is not available on this MATLAB installation.');
    end

    model = 'Chapter3_Lesson1_Simulink';
    if bdIsLoaded(model)
        close_system(model, 0);
    end
    if exist([model '.slx'], 'file')
        delete([model '.slx']);
    end

    new_system(model);
    open_system(model);

    add_block('simulink/Sources/Constant', [model '/One'], ...
        'Value', '1.0', 'Position', [40 45 80 75]);
    add_block('simulink/Sources/Sine Wave', [model '/SinTerm'], ...
        'Amplitude', '0.35', 'Frequency', '0.8', ...
        'Position', [40 105 90 135]);
    add_block('simulink/Sources/Sine Wave', [model '/CosTerm'], ...
        'Amplitude', '0.15', 'Frequency', '0.21', 'Phase', 'pi/2', ...
        'Position', [40 165 90 195]);
    add_block('simulink/Math Operations/Sum', [model '/a_of_t'], ...
        'Inputs', '+++', 'Position', [145 90 175 170]);
    add_block('simulink/Math Operations/Product', [model '/a_times_e'], ...
        'Position', [250 100 290 150]);
    add_block('simulink/Math Operations/Gain', [model '/Minus'], ...
        'Gain', '-1', 'Position', [340 105 390 145]);
    add_block('simulink/Continuous/Integrator', [model '/Integrator'], ...
        'InitialCondition', num2str(initialError), ...
        'Position', [440 105 470 145]);
    add_block('simulink/Math Operations/Math Function', [model '/Square'], ...
        'Operator', 'square', 'Position', [535 180 585 215]);
    add_block('simulink/Math Operations/Gain', [model '/Half'], ...
        'Gain', '0.5', 'Position', [625 180 675 215]);
    add_block('simulink/Sinks/To Workspace', [model '/e_out'], ...
        'VariableName', 'e_sim', 'SaveFormat', 'Timeseries', ...
        'Position', [535 95 625 125]);
    add_block('simulink/Sinks/To Workspace', [model '/a_out'], ...
        'VariableName', 'a_sim', 'SaveFormat', 'Timeseries', ...
        'Position', [250 45 340 75]);
    add_block('simulink/Sinks/To Workspace', [model '/V_out'], ...
        'VariableName', 'V_sim', 'SaveFormat', 'Timeseries', ...
        'Position', [715 180 805 210]);

    add_line(model, 'One/1', 'a_of_t/1');
    add_line(model, 'SinTerm/1', 'a_of_t/2');
    add_line(model, 'CosTerm/1', 'a_of_t/3');
    add_line(model, 'a_of_t/1', 'a_times_e/1');
    add_line(model, 'a_of_t/1', 'a_out/1');
    add_line(model, 'Integrator/1', 'a_times_e/2');
    add_line(model, 'a_times_e/1', 'Minus/1');
    add_line(model, 'Minus/1', 'Integrator/1');
    add_line(model, 'Integrator/1', 'e_out/1');
    add_line(model, 'Integrator/1', 'Square/1');
    add_line(model, 'Square/1', 'Half/1');
    add_line(model, 'Half/1', 'V_out/1');

    set_param(model, 'StopTime', num2str(finalTime), 'Solver', 'ode45');
    save_system(model);
    sim(model);
    fprintf('Created and simulated %s.slx\n', model);
end

Chapter3_Lesson1.nb

Notebook[{
  Cell["Chapter 3, Lesson 1: Lyapunov Stability for Adaptive Use", "Title"],

  Cell[
    "Verification of a Lyapunov inequality for the time-varying system \
e'(t) = -a(t)e(t), where a(t) is uniformly positive.",
    "Text"
  ],

  Cell[BoxData[ToBoxes[
    ClearAll["Global`*"];

    a[t_] := 1.0 + 0.35 Sin[0.8 t] + 0.15 Cos[0.21 t];
    aMin = 0.5; tFinal = 20.0; e0 = 2.5;

    sol = NDSolveValue[
      {e'[t] == -a[t] e[t], e[0] == e0},
      e, {t, 0, tFinal},
      MaxStepFraction -> 1/2000
    ];

    v[t_] := 0.5 sol[t]^2;
    vDot[t_] := -a[t] sol[t]^2;
    vBound[t_] := 0.5 e0^2 Exp[-2 aMin t];

    grid = Subdivide[0.0, tFinal, 10000];

    minA = Min[a /@ grid];
    maxVDot = Max[vDot /@ grid];
    derivativeMargin = Max[(vDot[#] + 2 aMin v[#]) & /@ grid];
    boundMargin = Max[(v[#] - vBound[#]) & /@ grid];

    Print["Minimum sampled a(t): ", N[minA, 10]];
    Print["Maximum sampled V_dot: ", ScientificForm[maxVDot, 6]];
    Print["Maximum inequality margin: ", ScientificForm[derivativeMargin, 6]];
    Print["Maximum bound margin: ", ScientificForm[boundMargin, 6]];
    Print["Final error: ", ScientificForm[sol[tFinal], 6]];

    data = Table[
      {t, a[t], sol[t], v[t], vDot[t], vBound[t]},
      {t, grid}
    ];

    Export[
      "Chapter3_Lesson1.csv",
      Prepend[data, {"time", "a", "error", "V", "V_dot", "V_bound"}]
    ];

    statePlot = Plot[
      sol[t], {t, 0, tFinal},
      Frame -> True, FrameLabel -> {"time", "e(t)"},
      PlotLabel -> "State Error", ImageSize -> Large
    ];

    valuePlot = LogPlot[
      {v[t], vBound[t]}, {t, 0, tFinal},
      Frame -> True, FrameLabel -> {"time", "Lyapunov value"},
      PlotLegends -> {"V(t)", "Guaranteed bound"},
      PlotLabel -> "Lyapunov Decay", ImageSize -> Large
    ];

    derivativePlot = Plot[
      vDot[t], {t, 0, tFinal},
      Frame -> True, FrameLabel -> {"time", "V_dot(t)"},
      PlotLabel -> "Lyapunov Derivative", ImageSize -> Large
    ];

    Column[{statePlot, valuePlot, derivativePlot}, Spacings -> 2]
  ]], "Input"]
},
WindowTitle -> "Chapter3_Lesson1",
Saveable -> True
]        

Problem 1: Asymptotic but Nonexponential Scalar Dynamics

Consider \( \dot x=-ax^{3} \) with \( a>0 \). Use \( V=\tfrac12x^{2} \) to prove global asymptotic stability. Then determine whether the decay is exponential.

Solution.

\[ \dot V=x\dot x=-ax^{4}<0\quad\forall x\ne0. \]

The candidate is positive definite and radially unbounded, and its derivative is negative definite. Hence the origin is globally asymptotically stable. The exact solution follows from separation of variables:

\[ x(t)=\frac{x(0)}{\sqrt{1+2ax^{2}(0)t}}. \]

This decay is algebraic, not exponential. Equivalently, no positive constant \( c \) satisfies \( ax^{4}\ge c x^{2} \) in a full neighborhood of the origin.

Problem 2: Solve a Continuous-Time Lyapunov Equation

For

\[ \mathbf{A}=\begin{bmatrix}0&1\\-2&-3\end{bmatrix}, \qquad \mathbf{Q}=\mathbf{I}, \]

solve \( \mathbf{A}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}=-\mathbf{I} \) for symmetric \( \mathbf{P} \) and verify positive definiteness.

Solution. Write

\[ \mathbf{P}=\begin{bmatrix}p_{11}&p_{12}\\p_{12}&p_{22}\end{bmatrix}. \]

Equating matrix entries gives

\[ -4p_{12}=-1,\qquad p_{11}-3p_{12}-2p_{22}=0,\qquad 2p_{12}-6p_{22}=-1. \]

Thus

\[ \mathbf{P}= \begin{bmatrix}5/4&1/4\\1/4&1/4\end{bmatrix}. \]

Sylvester's criterion gives \( p_{11}=5/4>0 \) and \( \det\mathbf{P}=1/4>0 \); hence \( \mathbf{P}>0 \).

Problem 3: Complete a LaSalle Argument

For the nonlinear oscillator in Section 6, explain why \( v=0 \) alone is not an invariant set except at the origin.

Solution.

On \( v=0 \), the first equation gives \( \dot q=0 \), but the second gives

\[ \dot v=-kq-\beta q^{3}=-q(k+\beta q^{2}). \]

Since \( k>0 \) and \( \beta>0 \), this derivative is zero only for \( q=0 \). Every point with \( q\ne0 \) immediately leaves the set \( v=0 \). Therefore the largest invariant subset of \( \dot V=0 \) is \( (q,v)=(0,0) \).

Problem 4: Diagnose an Inconclusive Adaptive Candidate

Let \( \dot e=-e+\tilde\theta\sin t \) and \( V=\tfrac12e^{2} \). Show that \( \dot V \) is not negative semidefinite for arbitrary \( \tilde\theta \). Does that establish instability?

Solution.

\[ \dot V=-e^{2}+e\tilde\theta\sin t. \]

At a time for which \( \sin t=1 \), choose \( e=0.1 \) and \( \tilde\theta=10 \). Then \( \dot V=-0.01+1=0.99>0 \). Thus the candidate does not provide a monotonic certificate. This is not an instability proof; another Lyapunov function, an additional bound, or an adaptive mechanism may still establish stability.

Problem 5: Common Lyapunov Function for an Uncertain Matrix Family

Consider

\[ \dot{\mathbf{x}}= \begin{bmatrix}-1&\rho\\-\rho&-2\end{bmatrix}\mathbf{x}, \qquad \rho\in\mathbb{R}. \]

Use \( V=\tfrac12(x_{1}^{2}+x_{2}^{2}) \) to prove uniform global exponential stability for every \( \rho \).

Solution.

\[ \begin{aligned} \dot V &=x_{1}(-x_{1}+\rho x_{2}) +x_{2}(-\rho x_{1}-2x_{2})\\ &=-x_{1}^{2}-2x_{2}^{2} \le-\|\mathbf{x}\|^{2}=-2V. \end{aligned} \]

The uncertain skew-symmetric cross terms cancel. Therefore \( V(t)\le e^{-2t}V(0) \) and \( \|\mathbf{x}(t)\|\le e^{-t}\|\mathbf{x}(0)\| \), with constants independent of \( \rho \).

Problem 6: Numerical Derivatives and False Violations

A simulation estimates \( \dot V(t_{k}) \) by a forward difference. Small positive values appear even though the analytical derivative is nonpositive. Explain why and give a better verification procedure.

Solution.

Forward differencing has truncation error of order \( O(\Delta t) \), and floating-point subtraction magnifies roundoff when successive values are close. Solver error also perturbs the sampled state. A better procedure evaluates the analytical orbital derivative

\[ \dot V(t_{k})=\nabla V(\mathbf{x}_{k})^{T} \mathbf{f}(t_{k},\mathbf{x}_{k}) \]

and compares it with the required inequality using a documented numerical tolerance. A centered difference can be used as a secondary consistency check, but not as the primary proof.

Lyapunov analysis replaces explicit trajectory solutions with scalar energy-like certificates. Positive definiteness links the certificate to distance from the equilibrium; a nonpositive derivative creates invariant sublevel sets; a negative-definite derivative establishes attraction; and quadratic comparison inequalities provide exponential rates. LaSalle's principle handles semidefinite derivatives by identifying the largest invariant set on which the derivative vanishes. For stable linear dynamics, the matrix Lyapunov equation constructs a quadratic certificate directly. In adaptive control, uniformity with respect to time and uncertainty is essential, and an indefinite mixed state–parameter term is a diagnostic that the state-only candidate is insufficient. Lesson 2 builds the corresponding combined state–parameter Lyapunov functions.