A numerical simulation can fail in at least three different ways: (i) the discretization is unstable for the chosen step size, (ii) the discretization is stable but inaccurate because truncation errors are too large, or (iii) the numerical solution is accurate for the wrong model. In this lesson we separate these issues explicitly.

Let \( \dot{\mathbf{x} } = \mathbf{f}(t,\mathbf{x}) \) with initial condition \( \mathbf{x}(t_0)=\mathbf{x}_0 \). A one-step method defines a map

\[ \mathbf{x}_{n+1} = \Phi_h(t_n,\mathbf{x}_n), \qquad t_{n+1}=t_n+h_n. \]

Reliability requires:

• Numerical stability: perturbations and roundoff should not amplify unphysically.
• Error control: local error indicators should keep global error within tolerance.
• Verification: the computed solution should satisfy expected order and consistency checks.

In control-oriented simulation, failure to enforce these conditions can produce false conclusions about damping, overshoot, settling time, or even stability of the physical system.

The canonical stability probe is the scalar linear test equation \( y'=\lambda y \), where \( \operatorname{Re}(\lambda) < 0 \). A one-step method induces

\[ y_{n+1} = R(z)\,y_n, \qquad z = h\lambda, \]

where \( R(z) \) is the method's stability function. The absolute stability region is

\[ \mathcal{S} = \{ z \in \mathbb{C} : |R(z)| \leq 1 \}. \]

For explicit Euler, \( R(z)=1+z \), so stability on the negative real axis requires

\[ |1+h\lambda| \leq 1. \]

If \( \lambda=-\alpha \) with \( \alpha > 0 \), then

\[ |1-\alpha h| \leq 1 \quad \Longrightarrow \quad 0 \leq h \leq \frac{2}{\alpha}. \]

This is a purely numerical restriction; the continuous system is asymptotically stable for all time, but the discrete scheme becomes unstable when \( h \) is too large.

flowchart TD
  A["Continuous model xdot = f(t,x)"] --> B["Choose integration method"]
  B --> C["Form test equation y' = lambda y"]
  C --> D["Compute stability function R(z)"]
  D --> E["Check |R(h lambda)| <= 1"]
  E -->|yes| F["Step size is numerically stable"]
  E -->|no| G["Reduce h or change solver"]
        

For a method of order \( p \), the local truncation error (LTE) satisfies

\[ \boldsymbol{\tau}_{n+1} = \frac{\mathbf{x}(t_{n+1}) - \Phi_h(t_n,\mathbf{x}(t_n))}{h} = \mathcal{O}(h^p), \]

equivalently the one-step defect is \( \mathcal{O}(h^{p+1}) \). Under Lipschitz continuity of \( \mathbf{f} \) and zero-stability, the global error \( \mathbf{e}_n = \mathbf{x}(t_n)-\mathbf{x}_n \) satisfies

\[ \|\mathbf{e}_n\| \leq C\left(e^{L(t_n-t_0)}-1\right) h^p, \]

for sufficiently small constant step \( h \), where \( L \) is a Lipschitz constant of \( \mathbf{f} \). Thus a stable order-\( p \) method gives global convergence of order \( p \).

In practice, adaptive solvers estimate the local error and choose \( h_n \) to keep a weighted norm below 1:

\[ \mathrm{err}_n = \sqrt{\frac{1}{m}\sum_{i=1}^m \left(\frac{\widehat{e}_{n,i} }{\mathrm{ATOL}_i + \mathrm{RTOL}_i\max(|x_{n,i}|,|x_{n+1,i}|)}\right)^2 }. \]

Accept the step when \( \mathrm{err}_n \leq 1 \); otherwise reject and retry with smaller \( h_n \).

Two common error estimators are:

1. Embedded pairs (e.g., RK45): compute two approximations of different orders from shared stages.
2. Step-doubling: compare one step of size \( h \) against two steps of size \( h/2 \).

For a method of order \( p \), denote by \( \mathbf{x}^{[h]}_{n+1} \) the one-step result and \( \mathbf{x}^{[h/2]}_{n+1} \) the two-half-step result. Then

\[ \mathbf{x}^{[h]}_{n+1} = \mathbf{x}(t_{n+1}) + C h^{p+1} + \mathcal{O}(h^{p+2}), \]

\[ \mathbf{x}^{[h/2]}_{n+1} = \mathbf{x}(t_{n+1}) + C \left(\frac{h}{2}\right)^{p+1}\cdot 2 + \mathcal{O}(h^{p+2}) = \mathbf{x}(t_{n+1}) + \frac{C}{2^p} h^{p+1} + \mathcal{O}(h^{p+2}). \]

Subtracting yields an asymptotically valid error estimator:

\[ \widehat{\mathbf{e} }_{n+1} \approx \frac{\mathbf{x}^{[h/2]}_{n+1} - \mathbf{x}^{[h]}_{n+1} }{2^p-1}. \]

A Richardson-improved accepted value is then

\[ \mathbf{x}^{\mathrm{rich} }_{n+1} = \mathbf{x}^{[h/2]}_{n+1} + \frac{\mathbf{x}^{[h/2]}_{n+1} - \mathbf{x}^{[h]}_{n+1} }{2^p-1}. \]

If the local error behaves like \( \|\widehat{\mathbf{e} }\| \approx K h^{p+1} \), then the asymptotically optimal next step size is

\[ h_{n+1} = h_n \,\gamma \left(\frac{1}{\mathrm{err}_n}\right)^{\frac{1}{p+1} }, \]

where \( \gamma \in (0,1) \) is a safety factor (commonly 0.8 to 0.95). In practice, bounds are imposed:

\[ h_{n+1} = h_n \cdot \min\!\Big(r_{\max}, \max\!\big(r_{\min}, \gamma\,\mathrm{err}_n^{-1/(p+1)}\big)\Big). \]

For control simulations, component-wise tolerances should reflect engineering scales: for example, a position state may require micrometer absolute tolerance, while a thermal state may tolerate larger absolute error but tighter relative error.

An important caveat: local error control does not guarantee model correctness. It only ensures numerical consistency relative to the chosen differential equation.

Verification answers: Are we solving the equations right? (as opposed to validation, which asks whether the equations represent reality). For system dynamics courses, a practical verification toolkit includes:

1. Order verification: estimate observed order \( p_{\mathrm{obs} } \) from successive step refinements.
2. Richardson extrapolation: estimate a higher-accuracy surrogate solution.
3. Invariant / dissipation checks: energy conservation or monotone energy decay when physically expected.
4. Residual checks: substitute the numerical trajectory into the ODE and inspect defect size.
5. Cross-solver consistency: compare explicit and stiff-capable solvers on the same problem.

For step sizes \( h, h/2, h/4 \), define errors (or proxies) \( E_h, E_{h/2}, E_{h/4} \). The observed order is estimated by

\[ p_{\mathrm{obs} } = \frac{\log(E_h/E_{h/2})}{\log 2}. \]

More generally, using three solutions at the same final time, \( \mathbf{x}_h, \mathbf{x}_{h/r}, \mathbf{x}_{h/r^2} \) with refinement ratio \( r \):

\[ p_{\mathrm{obs} } \approx \frac{\log \|\mathbf{x}_h-\mathbf{x}_{h/r}\| - \log \|\mathbf{x}_{h/r}-\mathbf{x}_{h/r^2}\|}{\log r}. \]

flowchart TD
  A["Build simulation model"] --> B["Run baseline solver"]
  B --> C["Refine step size (h, h/2, h/4)"]
  C --> D["Compute error proxy and p_obs"]
  D --> E["Check physics: energy, bounds, units"]
  E --> F["Cross-check with another solver"]
  F --> G["Verified numerical implementation"]
        

Consider the second-order damped oscillator \( \ddot{q} + 2\zeta\omega_n \dot{q} + \omega_n^2 q = 0 \). Defining state variables \( x_1=q \) and \( x_2=\dot{q} \) gives

\[ \dot{\mathbf{x} } = \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = \begin{bmatrix} x_2 \\ -2\zeta\omega_n x_2 - \omega_n^2 x_1 \end{bmatrix}. \]

For \( 0 \leq \zeta < 1 \), the exact displacement is

\[ q(t) = e^{-\zeta\omega_n t} \left( A \cos(\omega_d t) + B \sin(\omega_d t) \right), \qquad \omega_d = \omega_n\sqrt{1-\zeta^2}. \]

A useful physics-based diagnostic is the mechanical energy surrogate

\[ E(t) = \frac{1}{2}\dot{q}(t)^2 + \frac{1}{2}\omega_n^2 q(t)^2. \]

Differentiating along trajectories gives

\[ \dot{E}(t) = \dot{q}\ddot{q} + \omega_n^2 q\dot{q} = \dot{q}\left(\ddot{q}+\omega_n^2 q\right) = -2\zeta\omega_n \dot{q}^2 \leq 0. \]

Therefore, a numerically stable and sufficiently accurate simulation should produce a nonincreasing energy sequence (up to small tolerance and roundoff).

The following script demonstrates: (i) explicit Euler stability on the test equation, (ii) RK4 convergence with observed order, and (iii) adaptive RK4 with step-doubling error control and energy monotonicity checks for a damped oscillator.

Code title: Chapter15_Lesson5.py

# Chapter15_Lesson5.py
"""
Numerical Stability, Error Control, and Model Verification via Simulation
Python implementation for Chapter 15, Lesson 5 (System Dynamics)

Features
--------
1) Absolute-stability demonstration on the test equation y' = lambda y
2) Adaptive RK4 with step-doubling error control
3) Verification by Richardson extrapolation and observed order estimation
4) Physical sanity check on a damped oscillator (monotone energy decay)
"""

from __future__ import annotations
import math
from dataclasses import dataclass
from typing import Callable, List, Tuple

import numpy as np


Vector = np.ndarray


def rk4_step(f: Callable[[float, Vector], Vector], t: float, y: Vector, h: float) -> Vector:
    k1 = f(t, y)
    k2 = f(t + 0.5 * h, y + 0.5 * h * k1)
    k3 = f(t + 0.5 * h, y + 0.5 * h * k2)
    k4 = f(t + h, y + h * k3)
    return y + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)


@dataclass
class AdaptiveResult:
    t: np.ndarray
    y: np.ndarray
    accepted_steps: int
    rejected_steps: int


def adaptive_rk4_stepdoubling(
    f: Callable[[float, Vector], Vector],
    t0: float,
    y0: Vector,
    tf: float,
    h0: float = 0.05,
    rtol: float = 1e-6,
    atol: float = 1e-9,
    h_min: float = 1e-8,
    h_max: float = 0.5,
) -> AdaptiveResult:
    """
    Adaptive RK4 using step-doubling:
      y_big  = RK4(t, y, h)
      y_half = RK4(t+h/2, RK4(t, y, h/2), h/2)
    local error estimate: e ~= (y_half - y_big)/(2^p - 1), p=4
    """
    p = 4
    t = t0
    y = np.array(y0, dtype=float)
    h = h0

    ts: List[float] = [t]
    ys: List[Vector] = [y.copy()]
    accepted = 0
    rejected = 0

    while t < tf:
        if t + h > tf:
            h = tf - t

        y_big = rk4_step(f, t, y, h)
        y_half_1 = rk4_step(f, t, y, 0.5 * h)
        y_half_2 = rk4_step(f, t + 0.5 * h, y_half_1, 0.5 * h)

        err_est = (y_half_2 - y_big) / (2**p - 1)
        scale = atol + rtol * np.maximum(np.abs(y_half_2), np.abs(y))
        err_norm = float(np.sqrt(np.mean((err_est / scale) ** 2)))

        if err_norm <= 1.0:
            # Richardson-improved accepted value
            y = y_half_2 + err_est
            t = t + h
            ts.append(t)
            ys.append(y.copy())
            accepted += 1

            # PI-like update (simple variant)
            if err_norm == 0.0:
                factor = 2.0
            else:
                factor = 0.9 * (1.0 / err_norm) ** (1.0 / (p + 1))
            h = min(h_max, max(h_min, h * min(2.0, max(0.3, factor))))
        else:
            rejected += 1
            factor = 0.9 * (1.0 / max(err_norm, 1e-16)) ** (1.0 / (p + 1))
            h = max(h_min, h * max(0.1, min(0.5, factor)))
            if h <= h_min:
                raise RuntimeError("Step size underflow during adaptive integration.")

    return AdaptiveResult(
        t=np.array(ts),
        y=np.vstack(ys),
        accepted_steps=accepted,
        rejected_steps=rejected,
    )


def fixed_rk4(
    f: Callable[[float, Vector], Vector], t0: float, y0: Vector, tf: float, h: float
) -> Tuple[np.ndarray, np.ndarray]:
    n = int(round((tf - t0) / h))
    t = t0
    y = np.array(y0, dtype=float)
    ts = [t]
    ys = [y.copy()]
    for _ in range(n):
        y = rk4_step(f, t, y, h)
        t += h
        ts.append(t)
        ys.append(y.copy())
    return np.array(ts), np.vstack(ys)


def test_equation_stability_demo(lam: complex, h_values: List[float], n_steps: int = 50) -> None:
    """
    Demonstrate growth/decay on y' = lambda y with explicit Euler:
      y_{n+1} = (1 + h*lambda) y_n
    """
    print("\n=== Absolute stability demo (explicit Euler on y' = lambda y) ===")
    for h in h_values:
        amp = abs(1.0 + h * lam)
        y = 1.0 + 0j
        for _ in range(n_steps):
            y = y + h * lam * y
        print(f"h={h:8.4f}, |1+h*lambda|={amp:10.6f}, |y_N|={abs(y):12.6e}")


def damped_oscillator_rhs(omega_n: float, zeta: float) -> Callable[[float, Vector], Vector]:
    def f(_t: float, x: Vector) -> Vector:
        # x = [q, v]
        q, v = x
        dq = v
        dv = -2.0 * zeta * omega_n * v - (omega_n**2) * q
        return np.array([dq, dv], dtype=float)

    return f


def exact_damped_position(t: np.ndarray, omega_n: float, zeta: float, q0: float, v0: float) -> np.ndarray:
    # underdamped closed form
    if not (0.0 <= zeta < 1.0):
        raise ValueError("This exact formula is coded for 0 <= zeta < 1.")
    wd = omega_n * math.sqrt(1.0 - zeta * zeta)
    A = q0
    B = (v0 + zeta * omega_n * q0) / wd
    return np.exp(-zeta * omega_n * t) * (A * np.cos(wd * t) + B * np.sin(wd * t))


def energy(x: np.ndarray, omega_n: float) -> np.ndarray:
    q = x[:, 0]
    v = x[:, 1]
    return 0.5 * (v**2 + (omega_n**2) * q**2)


def convergence_verification() -> None:
    print("\n=== Verification via Richardson extrapolation (RK4) ===")
    omega_n = 4.0
    zeta = 0.1
    q0, v0 = 1.0, 0.0
    tf = 5.0
    f = damped_oscillator_rhs(omega_n, zeta)

    hs = [0.2, 0.1, 0.05, 0.025]
    errors = []
    for h in hs:
        t, x = fixed_rk4(f, 0.0, np.array([q0, v0]), tf, h)
        q_exact = exact_damped_position(t, omega_n, zeta, q0, v0)
        err_inf = float(np.max(np.abs(x[:, 0] - q_exact)))
        errors.append(err_inf)
        print(f"h={h:7.4f}, max|q-q_exact|={err_inf:12.6e}")

    print("\nObserved order estimates:")
    for i in range(len(errors) - 1):
        p_obs = math.log(errors[i] / errors[i + 1], 2.0)
        print(f"p_obs(h={hs[i]} -> {hs[i+1]}) = {p_obs:.4f}")

    # Richardson extrapolated final-state estimate using h and h/2
    h = 0.1
    _, xh = fixed_rk4(f, 0.0, np.array([q0, v0]), tf, h)
    _, xh2 = fixed_rk4(f, 0.0, np.array([q0, v0]), tf, h / 2.0)
    q_rich = xh2[-1, 0] + (xh2[-1, 0] - xh[-1, 0]) / (2**4 - 1)
    q_exact_tf = exact_damped_position(np.array([tf]), omega_n, zeta, q0, v0)[0]
    print(f"\nFinal-time q(T) exact          = {q_exact_tf:.10f}")
    print(f"Final-time q(T) RK4 h/2        = {xh2[-1,0]:.10f}")
    print(f"Final-time q(T) Richardson est = {q_rich:.10f}")
    print(f"Richardson abs error           = {abs(q_rich - q_exact_tf):.3e}")


def adaptive_demo() -> None:
    print("\n=== Adaptive RK4 + error control on damped oscillator ===")
    omega_n = 4.0
    zeta = 0.05
    q0, v0 = 1.0, 0.0
    tf = 12.0

    f = damped_oscillator_rhs(omega_n, zeta)
    res = adaptive_rk4_stepdoubling(
        f=f,
        t0=0.0,
        y0=np.array([q0, v0]),
        tf=tf,
        h0=0.1,
        rtol=1e-6,
        atol=1e-9,
        h_max=0.2,
    )

    E = energy(res.y, omega_n)
    energy_nonincreasing = bool(np.all(np.diff(E) <= 1e-8))

    print(f"Accepted steps: {res.accepted_steps}")
    print(f"Rejected steps: {res.rejected_steps}")
    print(f"Final state    : q={res.y[-1,0]: .6f}, v={res.y[-1,1]: .6f}")
    print(f"Energy check   : nonincreasing={energy_nonincreasing}")

    # Compare with exact displacement at the adaptive nodes
    q_exact = exact_damped_position(res.t, omega_n, zeta, q0, v0)
    err_inf = float(np.max(np.abs(res.y[:, 0] - q_exact)))
    print(f"max|q-q_exact| on adaptive mesh = {err_inf:.3e}")


def main() -> None:
    # Stability demo on a stiff-ish decay mode
    lam = -50.0 + 0.0j
    test_equation_stability_demo(lam=lam, h_values=[0.01, 0.03, 0.05, 0.06], n_steps=30)

    # Verification and adaptive control demos
    convergence_verification()
    adaptive_demo()


if __name__ == "__main__":
    main()

This C++17 implementation mirrors the Python workflow and is suitable for high-performance numerical simulation exercises in control labs.

Code title: Chapter15_Lesson5.cpp

// Chapter15_Lesson5.cpp
// Numerical Stability, Error Control, and Model Verification via Simulation
// C++17 implementation for Chapter 15, Lesson 5 (System Dynamics)

#include <array>
#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <utility>
#include <vector>

using State = std::array<double, 2>;
using RHS = std::function<State(double, const State&)>;

State add(const State& a, const State& b) {
    return {a[0] + b[0], a[1] + b[1]};
}

State scale(double s, const State& a) {
    return {s * a[0], s * a[1]};
}

State rk4_step(const RHS& f, double t, const State& y, double h) {
    State k1 = f(t, y);
    State k2 = f(t + 0.5 * h, add(y, scale(0.5 * h, k1)));
    State k3 = f(t + 0.5 * h, add(y, scale(0.5 * h, k2)));
    State k4 = f(t + h, add(y, scale(h, k3)));
    State out{};
    for (int i = 0; i < 2; ++i) {
        out[i] = y[i] + (h / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
    }
    return out;
}

struct AdaptiveResult {
    std::vector<double> t;
    std::vector<State> y;
    int accepted = 0;
    int rejected = 0;
};

AdaptiveResult adaptive_rk4_stepdoubling(
    const RHS& f, double t0, const State& y0, double tf,
    double h0, double rtol, double atol
) {
    const int p = 4;
    double t = t0;
    State y = y0;
    double h = h0;
    const double hMin = 1e-8;
    const double hMax = 0.2;

    AdaptiveResult out;
    out.t.push_back(t);
    out.y.push_back(y);

    while (t < tf) {
        if (t + h > tf) h = tf - t;

        State yBig = rk4_step(f, t, y, h);
        State yHalf1 = rk4_step(f, t, y, 0.5 * h);
        State yHalf2 = rk4_step(f, t + 0.5 * h, yHalf1, 0.5 * h);

        State err{};
        double accum = 0.0;
        for (int i = 0; i < 2; ++i) {
            err[i] = (yHalf2[i] - yBig[i]) / (std::pow(2.0, p) - 1.0);
            double sc = atol + rtol * std::max(std::abs(yHalf2[i]), std::abs(y[i]));
            double z = err[i] / sc;
            accum += z * z;
        }
        double errNorm = std::sqrt(accum / 2.0);

        if (errNorm <= 1.0) {
            for (int i = 0; i < 2; ++i) y[i] = yHalf2[i] + err[i];
            t += h;
            out.t.push_back(t);
            out.y.push_back(y);
            out.accepted++;

            double factor = (errNorm == 0.0) ? 2.0 : 0.9 * std::pow(1.0 / errNorm, 1.0 / (p + 1.0));
            factor = std::max(0.3, std::min(2.0, factor));
            h = std::max(hMin, std::min(hMax, h * factor));
        } else {
            out.rejected++;
            double factor = 0.9 * std::pow(1.0 / std::max(errNorm, 1e-16), 1.0 / (p + 1.0));
            factor = std::max(0.1, std::min(0.5, factor));
            h = std::max(hMin, h * factor);
            if (h <= hMin) throw std::runtime_error("Step size underflow.");
        }
    }
    return out;
}

RHS damped_oscillator_rhs(double omegaN, double zeta) {
    return [omegaN, zeta](double /*t*/, const State& x) -> State {
        double q = x[0], v = x[1];
        return {v, -2.0 * zeta * omegaN * v - omegaN * omegaN * q};
    };
}

double exact_q(double t, double omegaN, double zeta, double q0, double v0) {
    double wd = omegaN * std::sqrt(1.0 - zeta * zeta);
    double A = q0;
    double B = (v0 + zeta * omegaN * q0) / wd;
    return std::exp(-zeta * omegaN * t) * (A * std::cos(wd * t) + B * std::sin(wd * t));
}

double energy(const State& x, double omegaN) {
    return 0.5 * (x[1] * x[1] + omegaN * omegaN * x[0] * x[0]);
}

void euler_stability_demo(double lambda, const std::vector<double>& hs) {
    std::cout << "=== Explicit Euler absolute-stability demo ===\n";
    for (double h : hs) {
        double amp = std::abs(1.0 + h * lambda);
        double y = 1.0;
        for (int n = 0; n < 30; ++n) y = y + h * lambda * y;
        std::cout << "h=" << std::setw(6) << h
                  << ", |1+h*lambda|=" << std::setw(10) << amp
                  << ", |y_N|=" << std::scientific << y << std::defaultfloat << "\n";
    }
    std::cout << "\n";
}

void convergence_demo() {
    std::cout << "=== RK4 convergence verification ===\n";
    double omegaN = 4.0, zeta = 0.1, tf = 5.0;
    State y0{1.0, 0.0};
    RHS f = damped_oscillator_rhs(omegaN, zeta);

    std::vector<double> hs{0.2, 0.1, 0.05, 0.025};
    std::vector<double> errs;

    for (double h : hs) {
        int n = static_cast<int>(std::round(tf / h));
        double t = 0.0;
        State y = y0;
        double maxErr = 0.0;
        for (int k = 0; k <= n; ++k) {
            double qEx = exact_q(t, omegaN, zeta, y0[0], y0[1]);
            maxErr = std::max(maxErr, std::abs(y[0] - qEx));
            if (k < n) {
                y = rk4_step(f, t, y, h);
                t += h;
            }
        }
        errs.push_back(maxErr);
        std::cout << "h=" << h << ", max|q-q_exact|=" << std::scientific << maxErr << std::defaultfloat << "\n";
    }

    std::cout << "Observed orders:\n";
    for (size_t i = 0; i + 1 < errs.size(); ++i) {
        double pObs = std::log(errs[i] / errs[i + 1]) / std::log(2.0);
        std::cout << "p_obs(" << hs[i] << "->" << hs[i + 1] << ")=" << pObs << "\n";
    }
    std::cout << "\n";
}

void adaptive_demo() {
    std::cout << "=== Adaptive RK4 (step-doubling) demo ===\n";
    double omegaN = 4.0, zeta = 0.05, tf = 12.0;
    State y0{1.0, 0.0};
    RHS f = damped_oscillator_rhs(omegaN, zeta);

    AdaptiveResult r = adaptive_rk4_stepdoubling(f, 0.0, y0, tf, 0.1, 1e-6, 1e-9);

    bool monotoneEnergy = true;
    double prevE = energy(r.y.front(), omegaN);
    for (size_t i = 1; i < r.y.size(); ++i) {
        double Ei = energy(r.y[i], omegaN);
        if (Ei - prevE > 1e-8) monotoneEnergy = false;
        prevE = Ei;
    }

    double maxErr = 0.0;
    for (size_t i = 0; i < r.t.size(); ++i) {
        double qEx = exact_q(r.t[i], omegaN, zeta, y0[0], y0[1]);
        maxErr = std::max(maxErr, std::abs(r.y[i][0] - qEx));
    }

    std::cout << "Accepted=" << r.accepted << ", Rejected=" << r.rejected << "\n";
    std::cout << "Final state q=" << r.y.back()[0] << ", v=" << r.y.back()[1] << "\n";
    std::cout << "Energy nonincreasing = " << (monotoneEnergy ? "true" : "false") << "\n";
    std::cout << "max|q-q_exact| on adaptive mesh = " << std::scientific << maxErr << std::defaultfloat << "\n";
}

int main() {
    euler_stability_demo(-50.0, {0.01, 0.03, 0.05, 0.06});
    convergence_demo();
    adaptive_demo();
    return 0;
}

The Java version is useful for students implementing simulation pipelines in strongly typed object-oriented environments.

Code title: Chapter15_Lesson5.java

// Chapter15_Lesson5.java
// Numerical Stability, Error Control, and Model Verification via Simulation
// Java implementation for Chapter 15, Lesson 5 (System Dynamics)

import java.util.ArrayList;
import java.util.List;

public class Chapter15_Lesson5 {

    interface RHS {
        double[] eval(double t, double[] y);
    }

    static double[] rk4Step(RHS f, double t, double[] y, double h) {
        double[] k1 = f.eval(t, y);
        double[] y2 = new double[]{y[0] + 0.5*h*k1[0], y[1] + 0.5*h*k1[1]};
        double[] k2 = f.eval(t + 0.5*h, y2);
        double[] y3 = new double[]{y[0] + 0.5*h*k2[0], y[1] + 0.5*h*k2[1]};
        double[] k3 = f.eval(t + 0.5*h, y3);
        double[] y4 = new double[]{y[0] + h*k3[0], y[1] + h*k3[1]};
        double[] k4 = f.eval(t + h, y4);

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

    static class AdaptiveResult {
        List<Double> t = new ArrayList<>();
        List<double[]> y = new ArrayList<>();
        int accepted = 0;
        int rejected = 0;
    }

    static AdaptiveResult adaptiveRK4StepDoubling(
            RHS f, double t0, double[] y0, double tf, double h0, double rtol, double atol) {

        int p = 4;
        double t = t0;
        double[] y = new double[]{y0[0], y0[1]};
        double h = h0;
        double hMin = 1e-8;
        double hMax = 0.2;

        AdaptiveResult out = new AdaptiveResult();
        out.t.add(t);
        out.y.add(new double[]{y[0], y[1]});

        while (t < tf) {
            if (t + h > tf) h = tf - t;

            double[] yBig = rk4Step(f, t, y, h);
            double[] yHalf1 = rk4Step(f, t, y, 0.5*h);
            double[] yHalf2 = rk4Step(f, t + 0.5*h, yHalf1, 0.5*h);

            double[] err = new double[2];
            double accum = 0.0;
            for (int i = 0; i < 2; i++) {
                err[i] = (yHalf2[i] - yBig[i]) / (Math.pow(2.0, p) - 1.0);
                double sc = atol + rtol * Math.max(Math.abs(yHalf2[i]), Math.abs(y[i]));
                double z = err[i] / sc;
                accum += z*z;
            }
            double errNorm = Math.sqrt(accum / 2.0);

            if (errNorm <= 1.0) {
                for (int i = 0; i < 2; i++) y[i] = yHalf2[i] + err[i];
                t += h;
                out.t.add(t);
                out.y.add(new double[]{y[0], y[1]});
                out.accepted++;

                double factor = (errNorm == 0.0) ? 2.0 :
                        0.9 * Math.pow(1.0 / errNorm, 1.0 / (p + 1.0));
                factor = Math.max(0.3, Math.min(2.0, factor));
                h = Math.max(hMin, Math.min(hMax, h * factor));
            } else {
                out.rejected++;
                double factor = 0.9 * Math.pow(1.0 / Math.max(errNorm, 1e-16), 1.0 / (p + 1.0));
                factor = Math.max(0.1, Math.min(0.5, factor));
                h = Math.max(hMin, h * factor);
                if (h <= hMin) throw new RuntimeException("Step size underflow.");
            }
        }
        return out;
    }

    static RHS dampedOscillatorRHS(final double omegaN, final double zeta) {
        return (t, x) -> new double[]{
            x[1],
            -2.0 * zeta * omegaN * x[1] - omegaN * omegaN * x[0]
        };
    }

    static double exactQ(double t, double omegaN, double zeta, double q0, double v0) {
        double wd = omegaN * Math.sqrt(1.0 - zeta*zeta);
        double A = q0;
        double B = (v0 + zeta * omegaN * q0) / wd;
        return Math.exp(-zeta * omegaN * t) * (A * Math.cos(wd * t) + B * Math.sin(wd * t));
    }

    static double energy(double[] x, double omegaN) {
        return 0.5 * (x[1]*x[1] + omegaN*omegaN*x[0]*x[0]);
    }

    static void eulerStabilityDemo(double lambda) {
        System.out.println("=== Explicit Euler stability demo ===");
        double[] hs = {0.01, 0.03, 0.05, 0.06};
        for (double h : hs) {
            double amp = Math.abs(1.0 + h * lambda);
            double y = 1.0;
            for (int n = 0; n < 30; n++) y = y + h * lambda * y;
            System.out.printf("h=%.3f, |1+h*lambda|=%.6f, y_N=%e%n", h, amp, y);
        }
        System.out.println();
    }

    static void convergenceDemo() {
        System.out.println("=== RK4 convergence verification ===");
        double omegaN = 4.0, zeta = 0.1, tf = 5.0;
        double[] y0 = {1.0, 0.0};
        RHS f = dampedOscillatorRHS(omegaN, zeta);
        double[] hs = {0.2, 0.1, 0.05, 0.025};
        double[] errs = new double[hs.length];

        for (int i = 0; i < hs.length; i++) {
            double h = hs[i];
            int n = (int)Math.round(tf / h);
            double t = 0.0;
            double[] y = {y0[0], y0[1]};
            double maxErr = 0.0;
            for (int k = 0; k <= n; k++) {
                double qEx = exactQ(t, omegaN, zeta, y0[0], y0[1]);
                maxErr = Math.max(maxErr, Math.abs(y[0] - qEx));
                if (k < n) {
                    y = rk4Step(f, t, y, h);
                    t += h;
                }
            }
            errs[i] = maxErr;
            System.out.printf("h=%.4f, max|q-q_exact|=%e%n", h, maxErr);
        }

        for (int i = 0; i < errs.length - 1; i++) {
            double pObs = Math.log(errs[i] / errs[i + 1]) / Math.log(2.0);
            System.out.printf("p_obs(%.4f->%.4f)=%.4f%n", hs[i], hs[i + 1], pObs);
        }
        System.out.println();
    }

    static void adaptiveDemo() {
        System.out.println("=== Adaptive RK4 (step-doubling) ===");
        double omegaN = 4.0, zeta = 0.05, tf = 12.0;
        double[] y0 = {1.0, 0.0};
        RHS f = dampedOscillatorRHS(omegaN, zeta);

        AdaptiveResult res = adaptiveRK4StepDoubling(f, 0.0, y0, tf, 0.1, 1e-6, 1e-9);

        boolean monotoneEnergy = true;
        double prevE = energy(res.y.get(0), omegaN);
        double maxErr = 0.0;

        for (int i = 1; i < res.t.size(); i++) {
            double Ei = energy(res.y.get(i), omegaN);
            if (Ei - prevE > 1e-8) monotoneEnergy = false;
            prevE = Ei;
        }
        for (int i = 0; i < res.t.size(); i++) {
            double qEx = exactQ(res.t.get(i), omegaN, zeta, y0[0], y0[1]);
            maxErr = Math.max(maxErr, Math.abs(res.y.get(i)[0] - qEx));
        }

        double[] yf = res.y.get(res.y.size() - 1);
        System.out.printf("Accepted=%d, Rejected=%d%n", res.accepted, res.rejected);
        System.out.printf("Final state: q=%f, v=%f%n", yf[0], yf[1]);
        System.out.printf("Energy nonincreasing=%s%n", Boolean.toString(monotoneEnergy));
        System.out.printf("max|q-q_exact| on adaptive mesh=%e%n", maxErr);
    }

    public static void main(String[] args) {
        eulerStabilityDemo(-50.0);
        convergenceDemo();
        adaptiveDemo();
    }
}

This MATLAB script compares adaptive solvers, performs RK4 convergence verification, and includes a Simulink set_param block to configure variable-step simulation options.

Code title: Chapter15_Lesson5.m

% Chapter15_Lesson5.m
% Numerical Stability, Error Control, and Model Verification via Simulation
% MATLAB/Simulink implementation for Chapter 15, Lesson 5 (System Dynamics)
%
% Features
%   1) Explicit Euler stability test on y' = lambda y
%   2) ode45 vs ode15s comparison (error-control behavior)
%   3) RK4 convergence verification + observed order
%   4) Simulink solver configuration example (if a model exists)

clear; clc;

%% 1) Absolute-stability demo (explicit Euler)
lambda = -50;
hs = [0.01 0.03 0.05 0.06];
N = 30;
fprintf('=== Explicit Euler stability demo ===\n');
for h = hs
    amp = abs(1 + h*lambda);
    y = 1;
    for n = 1:N
        y = y + h*lambda*y;
    end
    fprintf('h = %.3f, |1+h*lambda| = %.6f, y_N = %.6e\n', h, amp, y);
end
fprintf('\n');

%% 2) Damped oscillator: adaptive solvers and verification target
omega_n = 4.0; zeta = 0.05;
f = @(t,x) [x(2); -2*zeta*omega_n*x(2) - omega_n^2*x(1)];
tspan = [0 12];
x0 = [1; 0];

% ode45 (nonstiff, adaptive)
opts45 = odeset('RelTol',1e-6,'AbsTol',1e-9);
[t45, x45] = ode45(f, tspan, x0, opts45);

% ode15s (stiff-capable, adaptive)
opts15s = odeset('RelTol',1e-6,'AbsTol',1e-9);
[t15s, x15s] = ode15s(f, tspan, x0, opts15s);

qExact45 = exactQ(t45, omega_n, zeta, x0(1), x0(2));
qExact15s = exactQ(t15s, omega_n, zeta, x0(1), x0(2));

err45 = max(abs(x45(:,1) - qExact45));
err15s = max(abs(x15s(:,1) - qExact15s));

E45 = 0.5*(x45(:,2).^2 + omega_n^2*x45(:,1).^2);
energyNonInc45 = all(diff(E45) <= 1e-8);

fprintf('=== Adaptive solver comparison ===\n');
fprintf('ode45 steps = %d, max|q-q_exact| = %.3e, energy nonincreasing = %d\n', ...
    numel(t45)-1, err45, energyNonInc45);
fprintf('ode15s steps = %d, max|q-q_exact| = %.3e\n', numel(t15s)-1, err15s);
fprintf('\n');

%% 3) Fixed-step RK4 convergence verification
fprintf('=== RK4 convergence verification ===\n');
hs = [0.2 0.1 0.05 0.025];
errs = zeros(size(hs));
tf = 5.0;
for i = 1:numel(hs)
    h = hs(i);
    [t, X] = rk4Fixed(f, [0 tf], x0, h);
    qEx = exactQ(t, 4.0, 0.1, x0(1), x0(2));
    errs(i) = max(abs(X(:,1) - qEx));
    fprintf('h = %.4f, max error = %.6e\n', h, errs(i));
end
for i = 1:numel(hs)-1
    pObs = log(errs(i)/errs(i+1))/log(2);
    fprintf('p_obs(%.4f -> %.4f) = %.4f\n', hs(i), hs(i+1), pObs);
end
fprintf('\n');

%% 4) Simulink solver configuration snippet (optional)
% If you have a Simulink model named 'Chapter15Lesson5Model', this block
% configures solver and tolerances programmatically.
mdl = 'Chapter15Lesson5Model';
if exist([mdl '.slx'],'file') || exist([mdl '.mdl'],'file')
    load_system(mdl);
    set_param(mdl, ...
        'SolverType','Variable-step', ...
        'Solver','ode45', ...
        'RelTol','1e-6', ...
        'AbsTol','1e-9', ...
        'MaxStep','0.2');
    fprintf('Configured Simulink model "%s" with variable-step ode45.\n', mdl);

    % For stiff behavior, you may switch to ode15s:
    % set_param(mdl, 'Solver', 'ode15s');

    % To simulate:
    % simOut = sim(mdl, 'StopTime', '12');
else
    fprintf('No Simulink model file found. Skipping set_param demo.\n');
end

%% Local functions
function [t, X] = rk4Fixed(f, tspan, x0, h)
    t0 = tspan(1); tf = tspan(2);
    N = round((tf - t0)/h);
    t = linspace(t0, tf, N+1)';
    X = zeros(N+1, numel(x0));
    X(1,:) = x0(:).';
    for k = 1:N
        tk = t(k);
        xk = X(k,:).';
        k1 = f(tk, xk);
        k2 = f(tk + h/2, xk + h*k1/2);
        k3 = f(tk + h/2, xk + h*k2/2);
        k4 = f(tk + h,   xk + h*k3);
        X(k+1,:) = (xk + h*(k1 + 2*k2 + 2*k3 + k4)/6).';
    end
end

function q = exactQ(t, omega_n, zeta, q0, v0)
    wd = omega_n * sqrt(1 - zeta^2);
    A = q0;
    B = (v0 + zeta*omega_n*q0)/wd;
    q = exp(-zeta*omega_n*t) .* (A*cos(wd*t) + B*sin(wd*t));
end

The Mathematica script uses NDSolveValue for adaptive integration and includes a custom RK4 verification routine for observed order estimation.

Code title: Chapter15_Lesson5.nb

(* Chapter15_Lesson5.nb *)
(* Numerical Stability, Error Control, and Model Verification via Simulation *)
(* Wolfram Mathematica implementation for Chapter 15, Lesson 5 (System Dynamics) *)

ClearAll["Global`*"];

(* 1) Explicit Euler absolute-stability demo for y' = lambda y *)
eulerAmplification[lambda_, h_] := Abs[1 + h lambda];

eulerMarch[lambda_, h_, n_] := NestList[# + h lambda # &, 1.0, n];

Print["=== Explicit Euler stability demo ==="];
lambda = -50.0;
Do[
  amp = eulerAmplification[lambda, h];
  yN = Last[eulerMarch[lambda, h, 30]];
  Print[
    "h = ", NumberForm[h, {4, 3}],
    ", |1 + h lambda| = ", NumberForm[amp, {8, 6}],
    ", y_N = ", ScientificForm[yN]
  ],
  {h, {0.01, 0.03, 0.05, 0.06} }
];

(* 2) Damped oscillator exact solution and energy *)
omegaN = 4.0;
zeta = 0.05;
q0 = 1.0; v0 = 0.0;

wd = omegaN Sqrt[1 - zeta^2];
qExact[t_] := Exp[-zeta omegaN t] (q0 Cos[wd t] + ((v0 + zeta omegaN q0)/wd) Sin[wd t]);

energy[q_, v_] := 1/2 (v^2 + omegaN^2 q^2);

(* 3) Adaptive NDSolve with error control *)
sol = NDSolveValue[
  {
    q'[t] == v[t],
    v'[t] == -2 zeta omegaN v[t] - omegaN^2 q[t],
    q[0] == q0,
    v[0] == v0
  },
  {q, v},
  {t, 0, 12},
  AccuracyGoal -> 9,
  PrecisionGoal -> 6,
  MaxStepSize -> 0.2
];

sample = Table[{tt, sol[[1]][tt], sol[[2]][tt]}, {tt, 0, 12, 0.02}];
qErrInf = Max[Abs[(sample[[All, 2]]) - (qExact /@ sample[[All, 1]])]];
energies = energy @@@ sample[[All, {2, 3}]];
energyNonInc = And @@ Thread[Differences[energies] <= 10^-8];

Print[" "];
Print["=== Adaptive NDSolve verification ==="];
Print["max|q - qExact| on sampled grid = ", ScientificForm[qErrInf]];
Print["Energy nonincreasing (sampled) = ", energyNonInc];

(* 4) Fixed-step RK4 verification and observed order *)
rk4Step[f_, t_, y_, h_] := Module[{k1, k2, k3, k4},
  k1 = f[t, y];
  k2 = f[t + h/2, y + (h/2) k1];
  k3 = f[t + h/2, y + (h/2) k2];
  k4 = f[t + h, y + h k3];
  y + (h/6) (k1 + 2 k2 + 2 k3 + k4)
];

fOsc[t_, y_] := {y[[2]], -2*0.1*4.0*y[[2]] - 4.0^2*y[[1]]};

rk4Fixed[t0_, tf_, y0_, h_] := Module[{n, t = t0, y = y0, pts},
  n = Round[(tf - t0)/h];
  pts = { {t, y[[1]], y[[2]]} };
  Do[
    y = rk4Step[fOsc, t, y, h];
    t = t + h;
    AppendTo[pts, {t, y[[1]], y[[2]]}],
    {n}
  ];
  pts
];

hs = {0.2, 0.1, 0.05, 0.025};
errs = Table[
   pts = rk4Fixed[0.0, 5.0, {1.0, 0.0}, h];
   Max[Abs[pts[[All, 2]] - (qExact /@ pts[[All, 1]])]],
   {h, hs}
];

Print[" "];
Print["=== RK4 convergence verification ==="];
Do[
  Print["h = ", hs[[i]], ", max error = ", ScientificForm[errs[[i]]]],
  {i, Length[hs]}
];

pObs = Table[Log[errs[[i]]/errs[[i + 1]]]/Log[2], {i, 1, Length[errs] - 1}];
Do[
  Print["p_obs(", hs[[i]], " -> ", hs[[i + 1]], ") = ", NumberForm[pObs[[i]], {5, 4}]],
  {i, Length[pObs]}
];

Problem 1 (Absolute stability of explicit Euler): For the test equation \( y'=\lambda y \) with \( \lambda=-\alpha \), \( \alpha > 0 \), derive the admissible step-size interval for explicit Euler.

Solution: Explicit Euler gives

\[ y_{n+1}=(1+h\lambda)y_n=(1-\alpha h)y_n. \]

Stability requires \( |1-\alpha h|\leq 1 \), hence

\[ -1 \leq 1-\alpha h \leq 1. \]

Subtract 1 and multiply by \( -1 \) (reversing inequalities):

\[ 0 \leq \alpha h \leq 2 \quad \Longrightarrow \quad 0 \leq h \leq \frac{2}{\alpha}. \]

Therefore the numerical scheme is stable only for sufficiently small \( h \), even though the continuous system is stable for all time.

Problem 2 (Step-doubling estimator): Assume a one-step method of order \( p \) and write the leading truncation-error term for one full step and two half steps. Derive the step-doubling error estimator.

Solution: Let the exact solution at \( t_{n+1} \) be \( x^\star \). Then

\[ x^{[h]} = x^\star + C h^{p+1} + \mathcal{O}(h^{p+2}), \qquad x^{[h/2]} = x^\star + \frac{C}{2^p}h^{p+1} + \mathcal{O}(h^{p+2}). \]

Subtracting:

\[ x^{[h/2]} - x^{[h]} = -\left(1-\frac{1}{2^p}\right) C h^{p+1} + \mathcal{O}(h^{p+2}). \]

Hence an asymptotically correct estimator for the local one-step error is

\[ \widehat{e} \approx \frac{x^{[h/2]} - x^{[h]} }{2^p-1}. \]

Problem 3 (Observed order from a convergence table): Suppose a fixed-step RK method yields final-time errors \( E_{0.2}=1.3157\times 10^{-2} \) and \( E_{0.1}=7.9849\times 10^{-4} \). Estimate the observed order using refinement ratio \( r=2 \).

Solution:

\[ p_{\mathrm{obs} }=\frac{\log(E_{0.2}/E_{0.1})}{\log 2}. \]

\[ p_{\mathrm{obs} } = \frac{\log\!\left((1.3157\times 10^{-2})/(7.9849\times 10^{-4})\right)}{\log 2} \approx 4.04. \]

This is consistent with fourth-order convergence (RK4).

Problem 4 (Energy dissipation verification): For the damped oscillator \( \ddot{q}+2\zeta\omega_n\dot{q}+\omega_n^2 q=0 \), prove that the energy surrogate \( E=\frac{1}{2}\dot{q}^2+\frac{1}{2}\omega_n^2 q^2 \) is nonincreasing.

Solution:

\[ \dot{E} = \dot{q}\ddot{q} + \omega_n^2 q\dot{q} = \dot{q}\left(\ddot{q}+\omega_n^2 q\right). \]

Using the equation of motion, \( \ddot{q}+\omega_n^2 q = -2\zeta\omega_n\dot{q} \), so

\[ \dot{E} = -2\zeta\omega_n \dot{q}^2 \leq 0. \]

Therefore energy must decay monotonically (or remain constant when \( \zeta=0 \)). If a numerical simulation violates this trend significantly, the step size may be too large, the method may be unstable, or the implementation may contain an error.

Problem 5 (Residual-based verification): Let a numerical trajectory \( \mathbf{x}_h(t) \) be reconstructed by interpolation. Define the residual \( \mathbf{r}_h(t)=\dot{\mathbf{x} }_h(t)-\mathbf{f}(t,\mathbf{x}_h(t)) \). Explain why \( \|\mathbf{r}_h\| \) should decrease under mesh refinement for a correct implementation.

Solution: For a consistent method, the interpolated numerical solution converges to the exact solution as \( h \to 0 \). Since the exact solution satisfies \( \dot{\mathbf{x} }-\mathbf{f}(t,\mathbf{x})=\mathbf{0} \), the reconstructed residual should tend to zero:

\[ \|\mathbf{r}_h\| = \mathcal{O}(h^p) \quad \text{(up to interpolation order effects)}. \]

If residual norms do not decrease when refining \( h \), likely causes include: coding mistakes, event-discontinuity mishandling, inconsistent units, or solver tolerances that are too loose.

We established the stability-function view of numerical integration, derived practical error-control formulas for adaptive time stepping, and formalized verification through convergence order, Richardson extrapolation, and physics-based checks. These tools are central to trustworthy system simulation in control engineering, especially before controller design decisions are made from simulated trajectories.