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

• form an augmented state containing tracking and parameter errors;
• construct a positive-definite block-diagonal quadratic Lyapunov function;
• derive an adaptive update law by cancelling the mixed term in the Lyapunov derivative;
• obtain explicit uniform bounds on both error components;
• distinguish bounded parameter error from parameter convergence; and
• verify the theoretical identity numerically without mistaking numerical monotonicity for a proof.

Students are assumed to know Lyapunov stability, Hurwitz matrices, quadratic forms, eigenvalue bounds, and the continuous-time Lyapunov equation from linear and nonlinear control. No model-reference adaptive-control structure is required yet; the analysis is intentionally stated for a generic adaptive error model.

Let \( \mathbf{e}(t)\in\mathbb{R}^{n} \) denote a state or tracking error. Let \( \boldsymbol{\theta}^{*}\in\mathbb{R}^{p} \) be a constant ideal parameter vector and \( \hat{\boldsymbol{\theta} }(t) \) its online estimate. We use the convention

\[ \tilde{\boldsymbol{\theta} }(t) =\hat{\boldsymbol{\theta} }(t)-\boldsymbol{\theta}^{*}. \]

The augmented error coordinate is therefore

\[ \boldsymbol{\xi}(t)= \begin{bmatrix} \mathbf{e}(t)\\ \tilde{\boldsymbol{\theta} }(t) \end{bmatrix}\in\mathbb{R}^{n+p}. \]

A basic linearly parameterized adaptive error model has the form

\[ \dot{\mathbf{e} }=\mathbf{A}\mathbf{e} +\mathbf{W}(t)\tilde{\boldsymbol{\theta} }, \]

where \( \mathbf{A}\in\mathbb{R}^{n\times n} \) is Hurwitz and \( \mathbf{W}(t)\in\mathbb{R}^{n\times p} \) is a known, locally bounded regressor matrix assembled from measurable signals. A common single-input specialization is

\[ \mathbf{W}(t)=\mathbf{b}\boldsymbol{\phi}^{T}(t),\qquad \dot{\mathbf{e} }=\mathbf{A}\mathbf{e} +\mathbf{b}\boldsymbol{\phi}^{T}(t)\tilde{\boldsymbol{\theta} }. \]

flowchart TD
  R["Known regressor W(t)"] --> E["Error dynamics"]
  T["Parameter error theta_tilde"] --> E
  E --> S["State error e"]
  S --> U["Update law: theta_tilde_dot = -Gamma W_T P e"]
  R --> U
  U --> T
  S --> V["Joint quadratic energy V"]
  T --> V
  V --> D["Energy rate: dV = -e_T Q e"]
        

Because \( \mathbf{A} \) is Hurwitz, for every selected symmetric matrix \( \mathbf{Q}=\mathbf{Q}^{T}\succ0 \) there exists a unique symmetric matrix \( \mathbf{P}=\mathbf{P}^{T}\succ0 \) satisfying

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

The matrix \( \mathbf{P} \) defines the energy metric for the state error. The choice of \( \mathbf{Q} \) changes the geometry of this metric but does not change the given matrix \( \mathbf{A} \). In adaptive proofs, the main role of this equation is to convert the nominal state contribution in \( \dot V \) into a negative quadratic term.

\[ \mathbf{e}^{T} \left(\mathbf{A}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}\right)\mathbf{e} =-\mathbf{e}^{T}\mathbf{Q}\mathbf{e}. \]

Select an adaptation-gain matrix \( \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}\succ0 \) and define

\[ V(\mathbf{e},\tilde{\boldsymbol{\theta} }) =\mathbf{e}^{T}\mathbf{P}\mathbf{e} +\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta} }. \]

Equivalently,

\[ V=\boldsymbol{\xi}^{T}\mathbf{H}\boldsymbol{\xi},\qquad \mathbf{H}=\begin{bmatrix} \mathbf{P} & \mathbf{0}\\ \mathbf{0} & \boldsymbol{\Gamma}^{-1} \end{bmatrix}\succ0. \]

The first term measures tracking energy and the second measures parameter mismatch. The inverse gain appears deliberately: it allows the factor \( \boldsymbol{\Gamma} \) in the update law to cancel when the derivative is formed.

Some texts place a factor of \( 1/2 \) in front of each quadratic term. Either convention is valid, but all derivative factors must then be handled consistently.

For symmetric positive-definite matrices, the Rayleigh-quotient bounds give

\[ \lambda_{\min}(\mathbf{P})\lVert\mathbf{e}\rVert^{2} \leq \mathbf{e}^{T}\mathbf{P}\mathbf{e} \leq \lambda_{\max}(\mathbf{P})\lVert\mathbf{e}\rVert^{2}, \]

\[ \lambda_{\min}(\boldsymbol{\Gamma}^{-1}) \lVert\tilde{\boldsymbol{\theta} }\rVert^{2} \leq \tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta} } \leq \lambda_{\max}(\boldsymbol{\Gamma}^{-1}) \lVert\tilde{\boldsymbol{\theta} }\rVert^{2}. \]

Consequently, with

\[ \alpha_{1}=\min\left\{ \lambda_{\min}(\mathbf{P}), \lambda_{\min}(\boldsymbol{\Gamma}^{-1})\right\},\qquad \alpha_{2}=\max\left\{ \lambda_{\max}(\mathbf{P}), \lambda_{\max}(\boldsymbol{\Gamma}^{-1})\right\}, \]

we obtain the augmented-state bounds

\[ \alpha_{1}\lVert\boldsymbol{\xi}\rVert^{2} \leq V(\boldsymbol{\xi}) \leq \alpha_{2}\lVert\boldsymbol{\xi}\rVert^{2}. \]

Thus \( V \) is positive definite and radially unbounded in the combined error coordinate.

Assume the ideal parameter is constant. Then \( \dot{\tilde{\boldsymbol{\theta} } } =\dot{\hat{\boldsymbol{\theta} } } \). Differentiating the candidate along the augmented dynamics gives

\[ \begin{aligned} \dot V &=2\mathbf{e}^{T}\mathbf{P}\dot{\mathbf{e} } +2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\tilde{\boldsymbol{\theta} } }\\ &=\mathbf{e}^{T} \left(\mathbf{A}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}\right)\mathbf{e} +2\mathbf{e}^{T}\mathbf{P}\mathbf{W} \tilde{\boldsymbol{\theta} } +2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\tilde{\boldsymbol{\theta} } }. \end{aligned} \]

Choose the parameter-error dynamics as

\[ \boxed{ \dot{\tilde{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\mathbf{W}^{T}(t)\mathbf{P}\mathbf{e} } \qquad \left(\dot{\hat{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\mathbf{W}^{T}(t)\mathbf{P}\mathbf{e}\right). \]

The parameter term becomes

\[ 2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\tilde{\boldsymbol{\theta} } } =-2\tilde{\boldsymbol{\theta} }^{T} \mathbf{W}^{T}\mathbf{P}\mathbf{e}. \]

Since both mixed expressions are scalars and one is the transpose of the other,

\[ \mathbf{e}^{T}\mathbf{P}\mathbf{W} \tilde{\boldsymbol{\theta} } =\tilde{\boldsymbol{\theta} }^{T} \mathbf{W}^{T}\mathbf{P}\mathbf{e}, \]

so they cancel exactly. Using the Lyapunov equation yields

\[ \boxed{\dot V=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e}\leq0}. \]

flowchart TD
  A["Verify A is Hurwitz"] --> B["Choose Q positive definite"]
  B --> C["Solve A_T P + P A = -Q"]
  C --> D["Form V = e_T P e + theta_T Gamma^-1 theta"]
  D --> E["Differentiate along combined dynamics"]
  E --> F["Identify the mixed state-parameter term"]
  F --> G["Choose update law for exact cancellation"]
  G --> H["Obtain dV = -e_T Q e"]
  H --> I["Conclude stability, boundedness, and finite error energy"]
        

Theorem. Consider

\[ \dot{\mathbf{e} }=\mathbf{A}\mathbf{e} +\mathbf{W}(t)\tilde{\boldsymbol{\theta} },\qquad \dot{\tilde{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\mathbf{W}^{T}(t)\mathbf{P}\mathbf{e}, \]

where \( \mathbf{A} \) is Hurwitz, \( \mathbf{Q}\succ0 \), \( \mathbf{P}\succ0 \) solves the Lyapunov equation, \( \boldsymbol{\Gamma}\succ0 \), and \( \mathbf{W}(t) \) is locally bounded. Then the augmented origin is Lyapunov stable, and both \( \mathbf{e}(t) \) and \( \tilde{\boldsymbol{\theta} }(t) \) remain bounded for all forward time for which the model is defined.

Proof. Positive definiteness gives

\[ V(t)\leq V(0), \]

and therefore

\[ \lVert\mathbf{e}(t)\rVert \leq\sqrt{\frac{V(0)}{\lambda_{\min}(\mathbf{P})} }, \]

\[ \lVert\tilde{\boldsymbol{\theta} }(t)\rVert \leq\sqrt{\frac{V(0)} {\lambda_{\min}(\boldsymbol{\Gamma}^{-1})} } =\sqrt{V(0)\lambda_{\max}(\boldsymbol{\Gamma})}. \]

For every prescribed neighbourhood of the origin, the quadratic bounds in Section 5 allow a sufficiently small initial level set to be selected inside that neighbourhood. Because \( V(t) \) cannot increase, the trajectory cannot leave that level set. This proves Lyapunov stability and the stated boundedness. \( \square \)

Integrating the exact derivative identity from zero to \( T \) gives

\[ \int_{0}^{T}\mathbf{e}^{T}(t)\mathbf{Q}\mathbf{e}(t)\,dt =V(0)-V(T)\leq V(0). \]

Because

\[ \mathbf{e}^{T}\mathbf{Q}\mathbf{e} \geq\lambda_{\min}(\mathbf{Q})\lVert\mathbf{e}\rVert^{2}, \]

we obtain

\[ \int_{0}^{\infty}\lVert\mathbf{e}(t)\rVert^{2}\,dt \leq\frac{V(0)}{\lambda_{\min}(\mathbf{Q})}. \]

Thus the tracking error has finite squared energy. This result is stronger than mere boundedness but is not, by itself, a proof that the tracking error converges pointwise to zero. The additional regularity argument needed for that conclusion is developed later in the chapter when Barbalat's lemma is introduced.

The derivative is negative definite in \( \mathbf{e} \) but only negative semidefinite in the full augmented coordinate:

\[ \dot V=0\quad\text{whenever}\quad\mathbf{e}=\mathbf{0}, \]

regardless of the instantaneous value of \( \tilde{\boldsymbol{\theta} } \). Therefore the proof does not force the parameter error to vanish. A parameter component that is not sufficiently represented in \( \mathbf{W}(t) \) can remain unidentified while the tracking error becomes small.

A simple extreme case is \( \mathbf{W}(t)=\mathbf{0} \) after some time \( T_{0} \). Then

\[ \dot{\tilde{\boldsymbol{\theta} } }=\mathbf{0},\qquad t\geq T_{0}, \]

so the parameter error freezes at its value at \( T_{0} \), while the stable nominal dynamics can continue to reduce \( \mathbf{e} \). Conditions for parameter convergence are postponed to the later chapter on persistent excitation.

The sublevel set

\[ \Omega_{c}=\left\{ (\mathbf{e},\tilde{\boldsymbol{\theta} }): \mathbf{e}^{T}\mathbf{P}\mathbf{e} +\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta} } \leq c\right\} \]

is an ellipsoid in the augmented space. Since \( \dot V\leq0 \), every such sublevel set is forward invariant under the assumed model.

Increasing an eigenvalue of \( \boldsymbol{\Gamma} \) decreases the corresponding weight in \( \boldsymbol{\Gamma}^{-1} \). Geometrically, the same level set permits a larger extent in that parameter-error direction. At the same time, the update law has a larger gain in that direction. This does not imply that a larger gain automatically gives a faster monotone decrease of \( V \), because after exact cancellation the expression \( \dot V=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e} \) contains no explicit adaptation gain.

Consider

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

Solving the Lyapunov equation gives

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

Its eigenvalues are \( (3-\sqrt{5})/4 \) and \( (3+\sqrt{5})/4 \), so it is positive definite. Let

\[ \mathbf{W}(t)=\mathbf{b}\boldsymbol{\phi}^{T}(t),\qquad \boldsymbol{\Gamma}=\operatorname{diag}(3,1.5). \]

The candidate becomes

\[ V=\frac{5}{4}e_{1}^{2} +\frac{1}{2}e_{1}e_{2} +\frac{1}{4}e_{2}^{2} +\frac{1}{3}\tilde{\theta}_{1}^{2} +\frac{2}{3}\tilde{\theta}_{2}^{2}. \]

The adaptive law is

\[ \dot{\tilde{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi}(t) \mathbf{b}^{T}\mathbf{P}\mathbf{e}, \]

and the derivative reduces exactly to

\[ \dot V=-(e_{1}^{2}+e_{2}^{2}). \]

The computational laboratories use the bounded decaying regressor

\[ \boldsymbol{\phi}(t)=e^{-0.15t} \begin{bmatrix}1\\0.5\end{bmatrix}. \]

It is deliberately not sufficiently rich to identify both independent parameter components. The numerical trajectory therefore demonstrates decreasing joint energy and state-error decay while a nonzero parameter residual can remain.

A more general quadratic candidate is

\[ V=\boldsymbol{\xi}^{T} \begin{bmatrix} \mathbf{P}&\mathbf{S}\\ \mathbf{S}^{T}&\mathbf{R} \end{bmatrix} \boldsymbol{\xi}. \]

If \( \mathbf{P}\succ0 \), the block matrix is positive definite if and only if its Schur complement is positive definite:

\[ \mathbf{R}-\mathbf{S}^{T}\mathbf{P}^{-1}\mathbf{S}\succ0. \]

Cross blocks can reduce conservatism in specialized analyses, but they also introduce additional derivative terms. The standard block-diagonal choice \( \mathbf{S}=\mathbf{0} \) is preferred at this stage because it exposes the cancellation mechanism and separates state and parameter metrics cleanly.

If the ideal parameter is time varying, the convention \( \tilde{\boldsymbol{\theta} } =\hat{\boldsymbol{\theta} }-\boldsymbol{\theta}^{*}(t) \) gives

\[ \dot{\tilde{\boldsymbol{\theta} } } =\dot{\hat{\boldsymbol{\theta} } } -\dot{\boldsymbol{\theta} }^{*}. \]

Using the same estimate update law produces

\[ \dot V=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e} -2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\boldsymbol{\theta} }^{*}. \]

The new term has no fixed sign. Thus the clean nonincrease proof depends on the constant-ideal-parameter assumption. Later robust modifications handle disturbances, drift, and modelling mismatch; they must not be silently omitted from the theorem stated here.

The implementations below perform the same steps:

1. define \( \mathbf{A},\mathbf{P},\mathbf{Q},\boldsymbol{\Gamma} \) and verify the Lyapunov-equation residual;
2. integrate the augmented continuous-time dynamics;
3. evaluate \( V(t) \) and the exact analytical rate;
4. check the largest sampled increment of \( V \); and
5. report the final parameter error to illustrate that stability does not imply identification.

A sampled increase caused by a coarse numerical solver would be a discretization artifact, not a contradiction of the continuous-time proof. Reducing the step size or using a higher-order solver should reduce that artifact.

Libraries: NumPy for linear algebra and Matplotlib for plots.

Chapter3_Lesson2.py

"""Chapter3_Lesson2.py
Quadratic Lyapunov function for combined state-parameter dynamics.

Requires: numpy, matplotlib
"""
from __future__ import annotations

import numpy as np
import matplotlib.pyplot as plt

A = np.array([[0.0, 1.0], [-2.0, -3.0]])
b = np.array([0.0, 1.0])
P = np.array([[1.25, 0.25], [0.25, 0.25]])
Q = np.eye(2)
GAMMA = np.diag([3.0, 1.5])
GAMMA_INV = np.linalg.inv(GAMMA)


def regressor(t: float) -> np.ndarray:
    """A bounded, decaying, non-persistently-exciting regressor."""
    return np.exp(-0.15 * t) * np.array([1.0, 0.5])


def dynamics(t: float, z: np.ndarray) -> np.ndarray:
    e = z[:2]
    theta_tilde = z[2:]
    phi = regressor(t)
    coupling = b * float(phi @ theta_tilde)
    e_dot = A @ e + coupling
    theta_tilde_dot = -GAMMA @ phi * float(b @ P @ e)
    return np.concatenate((e_dot, theta_tilde_dot))


def lyapunov(z: np.ndarray) -> float:
    e = z[:2]
    theta_tilde = z[2:]
    return float(e @ P @ e + theta_tilde @ GAMMA_INV @ theta_tilde)


def lyapunov_rate(z: np.ndarray) -> float:
    e = z[:2]
    return -float(e @ Q @ e)


def rk4_step(t: float, z: np.ndarray, dt: float) -> np.ndarray:
    k1 = dynamics(t, z)
    k2 = dynamics(t + 0.5 * dt, z + 0.5 * dt * k1)
    k3 = dynamics(t + 0.5 * dt, z + 0.5 * dt * k2)
    k4 = dynamics(t + dt, z + dt * k3)
    return z + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def main() -> None:
    dt = 0.002
    final_time = 30.0
    time = np.arange(0.0, final_time + dt, dt)
    z = np.zeros((time.size, 4))
    z[0] = np.array([1.0, -0.5, 1.2, -0.8])

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

    values = np.array([lyapunov(state) for state in z])
    rates = np.array([lyapunov_rate(state) for state in z])
    increments = np.diff(values)

    print(f"A^T P + P A + Q residual: {np.linalg.norm(A.T @ P + P @ A + Q):.3e}")
    print(f"V(0) = {values[0]:.9f}")
    print(f"V(T) = {values[-1]:.9f}")
    print(f"Largest sampled V increment = {increments.max():.3e}")
    print("Final parameter error:", z[-1, 2:])
    print("The nonzero residual illustrates that V_dot <= 0 alone does not prove parameter convergence.")

    np.savetxt(
        "Chapter3_Lesson2_results.csv",
        np.column_stack((time, z, values, rates)),
        delimiter=",",
        header="t,e1,e2,theta_tilde1,theta_tilde2,V,V_dot_exact",
        comments="",
    )

    plt.figure()
    plt.plot(time, z[:, 0], label="e1")
    plt.plot(time, z[:, 1], label="e2")
    plt.xlabel("Time")
    plt.ylabel("State error")
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter3_Lesson2_state_error.png", dpi=200)

    plt.figure()
    plt.plot(time, z[:, 2], label="theta_tilde1")
    plt.plot(time, z[:, 3], label="theta_tilde2")
    plt.xlabel("Time")
    plt.ylabel("Parameter error")
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter3_Lesson2_parameter_error.png", dpi=200)

    plt.figure()
    plt.plot(time, values, label="V")
    plt.xlabel("Time")
    plt.ylabel("Joint Lyapunov function")
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter3_Lesson2_lyapunov.png", dpi=200)

    plt.show()


if __name__ == "__main__":
    main()

This version uses only the C++17 standard library. Fixed-size arrays make every matrix operation visible rather than hiding the adaptive-law structure behind a numerical package.

Chapter3_Lesson2.cpp

// Chapter3_Lesson2.cpp
// Quadratic Lyapunov function for combined state-parameter dynamics.
// Build: g++ -std=c++17 -O2 Chapter3_Lesson2.cpp -o Chapter3_Lesson2

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

using State = std::array<double, 4>;

constexpr double DT = 0.002;
constexpr double FINAL_TIME = 30.0;

std::array<double, 2> regressor(double t) {
    const double scale = std::exp(-0.15 * t);
    return {scale, 0.5 * scale};
}

State add(const State& a, const State& b, double scale) {
    State result{};
    for (std::size_t i = 0; i < result.size(); ++i) {
        result[i] = a[i] + scale * b[i];
    }
    return result;
}

State dynamics(double t, const State& z) {
    const double e1 = z[0];
    const double e2 = z[1];
    const double th1 = z[2];
    const double th2 = z[3];
    const auto phi = regressor(t);

    const double phi_theta = phi[0] * th1 + phi[1] * th2;
    const double bTPe = 0.25 * e1 + 0.25 * e2;

    State dz{};
    dz[0] = e2;
    dz[1] = -2.0 * e1 - 3.0 * e2 + phi_theta;
    dz[2] = -3.0 * phi[0] * bTPe;
    dz[3] = -1.5 * phi[1] * bTPe;
    return dz;
}

State rk4Step(double t, const State& z) {
    const State k1 = dynamics(t, z);
    const State k2 = dynamics(t + 0.5 * DT, add(z, k1, 0.5 * DT));
    const State k3 = dynamics(t + 0.5 * DT, add(z, k2, 0.5 * DT));
    const State k4 = dynamics(t + DT, add(z, k3, DT));

    State next{};
    for (std::size_t i = 0; i < next.size(); ++i) {
        next[i] = z[i] + (DT / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
    }
    return next;
}

double lyapunov(const State& z) {
    const double e1 = z[0];
    const double e2 = z[1];
    const double th1 = z[2];
    const double th2 = z[3];
    const double statePart = 1.25 * e1 * e1 + 0.5 * e1 * e2 + 0.25 * e2 * e2;
    const double parameterPart = th1 * th1 / 3.0 + th2 * th2 / 1.5;
    return statePart + parameterPart;
}

double lyapunovRate(const State& z) {
    return -(z[0] * z[0] + z[1] * z[1]);
}

int main() {
    std::ofstream csv("Chapter3_Lesson2_results_cpp.csv");
    if (!csv) {
        throw std::runtime_error("Cannot create output CSV file.");
    }
    csv << "t,e1,e2,theta_tilde1,theta_tilde2,V,V_dot_exact\n";
    csv << std::setprecision(12);

    State z{1.0, -0.5, 1.2, -0.8};
    double previousV = lyapunov(z);
    const double initialV = previousV;
    double largestIncrement = -1.0e300;
    const int steps = static_cast<int>(std::llround(FINAL_TIME / DT));

    for (int k = 0; k <= steps; ++k) {
        const double t = k * DT;
        const double value = lyapunov(z);
        csv << t << ',' << z[0] << ',' << z[1] << ',' << z[2] << ',' << z[3]
            << ',' << value << ',' << lyapunovRate(z) << '\n';

        if (k > 0) {
            largestIncrement = std::max(largestIncrement, value - previousV);
        }
        previousV = value;
        if (k < steps) {
            z = rk4Step(t, z);
        }
    }

    std::cout << std::scientific << std::setprecision(6);
    std::cout << "V(0) = " << initialV << '\n';
    std::cout << "V(T) = " << lyapunov(z) << '\n';
    std::cout << "Largest sampled V increment = " << largestIncrement << '\n';
    std::cout << "Final parameter error = [" << z[2] << ", " << z[3] << "]\n";
    return 0;
}

This Java implementation uses the standard runtime only and writes the trajectory to a CSV file for external plotting.

Chapter3_Lesson2.java

// Chapter3_Lesson2.java
// Quadratic Lyapunov function for combined state-parameter dynamics.
// Build and run: javac Chapter3_Lesson2.java && java Chapter3_Lesson2

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

public final class Chapter3_Lesson2 {
    private static final double DT = 0.002;
    private static final double FINAL_TIME = 30.0;

    private Chapter3_Lesson2() {}

    private static double[] regressor(double t) {
        double scale = Math.exp(-0.15 * t);
        return new double[] {scale, 0.5 * scale};
    }

    private static double[] dynamics(double t, double[] z) {
        double e1 = z[0];
        double e2 = z[1];
        double th1 = z[2];
        double th2 = z[3];
        double[] phi = regressor(t);

        double phiTheta = phi[0] * th1 + phi[1] * th2;
        double bTPe = 0.25 * e1 + 0.25 * e2;

        return new double[] {
            e2,
            -2.0 * e1 - 3.0 * e2 + phiTheta,
            -3.0 * phi[0] * bTPe,
            -1.5 * phi[1] * bTPe
        };
    }

    private static double[] addScaled(double[] a, double[] b, double scale) {
        double[] result = new double[a.length];
        for (int i = 0; i < a.length; i++) {
            result[i] = a[i] + scale * b[i];
        }
        return result;
    }

    private static double[] rk4Step(double t, double[] z) {
        double[] k1 = dynamics(t, z);
        double[] k2 = dynamics(t + 0.5 * DT, addScaled(z, k1, 0.5 * DT));
        double[] k3 = dynamics(t + 0.5 * DT, addScaled(z, k2, 0.5 * DT));
        double[] k4 = dynamics(t + DT, addScaled(z, k3, DT));

        double[] next = new double[z.length];
        for (int i = 0; i < z.length; i++) {
            next[i] = z[i] + (DT / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        }
        return next;
    }

    private static double lyapunov(double[] z) {
        double e1 = z[0];
        double e2 = z[1];
        double th1 = z[2];
        double th2 = z[3];
        double statePart = 1.25 * e1 * e1 + 0.5 * e1 * e2 + 0.25 * e2 * e2;
        double parameterPart = th1 * th1 / 3.0 + th2 * th2 / 1.5;
        return statePart + parameterPart;
    }

    private static double lyapunovRate(double[] z) {
        return -(z[0] * z[0] + z[1] * z[1]);
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.ROOT);
        double[] z = {1.0, -0.5, 1.2, -0.8};
        int steps = (int) Math.round(FINAL_TIME / DT);
        double initialV = lyapunov(z);
        double previousV = initialV;
        double largestIncrement = -Double.MAX_VALUE;

        try (BufferedWriter writer = Files.newBufferedWriter(Path.of("Chapter3_Lesson2_results_java.csv"))) {
            writer.write("t,e1,e2,theta_tilde1,theta_tilde2,V,V_dot_exact\n");
            for (int k = 0; k <= steps; k++) {
                double t = k * DT;
                double value = lyapunov(z);
                writer.write(String.format(
                    Locale.ROOT,
                    "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                    t, z[0], z[1], z[2], z[3], value, lyapunovRate(z)
                ));

                if (k > 0) {
                    largestIncrement = Math.max(largestIncrement, value - previousV);
                }
                previousV = value;
                if (k < steps) {
                    z = rk4Step(t, z);
                }
            }
        }

        System.out.printf(Locale.ROOT, "V(0) = %.9f%n", initialV);
        System.out.printf(Locale.ROOT, "V(T) = %.9f%n", lyapunov(z));
        System.out.printf(Locale.ROOT, "Largest sampled V increment = %.3e%n", largestIncrement);
        System.out.printf(Locale.ROOT, "Final parameter error = [%.9f, %.9f]%n", z[2], z[3]);
    }
}

The MATLAB script uses ode45, verifies the matrix identity, plots all error variables, and exports a table.

Chapter3_Lesson2.m

% Chapter3_Lesson2.m
% Quadratic Lyapunov function for combined state-parameter dynamics.

clear; clc; close all;

A = [0 1; -2 -3];
b = [0; 1];
P = [5/4 1/4; 1/4 1/4];
Q = eye(2);
Gamma = diag([3, 1.5]);
GammaInv = inv(Gamma);

residual = norm(A' * P + P * A + Q, 'fro');
fprintf('Lyapunov-equation residual: %.3e\n', residual);

z0 = [1; -0.5; 1.2; -0.8];
tspan = [0 30];
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-11);
[t, z] = ode45(@augmentedDynamics, tspan, z0, options, A, b, P, Gamma);

V = zeros(size(t));
VdotExact = zeros(size(t));
for k = 1:numel(t)
    e = z(k, 1:2)';
    thetaTilde = z(k, 3:4)';
    V(k) = e' * P * e + thetaTilde' * GammaInv * thetaTilde;
    VdotExact(k) = -e' * Q * e;
end

fprintf('V(0) = %.9f\n', V(1));
fprintf('V(T) = %.9f\n', V(end));
fprintf('Largest sampled V increment = %.3e\n', max(diff(V)));
fprintf('Final parameter error = [%.9f, %.9f]\n', z(end, 3), z(end, 4));

results = table(t, z(:,1), z(:,2), z(:,3), z(:,4), V, VdotExact, ...
    'VariableNames', {'t','e1','e2','thetaTilde1','thetaTilde2','V','VdotExact'});
writetable(results, 'Chapter3_Lesson2_results_matlab.csv');

figure;
plot(t, z(:,1), t, z(:,2), 'LineWidth', 1.2);
xlabel('Time'); ylabel('State error'); legend('e_1','e_2'); grid on;

figure;
plot(t, z(:,3), t, z(:,4), 'LineWidth', 1.2);
xlabel('Time'); ylabel('Parameter error');
legend('theta-tilde_1','theta-tilde_2'); grid on;

figure;
plot(t, V, 'LineWidth', 1.2);
xlabel('Time'); ylabel('Joint Lyapunov function V'); grid on;

function dz = augmentedDynamics(t, z, A, b, P, Gamma)
    e = z(1:2);
    thetaTilde = z(3:4);
    phi = exp(-0.15 * t) * [1; 0.5];
    eDot = A * e + b * (phi' * thetaTilde);
    thetaTildeDot = -Gamma * phi * (b' * P * e);
    dz = [eDot; thetaTildeDot];
end

Running the following script creates Chapter3_Lesson2_Model.slx. It places the four augmented states in one vector Integrator block and implements the differential equations in a MATLAB Function block. Simulink is required.

Chapter3_Lesson2_Simulink.m

% Chapter3_Lesson2_Simulink.m
% Programmatically builds and runs a Simulink realization of the augmented ODE.
% Requires Simulink. The generated model is Chapter3_Lesson2_Model.slx.

clear; clc;
model = 'Chapter3_Lesson2_Model';

if bdIsLoaded(model)
    close_system(model, 0);
end
if isfile([model '.slx'])
    delete([model '.slx']);
end

new_system(model);
open_system(model);

add_block('simulink/Sources/Clock', [model '/Clock'], ...
    'Position', [40 55 70 85]);
add_block('simulink/Continuous/Integrator', [model '/AugmentedState'], ...
    'InitialCondition', '[1;-0.5;1.2;-0.8]', ...
    'Position', [365 110 405 160]);
add_block('simulink/User-Defined Functions/MATLAB Function', ...
    [model '/AdaptiveDynamics'], 'Position', [145 75 300 195]);
add_block('simulink/Sinks/To Workspace', [model '/StateLog'], ...
    'VariableName', 'z_log', 'SaveFormat', 'Structure With Time', ...
    'Position', [470 95 560 125]);
add_block('simulink/Sinks/To Workspace', [model '/VLog'], ...
    'VariableName', 'V_log', 'SaveFormat', 'Structure With Time', ...
    'Position', [470 145 560 175]);
add_block('simulink/Sinks/To Workspace', [model '/VdotLog'], ...
    'VariableName', 'Vdot_log', 'SaveFormat', 'Structure With Time', ...
    'Position', [470 195 560 225]);

root = sfroot;
chart = find(root, '-isa', 'Stateflow.EMChart', ...
    'Path', [model '/AdaptiveDynamics']);
chart.Script = sprintf([ ...
    'function [dz,V,Vdot] = fcn(t,z)\n' ...
    'A = [0 1; -2 -3];\n' ...
    'b = [0;1];\n' ...
    'P = [5/4 1/4; 1/4 1/4];\n' ...
    'Gamma = diag([3 1.5]);\n' ...
    'GammaInv = inv(Gamma);\n' ...
    'Q = eye(2);\n' ...
    'e = z(1:2);\n' ...
    'thetaTilde = z(3:4);\n' ...
    'phi = exp(-0.15*t)*[1;0.5];\n' ...
    'eDot = A*e + b*(phi''*thetaTilde);\n' ...
    'thetaTildeDot = -Gamma*phi*(b''*P*e);\n' ...
    'dz = [eDot;thetaTildeDot];\n' ...
    'V = e''*P*e + thetaTilde''*GammaInv*thetaTilde;\n' ...
    'Vdot = -e''*Q*e;\n' ...
    'end\n']);

add_line(model, 'Clock/1', 'AdaptiveDynamics/1', 'autorouting', 'on');
add_line(model, 'AugmentedState/1', 'AdaptiveDynamics/2', 'autorouting', 'on');
add_line(model, 'AdaptiveDynamics/1', 'AugmentedState/1', 'autorouting', 'on');
add_line(model, 'AugmentedState/1', 'StateLog/1', 'autorouting', 'on');
add_line(model, 'AdaptiveDynamics/2', 'VLog/1', 'autorouting', 'on');
add_line(model, 'AdaptiveDynamics/3', 'VdotLog/1', 'autorouting', 'on');

set_param(model, 'Solver', 'ode45', 'StopTime', '30');
save_system(model, [model '.slx']);
set_param(model, 'SimulationCommand', 'update');
sim(model);

fprintf('Generated and simulated %s.slx\n', model);
open_system(model);

The notebook uses NDSolveValue for integration, evaluates the quadratic form, exports a CSV file, and generates three plots.

Chapter3_Lesson2.nb

Notebook[{
 Cell["Chapter 3, Lesson 2: Quadratic Lyapunov Functions", "Title"],
 Cell["Combined State–Parameter Dynamics", "Subtitle"],
 Cell["Simulation of a coupled state-error and parameter-error system using a block-diagonal quadratic Lyapunov function.", "Text"],
 Cell[BoxData @ ToBoxes @ Defer[
  ClearAll["Global`*"];
  initialTime = 0.; finalTime = 30.; numberOfSamplingIntervals = 15000;
  initialStateError = {1., -0.5}; initialParameterError = {1.2, -0.8}; initialCombinedState = Join[initialStateError, initialParameterError];
  A = { {0., 1.}, {-2., -3.} }; b = {0., 1.}; Q = IdentityMatrix[2];
  P = { {5./4., 1./4.}, {1./4., 1./4.} }; pSymmetryResidual = Norm[P - Transpose[P]]; pEigenvalues = Eigenvalues[P];
  Gamma = DiagonalMatrix[{3., 1.5}]; GammaInverse = Inverse[Gamma]; gammaEigenvalues = Eigenvalues[Gamma];
  zeroBlock = ConstantArray[0., {2, 2}]; PCombined = ArrayFlatten[{ {P, zeroBlock}, {zeroBlock, GammaInverse} }]; combinedMatrixEigenvalues = Eigenvalues[PCombined];
  regressorDecayRate = 0.15; regressorDirection = {1., 0.5}; phi[t_?NumericQ] := Exp[-regressorDecayRate t] regressorDirection;
  combinedDynamics[t_?NumericQ, z_?VectorQ] := Module[{e = z[[1 ;; 2]], th = z[[3 ;; 4]], ph = phi[t], r, s}, r = ph.th; s = b.P.e; Join[A.e + b r, -(Gamma.ph) s]];
  combinedSolution = NDSolveValue[{z'[t] == combinedDynamics[t, z[t]], z[initialTime] == initialCombinedState}, z, {t, initialTime, finalTime},
    Method -> {"TimeIntegration" -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4} }, MaxSteps -> Infinity, AccuracyGoal -> 10, PrecisionGoal -> 8];
  combinedVector[t_?NumericQ] := combinedSolution[t]; stateErrorVector[t_?NumericQ] := combinedSolution[t][[1 ;; 2]];
  parameterErrorVector[t_?NumericQ] := combinedSolution[t][[3 ;; 4]]; stateError1[t_?NumericQ] := stateErrorVector[t][[1]]; stateError2[t_?NumericQ] := stateErrorVector[t][[2]];
  parameterError1[t_?NumericQ] := parameterErrorVector[t][[1]]; parameterError2[t_?NumericQ] := parameterErrorVector[t][[2]];
  stateLyapunovTerm[t_?NumericQ] := With[{e = stateErrorVector[t]}, e.P.e]; parameterLyapunovTerm[t_?NumericQ] := With[{th = parameterErrorVector[t]}, th.GammaInverse.th];
  lyapunovFunction[t_?NumericQ] := stateLyapunovTerm[t] + parameterLyapunovTerm[t]; combinedMatrixLyapunovFunction[t_?NumericQ] := With[{x = combinedVector[t]}, x.PCombined.x];
  exactLyapunovDerivative[t_?NumericQ] := With[{e = stateErrorVector[t]}, -e.Q.e];
  directLyapunovDerivative[t_?NumericQ] := Module[{x = combinedVector[t], dx, e, th, ed, thd}, dx = combinedDynamics[t, x]; e = x[[1 ;; 2]]; th = x[[3 ;; 4]];
    ed = dx[[1 ;; 2]]; thd = dx[[3 ;; 4]]; 2. e.P.ed + 2. th.GammaInverse.thd];
  lyapunovEquationResidualMatrix = Transpose[A].P + P.A + Q; lyapunovEquationResidualNorm = Norm[lyapunovEquationResidualMatrix];
  samplingGrid = Subdivide[initialTime, finalTime, numberOfSamplingIntervals]; numberOfSamplingPoints = Length[samplingGrid];
  sampledCombinedVectors = combinedVector /@ samplingGrid; sampledStateErrors = stateErrorVector /@ samplingGrid; sampledParameterErrors = parameterErrorVector /@ samplingGrid;
  sampledStateLyapunovTerms = stateLyapunovTerm /@ samplingGrid; sampledParameterLyapunovTerms = parameterLyapunovTerm /@ samplingGrid; sampledLyapunovValues = lyapunovFunction /@ samplingGrid;
  sampledCombinedMatrixLyapunovValues = combinedMatrixLyapunovFunction /@ samplingGrid; sampledExactLyapunovDerivatives = exactLyapunovDerivative /@ samplingGrid; sampledDirectLyapunovDerivatives = directLyapunovDerivative /@ samplingGrid;
  maximumLyapunovRepresentationError = Max@Abs[sampledLyapunovValues - sampledCombinedMatrixLyapunovValues]; maximumLyapunovDerivativeError = Max@Abs[sampledDirectLyapunovDerivatives - sampledExactLyapunovDerivatives];
  sampledLyapunovIncrements = Differences[sampledLyapunovValues]; largestSampledLyapunovIncrement = Max[sampledLyapunovIncrements]; monotonicityTolerance = 10.^-8;
  numberOfSignificantPositiveIncrements = Count[sampledLyapunovIncrements, x_ /; x > monotonicityTolerance]; largestSampledExactDerivative = Max[sampledExactLyapunovDerivatives];
  initialStateErrorValue = stateErrorVector[initialTime]; finalStateErrorValue = stateErrorVector[finalTime]; initialParameterErrorValue = parameterErrorVector[initialTime]; finalParameterErrorValue = parameterErrorVector[finalTime];
  initialLyapunovValue = lyapunovFunction[initialTime]; finalLyapunovValue = lyapunovFunction[finalTime]; initialStateErrorNorm = Norm[initialStateErrorValue]; finalStateErrorNorm = Norm[finalStateErrorValue];
  initialParameterErrorNorm = Norm[initialParameterErrorValue]; finalParameterErrorNorm = Norm[finalParameterErrorValue];
  Scan[Print, {"------------------------------------------------------------", "Chapter 3, Lesson 2: Mathematica Simulation Results",
    "P symmetry residual = " <> ToString[N[pSymmetryResidual, 12]], "Eigenvalues of P = " <> ToString[N[pEigenvalues, 12]], "Eigenvalues of Gamma = " <> ToString[N[gammaEigenvalues, 12]],
    "Combined-matrix eigenvalues = " <> ToString[N[combinedMatrixEigenvalues, 12]], "Lyapunov residual norm = " <> ToString[N[lyapunovEquationResidualNorm, 12]], "Sampled points = " <> ToString[numberOfSamplingPoints],
    "Initial state error = " <> ToString[N[initialStateErrorValue, 10]], "Final state error = " <> ToString[N[finalStateErrorValue, 10]], "Initial parameter error = " <> ToString[N[initialParameterErrorValue, 10]],
    "Final parameter error = " <> ToString[N[finalParameterErrorValue, 10]], "Initial state-error norm = " <> ToString[N[initialStateErrorNorm, 10]], "Final state-error norm = " <> ToString[N[finalStateErrorNorm, 10]],
    "Initial parameter-error norm = " <> ToString[N[initialParameterErrorNorm, 10]], "Final parameter-error norm = " <> ToString[N[finalParameterErrorNorm, 10]], "V(0) = " <> ToString[N[initialLyapunovValue, 12]],
    "V(30) = " <> ToString[N[finalLyapunovValue, 12]], "Largest V increment = " <> ToString[ScientificForm[largestSampledLyapunovIncrement, 8]], "Positive V increments = " <> ToString[numberOfSignificantPositiveIncrements],
    "Largest exact VDot = " <> ToString[ScientificForm[largestSampledExactDerivative, 8]], "Maximum V representation error = " <> ToString[ScientificForm[maximumLyapunovRepresentationError, 8]],
    "Maximum VDot error = " <> ToString[ScientificForm[maximumLyapunovDerivativeError, 8]], "------------------------------------------------------------"}];
  exportHeader = {"t", "e1", "e2", "thetaTilde1", "thetaTilde2", "stateLyapunovTerm", "parameterLyapunovTerm", "V", "VdotExact", "VdotDirect"};
  exportRows = MapThread[Join[{ #1}, #2, #3, { #4, #5, #6, #7, #8}] &, {samplingGrid, sampledStateErrors, sampledParameterErrors, sampledStateLyapunovTerms, sampledParameterLyapunovTerms, sampledLyapunovValues, sampledExactLyapunovDerivatives, sampledDirectLyapunovDerivatives}];
  csvOutputFile = FileNameJoin[{NotebookDirectory[], "Chapter3_Lesson2_results_mathematica.csv"}]; Export[csvOutputFile, Prepend[exportRows, exportHeader], "CSV"]; Print["Simulation data exported to: ", csvOutputFile];
  stateErrorPlot = Plot[{stateError1[t], stateError2[t]}, {t, initialTime, finalTime}, PlotRange -> All, PlotLegends -> {"e1(t)", "e2(t)"}, AxesLabel -> {"Time", "State error"}, PlotLabel -> "Evolution of State-Error Components", GridLines -> Automatic, ImageSize -> Large];
  parameterErrorPlot = Plot[{parameterError1[t], parameterError2[t]}, {t, initialTime, finalTime}, PlotRange -> All, PlotLegends -> {"thetaTilde1(t)", "thetaTilde2(t)"}, AxesLabel -> {"Time", "Parameter error"}, PlotLabel -> "Evolution of Parameter-Error Components", GridLines -> Automatic, ImageSize -> Large];
  lyapunovFunctionPlot = Plot[lyapunovFunction[t], {t, initialTime, finalTime}, PlotRange -> All, AxesLabel -> {"Time", "V(t)"}, PlotLabel -> "Combined Quadratic Lyapunov Function", GridLines -> Automatic, ImageSize -> Large];
  lyapunovComponentsPlot = Plot[{stateLyapunovTerm[t], parameterLyapunovTerm[t], lyapunovFunction[t]}, {t, initialTime, finalTime}, PlotRange -> All, PlotLegends -> {"State contribution", "Parameter contribution", "Total V"}, AxesLabel -> {"Time", "Lyapunov energy"}, PlotLabel -> "State and Parameter Contributions to V", GridLines -> Automatic, ImageSize -> Large];
  lyapunovDerivativePlot = Plot[{exactLyapunovDerivative[t], directLyapunovDerivative[t]}, {t, initialTime, finalTime}, PlotRange -> All, PlotLegends -> {"Exact: -e^T Q e", "Direct differentiation"}, AxesLabel -> {"Time", "VDot(t)"}, PlotLabel -> "Analytical and Direct Lyapunov Derivatives", GridLines -> Automatic, ImageSize -> Large];
  stateErrorPhasePlot = ParametricPlot[{stateError1[t], stateError2[t]}, {t, initialTime, finalTime}, PlotRange -> All, AxesLabel -> {"e1", "e2"}, PlotLabel -> "State-Error Phase Trajectory", GridLines -> Automatic, ImageSize -> Large];
  regressorPlot = Plot[{phi[t][[1]], phi[t][[2]]}, {t, initialTime, finalTime}, PlotRange -> All, PlotLegends -> {"phi1(t)", "phi2(t)"}, AxesLabel -> {"Time", "Regressor"}, PlotLabel -> "Exponentially Decaying Regressor Components", GridLines -> Automatic, ImageSize -> Large];
  simulationGraphics = Column[{stateErrorPlot, parameterErrorPlot, lyapunovFunctionPlot, lyapunovComponentsPlot, lyapunovDerivativePlot, stateErrorPhasePlot, regressorPlot}, Spacings -> 2];
  Export[FileNameJoin[{NotebookDirectory[], "Chapter3_Lesson2_state_error.png"}], stateErrorPlot, ImageResolution -> 200]; Export[FileNameJoin[{NotebookDirectory[], "Chapter3_Lesson2_parameter_error.png"}], parameterErrorPlot, ImageResolution -> 200];
  Export[FileNameJoin[{NotebookDirectory[], "Chapter3_Lesson2_lyapunov_function.png"}], lyapunovFunctionPlot, ImageResolution -> 200]; Export[FileNameJoin[{NotebookDirectory[], "Chapter3_Lesson2_lyapunov_derivative.png"}], lyapunovDerivativePlot, ImageResolution -> 200];
  stabilityInterpretation = Column[{Style["Interpretation", Bold, 16], "1. P is symmetric and positive definite.", "2. Gamma is positive definite.", "3. The combined Lyapunov matrix is positive definite.",
    "4. The adaptive law cancels the state–parameter cross term.", "5. VDot = -e^T Q e is nonpositive.", "6. State and parameter errors remain bounded.", "7. The state error approaches zero.",
    "8. Parameter convergence is not guaranteed because the decaying regressor is not persistently exciting."}, Spacings -> 1];
  Column[{simulationGraphics, stabilityInterpretation}]
 ], "Input"]
}, WindowTitle -> "Chapter3_Lesson2", Saveable -> True, Editable -> True, Magnification -> 1., StyleDefinitions -> "Default.nb"]        

• Wrong parameter-error sign: changing the definition of \( \tilde{\boldsymbol{\theta} } \) without changing the adaptive-law sign destroys cancellation.
• Missing transpose: the dimensionally valid update is \( -\boldsymbol{\Gamma}\mathbf{W}^{T}\mathbf{P}\mathbf{e} \).
• Using a non-Hurwitz state matrix: a positive-definite solution of the required Lyapunov equation is then not guaranteed.
• Claiming parameter convergence: \( \dot V\leq0 \) proves bounded parameter error, not generally zero parameter error.
• Ignoring time variation of ideal parameters: this introduces an uncancelled derivative term.
• Confusing continuous and discrete behavior: forward Euler can produce sampled increases in \( V \) even when the continuous derivative is nonpositive.

Problem 1 (Quadratic bounds): Prove that the joint candidate satisfies bounds of the form \( \alpha_{1}\lVert\boldsymbol{\xi}\rVert^{2} \leq V\leq\alpha_{2}\lVert\boldsymbol{\xi}\rVert^{2} \) and state suitable values of \( \alpha_{1} \) and \( \alpha_{2} \).

Solution: Apply the Rayleigh-quotient inequality separately to the state and parameter terms:

\[ \begin{aligned} V&\geq\lambda_{\min}(\mathbf{P})\lVert\mathbf{e}\rVert^{2} +\lambda_{\min}(\boldsymbol{\Gamma}^{-1}) \lVert\tilde{\boldsymbol{\theta} }\rVert^{2}\\ &\geq\min\left\{\lambda_{\min}(\mathbf{P}), \lambda_{\min}(\boldsymbol{\Gamma}^{-1})\right\} \left(\lVert\mathbf{e}\rVert^{2} +\lVert\tilde{\boldsymbol{\theta} }\rVert^{2}\right). \end{aligned} \]

The upper bound follows similarly using the maximum eigenvalues. Therefore the values in Section 5 are valid.

Problem 2 (Derive the adaptive law): For \( \dot{\mathbf{e} }=\mathbf{A}\mathbf{e} +\mathbf{b}\boldsymbol{\phi}^{T}\tilde{\boldsymbol{\theta} } \) and

\[ V=\frac{1}{2}\mathbf{e}^{T}\mathbf{P}\mathbf{e} +\frac{1}{2}\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta} }, \]

derive an update law that eliminates the mixed term.

Solution: Differentiation gives

\[ \dot V=-\frac{1}{2}\mathbf{e}^{T}\mathbf{Q}\mathbf{e} +\mathbf{e}^{T}\mathbf{P}\mathbf{b} \boldsymbol{\phi}^{T}\tilde{\boldsymbol{\theta} } +\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\tilde{\boldsymbol{\theta} } }. \]

Choose

\[ \dot{\tilde{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi} \mathbf{b}^{T}\mathbf{P}\mathbf{e}. \]

The last two scalar terms cancel, leaving \( \dot V=-(1/2)\mathbf{e}^{T}\mathbf{Q}\mathbf{e} \).

Problem 3 (Solve the Lyapunov equation): For the matrix \( \mathbf{A} \) in Section 11 and \( \mathbf{Q}=\mathbf{I}_{2} \), let

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

Find \( \mathbf{P} \) and verify positive definiteness.

Solution: Substitution into the Lyapunov equation gives

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

Hence

\[ p_{12}=\frac{1}{4},\qquad p_{22}=\frac{1}{4},\qquad p_{11}=\frac{5}{4}. \]

The leading principal minors are \( 5/4 \) and \( (5/4)(1/4)-(1/4)^{2}=1/4 \), both positive; therefore \( \mathbf{P}\succ0 \).

Problem 4 (No parameter-convergence conclusion): Suppose \( \mathbf{W}(t)=\mathbf{0} \) for all \( t\geq T_{0} \). Determine the parameter-error behavior after \( T_{0} \) and explain why the Lyapunov result remains valid.

Solution: The adaptive law becomes

\[ \dot{\tilde{\boldsymbol{\theta} } }(t)=\mathbf{0},\qquad t\geq T_{0}. \]

Thus \( \tilde{\boldsymbol{\theta} }(t) =\tilde{\boldsymbol{\theta} }(T_{0}) \). Meanwhile the state error satisfies \( \dot{\mathbf{e} }=\mathbf{A}\mathbf{e} \), so its nominal stable dynamics continue. The parameter contribution to \( V \) becomes constant and the state contribution decreases; hence nonincrease of \( V \) is fully consistent with a nonzero frozen parameter error.

Problem 5 (Time-varying ideal parameter): Let \( \boldsymbol{\theta}^{*}=\boldsymbol{\theta}^{*}(t) \). Derive the uncancelled term in \( \dot V \) when the same estimate update law is retained.

Solution: Since

\[ \dot{\tilde{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\mathbf{W}^{T}\mathbf{P}\mathbf{e} -\dot{\boldsymbol{\theta} }^{*}, \]

substitution gives

\[ \dot V=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e} -2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\boldsymbol{\theta} }^{*}. \]

The added term can be positive or negative, so the original theorem cannot be applied without an additional bound or a modified adaptive design.

Problem 6 (Schur-complement condition): Assume \( \mathbf{P}\succ0 \). State and justify the condition on \( \mathbf{R} \) and \( \mathbf{S} \) that makes \( \begin{bmatrix}\mathbf{P}&\mathbf{S}\\ \mathbf{S}^{T}&\mathbf{R}\end{bmatrix} \) positive definite.

Solution: Complete the square:

\[ \begin{aligned} \begin{bmatrix}\mathbf{x}\\\mathbf{y}\end{bmatrix}^{T} \begin{bmatrix}\mathbf{P}&\mathbf{S}\\ \mathbf{S}^{T}&\mathbf{R}\end{bmatrix} \begin{bmatrix}\mathbf{x}\\\mathbf{y}\end{bmatrix} &=\left(\mathbf{x}+\mathbf{P}^{-1}\mathbf{S}\mathbf{y}\right)^{T} \mathbf{P} \left(\mathbf{x}+\mathbf{P}^{-1}\mathbf{S}\mathbf{y}\right)\\ &\quad+\mathbf{y}^{T} \left(\mathbf{R}-\mathbf{S}^{T}\mathbf{P}^{-1}\mathbf{S}\right) \mathbf{y}. \end{aligned} \]

Therefore the full matrix is positive definite exactly when \( \mathbf{R}-\mathbf{S}^{T}\mathbf{P}^{-1}\mathbf{S}\succ0 \).

A combined adaptive error system is naturally analyzed with a quadratic function containing a state-error term and a weighted parameter-error term. The state metric comes from a continuous-time Lyapunov equation, while the inverse adaptation gain supplies the parameter metric. Differentiation exposes a mixed term, and the adaptive update law is selected to cancel that term exactly. The resulting identity proves stability, boundedness, forward invariance of quadratic level sets, and finite tracking-error energy. It does not, without further excitation and convergence arguments, prove that the parameter estimate approaches the ideal parameter.