After completing this lesson, students should be able to:

• state the classical Nussbaum-function property in integral-average form;
• prove the Nussbaum boundedness lemma by contradiction;
• construct a Lyapunov function for a reduced unknown-control-direction error model;
• show boundedness of the adaptive variable, gain, error, and control signal;
• apply Barbalat's lemma to conclude asymptotic regulation;
• distinguish constant-sign proofs from stronger time-varying-gain results.

The lesson uses Lyapunov stability, positive-definite functions, square-integrability, uniform continuity, and Barbalat's lemma from earlier chapters. Lessons 1 and 2 of this chapter introduced the unknown-control-direction problem and the conceptual role of Nussbaum-type functions.

Consider an error variable \( e(t) \in \mathbb{R} \) whose nominally compensated dynamics have been reduced to

\[ \dot e = -\lambda e + b u, \qquad \lambda > 0, \qquad b \in \mathbb{R}\setminus\{0\}. \]

The magnitude and sign of \( b \) are unknown, but \( b \) is assumed constant and nonzero. If its sign were known, the feedback \(u=-\operatorname{sgn}(b)k e\) with \(k>0\) would add damping. With an unknown sign, one fixed feedback polarity stabilizes one case and destabilizes the other.

The Nussbaum design used throughout this lesson is

\[ u=N(\kappa)e, \qquad \dot\kappa=e^2, \qquad \kappa(0)=\kappa_0. \]

Because \(\dot\kappa=e^2\geq 0\), the adaptive variable is continuous and nondecreasing. The gain \(N(\kappa)\) explores alternating effective feedback directions while its integral average grows without bound in both signs.

flowchart TD
  R["Desired value: zero"] --> S["Error e"]
  S --> C["Controller u = N(kappa) e"]
  C --> P["Reduced channel: e_dot = -lambda e + b u"]
  P --> S
  S --> A["Adaptive law: kappa_dot = e squared"]
  A --> C
  U["Unknown constant b: \npositive or negative"] --> P
        

A continuous function \(N:[0,\infty)\to\mathbb{R}\) is a classical Nussbaum function when

\[ \limsup_{K\to\infty}\frac{1}{K}\int_0^K N(s)\,ds=+\infty, \qquad \liminf_{K\to\infty}\frac{1}{K}\int_0^K N(s)\,ds=-\infty. \]

The pointwise sign changes of \(N\) are not enough. Stability uses the unbounded positive and negative excursions of its integral average. A standard polynomial-growth choice is

\[ N(\kappa)=\kappa^2\cos\kappa. \]

Its integral is available in closed form:

\[ \int_0^K s^2\cos s\,ds =K^2\sin K+2K\cos K-2\sin K. \]

Therefore the integral average is

\[ A(K)=\frac{1}{K}\int_0^K s^2\cos s\,ds =K\sin K+2\cos K-\frac{2\sin K}{K}. \]

Along \(K_n^+=\frac{\pi}{2}+2\pi n\), \(A(K_n^+)\to+\infty\). Along \(K_n^-=\frac{3\pi}{2}+2\pi n\), \(A(K_n^-)\to-\infty\). Thus \(\kappa^2\cos\kappa\) satisfies the classical Nussbaum property.

Other commonly used candidates include

\[ N_1(\kappa)=\kappa^2\sin\kappa, \qquad N_2(\kappa)=\exp(\kappa^2)\cos\left(\frac{\pi\kappa}{2}\right). \]

The exponential candidate can create severe numerical transients and overflow. Polynomial-growth functions are usually preferable for a first implementation, although numerical convenience does not replace the analytical Nussbaum conditions.

Theorem 1 (monotone integral form). Suppose:

1. \(V(t)\geq0\) is continuous;
2. \(\kappa(t)\) is continuously differentiable and nondecreasing;
3. \(N\) is a classical Nussbaum function;
4. \(b\neq0\) is an unknown constant;
5. for constants \(c_0\geq0\) and \(c_1\in\mathbb{R}\),

\[ V(t)\leq c_0+c_1\bigl(\kappa(t)-\kappa_0\bigr) +b\int_{\kappa_0}^{\kappa(t)}N(s)\,ds. \]

Then \(\kappa(t)\) and \(V(t)\) are bounded on their interval of existence.

Proof. Assume, for contradiction, that \(\kappa(t)\) is unbounded. Since it is continuous and nondecreasing, it reaches arbitrarily large values. Define

\[ F(K)=c_1(K-\kappa_0)+b\int_{\kappa_0}^{K}N(s)\,ds. \]

Dividing by \(K\) gives

\[ \frac{F(K)}{K} =c_1\left(1-\frac{\kappa_0}{K}\right) +b\frac{1}{K}\int_{\kappa_0}^{K}N(s)\,ds. \]

The finite lower integration interval only contributes a term that vanishes after division by \(K\). If \(b>0\), the negative Nussbaum excursion implies

\[ \liminf_{K\to\infty}\frac{F(K)}{K}=-\infty. \]

If \(b<0\), multiplying the positive Nussbaum excursion by \(b\) gives the same conclusion. Hence there exists a sequence \(K_j\to\infty\) for which \(c_0+F(K_j)<0\). Because an unbounded continuous nondecreasing \(\kappa\) reaches each \(K_j\), there is a time \(t_j\) such that \(\kappa(t_j)=K_j\). The assumed inequality then yields \(V(t_j)<0\), contradicting \(V(t)\geq0\). Therefore \(\kappa\) is bounded. Continuity of the right-hand side on the resulting compact interval then implies boundedness of \(V\). \(\square\)

The lemma is a contradiction mechanism, not an ordinary sign-definite Lyapunov argument. The derivative of the Lyapunov function may become positive during portions of the trajectory. Stability follows because an unbounded adaptive excursion would force the integrated upper bound for a nonnegative Lyapunov function to become negative.

Substituting the Nussbaum controller into the error channel gives

\[ \dot e=\bigl[-\lambda+bN(\kappa)\bigr]e, \qquad \dot\kappa=e^2. \]

Choose the quadratic Lyapunov function

\[ V(e)=\frac{1}{2}e^2. \]

Its derivative is

\[ \dot V=e\dot e =\bigl[-\lambda+bN(\kappa)\bigr]e^2 =\bigl[-\lambda+bN(\kappa)\bigr]\dot\kappa. \]

Integrating from \(0\) to \(t\) yields

\[ V(t)=V(0)-\lambda\bigl(\kappa(t)-\kappa_0\bigr) +b\int_{\kappa_0}^{\kappa(t)}N(s)\,ds. \]

This is exactly Theorem 1 with \(c_0=V(0)\) and \(c_1=-\lambda\). Consequently, \(V(t)\) and \(\kappa(t)\) are bounded.

Theorem 2. Consider

\[ \dot e=-\lambda e+bu, \qquad u=N(\kappa)e, \qquad \dot\kappa=e^2, \]

where \(\lambda>0\), \(b\neq0\) is constant, and \(N\) is a classical Nussbaum function. Then all closed-loop signals are bounded and \(\lim_{t\to\infty}e(t)=0\).

Proof.

Step 1: Boundedness of \(V\) and \(\kappa\). The identity in Section 5 and Theorem 1 imply

\[ 0\leq V(t)\leq \bar V, \qquad \kappa_0\leq\kappa(t)\leq\bar\kappa. \]

Step 2: Boundedness of the error. Since \(V=e^2/2\),

\[ |e(t)|\leq\sqrt{2\bar V}. \]

Step 3: Square-integrability. From \(\dot\kappa=e^2\),

\[ \int_0^t e^2(\tau)\,d\tau=\kappa(t)-\kappa_0 \leq\bar\kappa-\kappa_0. \]

Therefore \(e\in L_2\).

Step 4: Boundedness of the gain and control. Since \(\kappa(t)\) stays in a compact interval and \(N\) is continuous,

\[ |N(\kappa(t))|\leq\bar N. \]

Thus

\[ |u(t)|=|N(\kappa(t))e(t)| \leq\bar N\sqrt{2\bar V}. \]

Step 5: Uniform continuity. The error derivative satisfies

\[ |\dot e| \leq\bigl(\lambda+|b|\bar N\bigr)|e|, \]

so \(\dot e\in L_\infty\). Hence \(e(t)\) is uniformly continuous.

Step 6: Asymptotic convergence. Since \(e\in L_2\) and \(e\) is uniformly continuous, Barbalat's lemma gives

\[ \lim_{t\to\infty}e(t)=0. \]

Finally, bounded monotone \(\kappa(t)\) converges to a finite limit \(\kappa_\infty\). This limit need not identify \(b\); the mechanism guarantees regulation, not parameter identification. \(\square\)

flowchart TD
  A["Choose V = one half e squared"] --> B["Obtain V_dot = (-lambda + b N(kappa)) kappa_dot"]
  B --> C["Integrate with respect to kappa"]
  C --> D["Assume kappa is unbounded"]
  D --> E["Nussbaum average creates a negative subsequence"]
  E --> F["Integrated upper bound would give V below zero"]
  F --> G["Contradiction: kappa and V are bounded"]
  G --> H["kappa_dot = e squared gives e in L2"]
  H --> I["Bounded e and kappa give bounded e_dot"]
  I --> J["Barbalat lemma gives e converges to zero"]
        

In higher-order adaptive designs, the exact scalar equality is usually replaced by an inequality of the form

\[ \dot V\leq-W(z)+\bigl[\alpha+bN(\kappa)\bigr]\dot\kappa, \qquad \dot\kappa\geq0, \]

where \(V(z)\) is positive definite, \(W(z)\geq0\), and \(\alpha\) is a known finite constant. Integration gives

\[ V(t)+\int_0^t W(z(\tau))\,d\tau \leq V(0)+\alpha\bigl(\kappa(t)-\kappa_0\bigr) +b\int_{\kappa_0}^{\kappa(t)}N(s)\,ds. \]

The Nussbaum lemma bounds \(V\) and \(\kappa\). It also yields

\[ \int_0^\infty W(z(\tau))\,d\tau<\infty. \]

If \(W(z(t))\) is uniformly continuous, Barbalat's lemma gives \(W(z(t))\to0\). This template is the reusable proof pattern behind many adaptive backstepping and output-feedback constructions with unknown control direction.

For a time-varying unknown coefficient, a classical Nussbaum function is not automatically sufficient. Stronger function classes and modified lemmas are required. Similarly, actuator constraints require a separate well-posedness and bounded-command analysis; these topics belong to the next lesson and to the later chapter on constraints.

The implementations simulate

\[ \dot e=\left[-1.5+b\kappa^2\cos\kappa\right]e, \qquad \dot\kappa=e^2, \qquad e(0)=3, \quad \kappa(0)=0, \]

for \(b=+1\) and \(b=-1\). The expected qualitative observations are:

• the two cases can have very different transient excursions;
• \(\kappa(t)\) is nondecreasing because its derivative is \(e^2\);
• the adaptive search stops at a finite value once the accumulated error energy is finite;
• both error trajectories converge toward zero despite opposite values of \(b\);
• the negative-direction case may exhibit a larger control peak, illustrating that asymptotic stability does not imply mild transients.

A small integration step is important because the effective coefficient \(-\lambda+bN(\kappa)\) can change rapidly. The Python, C++, and Java files use a from-scratch fourth-order Runge–Kutta method. MATLAB uses ode45, and Wolfram Mathematica uses NDSolveValue.

This implementation uses NumPy for arrays and CSV output, Matplotlib for plots, and a transparent from-scratch RK4 integrator. For larger adaptive systems, scipy.integrate.solve_ivp is an appropriate alternative with adaptive step-size control.

"""Chapter16_Lesson3.py

Stability experiment for a scalar adaptive system with unknown control direction.
The reduced error model is
    e_dot = -lambda*e + b*u,
with unknown nonzero constant b. The Nussbaum-gain controller is
    u = N(kappa)*e,
    kappa_dot = e**2,
    N(kappa) = kappa**2*cos(kappa).

Dependencies:
    numpy, matplotlib
"""

from __future__ import annotations

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

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class SimulationConfig:
    b: float
    lambda_gain: float = 1.5
    e0: float = 3.0
    kappa0: float = 0.0
    final_time: float = 12.0
    step: float = 1.0e-3


def nussbaum(kappa: float) -> float:
    """Classical polynomial Nussbaum function N(kappa)=kappa^2 cos(kappa)."""
    return kappa * kappa * np.cos(kappa)


def dynamics(state: np.ndarray, cfg: SimulationConfig) -> np.ndarray:
    e, kappa = state
    control = nussbaum(kappa) * e
    return np.array(
        [
            -cfg.lambda_gain * e + cfg.b * control,
            e * e,
        ],
        dtype=float,
    )


def rk4_step(
    state: np.ndarray,
    step: float,
    rhs: Callable[[np.ndarray], np.ndarray],
) -> np.ndarray:
    k1 = rhs(state)
    k2 = rhs(state + 0.5 * step * k1)
    k3 = rhs(state + 0.5 * step * k2)
    k4 = rhs(state + step * k3)
    return state + (step / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(cfg: SimulationConfig) -> dict[str, np.ndarray]:
    sample_count = int(round(cfg.final_time / cfg.step)) + 1
    time = np.linspace(0.0, cfg.final_time, sample_count)
    state = np.empty((sample_count, 2), dtype=float)
    state[0] = np.array([cfg.e0, cfg.kappa0], dtype=float)

    for index in range(sample_count - 1):
        state[index + 1] = rk4_step(
            state[index], cfg.step, lambda x: dynamics(x, cfg)
        )
        if not np.all(np.isfinite(state[index + 1])):
            raise FloatingPointError(
                "Simulation diverged numerically. Reduce the integration step."
            )

    error = state[:, 0]
    kappa = state[:, 1]
    gain = nussbaum(kappa)
    control = gain * error
    lyapunov = 0.5 * error * error

    return {
        "time": time,
        "error": error,
        "kappa": kappa,
        "gain": gain,
        "control": control,
        "lyapunov": lyapunov,
    }


def save_csv(path: Path, data: dict[str, np.ndarray], b: float) -> None:
    matrix = np.column_stack(
        [
            data["time"],
            data["error"],
            data["kappa"],
            data["gain"],
            data["control"],
            data["lyapunov"],
        ]
    )
    header = "time,error,kappa,nussbaum_gain,control,lyapunov"
    np.savetxt(path, matrix, delimiter=",", header=header, comments="")
    print(f"Saved {path.name} for b={b:+g}")


def print_metrics(data: dict[str, np.ndarray], b: float) -> None:
    print(
        f"b={b:+g}: final |e|={abs(data['error'][-1]):.6e}, "
        f"max |e|={np.max(np.abs(data['error'])):.6f}, "
        f"final kappa={data['kappa'][-1]:.6f}, "
        f"max |u|={np.max(np.abs(data['control'])):.6f}"
    )


def main() -> None:
    output_dir = Path(__file__).resolve().parent
    results: dict[float, dict[str, np.ndarray]] = {}

    for b in (1.0, -1.0):
        cfg = SimulationConfig(b=b)
        data = simulate(cfg)
        results[b] = data
        suffix = "positive_b" if b > 0.0 else "negative_b"
        save_csv(output_dir / f"Chapter16_Lesson3_{suffix}.csv", data, b)
        print_metrics(data, b)

    plt.figure(figsize=(9, 5))
    for b, data in results.items():
        plt.plot(data["time"], data["error"], label=f"b={b:+g}")
    plt.xlabel("Time")
    plt.ylabel("Tracking error e(t)")
    plt.title("Nussbaum-Gain Stabilization for Both Control Directions")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter16_Lesson3_error.png", dpi=180)

    plt.figure(figsize=(9, 5))
    for b, data in results.items():
        plt.plot(data["time"], data["kappa"], label=f"b={b:+g}")
    plt.xlabel("Time")
    plt.ylabel("Adaptive variable kappa(t)")
    plt.title("Bounded Monotone Adaptive Variable")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter16_Lesson3_kappa.png", dpi=180)

    plt.show()


if __name__ == "__main__":
    main()

The C++17 program uses only the standard library and writes CSV files. In larger control projects, Eigen can be used for vector and matrix operations, while Boost.Odeint provides reusable numerical integrators.

// Chapter16_Lesson3.cpp
// From-scratch RK4 simulation of Nussbaum-gain stabilization.
// Build: g++ -std=c++17 -O2 Chapter16_Lesson3.cpp -o Chapter16_Lesson3

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

struct State {
    double error;
    double kappa;
};

struct Config {
    double b;
    double lambda_gain{1.5};
    double error0{3.0};
    double kappa0{0.0};
    double final_time{12.0};
    double step{1.0e-3};
};

double nussbaum(double kappa) {
    return kappa * kappa * std::cos(kappa);
}

State add(const State& a, const State& b, double scale = 1.0) {
    return {a.error + scale * b.error, a.kappa + scale * b.kappa};
}

State dynamics(const State& state, const Config& cfg) {
    const double control = nussbaum(state.kappa) * state.error;
    return {
        -cfg.lambda_gain * state.error + cfg.b * control,
        state.error * state.error
    };
}

State rk4Step(const State& state, const Config& cfg) {
    const State k1 = dynamics(state, cfg);
    const State k2 = dynamics(add(state, k1, 0.5 * cfg.step), cfg);
    const State k3 = dynamics(add(state, k2, 0.5 * cfg.step), cfg);
    const State k4 = dynamics(add(state, k3, cfg.step), cfg);

    return {
        state.error + (cfg.step / 6.0) *
            (k1.error + 2.0 * k2.error + 2.0 * k3.error + k4.error),
        state.kappa + (cfg.step / 6.0) *
            (k1.kappa + 2.0 * k2.kappa + 2.0 * k3.kappa + k4.kappa)
    };
}

void simulate(const Config& cfg, const std::string& filename) {
    std::ofstream output(filename);
    if (!output) {
        throw std::runtime_error("Cannot open output file: " + filename);
    }

    output << "time,error,kappa,nussbaum_gain,control,lyapunov\n";
    output << std::setprecision(15);

    State state{cfg.error0, cfg.kappa0};
    const long steps = static_cast<long>(std::llround(cfg.final_time / cfg.step));
    double max_error = std::abs(state.error);
    double max_control = 0.0;

    for (long i = 0; i <= steps; ++i) {
        const double time = i * cfg.step;
        const double gain = nussbaum(state.kappa);
        const double control = gain * state.error;
        const double lyapunov = 0.5 * state.error * state.error;

        output << time << ',' << state.error << ',' << state.kappa << ','
               << gain << ',' << control << ',' << lyapunov << '\n';

        max_error = std::max(max_error, std::abs(state.error));
        max_control = std::max(max_control, std::abs(control));

        if (i < steps) {
            state = rk4Step(state, cfg);
            if (!std::isfinite(state.error) || !std::isfinite(state.kappa)) {
                throw std::runtime_error(
                    "Numerical divergence. Reduce the integration step."
                );
            }
        }
    }

    std::cout << "b=" << std::showpos << cfg.b << std::noshowpos
              << ": final |e|=" << std::abs(state.error)
              << ", max |e|=" << max_error
              << ", final kappa=" << state.kappa
              << ", max |u|=" << max_control << '\n';
}

int main() {
    try {
        simulate(Config{1.0}, "Chapter16_Lesson3_positive_b_cpp.csv");
        simulate(Config{-1.0}, "Chapter16_Lesson3_negative_b_cpp.csv");
    } catch (const std::exception& exception) {
        std::cerr << "Error: " << exception.what() << '\n';
        return 1;
    }
    return 0;
}

The Java program uses records for immutable states and configurations, a from-scratch RK4 step, and standard NIO file output. Apache Commons Math provides production-grade ODE solvers when a larger numerical framework is needed.

// Chapter16_Lesson3.java
// From-scratch RK4 simulation of Nussbaum-gain stabilization.
// Build and run: javac Chapter16_Lesson3.java && java Chapter16_Lesson3

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 Chapter16_Lesson3 {
    private Chapter16_Lesson3() {}

    private record State(double error, double kappa) {}

    private record Config(
            double b,
            double lambdaGain,
            double error0,
            double kappa0,
            double finalTime,
            double step) {
        static Config forControlDirection(double b) {
            return new Config(b, 1.5, 3.0, 0.0, 12.0, 1.0e-3);
        }
    }

    private static double nussbaum(double kappa) {
        return kappa * kappa * Math.cos(kappa);
    }

    private static State add(State a, State b, double scale) {
        return new State(
                a.error() + scale * b.error(),
                a.kappa() + scale * b.kappa());
    }

    private static State dynamics(State state, Config cfg) {
        double control = nussbaum(state.kappa()) * state.error();
        return new State(
                -cfg.lambdaGain() * state.error() + cfg.b() * control,
                state.error() * state.error());
    }

    private static State rk4Step(State state, Config cfg) {
        State k1 = dynamics(state, cfg);
        State k2 = dynamics(add(state, k1, 0.5 * cfg.step()), cfg);
        State k3 = dynamics(add(state, k2, 0.5 * cfg.step()), cfg);
        State k4 = dynamics(add(state, k3, cfg.step()), cfg);

        return new State(
                state.error() + (cfg.step() / 6.0)
                        * (k1.error() + 2.0 * k2.error()
                        + 2.0 * k3.error() + k4.error()),
                state.kappa() + (cfg.step() / 6.0)
                        * (k1.kappa() + 2.0 * k2.kappa()
                        + 2.0 * k3.kappa() + k4.kappa()));
    }

    private static void simulate(Config cfg, Path outputPath) throws IOException {
        State state = new State(cfg.error0(), cfg.kappa0());
        long steps = Math.round(cfg.finalTime() / cfg.step());
        double maxError = Math.abs(state.error());
        double maxControl = 0.0;

        try (BufferedWriter writer = Files.newBufferedWriter(outputPath)) {
            writer.write("time,error,kappa,nussbaum_gain,control,lyapunov\n");

            for (long i = 0; i <= steps; i++) {
                double time = i * cfg.step();
                double gain = nussbaum(state.kappa());
                double control = gain * state.error();
                double lyapunov = 0.5 * state.error() * state.error();

                writer.write(String.format(
                        Locale.US,
                        "%.12f,%.12e,%.12e,%.12e,%.12e,%.12e%n",
                        time, state.error(), state.kappa(), gain, control, lyapunov));

                maxError = Math.max(maxError, Math.abs(state.error()));
                maxControl = Math.max(maxControl, Math.abs(control));

                if (i < steps) {
                    state = rk4Step(state, cfg);
                    if (!Double.isFinite(state.error())
                            || !Double.isFinite(state.kappa())) {
                        throw new ArithmeticException(
                                "Numerical divergence. Reduce the integration step.");
                    }
                }
            }
        }

        System.out.printf(
                Locale.US,
                "b=%+.0f: final |e|=%.6e, max |e|=%.6f, "
                        + "final kappa=%.6f, max |u|=%.6f%n",
                cfg.b(), Math.abs(state.error()), maxError,
                state.kappa(), maxControl);
    }

    public static void main(String[] args) {
        try {
            simulate(
                    Config.forControlDirection(1.0),
                    Path.of("Chapter16_Lesson3_positive_b_java.csv"));
            simulate(
                    Config.forControlDirection(-1.0),
                    Path.of("Chapter16_Lesson3_negative_b_java.csv"));
        } catch (IOException | ArithmeticException exception) {
            System.err.println("Error: " + exception.getMessage());
            System.exit(1);
        }
    }
}

MATLAB's ode45 integrates the two-state nonlinear closed loop. The same equations can be placed in a Simulink MATLAB Function block or implemented using Integrator blocks for \(e\) and \(\kappa\), a Product block for \(e^2\), and a MATLAB Function block for \(N(\kappa)\).

% Chapter16_Lesson3.m
% Nussbaum-gain stabilization for both unknown control directions.
% MATLAB uses ode45 with tight tolerances and exports numerical results.

clear; close all; clc;

lambda = 1.5;
e0 = 3.0;
kappa0 = 0.0;
tspan = [0, 12];
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-11, 'MaxStep', 1e-2);
controlDirections = [1, -1];
results = cell(size(controlDirections));

for index = 1:numel(controlDirections)
    b = controlDirections(index);
    rhs = @(t, z) chapter16Lesson3Dynamics(t, z, b, lambda);
    [time, state] = ode45(rhs, tspan, [e0; kappa0], options);

    errorSignal = state(:, 1);
    kappa = state(:, 2);
    gain = kappa.^2 .* cos(kappa);
    control = gain .* errorSignal;
    lyapunov = 0.5 .* errorSignal.^2;

    results{index} = struct( ...
        'b', b, ...
        'time', time, ...
        'error', errorSignal, ...
        'kappa', kappa, ...
        'gain', gain, ...
        'control', control, ...
        'lyapunov', lyapunov);

    outputTable = table(time, errorSignal, kappa, gain, control, lyapunov, ...
        'VariableNames', {'time', 'error', 'kappa', 'nussbaum_gain', ...
        'control', 'lyapunov'});

    if b > 0
        filename = 'Chapter16_Lesson3_positive_b_matlab.csv';
    else
        filename = 'Chapter16_Lesson3_negative_b_matlab.csv';
    end
    writetable(outputTable, filename);

    fprintf(['b=%+g: final |e|=%.6e, max |e|=%.6f, ', ...
        'final kappa=%.6f, max |u|=%.6f\n'], ...
        b, abs(errorSignal(end)), max(abs(errorSignal)), ...
        kappa(end), max(abs(control)));
end

figure('Name', 'Nussbaum-Gain Error');
hold on;
for index = 1:numel(results)
    plot(results{index}.time, results{index}.error, 'LineWidth', 1.4, ...
        'DisplayName', sprintf('b=%+g', results{index}.b));
end
xlabel('Time'); ylabel('Tracking error e(t)');
title('Nussbaum-Gain Stabilization for Both Control Directions');
grid on; legend('Location', 'best');
exportgraphics(gcf, 'Chapter16_Lesson3_error_matlab.png', 'Resolution', 180);

figure('Name', 'Nussbaum Adaptive Variable');
hold on;
for index = 1:numel(results)
    plot(results{index}.time, results{index}.kappa, 'LineWidth', 1.4, ...
        'DisplayName', sprintf('b=%+g', results{index}.b));
end
xlabel('Time'); ylabel('Adaptive variable kappa(t)');
title('Bounded Monotone Adaptive Variable');
grid on; legend('Location', 'best');
exportgraphics(gcf, 'Chapter16_Lesson3_kappa_matlab.png', 'Resolution', 180);

function dz = chapter16Lesson3Dynamics(~, z, b, lambda)
    errorSignal = z(1);
    kappa = z(2);
    nussbaumGain = kappa^2 * cos(kappa);
    control = nussbaumGain * errorSignal;

    dz = [ ...
        -lambda * errorSignal + b * control; ...
        errorSignal^2 ...
    ];
end

The Wolfram Language implementation uses NDSolveValue, exports CSV data, reports stability metrics, and creates error and adaptive-state plots. The source can be pasted directly into a notebook and saved as Chapter16_Lesson3.nb.

(* Chapter16_Lesson3.nb
   Wolfram Language source that can be pasted into a Mathematica notebook.
   It solves the Nussbaum-gain closed loop for b=+1 and b=-1. *)

ClearAll["Global`*"];

lambda = 1.5;
e0 = 3.0;
kappa0 = 0.0;
tFinal = 12.0;

nussbaum[k_] := k^2 Cos[k];

simulate[b_?NumericQ] := Module[
  {solution, errorFunction, kappaFunction, gainFunction, controlFunction,
   sampleTimes, data},

  solution = NDSolveValue[
    {
      e'[t] == -lambda e[t] + b nussbaum[kappa[t]] e[t],
      kappa'[t] == e[t]^2,
      e[0] == e0,
      kappa[0] == kappa0
    },
    {e, kappa},
    {t, 0, tFinal},
    Method -> {"TimeIntegration" -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5}},
    AccuracyGoal -> 10,
    PrecisionGoal -> 10,
    MaxStepFraction -> 1/2000
  ];

  errorFunction[t_] := solution[[1]][t];
  kappaFunction[t_] := solution[[2]][t];
  gainFunction[t_] := nussbaum[kappaFunction[t]];
  controlFunction[t_] := gainFunction[t] errorFunction[t];

  sampleTimes = Subdivide[0, tFinal, 12000];
  data = Table[
    {
      time,
      errorFunction[time],
      kappaFunction[time],
      gainFunction[time],
      controlFunction[time],
      1/2 errorFunction[time]^2
    },
    {time, sampleTimes}
  ];

  <|
    "b" -> b,
    "Solution" -> solution,
    "Data" -> data,
    "FinalError" -> Abs[errorFunction[tFinal]],
    "MaximumError" -> Max[Abs[data[[All, 2]]]],
    "FinalKappa" -> kappaFunction[tFinal],
    "MaximumControl" -> Max[Abs[data[[All, 5]]]]
  |>
];

positiveResult = simulate[1.0];
negativeResult = simulate[-1.0];
results = {positiveResult, negativeResult};

Scan[
  Function[result,
    Print@Row[{
      "b=", result["b"],
      ": final |e|=", ScientificForm[result["FinalError"], 6],
      ", max |e|=", NumberForm[result["MaximumError"], 6],
      ", final kappa=", NumberForm[result["FinalKappa"], 6],
      ", max |u|=", NumberForm[result["MaximumControl"], 6]
    }]
  ],
  results
];

Export[
  "Chapter16_Lesson3_positive_b_mathematica.csv",
  Prepend[positiveResult["Data"],
    {"time", "error", "kappa", "nussbaum_gain", "control", "lyapunov"}]
];

Export[
  "Chapter16_Lesson3_negative_b_mathematica.csv",
  Prepend[negativeResult["Data"],
    {"time", "error", "kappa", "nussbaum_gain", "control", "lyapunov"}]
];

errorPlot = Plot[
  Evaluate@Table[result["Solution"][[1]][t], {result, results}],
  {t, 0, tFinal},
  PlotLegends -> {"b=+1", "b=-1"},
  AxesLabel -> {"time", "e(t)"},
  PlotLabel -> "Nussbaum-Gain Stabilization for Both Control Directions",
  PlotRange -> All,
  ImageSize -> Large
];
Export["Chapter16_Lesson3_error_mathematica.png", errorPlot, ImageResolution -> 180];

kappaPlot = Plot[
  Evaluate@Table[result["Solution"][[2]][t], {result, results}],
  {t, 0, tFinal},
  PlotLegends -> {"b=+1", "b=-1"},
  AxesLabel -> {"time", "kappa(t)"},
  PlotLabel -> "Bounded Monotone Adaptive Variable",
  PlotRange -> All,
  ImageSize -> Large
];
Export["Chapter16_Lesson3_kappa_mathematica.png", kappaPlot, ImageResolution -> 180];

A numerical run should not be accepted merely because the final error is small. Verify the analytical invariants and convergence indicators:

1. Monotonicity of the adaptive state:

   \[ \kappa(t_{i+1})-\kappa(t_i)\geq-\varepsilon_{\mathrm{num}}. \]

2. Integral identity:

   \[ \kappa(t)-\kappa_0\approx\int_0^t e^2(\tau)\,d\tau. \]

3. Lyapunov identity residual:

   \[ r_V(t)=V(t)-V(0)+\lambda\bigl(\kappa(t)-\kappa_0\bigr) -b\int_{\kappa_0}^{\kappa(t)}N(s)\,ds. \]

   A converged integration should keep \(|r_V|\) small.
4. Tail error: compare the maximum value of \(|e(t)|\) on the final ten percent of the simulation with the initial error.
5. Step-size refinement: halve the integration step and confirm that final values and peak values change only slightly.

Problem 1 (Verifying a Nussbaum Function). Prove that \(N(\kappa)=\kappa^2\sin\kappa\) is a classical Nussbaum function.

Solution. Integration by parts gives

\[ \int_0^K s^2\sin s\,ds =-K^2\cos K+2K\sin K+2\cos K-2. \]

Hence

\[ \frac{1}{K}\int_0^K s^2\sin s\,ds =-K\cos K+2\sin K+\frac{2\cos K-2}{K}. \]

For \(K_n^+=(2n+1)\pi\), the leading term is \(+K_n^+\), so the average tends to \(+\infty\). For \(K_n^-=2n\pi\), the leading term is \(-K_n^-\), so the average tends to \(-\infty\). Both Nussbaum limits hold.

Problem 2 (Why Oscillation Alone Is Insufficient). Determine whether \(N(\kappa)=\sin\kappa\) is a Nussbaum function.

Solution.

\[ \frac{1}{K}\int_0^K\sin s\,ds =\frac{1-\cos K}{K}. \]

The numerator is bounded and the denominator grows without bound, so the expression tends to zero. Therefore its limsup and liminf are both zero. The function changes sign infinitely often, but it does not have the required unbounded integral-average excursions.

Problem 3 (Boundedness Before Convergence). For the closed loop in Theorem 2, prove that \(u\in L_\infty\) before invoking Barbalat's lemma.

Solution. The Nussbaum lemma first establishes boundedness of \(\kappa\) and \(V=e^2/2\). Therefore \(e\in L_\infty\). Continuity of \(N\) on the compact interval containing the trajectory of \(\kappa\) implies \(N(\kappa)\in L_\infty\). Their product \(u=N(\kappa)e\) is consequently bounded. This ordering is essential because bounded \(u\) is used to establish bounded \(\dot e\), which supplies uniform continuity.

Problem 4 (A General Dissipation Term). Suppose

\[ \dot V\leq-c|z|^2+\bigl[\alpha+bN(\kappa)\bigr]\dot\kappa, \qquad c>0, \qquad \dot\kappa\geq0. \]

Show that \(z\in L_2\) once the Nussbaum lemma has established boundedness of \(V\) and \(\kappa\).

Solution. Integrating gives

\[ V(t)+c\int_0^t|z(\tau)|^2\,d\tau \leq V(0)+\alpha(\kappa(t)-\kappa_0) +b\int_{\kappa_0}^{\kappa(t)}N(s)\,ds. \]

The right-hand side is bounded because \(\kappa\) is bounded and \(N\) is continuous. Since \(V(t)\geq0\), the integral of \(|z|^2\) is bounded uniformly in time. Thus \(z\in L_2\).

Problem 5 (Failure Under Sign Switching). Explain why Theorem 2 cannot be directly applied when \(b=b(t)\) changes sign.

Solution. The proof converts the term \(bN(\kappa)\dot\kappa\) into the state integral

\[ b\int_{\kappa_0}^{\kappa(t)}N(s)\,ds, \]

which requires \(b\) to be constant. When \(b=b(t)\), the integral is instead \(\int_0^t b(\tau)N(\kappa(\tau))\dot\kappa(\tau)d\tau\), and it cannot generally be rewritten as a fixed scalar times an integral in \(\kappa\). Arbitrary sign switching can cancel the favorable Nussbaum excursions. Additional assumptions and stronger Nussbaum-function constructions are required.

Problem 6 (Programming Exercise). Create Chapter16_Lesson3_Ex1.py by adding three automated checks to the Python program: monotonicity of \(\kappa\), the numerical identity \(\kappa(t)-\kappa_0\approx\int_0^t e^2d\tau\), and a step-refinement comparison.

Solution outline. Use np.diff(kappa) for monotonicity, a trapezoidal quadrature for accumulated error energy, and two simulations with steps 1e-3 and 5e-4. Report tolerances rather than testing exact equality, because all three checks are numerical.

Nussbaum-gain stability is obtained through an integral contradiction, not by forcing the Lyapunov derivative to remain negative. For the reduced unknown-direction channel, the adaptive law \(\dot\kappa=e^2\) makes the Lyapunov derivative integrable with respect to \(\kappa\). The unbounded positive and negative averages of \(N\) prevent \(\kappa\) from escaping to infinity for either sign of the unknown constant gain. Boundedness then gives \(e\in L_2\cap L_\infty\), bounded \(\dot e\), and finally \(e(t)\to0\) by Barbalat's lemma. The proof does not imply gain identification, exponential convergence, tolerance of arbitrary sign switching, or automatic compatibility with saturation.