After completing this lesson, students should be able to:

• state the common tracking objective of continuous-time and discrete-time MRAC;
• distinguish differential adaptation laws from recursive sample-to-sample updates;
• derive continuous-time and discrete-time tracking-error equations under matching;
• explain why continuous-time Lyapunov cancellation does not automatically survive discretization;
• map continuous poles to discrete poles under exact zero-order-hold sampling;
• implement a first-order comparison in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The lesson uses only the reference-model, tracking-error, parameter-error, and direct/indirect MRAC concepts introduced in the preceding lessons. Formal parameter-convergence conditions are deferred to Chapter 10.

In both settings, the plant output or state should follow a stable reference model despite unknown constant parameters. The difference is the mathematical clock:

\[ \begin{aligned} \text{continuous time:}\quad & t\in\mathbb{R}_{\geq 0},\qquad \dot{x}(t)=f(x(t),u(t),\theta), \\ \text{discrete time:}\quad & k\in\mathbb{Z}_{\geq 0},\qquad x[k+1]=F(x[k],u[k],\theta). \end{aligned} \]

A continuous adaptive parameter \( \hat{\theta}(t) \) evolves through a differential equation. A discrete adaptive parameter \( \hat{\theta}[k] \) changes through a recursion. The associated stability tools are therefore different:

\[ \dot{V}(t)\leq 0 \qquad\text{versus}\qquad \Delta V[k]=V[k+1]-V[k]\leq 0. \]

Both formulations contain a reference model, an adjustable controller, a tracking error, and an adaptation mechanism. Their signal timing is different, but their conceptual structure is the same.

flowchart TD
  R["Command r"] --> RM["Reference model"]
  R --> C["Adaptive controller"]
  P["Plant output or state"] --> C
  C --> U["Control input u"]
  U --> PL["Uncertain plant"]
  PL --> P
  RM --> E["Tracking error"]
  P --> E
  E --> A["Adaptation mechanism"]
  P --> A
  R --> A
  A --> C
  CT["Continuous case: \ndifferential update"] -.-> A
  DT["Discrete case: \nrecursive update \nafter sampling"] -.-> A
        

The diagram is intentionally architecture-level. The exact error signal supplied to the discrete adaptation mechanism may be filtered or augmented so that it is computable and compatible with a discrete stability proof.

Consider the state-space plant and reference model

\[ \dot{x}=Ax+b\,u,\qquad \dot{x}_m=A_mx_m+b_mr, \]

where \( A_m \) is Hurwitz. Choose a direct adaptive control law

\[ u=\hat{\theta}^T\omega,\qquad \omega=\begin{bmatrix}x^T & r\end{bmatrix}^T. \]

Suppose ideal parameters \( \theta^* \) satisfy the matching equations. With \( \tilde{\theta}=\hat{\theta}-\theta^* \) and \( e=x-x_m \), the closed-loop tracking-error model can be written as

\[ \dot{e}=A_me+b\,\tilde{\theta}^T\omega. \]

Let \( Q=Q^T>0 \) and let \( P=P^T>0 \) solve

\[ A_m^TP+PA_m=-Q. \]

A standard Lyapunov candidate is

\[ V=e^TPe+\tilde{\theta}^T\Gamma^{-1}\tilde{\theta}, \qquad \Gamma=\Gamma^T>0. \]

Differentiation gives

\[ \dot{V} =-e^TQe +2e^TPb\,\tilde{\theta}^T\omega +2\tilde{\theta}^T\Gamma^{-1}\dot{\hat{\theta}}. \]

Selecting the continuous adaptation law

\[ \dot{\hat{\theta}} =-\Gamma\omega\,e^TPb \]

cancels the mixed term exactly, yielding

\[ \dot{V}=-e^TQe\leq 0. \]

This cancellation is algebraically tied to the product rule for derivatives. That fact becomes important when moving to discrete time.

For a sampled plant and a discrete reference model, write

\[ x[k+1]=A_dx[k]+b_du[k],\qquad x_m[k+1]=A_{m,d}x_m[k]+b_{m,d}r[k]. \]

With \( u[k]=\hat{\theta}^T[k]\omega[k] \) and exact discrete matching, the tracking-error equation becomes

\[ e[k+1] =A_{m,d}e[k]+b_d\tilde{\theta}^T[k]\omega[k]. \]

The reference matrix must be Schur stable:

\[ \rho(A_{m,d})<1, \]

where \( \rho(\cdot) \) is the spectral radius. The corresponding discrete Lyapunov equation is

\[ A_{m,d}^TPA_{m,d}-P=-Q, \qquad P=P^T>0,\quad Q=Q^T>0. \]

Stability analysis must evaluate the complete difference \( V[k+1]-V[k] \). There is no differential product rule that removes all second-order increments automatically.

If a continuous plant is sampled every \( h>0 \) seconds and the control input is held constant between samples, exact zero-order-hold discretization gives

\[ A_d=e^{Ah},\qquad b_d=\int_0^h e^{A\tau}b\,d\tau. \]

A continuous eigenvalue \( \lambda \) maps to the discrete eigenvalue

\[ z=e^{\lambda h}. \]

Therefore \( \operatorname{Re}(\lambda)<0 \) implies \( |z|<1 \) under exact sampling. By contrast, forward Euler uses

\[ z_{\mathrm{FE}}=1+h\lambda. \]

For a real pole \( \lambda<0 \), forward-Euler stability requires

\[ |1+h\lambda|<1 \quad\Longleftrightarrow\quad 0<h<-\frac{2}{\lambda}. \]

Thus a numerically convenient step size is also a control-design variable: it affects model poles, sample delay, computational burden, and the effective speed of parameter updates.

Suppose a continuous update is approximated by forward Euler:

\[ \hat{\theta}[k+1] =\hat{\theta}[k] -h\Gamma\omega[k]e^T[k]Pb. \]

Even if this recursion resembles the continuous law, its Lyapunov difference contains quadratic increment terms. For a generic vector \( z[k+1]=z[k]+\delta z[k] \),

\[ z^T[k+1]Mz[k+1]-z^T[k]Mz[k] =2z^T[k]M\delta z[k]+\delta z^T[k]M\delta z[k]. \]

The final term is nonnegative when \( M>0 \). It has no direct continuous-time analogue because it is second order in the increment. Consequently, cancellation of the first-order cross term does not by itself prove \( \Delta V[k]\leq 0 \).

A valid discrete design therefore uses a recursion, normalization, error filtering, or augmented error specifically constructed for the sample-to-sample dynamics.

To expose the discrete mechanism without introducing the full higher-order MRAC machinery, consider a computable one-step regression error

\[ \varepsilon[k+1] =\tilde{\theta}^T[k]\phi[k]. \]

Use the normalized update

\[ \hat{\theta}[k+1] =\hat{\theta}[k] -\frac{\Gamma\phi[k]\varepsilon[k+1]} {1+\phi^T[k]\Gamma\phi[k]}. \]

Define \( s[k]=\phi^T[k]\Gamma\phi[k]\geq 0 \) and the parameter energy

\[ W[k]=\tilde{\theta}^T[k]\Gamma^{-1}\tilde{\theta}[k]. \]

Substituting the recursion gives the exact difference

\[ \begin{aligned} W[k+1]-W[k] &=-\frac{2\varepsilon^2[k+1]}{1+s[k]} +\frac{s[k]\varepsilon^2[k+1]}{(1+s[k])^2} \\ &=-\frac{(2+s[k])\varepsilon^2[k+1]}{(1+s[k])^2} \leq 0. \end{aligned} \]

The normalization denominator is not merely a numerical convenience. It limits the update when the regressor is large and makes the descent inequality hold without requiring an arbitrarily small fixed step.

Consider

\[ \dot{y}=-ay+bu,\qquad \dot{y}_m=-a_my_m+b_mr, \]

with controller

\[ u=\hat{\theta}_y y+\hat{\theta}_r r. \]

The continuous-time ideal gains are

\[ \theta_{y,c}^*=\frac{a-a_m}{b}, \qquad \theta_{r,c}^*=\frac{b_m}{b}. \]

For \( a_m>0 \), a scalar Lyapunov choice leads to

\[ \dot{\hat{\theta}} =-\gamma\,\operatorname{sgn}(b) \begin{bmatrix}y\\r\end{bmatrix}e, \qquad e=y-y_m. \]

Under exact zero-order-hold sampling,

\[ a_d=e^{-ah},\qquad b_d=\frac{b}{a}\left(1-e^{-ah}\right), \]

\[ a_{m,d}=e^{-a_mh},\qquad b_{m,d}=\frac{b_m}{a_m} \left(1-e^{-a_mh}\right). \]

The exact discrete matching gains are

\[ \theta_{y,d}^* =\frac{a_{m,d}-a_d}{b_d}, \qquad \theta_{r,d}^* =\frac{b_{m,d}}{b_d}. \]

These discrete ideal gains generally differ from the continuous ideal gains for finite \( h \), although they approach them as \( h→0 \). The sample-to-sample error equation is

\[ e[k+1] =a_{m,d}e[k] +b_d\tilde{\theta}^T[k]\phi[k], \qquad \phi[k]=\begin{bmatrix}y[k]\\r[k]\end{bmatrix}. \]

Therefore the one-step regression error is directly recovered as

\[ \varepsilon[k+1] =\frac{e[k+1]-a_{m,d}e[k]}{b_d} =\tilde{\theta}^T[k]\phi[k]. \]

flowchart TD
  S0["At sample k: read y[k] and command r[k]"] --> S1["Compute regressor phi[k]"]
  S1 --> S2["Compute control u[k] from current parameters"]
  S2 --> S3["Hold u[k] over one sampling interval"]
  S3 --> S4["Plant and reference model advance to k+1"]
  S4 --> S5["Read y[k+1] and compute e[k+1]"]
  S5 --> S6["Form one-step error epsilon[k+1]"]
  S6 --> S7["Normalize and update parameters"]
  S7 --> S8["Use updated parameters at sample k+1"]
        

The update uses information that becomes available after the plant has evolved through one sample. This timing prevents the controller from using a future measurement. In an embedded implementation, computation time must fit inside the sample period or be represented as an additional delay.

1. Choose a stable reference model in the same time domain as the proof.
2. Verify continuous or discrete matching conditions before defining parameter errors.
3. Use exact sampled coefficients when the zero-order-hold model is known.
4. Define the order in which measurement, control, model propagation, and adaptation occur.
5. Do not reuse a continuous Lyapunov proof after discretization without bounding the extra increment terms.
6. Scale or normalize discrete updates when large regressors can cause large parameter jumps.
7. Compare the sample period with plant and reference-model time constants.
8. Log tracking error, control effort, parameter estimates, and the one-step adaptation error.

The Python program simulates the continuous law with a small integration step and the discrete law with exact zero-order-hold coefficients. It writes CSV files and plots both trajectories.

Chapter4_Lesson4.py

# Chapter4_Lesson4.py
"""
Continuous-time and discrete-time MRAC comparison for a scalar first-order plant.

Continuous plant:
    y_dot = -a*y + b*u
Reference model:
    ym_dot = -am*ym + bm*r
Controller:
    u = theta_y*y + theta_r*r

The continuous adaptive law is Lyapunov-based. The discrete implementation uses
the exact zero-order-hold plant/model coefficients and a normalized one-step
error update. This is an educational comparison, not a robustness-certified
implementation.
"""

from __future__ import annotations

import csv
import math
from dataclasses import dataclass
from pathlib import Path
from typing import Callable

import matplotlib.pyplot as plt


@dataclass(frozen=True)
class Parameters:
    a: float = 1.6
    b: float = 1.2
    am: float = 1.0
    bm: float = 1.0
    gamma_ct: float = 4.0
    gamma_dt: float = 0.35
    t_final: float = 20.0
    dt: float = 0.001
    sample_time: float = 0.05


def reference(t: float) -> float:
    """Bounded command containing steps and a low-frequency sinusoid."""
    if t < 5.0:
        return 1.0
    if t < 10.0:
        return -0.5
    return 0.75 + 0.25 * math.sin(0.7 * t)


def sign_nonzero(value: float) -> float:
    if value == 0.0:
        raise ValueError("The input gain b must be nonzero.")
    return 1.0 if value > 0.0 else -1.0


def simulate_continuous(p: Parameters) -> dict[str, list[float]]:
    """Forward-Euler integration with a small step for transparent pedagogy."""
    n = int(round(p.t_final / p.dt))
    y = 0.0
    ym = 0.0
    theta_y = 0.0
    theta_r = 0.0
    sign_b = sign_nonzero(p.b)

    data = {key: [] for key in ("t", "r", "y", "ym", "e", "u", "theta_y", "theta_r")}

    for i in range(n + 1):
        t = i * p.dt
        r = reference(t)
        e = y - ym
        u = theta_y * y + theta_r * r

        for key, value in (
            ("t", t),
            ("r", r),
            ("y", y),
            ("ym", ym),
            ("e", e),
            ("u", u),
            ("theta_y", theta_y),
            ("theta_r", theta_r),
        ):
            data[key].append(value)

        y_dot = -p.a * y + p.b * u
        ym_dot = -p.am * ym + p.bm * r
        theta_y_dot = -p.gamma_ct * sign_b * y * e
        theta_r_dot = -p.gamma_ct * sign_b * r * e

        y += p.dt * y_dot
        ym += p.dt * ym_dot
        theta_y += p.dt * theta_y_dot
        theta_r += p.dt * theta_r_dot

    return data


def exact_scalar_zoh(decay_rate: float, input_gain: float, h: float) -> tuple[float, float]:
    """Return x[k+1] = ad*x[k] + bd*u[k] for x_dot=-decay_rate*x+input_gain*u."""
    if decay_rate <= 0.0 or h <= 0.0:
        raise ValueError("decay_rate and sample time must be positive.")
    ad = math.exp(-decay_rate * h)
    bd = input_gain * (1.0 - ad) / decay_rate
    return ad, bd


def simulate_discrete(p: Parameters) -> dict[str, list[float]]:
    """
    Exact-ZOH scalar MRAC with a normalized one-step update.

    Let e[k]=y[k]-ym[k] and phi[k]=[y[k], r[k]]^T. Under matching,
        (e[k+1]-amd*e[k]) / bd = tilde_theta[k]^T phi[k].
    The measured quantity on the left is used after y[k+1] becomes available.
    """
    ad, bd = exact_scalar_zoh(p.a, p.b, p.sample_time)
    amd, bmd = exact_scalar_zoh(p.am, p.bm, p.sample_time)

    n = int(round(p.t_final / p.sample_time))
    y = 0.0
    ym = 0.0
    theta_y = 0.0
    theta_r = 0.0

    data = {
        key: []
        for key in ("t", "r", "y", "ym", "e", "u", "theta_y", "theta_r", "epsilon")
    }

    for k in range(n + 1):
        t = k * p.sample_time
        r = reference(t)
        e = y - ym
        u = theta_y * y + theta_r * r

        epsilon = 0.0
        if k < n:
            y_next = ad * y + bd * u
            ym_next = amd * ym + bmd * r
            e_next = y_next - ym_next

            epsilon = (e_next - amd * e) / bd
            phi_norm_sq = y * y + r * r
            normalizer = 1.0 + p.gamma_dt * phi_norm_sq

            theta_y_next = theta_y - p.gamma_dt * y * epsilon / normalizer
            theta_r_next = theta_r - p.gamma_dt * r * epsilon / normalizer
        else:
            y_next = y
            ym_next = ym
            theta_y_next = theta_y
            theta_r_next = theta_r

        for key, value in (
            ("t", t),
            ("r", r),
            ("y", y),
            ("ym", ym),
            ("e", e),
            ("u", u),
            ("theta_y", theta_y),
            ("theta_r", theta_r),
            ("epsilon", epsilon),
        ):
            data[key].append(value)

        y, ym = y_next, ym_next
        theta_y, theta_r = theta_y_next, theta_r_next

    return data


def write_csv(path: Path, data: dict[str, list[float]]) -> None:
    keys = list(data)
    rows = zip(*(data[key] for key in keys))
    with path.open("w", newline="", encoding="utf-8") as handle:
        writer = csv.writer(handle)
        writer.writerow(keys)
        writer.writerows(rows)


def main() -> None:
    p = Parameters()
    ct = simulate_continuous(p)
    dt = simulate_discrete(p)

    output_dir = Path(__file__).resolve().parent
    write_csv(output_dir / "Chapter4_Lesson4_continuous.csv", ct)
    write_csv(output_dir / "Chapter4_Lesson4_discrete.csv", dt)

    ad, bd = exact_scalar_zoh(p.a, p.b, p.sample_time)
    amd, bmd = exact_scalar_zoh(p.am, p.bm, p.sample_time)
    ideal_theta_y_ct = (p.a - p.am) / p.b
    ideal_theta_r_ct = p.bm / p.b
    ideal_theta_y_dt = (amd - ad) / bd
    ideal_theta_r_dt = bmd / bd
    print(
        f"Continuous ideal parameters: theta_y*={ideal_theta_y_ct:.6f}, "
        f"theta_r*={ideal_theta_r_ct:.6f}"
    )
    print(
        f"Discrete ideal parameters: theta_y*={ideal_theta_y_dt:.6f}, "
        f"theta_r*={ideal_theta_r_dt:.6f}"
    )
    print(
        "Final CT parameters:",
        f"theta_y={ct['theta_y'][-1]:.6f}, theta_r={ct['theta_r'][-1]:.6f}",
    )
    print(
        "Final DT parameters:",
        f"theta_y={dt['theta_y'][-1]:.6f}, theta_r={dt['theta_r'][-1]:.6f}",
    )

    plt.figure()
    plt.plot(ct["t"], ct["ym"], label="continuous reference model")
    plt.plot(ct["t"], ct["y"], label="continuous plant")
    plt.plot(dt["t"], dt["ym"], "o", markevery=12, label="discrete reference model")
    plt.plot(dt["t"], dt["y"], ".", markevery=12, label="discrete plant")
    plt.xlabel("time (s)")
    plt.ylabel("output")
    plt.title("Continuous-time and discrete-time MRAC comparison")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

This standard C++17 implementation has no external numerical dependencies. It exports the continuous and discrete simulations as CSV.

Chapter4_Lesson4.cpp

// Chapter4_Lesson4.cpp
// Continuous-time and discrete-time MRAC comparison for a scalar first-order plant.
// Build: g++ -std=c++17 -O2 Chapter4_Lesson4.cpp -o Chapter4_Lesson4

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

struct Parameters {
    double a = 1.6;
    double b = 1.2;
    double am = 1.0;
    double bm = 1.0;
    double gamma_ct = 4.0;
    double gamma_dt = 0.35;
    double t_final = 20.0;
    double dt = 0.001;
    double sample_time = 0.05;
};

struct Record {
    double t{};
    double r{};
    double y{};
    double ym{};
    double e{};
    double u{};
    double theta_y{};
    double theta_r{};
    double epsilon{};
};

double reference_signal(double t) {
    if (t < 5.0) {
        return 1.0;
    }
    if (t < 10.0) {
        return -0.5;
    }
    return 0.75 + 0.25 * std::sin(0.7 * t);
}

double sign_nonzero(double value) {
    if (value == 0.0) {
        throw std::invalid_argument("The input gain b must be nonzero.");
    }
    return value > 0.0 ? 1.0 : -1.0;
}

std::pair<double, double> exact_scalar_zoh(
    double decay_rate,
    double input_gain,
    double h
) {
    if (decay_rate <= 0.0 || h <= 0.0) {
        throw std::invalid_argument("decay_rate and sample time must be positive.");
    }
    const double ad = std::exp(-decay_rate * h);
    const double bd = input_gain * (1.0 - ad) / decay_rate;
    return {ad, bd};
}

std::vector<Record> simulate_continuous(const Parameters& p) {
    const int n = static_cast<int>(std::llround(p.t_final / p.dt));
    const double sign_b = sign_nonzero(p.b);

    double y = 0.0;
    double ym = 0.0;
    double theta_y = 0.0;
    double theta_r = 0.0;

    std::vector<Record> data;
    data.reserve(static_cast<std::size_t>(n + 1));

    for (int i = 0; i <= n; ++i) {
        const double t = i * p.dt;
        const double r = reference_signal(t);
        const double e = y - ym;
        const double u = theta_y * y + theta_r * r;

        data.push_back({t, r, y, ym, e, u, theta_y, theta_r, 0.0});

        const double y_dot = -p.a * y + p.b * u;
        const double ym_dot = -p.am * ym + p.bm * r;
        const double theta_y_dot = -p.gamma_ct * sign_b * y * e;
        const double theta_r_dot = -p.gamma_ct * sign_b * r * e;

        y += p.dt * y_dot;
        ym += p.dt * ym_dot;
        theta_y += p.dt * theta_y_dot;
        theta_r += p.dt * theta_r_dot;
    }

    return data;
}

std::vector<Record> simulate_discrete(const Parameters& p) {
    const auto [ad, bd] = exact_scalar_zoh(p.a, p.b, p.sample_time);
    const auto [amd, bmd] = exact_scalar_zoh(p.am, p.bm, p.sample_time);
    const int n = static_cast<int>(std::llround(p.t_final / p.sample_time));

    double y = 0.0;
    double ym = 0.0;
    double theta_y = 0.0;
    double theta_r = 0.0;

    std::vector<Record> data;
    data.reserve(static_cast<std::size_t>(n + 1));

    for (int k = 0; k <= n; ++k) {
        const double t = k * p.sample_time;
        const double r = reference_signal(t);
        const double e = y - ym;
        const double u = theta_y * y + theta_r * r;

        double epsilon = 0.0;
        double y_next = y;
        double ym_next = ym;
        double theta_y_next = theta_y;
        double theta_r_next = theta_r;

        if (k < n) {
            y_next = ad * y + bd * u;
            ym_next = amd * ym + bmd * r;
            const double e_next = y_next - ym_next;

            // Under exact matching:
            // epsilon[k+1] = (e[k+1] - amd*e[k]) / bd = tilde_theta[k]^T*phi[k].
            epsilon = (e_next - amd * e) / bd;
            const double phi_norm_sq = y * y + r * r;
            const double normalizer = 1.0 + p.gamma_dt * phi_norm_sq;

            theta_y_next = theta_y - p.gamma_dt * y * epsilon / normalizer;
            theta_r_next = theta_r - p.gamma_dt * r * epsilon / normalizer;
        }

        data.push_back({t, r, y, ym, e, u, theta_y, theta_r, epsilon});

        y = y_next;
        ym = ym_next;
        theta_y = theta_y_next;
        theta_r = theta_r_next;
    }

    return data;
}

void write_csv(const std::string& filename, const std::vector<Record>& data) {
    std::ofstream file(filename);
    if (!file) {
        throw std::runtime_error("Cannot open output file: " + filename);
    }

    file << "t,r,y,ym,e,u,theta_y,theta_r,epsilon\n";
    file << std::setprecision(12);

    for (const auto& row : data) {
        file << row.t << ',' << row.r << ',' << row.y << ',' << row.ym << ','
             << row.e << ',' << row.u << ',' << row.theta_y << ','
             << row.theta_r << ',' << row.epsilon << '\n';
    }
}

int main() {
    try {
        const Parameters p;
        const auto continuous = simulate_continuous(p);
        const auto discrete = simulate_discrete(p);

        write_csv("Chapter4_Lesson4_continuous_cpp.csv", continuous);
        write_csv("Chapter4_Lesson4_discrete_cpp.csv", discrete);

        const auto [ad, bd] = exact_scalar_zoh(p.a, p.b, p.sample_time);
        const auto [amd, bmd] = exact_scalar_zoh(p.am, p.bm, p.sample_time);
        const double ideal_theta_y_ct = (p.a - p.am) / p.b;
        const double ideal_theta_r_ct = p.bm / p.b;
        const double ideal_theta_y_dt = (amd - ad) / bd;
        const double ideal_theta_r_dt = bmd / bd;

        std::cout << std::fixed << std::setprecision(6);
        std::cout << "Continuous ideal parameters: theta_y*=" << ideal_theta_y_ct
                  << ", theta_r*=" << ideal_theta_r_ct << '\n';
        std::cout << "Discrete ideal parameters: theta_y*=" << ideal_theta_y_dt
                  << ", theta_r*=" << ideal_theta_r_dt << '\n';
        std::cout << "Final CT parameters: theta_y="
                  << continuous.back().theta_y
                  << ", theta_r=" << continuous.back().theta_r << '\n';
        std::cout << "Final DT parameters: theta_y="
                  << discrete.back().theta_y
                  << ", theta_r=" << discrete.back().theta_r << '\n';

        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
}

The Java implementation uses only the standard library and follows the same update ordering as the C++ program.

Chapter4_Lesson4.java

// Chapter4_Lesson4.java
// Continuous-time and discrete-time MRAC comparison for a scalar first-order plant.
// Build: javac Chapter4_Lesson4.java
// Run:   java Chapter4_Lesson4

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.ArrayList;
import java.util.List;
import java.util.Locale;

public final class Chapter4_Lesson4 {
    private Chapter4_Lesson4() {
    }

    private static final class Parameters {
        final double a = 1.6;
        final double b = 1.2;
        final double am = 1.0;
        final double bm = 1.0;
        final double gammaCt = 4.0;
        final double gammaDt = 0.35;
        final double tFinal = 20.0;
        final double dt = 0.001;
        final double sampleTime = 0.05;
    }

    private static final class Record {
        final double t;
        final double r;
        final double y;
        final double ym;
        final double e;
        final double u;
        final double thetaY;
        final double thetaR;
        final double epsilon;

        Record(
            double t,
            double r,
            double y,
            double ym,
            double e,
            double u,
            double thetaY,
            double thetaR,
            double epsilon
        ) {
            this.t = t;
            this.r = r;
            this.y = y;
            this.ym = ym;
            this.e = e;
            this.u = u;
            this.thetaY = thetaY;
            this.thetaR = thetaR;
            this.epsilon = epsilon;
        }
    }

    private static double referenceSignal(double t) {
        if (t < 5.0) {
            return 1.0;
        }
        if (t < 10.0) {
            return -0.5;
        }
        return 0.75 + 0.25 * Math.sin(0.7 * t);
    }

    private static double signNonzero(double value) {
        if (value == 0.0) {
            throw new IllegalArgumentException("The input gain b must be nonzero.");
        }
        return value > 0.0 ? 1.0 : -1.0;
    }

    private static double[] exactScalarZoh(
        double decayRate,
        double inputGain,
        double h
    ) {
        if (decayRate <= 0.0 || h <= 0.0) {
            throw new IllegalArgumentException(
                "decayRate and sample time must be positive."
            );
        }
        double ad = Math.exp(-decayRate * h);
        double bd = inputGain * (1.0 - ad) / decayRate;
        return new double[] {ad, bd};
    }

    private static List<Record> simulateContinuous(Parameters p) {
        int n = (int) Math.round(p.tFinal / p.dt);
        double signB = signNonzero(p.b);

        double y = 0.0;
        double ym = 0.0;
        double thetaY = 0.0;
        double thetaR = 0.0;

        List<Record> data = new ArrayList<>(n + 1);

        for (int i = 0; i <= n; i++) {
            double t = i * p.dt;
            double r = referenceSignal(t);
            double e = y - ym;
            double u = thetaY * y + thetaR * r;

            data.add(new Record(t, r, y, ym, e, u, thetaY, thetaR, 0.0));

            double yDot = -p.a * y + p.b * u;
            double ymDot = -p.am * ym + p.bm * r;
            double thetaYDot = -p.gammaCt * signB * y * e;
            double thetaRDot = -p.gammaCt * signB * r * e;

            y += p.dt * yDot;
            ym += p.dt * ymDot;
            thetaY += p.dt * thetaYDot;
            thetaR += p.dt * thetaRDot;
        }

        return data;
    }

    private static List<Record> simulateDiscrete(Parameters p) {
        double[] plant = exactScalarZoh(p.a, p.b, p.sampleTime);
        double[] model = exactScalarZoh(p.am, p.bm, p.sampleTime);
        double ad = plant[0];
        double bd = plant[1];
        double amd = model[0];
        double bmd = model[1];

        int n = (int) Math.round(p.tFinal / p.sampleTime);

        double y = 0.0;
        double ym = 0.0;
        double thetaY = 0.0;
        double thetaR = 0.0;

        List<Record> data = new ArrayList<>(n + 1);

        for (int k = 0; k <= n; k++) {
            double t = k * p.sampleTime;
            double r = referenceSignal(t);
            double e = y - ym;
            double u = thetaY * y + thetaR * r;

            double epsilon = 0.0;
            double yNext = y;
            double ymNext = ym;
            double thetaYNext = thetaY;
            double thetaRNext = thetaR;

            if (k < n) {
                yNext = ad * y + bd * u;
                ymNext = amd * ym + bmd * r;
                double eNext = yNext - ymNext;

                // One-step filtered error:
                // epsilon[k+1] = (e[k+1] - amd*e[k]) / bd.
                epsilon = (eNext - amd * e) / bd;
                double phiNormSquared = y * y + r * r;
                double normalizer = 1.0 + p.gammaDt * phiNormSquared;

                thetaYNext = thetaY - p.gammaDt * y * epsilon / normalizer;
                thetaRNext = thetaR - p.gammaDt * r * epsilon / normalizer;
            }

            data.add(
                new Record(t, r, y, ym, e, u, thetaY, thetaR, epsilon)
            );

            y = yNext;
            ym = ymNext;
            thetaY = thetaYNext;
            thetaR = thetaRNext;
        }

        return data;
    }

    private static void writeCsv(Path path, List<Record> data) throws IOException {
        try (
            BufferedWriter writer = Files.newBufferedWriter(
                path,
                StandardCharsets.UTF_8
            )
        ) {
            writer.write("t,r,y,ym,e,u,theta_y,theta_r,epsilon");
            writer.newLine();

            for (Record row : data) {
                writer.write(
                    String.format(
                        Locale.US,
                        "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g",
                        row.t,
                        row.r,
                        row.y,
                        row.ym,
                        row.e,
                        row.u,
                        row.thetaY,
                        row.thetaR,
                        row.epsilon
                    )
                );
                writer.newLine();
            }
        }
    }

    public static void main(String[] args) {
        try {
            Parameters p = new Parameters();
            List<Record> continuous = simulateContinuous(p);
            List<Record> discrete = simulateDiscrete(p);

            writeCsv(Path.of("Chapter4_Lesson4_continuous_java.csv"), continuous);
            writeCsv(Path.of("Chapter4_Lesson4_discrete_java.csv"), discrete);

            double[] plant = exactScalarZoh(p.a, p.b, p.sampleTime);
            double[] model = exactScalarZoh(p.am, p.bm, p.sampleTime);
            double idealThetaYCt = (p.a - p.am) / p.b;
            double idealThetaRCt = p.bm / p.b;
            double idealThetaYDt = (model[0] - plant[0]) / plant[1];
            double idealThetaRDt = model[1] / plant[1];
            Record ctFinal = continuous.get(continuous.size() - 1);
            Record dtFinal = discrete.get(discrete.size() - 1);

            System.out.printf(
                Locale.US,
                "Continuous ideal parameters: theta_y*=%.6f, theta_r*=%.6f%n",
                idealThetaYCt,
                idealThetaRCt
            );
            System.out.printf(
                Locale.US,
                "Discrete ideal parameters: theta_y*=%.6f, theta_r*=%.6f%n",
                idealThetaYDt,
                idealThetaRDt
            );
            System.out.printf(
                Locale.US,
                "Final CT parameters: theta_y=%.6f, theta_r=%.6f%n",
                ctFinal.thetaY,
                ctFinal.thetaR
            );
            System.out.printf(
                Locale.US,
                "Final DT parameters: theta_y=%.6f, theta_r=%.6f%n",
                dtFinal.thetaY,
                dtFinal.thetaR
            );
        } catch (IOException | RuntimeException error) {
            System.err.println("Error: " + error.getMessage());
            System.exit(1);
        }
    }
}

The MATLAB script performs both simulations, plots the results, exports tables, and programmatically constructs an executable continuous-time Simulink model when Simulink is available.

Chapter4_Lesson4.m

% Chapter4_Lesson4.m
% Continuous-time and discrete-time MRAC comparison for a scalar plant.
% The script also builds an executable continuous-time Simulink model when
% Simulink is installed.

clear;
clc;
close all;

p.a = 1.6;
p.b = 1.2;
p.am = 1.0;
p.bm = 1.0;
p.gammaCT = 4.0;
p.gammaDT = 0.35;
p.tFinal = 20.0;
p.dt = 0.001;
p.Ts = 0.05;

if p.b == 0
    error('The input gain b must be nonzero.');
end

%% Continuous-time pedagogical simulation
nCT = round(p.tFinal / p.dt);
tCT = (0:nCT)' * p.dt;

yCT = zeros(nCT + 1, 1);
ymCT = zeros(nCT + 1, 1);
eCT = zeros(nCT + 1, 1);
uCT = zeros(nCT + 1, 1);
thetaYCT = zeros(nCT + 1, 1);
thetaRCT = zeros(nCT + 1, 1);
rCT = zeros(nCT + 1, 1);

for i = 1:nCT
    t = tCT(i);
    rCT(i) = referenceSignal(t);
    eCT(i) = yCT(i) - ymCT(i);
    uCT(i) = thetaYCT(i) * yCT(i) + thetaRCT(i) * rCT(i);

    yDot = -p.a * yCT(i) + p.b * uCT(i);
    ymDot = -p.am * ymCT(i) + p.bm * rCT(i);
    thetaYDot = -p.gammaCT * sign(p.b) * yCT(i) * eCT(i);
    thetaRDot = -p.gammaCT * sign(p.b) * rCT(i) * eCT(i);

    yCT(i + 1) = yCT(i) + p.dt * yDot;
    ymCT(i + 1) = ymCT(i) + p.dt * ymDot;
    thetaYCT(i + 1) = thetaYCT(i) + p.dt * thetaYDot;
    thetaRCT(i + 1) = thetaRCT(i) + p.dt * thetaRDot;
end

rCT(end) = referenceSignal(tCT(end));
eCT(end) = yCT(end) - ymCT(end);
uCT(end) = thetaYCT(end) * yCT(end) + thetaRCT(end) * rCT(end);

%% Exact-ZOH discrete-time simulation
[ad, bd] = exactScalarZOH(p.a, p.b, p.Ts);
[amd, bmd] = exactScalarZOH(p.am, p.bm, p.Ts);

nDT = round(p.tFinal / p.Ts);
tDT = (0:nDT)' * p.Ts;

yDT = zeros(nDT + 1, 1);
ymDT = zeros(nDT + 1, 1);
eDT = zeros(nDT + 1, 1);
uDT = zeros(nDT + 1, 1);
thetaYDT = zeros(nDT + 1, 1);
thetaRDT = zeros(nDT + 1, 1);
epsilonDT = zeros(nDT + 1, 1);
rDT = zeros(nDT + 1, 1);

for k = 1:nDT
    t = tDT(k);
    rDT(k) = referenceSignal(t);
    eDT(k) = yDT(k) - ymDT(k);
    uDT(k) = thetaYDT(k) * yDT(k) + thetaRDT(k) * rDT(k);

    yDT(k + 1) = ad * yDT(k) + bd * uDT(k);
    ymDT(k + 1) = amd * ymDT(k) + bmd * rDT(k);
    eNext = yDT(k + 1) - ymDT(k + 1);

    % Under matching:
    % epsilon[k+1] = (e[k+1] - amd*e[k]) / bd
    %              = tildeTheta[k]'*phi[k].
    epsilonDT(k) = (eNext - amd * eDT(k)) / bd;

    normalizer = 1 + p.gammaDT * (yDT(k)^2 + rDT(k)^2);
    thetaYDT(k + 1) = thetaYDT(k) ...
        - p.gammaDT * yDT(k) * epsilonDT(k) / normalizer;
    thetaRDT(k + 1) = thetaRDT(k) ...
        - p.gammaDT * rDT(k) * epsilonDT(k) / normalizer;
end

rDT(end) = referenceSignal(tDT(end));
eDT(end) = yDT(end) - ymDT(end);
uDT(end) = thetaYDT(end) * yDT(end) + thetaRDT(end) * rDT(end);

%% Export and report
continuousTable = table( ...
    tCT, rCT, yCT, ymCT, eCT, uCT, thetaYCT, thetaRCT, ...
    'VariableNames', {'t','r','y','ym','e','u','thetaY','thetaR'} ...
);
discreteTable = table( ...
    tDT, rDT, yDT, ymDT, eDT, uDT, thetaYDT, thetaRDT, epsilonDT, ...
    'VariableNames', {'t','r','y','ym','e','u','thetaY','thetaR','epsilon'} ...
);

writetable(continuousTable, 'Chapter4_Lesson4_continuous_matlab.csv');
writetable(discreteTable, 'Chapter4_Lesson4_discrete_matlab.csv');

idealThetaYCT = (p.a - p.am) / p.b;
idealThetaRCT = p.bm / p.b;
idealThetaYDT = (amd - ad) / bd;
idealThetaRDT = bmd / bd;

fprintf('Continuous ideal parameters: theta_y*=%.6f, theta_r*=%.6f\n', ...
    idealThetaYCT, idealThetaRCT);
fprintf('Discrete ideal parameters: theta_y*=%.6f, theta_r*=%.6f\n', ...
    idealThetaYDT, idealThetaRDT);
fprintf('Final CT parameters: theta_y=%.6f, theta_r=%.6f\n', ...
    thetaYCT(end), thetaRCT(end));
fprintf('Final DT parameters: theta_y=%.6f, theta_r=%.6f\n', ...
    thetaYDT(end), thetaRDT(end));

figure;
plot(tCT, ymCT, 'LineWidth', 1.4);
hold on;
plot(tCT, yCT, 'LineWidth', 1.1);
stairs(tDT, ymDT, 'LineWidth', 1.2);
stairs(tDT, yDT, 'LineWidth', 1.0);
grid on;
xlabel('Time (s)');
ylabel('Output');
title('Continuous-Time and Discrete-Time MRAC');
legend( ...
    'CT reference model', 'CT plant', ...
    'DT reference model', 'DT plant', ...
    'Location', 'best' ...
);

figure;
plot(tCT, thetaYCT, 'LineWidth', 1.2);
hold on;
plot(tCT, thetaRCT, 'LineWidth', 1.2);
stairs(tDT, thetaYDT, 'LineWidth', 1.2);
stairs(tDT, thetaRDT, 'LineWidth', 1.2);
yline(idealThetaYCT, '--');
yline(idealThetaRCT, '--');
yline(idealThetaYDT, ':');
yline(idealThetaRDT, ':');
grid on;
xlabel('Time (s)');
ylabel('Adaptive parameter');
title('Adaptive-Parameter Evolution');
legend( ...
    'CT theta_y', 'CT theta_r', ...
    'DT theta_y', 'DT theta_r', ...
    'CT ideal theta_y', 'CT ideal theta_r', ...
    'DT ideal theta_y', 'DT ideal theta_r', ...
    'Location', 'best' ...
);

%% Optional Simulink construction
if license('test', 'Simulink')
    modelName = 'Chapter4_Lesson4_Simulink';
    buildContinuousMRACSimulinkModel(modelName, p);
    fprintf('Created %s.slx\n', modelName);
else
    fprintf('Simulink is unavailable; skipped .slx model generation.\n');
end

function value = referenceSignal(t)
    if t < 5
        value = 1.0;
    elseif t < 10
        value = -0.5;
    else
        value = 0.75 + 0.25 * sin(0.7 * t);
    end
end

function [ad, bd] = exactScalarZOH(decayRate, inputGain, sampleTime)
    if decayRate <= 0 || sampleTime <= 0
        error('decayRate and sampleTime must be positive.');
    end
    ad = exp(-decayRate * sampleTime);
    bd = inputGain * (1 - ad) / decayRate;
end

function buildContinuousMRACSimulinkModel(modelName, p)
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end

    new_system(modelName);
    open_system(modelName);
    set_param(modelName, 'StopTime', num2str(p.tFinal));

    % Reference input
    add_block('simulink/Sources/Step', [modelName '/Reference'], ...
        'Time', '0', 'Before', '0', 'After', '1', ...
        'Position', [30 120 60 150]);

    % Reference model: ym_dot = -am*ym + bm*r
    add_block('simulink/Math Operations/Gain', [modelName '/bm'], ...
        'Gain', num2str(p.bm), 'Position', [100 105 150 135]);
    add_block('simulink/Math Operations/Sum', [modelName '/ReferenceSum'], ...
        'Inputs', '+-', 'Position', [190 100 215 145]);
    add_block('simulink/Continuous/Integrator', [modelName '/ym'], ...
        'Position', [255 105 285 135]);
    add_block('simulink/Math Operations/Gain', [modelName '/am'], ...
        'Gain', num2str(p.am), 'Position', [255 165 305 195]);

    % Error e = y - ym
    add_block('simulink/Math Operations/Sum', [modelName '/TrackingError'], ...
        'Inputs', '+-', 'Position', [610 195 635 240]);

    % Adaptive controller u = theta_y*y + theta_r*r
    add_block('simulink/Math Operations/Product', [modelName '/thetaY_times_y'], ...
        'Position', [420 300 455 335]);
    add_block('simulink/Math Operations/Product', [modelName '/thetaR_times_r'], ...
        'Position', [420 370 455 405]);
    add_block('simulink/Math Operations/Sum', [modelName '/ControlSum'], ...
        'Inputs', '++', 'Position', [500 325 525 380]);

    % Plant: y_dot = -a*y + b*u
    add_block('simulink/Math Operations/Gain', [modelName '/b'], ...
        'Gain', num2str(p.b), 'Position', [565 325 615 355]);
    add_block('simulink/Math Operations/Sum', [modelName '/PlantSum'], ...
        'Inputs', '+-', 'Position', [655 315 680 365]);
    add_block('simulink/Continuous/Integrator', [modelName '/y'], ...
        'Position', [720 325 750 355]);
    add_block('simulink/Math Operations/Gain', [modelName '/a'], ...
        'Gain', num2str(p.a), 'Position', [720 390 770 420]);

    % theta_y_dot = -gamma*y*e
    add_block('simulink/Math Operations/Product', [modelName '/y_times_e'], ...
        'Position', [650 500 685 535]);
    add_block('simulink/Math Operations/Gain', [modelName '/minusGammaY'], ...
        'Gain', num2str(-p.gammaCT * sign(p.b)), ...
        'Position', [720 495 790 535]);
    add_block('simulink/Continuous/Integrator', [modelName '/thetaY'], ...
        'Position', [830 500 860 530]);

    % theta_r_dot = -gamma*r*e
    add_block('simulink/Math Operations/Product', [modelName '/r_times_e'], ...
        'Position', [650 590 685 625]);
    add_block('simulink/Math Operations/Gain', [modelName '/minusGammaR'], ...
        'Gain', num2str(-p.gammaCT * sign(p.b)), ...
        'Position', [720 585 790 625]);
    add_block('simulink/Continuous/Integrator', [modelName '/thetaR'], ...
        'Position', [830 590 860 620]);

    % Visualization
    add_block('simulink/Signal Routing/Mux', [modelName '/OutputMux'], ...
        'Inputs', '3', 'Position', [820 170 825 250]);
    add_block('simulink/Sinks/Scope', [modelName '/Scope'], ...
        'Position', [880 185 930 235]);

    % Reference-model connections
    add_line(modelName, 'Reference/1', 'bm/1', 'autorouting', 'on');
    add_line(modelName, 'bm/1', 'ReferenceSum/1', 'autorouting', 'on');
    add_line(modelName, 'ReferenceSum/1', 'ym/1', 'autorouting', 'on');
    add_line(modelName, 'ym/1', 'am/1', 'autorouting', 'on');
    add_line(modelName, 'am/1', 'ReferenceSum/2', 'autorouting', 'on');

    % Controller and plant connections
    add_line(modelName, 'thetaY/1', 'thetaY_times_y/1', 'autorouting', 'on');
    add_line(modelName, 'y/1', 'thetaY_times_y/2', 'autorouting', 'on');
    add_line(modelName, 'thetaR/1', 'thetaR_times_r/1', 'autorouting', 'on');
    add_line(modelName, 'Reference/1', 'thetaR_times_r/2', 'autorouting', 'on');
    add_line(modelName, 'thetaY_times_y/1', 'ControlSum/1', 'autorouting', 'on');
    add_line(modelName, 'thetaR_times_r/1', 'ControlSum/2', 'autorouting', 'on');
    add_line(modelName, 'ControlSum/1', 'b/1', 'autorouting', 'on');
    add_line(modelName, 'b/1', 'PlantSum/1', 'autorouting', 'on');
    add_line(modelName, 'PlantSum/1', 'y/1', 'autorouting', 'on');
    add_line(modelName, 'y/1', 'a/1', 'autorouting', 'on');
    add_line(modelName, 'a/1', 'PlantSum/2', 'autorouting', 'on');

    % Error and adaptation connections
    add_line(modelName, 'y/1', 'TrackingError/1', 'autorouting', 'on');
    add_line(modelName, 'ym/1', 'TrackingError/2', 'autorouting', 'on');
    add_line(modelName, 'y/1', 'y_times_e/1', 'autorouting', 'on');
    add_line(modelName, 'TrackingError/1', 'y_times_e/2', 'autorouting', 'on');
    add_line(modelName, 'y_times_e/1', 'minusGammaY/1', 'autorouting', 'on');
    add_line(modelName, 'minusGammaY/1', 'thetaY/1', 'autorouting', 'on');
    add_line(modelName, 'Reference/1', 'r_times_e/1', 'autorouting', 'on');
    add_line(modelName, 'TrackingError/1', 'r_times_e/2', 'autorouting', 'on');
    add_line(modelName, 'r_times_e/1', 'minusGammaR/1', 'autorouting', 'on');
    add_line(modelName, 'minusGammaR/1', 'thetaR/1', 'autorouting', 'on');

    % Scope: y, ym, e
    add_line(modelName, 'y/1', 'OutputMux/1', 'autorouting', 'on');
    add_line(modelName, 'ym/1', 'OutputMux/2', 'autorouting', 'on');
    add_line(modelName, 'TrackingError/1', 'OutputMux/3', 'autorouting', 'on');
    add_line(modelName, 'OutputMux/1', 'Scope/1', 'autorouting', 'on');

    save_system(modelName);
end

The notebook uses NDSolveValue for continuous dynamics and FoldList for the discrete recursion.

Chapter4_Lesson4.nb

(* Chapter4_Lesson4.nb
   Continuous-time and discrete-time MRAC comparison for a scalar plant. *)

ClearAll["Global`*"];

a = 1.6;
b = 1.2;
am = 1.0;
bm = 1.0;
gammaCT = 4.0;
gammaDT = 0.35;
tFinal = 20.0;
sampleTime = 0.05;

If[b == 0.0, Print["The input gain b must be nonzero."]; Abort[]];

reference[t_?NumericQ] := Piecewise[
  {
    {1.0, t < 5.0},
    {-0.5, t < 10.0}
  },
  0.75 + 0.25 Sin[0.7 t]
];

(* Continuous-time Lyapunov-based MRAC. *)
continuousSolution = NDSolveValue[
  {
    y'[t] == -a y[t] + b (thetaY[t] y[t] + thetaR[t] reference[t]),
    ym'[t] == -am ym[t] + bm reference[t],
    thetaY'[t] == -gammaCT Sign[b] y[t] (y[t] - ym[t]),
    thetaR'[t] == -gammaCT Sign[b] reference[t] (y[t] - ym[t]),
    y[0] == 0.0,
    ym[0] == 0.0,
    thetaY[0] == 0.0,
    thetaR[0] == 0.0
  },
  {y, ym, thetaY, thetaR},
  {t, 0.0, tFinal},
  MaxStepFraction -> 1/1000
];

{yCT, ymCT, thetaYCT, thetaRCT} = continuousSolution;
tCT = Range[0.0, tFinal, 0.01];

continuousRows = Table[
  With[
    {
      r = reference[t0],
      yp = yCT[t0],
      ymp = ymCT[t0],
      typ = thetaYCT[t0],
      trp = thetaRCT[t0]
    },
    {
      t0,
      r,
      yp,
      ymp,
      yp - ymp,
      typ yp + trp r,
      typ,
      trp
    }
  ],
  {t0, tCT}
];

(* Exact zero-order-hold coefficients for x_dot=-rate*x+gain*u. *)
exactScalarZOH[rate_?NumericQ, gain_?NumericQ, h_?NumericQ] := Module[
  {ad, bd},
  If[rate <= 0.0 || h <= 0.0,
    Print["rate and sample time must be positive."];
    Abort[]
  ];
  ad = Exp[-rate h];
  bd = gain (1.0 - ad)/rate;
  {ad, bd}
];

{ad, bd} = exactScalarZOH[a, b, sampleTime];
{amd, bmd} = exactScalarZOH[am, bm, sampleTime];
nDT = Round[tFinal/sampleTime];

discreteStep[state_List, k_Integer] := Module[
  {
    yk = state[[1]],
    ymk = state[[2]],
    thetaYk = state[[3]],
    thetaRk = state[[4]],
    tk,
    rk,
    ek,
    uk,
    yNext,
    ymNext,
    eNext,
    epsilon,
    normalizer
  },
  tk = k sampleTime;
  rk = reference[tk];
  ek = yk - ymk;
  uk = thetaYk yk + thetaRk rk;

  yNext = ad yk + bd uk;
  ymNext = amd ymk + bmd rk;
  eNext = yNext - ymNext;

  epsilon = (eNext - amd ek)/bd;
  normalizer = 1.0 + gammaDT (yk^2 + rk^2);

  {
    yNext,
    ymNext,
    thetaYk - gammaDT yk epsilon/normalizer,
    thetaRk - gammaDT rk epsilon/normalizer
  }
];

discreteStates = FoldList[
  discreteStep,
  {0.0, 0.0, 0.0, 0.0},
  Range[0, nDT - 1]
];

tDT = Range[0, nDT] sampleTime;

discreteRows = Table[
  Module[
    {
      state = discreteStates[[k + 1]],
      tk = tDT[[k + 1]],
      yk,
      ymk,
      thetaYk,
      thetaRk,
      rk,
      ek,
      uk,
      epsilon = 0.0,
      yNext,
      ymNext,
      eNext
    },
    {yk, ymk, thetaYk, thetaRk} = state;
    rk = reference[tk];
    ek = yk - ymk;
    uk = thetaYk yk + thetaRk rk;

    If[k < nDT,
      yNext = ad yk + bd uk;
      ymNext = amd ymk + bmd rk;
      eNext = yNext - ymNext;
      epsilon = (eNext - amd ek)/bd;
    ];

    {tk, rk, yk, ymk, ek, uk, thetaYk, thetaRk, epsilon}
  ],
  {k, 0, nDT}
];

idealThetaYCT = (a - am)/b;
idealThetaRCT = bm/b;
idealThetaYDT = (amd - ad)/bd;
idealThetaRDT = bmd/bd;

Print[
  "Continuous ideal parameters: theta_y*=",
  NumberForm[idealThetaYCT, {8, 6}],
  ", theta_r*=", NumberForm[idealThetaRCT, {8, 6}]
];
Print[
  "Discrete ideal parameters: theta_y*=",
  NumberForm[idealThetaYDT, {8, 6}],
  ", theta_r*=", NumberForm[idealThetaRDT, {8, 6}]
];
Print[
  "Final CT parameters: theta_y=", NumberForm[thetaYCT[tFinal], {8, 6}],
  ", theta_r=", NumberForm[thetaRCT[tFinal], {8, 6}]
];
Print[
  "Final DT parameters: theta_y=",
  NumberForm[discreteStates[[-1, 3]], {8, 6}],
  ", theta_r=",
  NumberForm[discreteStates[[-1, 4]], {8, 6}]
];

Export[
  FileNameJoin[{NotebookDirectory[], "Chapter4_Lesson4_continuous_wolfram.csv"}],
  Prepend[
    continuousRows,
    {"t", "r", "y", "ym", "e", "u", "thetaY", "thetaR"}
  ]
];

Export[
  FileNameJoin[{NotebookDirectory[], "Chapter4_Lesson4_discrete_wolfram.csv"}],
  Prepend[
    discreteRows,
    {"t", "r", "y", "ym", "e", "u", "thetaY", "thetaR", "epsilon"}
  ]
];

trackingPlot = Show[
  ListLinePlot[
    {
      continuousRows[[All, {1, 3}]],
      continuousRows[[All, {1, 4}]]
    },
    PlotLegends -> {"CT plant", "CT reference model"},
    PlotRange -> All
  ],
  ListStepPlot[
    {
      discreteRows[[All, {1, 3}]],
      discreteRows[[All, {1, 4}]]
    },
    PlotLegends -> {"DT plant", "DT reference model"},
    PlotRange -> All
  ],
  Frame -> True,
  FrameLabel -> {"Time (s)", "Output"},
  PlotLabel -> "Continuous-Time and Discrete-Time MRAC"
];

parameterPlot = Show[
  Plot[
    {thetaYCT[t], thetaRCT[t]},
    {t, 0.0, tFinal},
    PlotLegends -> {"CT theta_y", "CT theta_r"}
  ],
  ListStepPlot[
    {
      discreteRows[[All, {1, 7}]],
      discreteRows[[All, {1, 8}]]
    },
    PlotLegends -> {"DT theta_y", "DT theta_r"},
    PlotRange -> All
  ],
  Frame -> True,
  FrameLabel -> {"Time (s)", "Adaptive parameter"},
  PlotLabel -> "Adaptive-Parameter Evolution"
];

Column[{trackingPlot, parameterPlot}]

The programs use the same physical plant and desired continuous model, but the discrete controller is matched to the exact sampled model. Consequently, its ideal gains are computed from \( a_d,b_d,a_{m,d},b_{m,d} \), not copied directly from the continuous design.

Similar tracking curves do not imply mathematically identical adaptive systems. The continuous implementation updates parameters at every integration step; the discrete implementation applies one parameter vector during an entire sample interval and updates only after the next output is available.

Students should vary \( h \) and \( \gamma \) separately. A smaller sample period generally improves the sampled approximation but increases computation and data rate. A larger adaptation gain can speed transient adjustment but can also create large parameter motion and sensitivity to sample-level effects.

Problem 1 (Pole Mapping): A continuous reference model has pole \( \lambda=-4 \). Find its exact discrete pole for \( h=0.1 \) s and compare it with the forward-Euler pole.

Solution:

\[ z=e^{\lambda h}=e^{-0.4}\approx 0.6703, \qquad z_{\mathrm{FE}}=1+h\lambda=0.6. \]

Both are inside the unit circle, but they are not equal. Exact sampling preserves the exact sample values of the continuous mode under a held input, whereas forward Euler introduces integration error.

Problem 2 (Forward-Euler Stability Limit): For the scalar mode \( \dot{x}=-12x \), determine the interval of forward-Euler step sizes that keeps the numerical mode stable.

Solution:

\[ |1-12h|<1 \quad\Longleftrightarrow\quad -1<1-12h<1. \]

\[ 0<h<\frac{1}{6}\approx 0.1667\ \text{s}. \]

Problem 3 (Continuous Lyapunov Cancellation): Starting from \( \dot{e}=A_me+b\tilde{\theta}^T\omega \) and \( V=e^TPe+\tilde{\theta}^T\Gamma^{-1}\tilde{\theta} \), derive an update law that removes the mixed term.

Solution:

\[ \dot{V} =-e^TQe +2e^TPb\,\tilde{\theta}^T\omega +2\tilde{\theta}^T\Gamma^{-1}\dot{\hat{\theta}}. \]

Choose

\[ \dot{\hat{\theta}} =-\Gamma\omega e^TPb. \]

Then \( 2\tilde{\theta}^T\Gamma^{-1}\dot{\hat{\theta}} =-2e^TPb\,\tilde{\theta}^T\omega \), so \( \dot{V}=-e^TQe\leq 0 \).

Problem 4 (Discrete Quadratic Increment): Let \( \theta^+=\theta+\delta\theta \) and \( W=\theta^TM\theta \), with \( M=M^T>0 \). Derive \( W^+-W \).

Solution:

\[ \begin{aligned} W^+-W &=(\theta+\delta\theta)^TM(\theta+\delta\theta)-\theta^TM\theta \\ &=2\theta^TM\delta\theta+\delta\theta^TM\delta\theta. \end{aligned} \]

The last term is nonnegative and must be included in a discrete stability proof.

Problem 5 (Exact Discrete Ideal Gains): For \( a=1.6 \), \( b=1.2 \), \( a_m=1 \), \( b_m=1 \), and \( h=0.05 \) s, compute the continuous and exact discrete ideal gains.

Solution:

\[ \theta_{y,c}^*=\frac{1.6-1}{1.2}=0.5, \qquad \theta_{r,c}^*=\frac{1}{1.2}\approx 0.8333. \]

\[ a_d=e^{-0.08}\approx 0.9231,\qquad b_d=\frac{1.2}{1.6}(1-e^{-0.08})\approx 0.0577, \]

\[ a_{m,d}=e^{-0.05}\approx 0.9512,\qquad b_{m,d}=1-e^{-0.05}\approx 0.0488. \]

\[ \theta_{y,d}^* =\frac{0.9512-0.9231}{0.0577}\approx 0.487, \qquad \theta_{r,d}^* =\frac{0.0488}{0.0577}\approx 0.846. \]

The finite-sample ideal gains are close to, but not identical with, the continuous ideal gains.

Problem 6 (Normalized Descent): Prove that the normalized update in Section 8 makes the parameter energy nonincreasing.

Solution:

Let \( m=1+s \), where \( s=\phi^T\Gamma\phi \), and let \( \varepsilon=\tilde{\theta}^T\phi \). Then

\[ \tilde{\theta}^+ =\tilde{\theta}-\frac{\Gamma\phi\varepsilon}{m}. \]

Expanding the quadratic energy gives

\[ \begin{aligned} W^+-W &=-\frac{2\varepsilon^2}{m} +\frac{s\varepsilon^2}{m^2} \\ &=-\frac{(2+s)\varepsilon^2}{(1+s)^2}\leq 0. \end{aligned} \]

Problem 7 (Timing): Why can the scalar one-step error \( \varepsilon[k+1] \) not generally be used to compute \( u[k] \)?

Solution:

It depends on \( y[k+1] \), which is measured only after \( u[k] \) has been applied over the current sampling interval. It may update the parameters used for \( u[k+1] \), but using it for \( u[k] \) would require a future measurement.

Continuous-time and discrete-time MRAC share the model-following objective and matching-based error structure, but they are distinct dynamical systems. Continuous design studies \( \dot{V} \) and differential parameter laws. Discrete design studies \( V[k+1]-V[k] \), update timing, exact sampled models, and quadratic increment terms. Exact zero-order-hold discretization maps continuous poles by \( z=e^{\lambda h} \), while the adaptive recursion must be designed directly for sample-to-sample dynamics. A normalized one-step update provides a simple example of an exact discrete descent argument.