In many engineering systems, the state trajectory is not governed by a single smooth ODE for all time. Instead, the flow is interrupted by \( \text{events} \) such as impacts, relay switching, contact loss, valve opening/closing, saturation entry, and sample-triggered updates. These systems are hybrid because they combine:

• Continuous dynamics: \( \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, t, q) \)
• Discrete mode updates: \( q^+ = \delta(q^-, \mathbf{x}^-, t_e) \)
• State reset maps: \( \mathbf{x}^+ = \mathbf{R}_{q^- → q^+}(\mathbf{x}^-, t_e) \)

Here \( \mathbf{x}(t) \in \mathbb{R}^n \) is the continuous state and \( q \in \mathcal{Q} \) is a discrete mode (e.g., heater ON/OFF, contact/no-contact). The superscripts \( - \) and \( + \) denote values immediately before and after an event time \( t_e \).

flowchart TD
  A["Initialize t0, x0, mode q0"] --> B["Integrate ODE in current mode"]
  B --> C["Monitor event functions g_j(x,t)"]
  C --> D["Zero crossing detected?"]
  D -->|No| B
  D -->|Yes| E["Refine event time t_e"]
  E --> F["Evaluate x(t_e-)"]
  F --> G["Apply reset x+ = R(x-)"]
  G --> H["Update mode q+ = delta(q-, x-, t_e)"]
  H --> I["Restart integration"]
  I --> B
        

A control-engineering hybrid model can be written as a family of mode-dependent vector fields. For each mode \( q \in \mathcal{Q} \), define a \( \text{flow set} \) \( \mathcal{C}_q \subseteq \mathbb{R}^n \times \mathbb{R} \) and a \( \text{jump set} \) \( \mathcal{D}_q \subseteq \mathbb{R}^n \times \mathbb{R} \).

\[ \begin{aligned} (\mathbf{x}, t) \in \mathcal{C}_q &\Rightarrow \dot{\mathbf{x}} = \mathbf{f}_q(\mathbf{x}, t), \\ (\mathbf{x}, t) \in \mathcal{D}_q &\Rightarrow \begin{cases} \mathbf{x}^+ = \mathbf{R}_q(\mathbf{x}, t), \\ q^+ = \delta_q(\mathbf{x}, t). \end{cases} \end{aligned} \]

A common practical representation uses \( m \) event functions \( g_j(\mathbf{x}, t) \). Event \( j \) is triggered when \( g_j(\mathbf{x}, t) = 0 \) with a required crossing direction. For example, descending-only detection for a ground impact uses \( g(\mathbf{x}, t) = h \) and requires \( \dot{g} = \nabla g^\top \mathbf{f}_q < 0 \).

The hybrid trajectory is therefore piecewise smooth on intervals \( [t_k, t_{k+1}) \) and may have jump discontinuities \( \mathbf{x}(t_k^+) \neq \mathbf{x}(t_k^-) \). Numerical solvers must preserve the event ordering \( t_0 < t_1 < t_2 < \cdots \) and apply the correct reset at each event.

Let \( g(\mathbf{x}, t) \) be a smooth scalar event function and suppose the exact event occurs at \( t_e \) such that \( g(\mathbf{x}(t_e), t_e) = 0 \). A key condition is \( \text{transversality} \):

\[ \frac{d}{dt}g(\mathbf{x}(t), t)\bigg|_{t=t_e} = \nabla_{\mathbf{x}} g(\mathbf{x}(t_e), t_e)^\top \mathbf{f}_q(\mathbf{x}(t_e), t_e) + \frac{\partial g}{\partial t}(\mathbf{x}(t_e), t_e) \neq 0. \]

If this derivative is nonzero, the implicit function theorem implies a locally unique event time. This avoids ambiguous grazing/contact tangency situations in which \( \frac{d}{dt}g = 0 \).

Theorem (first-order event-time correction from interpolation): Suppose an integration step from \( t_n \) to \( t_{n+1}=t_n+h \) brackets a unique transverse event and let \( \hat{\mathbf{x}}(t) \) be a continuous extension of order \( p \) on the step, i.e., \( \|\hat{\mathbf{x}}(t)-\mathbf{x}(t)\| = \mathcal{O}(h^{p+1}) \). Then the root \( \hat{t}_e \) of \( g(\hat{\mathbf{x}}(t), t)=0 \) satisfies \( |\hat{t}_e - t_e| = \mathcal{O}(h^{p+1}) \).

Proof sketch: Define \( \phi(t)=g(\mathbf{x}(t),t) \) and \( \hat{\phi}(t)=g(\hat{\mathbf{x}}(t),t) \). By smoothness of \( g \), the interpolation error implies \( \hat{\phi}(t)-\phi(t)=\mathcal{O}(h^{p+1}) \). Since \( \phi'(t_e)\neq 0 \) (transversality), a Taylor expansion around \( t_e \) gives \( \phi(\hat{t}_e)=\phi'(t_e)(\hat{t}_e-t_e)+\mathcal{O}((\hat{t}_e-t_e)^2) \). But \( \hat{\phi}(\hat{t}_e)=0 \), so \( \phi(\hat{t}_e)=\phi(\hat{t}_e)-\hat{\phi}(\hat{t}_e)=\mathcal{O}(h^{p+1}) \). Dividing by \( \phi'(t_e) \) yields \( \hat{t}_e-t_e=\mathcal{O}(h^{p+1}) \).

In practice, solvers use sign-change detection \( g_n g_{n+1} \le 0 \) plus root refinement (bisection, secant, or safeguarded Newton) on a dense output polynomial.

A switching system without state jumps has \( \mathbf{x}^+ = \mathbf{x}^- \) but changes the active vector field. A reset system (impact, sampler, relay with memory) may change both state and mode. Formally, for event index \( j \):

\[ (q^-, \mathbf{x}^-) \rightarrow (q^+, \mathbf{x}^+) = \left(\delta_j(q^-, \mathbf{x}^-, t_e), \; \mathbf{R}_j(\mathbf{x}^-, t_e)\right). \]

Example 1 (impact reset): Bouncing ball state \( \mathbf{x}=[h, v]^\top \) with ground event \( g(\mathbf{x})=h \) and restitution \( e\in(0,1) \):

\[ \dot{h}=v,\qquad \dot{v}=-g_0,\qquad h=0,\; v^- < 0 \Rightarrow \begin{cases} h^+ = 0, \\ v^+ = -e\,v^-. \end{cases} \]

Example 2 (thermostat with hysteresis): Let mode \( q \in \{\text{ON}, \text{OFF}\} \) and temperature \( T \) satisfy

\[ \dot{T} = -\frac{1}{\tau}(T-T_{\text{env}}) + \frac{1}{\tau}Q\,\mathbf{1}_{\{q=\text{ON}\}}. \]

To avoid chattering, use two thresholds \( T_{\text{low}} < T_{\text{high}} \):

\[ q^+ = \begin{cases} \text{OFF}, & q^-=\text{ON} \text{ and } T^- \ge T_{\text{high}}, \\ \text{ON}, & q^-=\text{OFF} \text{ and } T^- \le T_{\text{low}}, \\ q^-, & \text{otherwise.} \end{cases} \]

flowchart TD
  ON["Mode: Heater ON"] -->|T reaches T_high| OFF["Mode: Heater OFF"]
  OFF -->|T reaches T_low| ON
  ON --> A["Flow: dT/dt = -(T-Tenv)/tau + Q/tau"]
  OFF --> B["Flow: dT/dt = -(T-Tenv)/tau"]
        

A consistent hybrid integrator must separate three numerical tasks:

1. Continuous propagation within a mode
2. Event detection and time refinement
3. Jump/reset application and mode update

For a one-step method of order \( p \), the global state error in smooth regions is typically \( \mathcal{O}(h^p) \). However, event timing errors contribute an extra term. If \( \hat{t}_e-t_e = \mathcal{O}(h^{p+1}) \) and the pre-event flow is smooth, then the induced pre-jump state perturbation is

\[ \mathbf{x}(\hat{t}_e^-) - \mathbf{x}(t_e^-) = \mathbf{f}_{q^-}(\mathbf{x}(t_e^-), t_e)(\hat{t}_e-t_e) + \mathcal{O}((\hat{t}_e-t_e)^2). \]

After the jump, linearization of the reset map gives

\[ \mathbf{x}(\hat{t}_e^+) - \mathbf{x}(t_e^+) = D\mathbf{R}(\mathbf{x}(t_e^-), t_e)\left[\mathbf{x}(\hat{t}_e^-) - \mathbf{x}(t_e^-)\right] + \mathcal{O}(\|\mathbf{x}(\hat{t}_e^-)-\mathbf{x}(t_e^-)\|^2). \]

Thus, accurate event localization is essential for preserving the expected order of accuracy.

Saltation matrix (advanced local sensitivity across a jump): For trajectory sensitivity and variational analysis across a transverse event surface \( g(\mathbf{x}, t)=0 \), the perturbation may be mapped by the saltation matrix \( \mathbf{S} \):

\[ \delta \mathbf{x}^+ = \mathbf{S}\,\delta \mathbf{x}^-, \qquad \mathbf{S} = D\mathbf{R} + \frac{\left(\mathbf{f}^+ - D\mathbf{R}\,\mathbf{f}^- - \partial_t\mathbf{R}\right)\nabla g^\top} {\nabla g^\top \mathbf{f}^- + \partial_t g}. \]

This expression is important in hybrid stability analysis, shooting methods, and sensitivity-based design.

Numerical pathology: if switching is implemented without hysteresis near a threshold, repeated crossings may cause chattering, producing an excessive number of events. In some models, infinitely many events can accumulate in finite time (Zeno behavior). Practical simulators use hysteresis, minimum dwell time, or regularization to avoid solver failure.

The Python implementation below simulates a bouncing ball with exact event handling logic: continuous integration between impacts, root detection at \( h=0 \), and reset \( v^+ = -e v^- \). This is a canonical hybrid benchmark for event-driven solvers.

Chapter15_Lesson4.py

# Chapter15_Lesson4.py
"""
Event handling, switching, and hybrid continuous-discrete dynamics
Example: Bouncing ball with restitution using solve_ivp event detection.
"""

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt


g = 9.81           # gravitational acceleration (m/s^2)
e = 0.82           # coefficient of restitution
t0, tf = 0.0, 8.0  # simulation horizon
y0 = np.array([1.5, 0.0])  # [height, velocity]
min_bounce_speed = 0.05     # terminate when impacts become negligible


def flow(t, y):
    h, v = y
    return [v, -g]


def ground_event(t, y):
    # Event function g(y) = h. Trigger at h = 0 while descending.
    return y[0]


ground_event.terminal = True
ground_event.direction = -1.0


def jump_map(y_minus):
    h_minus, v_minus = y_minus
    # Reset state at impact: h^+ = 0, v^+ = -e v^-
    return np.array([0.0, -e * v_minus])


def simulate_hybrid_bouncing_ball():
    t_segments = []
    y_segments = []
    impact_times = []

    t_curr = t0
    y_curr = y0.copy()

    while t_curr < tf:
        sol = solve_ivp(
            flow,
            (t_curr, tf),
            y_curr,
            events=ground_event,
            dense_output=True,
            max_step=0.03,
            rtol=1e-8,
            atol=1e-10,
        )

        t_segments.append(sol.t)
        y_segments.append(sol.y.T)

        if sol.t_events[0].size == 0:
            break

        te = float(sol.t_events[0][0])
        ye_minus = sol.sol(te)
        impact_times.append(te)

        if abs(ye_minus[1]) < min_bounce_speed:
            # Stop after final small impact.
            t_segments.append(np.array([te]))
            y_segments.append(np.array([[0.0, 0.0]]))
            break

        y_plus = jump_map(ye_minus)
        t_curr = te
        # Restart slightly after the event to avoid immediate re-detection.
        y_curr = y_plus + np.array([0.0, 1e-12])

    return t_segments, y_segments, impact_times


if __name__ == "__main__":
    ts, ys, impacts = simulate_hybrid_bouncing_ball()

    t_plot = np.concatenate(ts)
    y_plot = np.vstack(ys)

    plt.figure(figsize=(8, 4))
    plt.plot(t_plot, y_plot[:, 0], label="height h(t)")
    for te in impacts:
        plt.axvline(te, linestyle="--", linewidth=0.8)
    plt.xlabel("Time (s)")
    plt.ylabel("Height (m)")
    plt.title("Hybrid Simulation: Bouncing Ball with Events")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

    print("Impact times:", impacts)

Library note: scipy.integrate.solve_ivp supports event functions, directional zero-crossing control, and dense output for accurate event localization.

This implementation is intentionally low-level and suitable for understanding how simulation engines are built internally. We use RK4 for continuous propagation and a bisection routine to refine the event time inside a step when \( g_n g_{n+1} \le 0 \).

Chapter15_Lesson4.cpp

// Chapter15_Lesson4.cpp
// Event handling and hybrid simulation from scratch (RK4 + bisection)
// Example: Bouncing ball with restitution

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

struct State {
    double h; // height
    double v; // velocity
};

static constexpr double g = 9.81;
static constexpr double e = 0.82;

State flow(const State& s) {
    return {s.v, -g};
}

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

State rk4Step(const State& s, double dt) {
    State k1 = flow(s);
    State k2 = flow(add(s, k1, 0.5 * dt));
    State k3 = flow(add(s, k2, 0.5 * dt));
    State k4 = flow(add(s, k3, dt));
    State out;
    out.h = s.h + (dt / 6.0) * (k1.h + 2.0 * k2.h + 2.0 * k3.h + k4.h);
    out.v = s.v + (dt / 6.0) * (k1.v + 2.0 * k2.v + 2.0 * k3.v + k4.v);
    return out;
}

double eventFunction(const State& s) {
    return s.h; // ground event when h = 0
}

State jumpMap(const State& sMinus) {
    return {0.0, -e * sMinus.v};
}

// Refine event time in [tL, tR] assuming sign change in h and monotone descent near event.
void refineEventBisection(double tL, const State& sL,
                         double tR, const State& sR,
                         double tol,
                         double& tEvent,
                         State& sEvent) {
    double a = tL, b = tR;
    State sa = sL, sb = sR;

    for (int iter = 0; iter < 60; ++iter) {
        double m = 0.5 * (a + b);
        double dt = m - a;
        State sm = rk4Step(sa, dt);

        if (std::abs(eventFunction(sm)) < tol || (b - a) < tol) {
            tEvent = m;
            sEvent = sm;
            return;
        }

        if (eventFunction(sa) * eventFunction(sm) <= 0.0) {
            b = m;
            sb = sm;
        } else {
            a = m;
            sa = sm;
        }
    }

    tEvent = 0.5 * (a + b);
    sEvent = rk4Step(sa, tEvent - a);
}

int main() {
    const double tFinal = 8.0;
    const double dt = 0.01;
    const double eventTol = 1e-10;
    const double minBounceSpeed = 0.05;

    double t = 0.0;
    State s{1.5, 0.0};

    std::ofstream csv("Chapter15_Lesson4_cpp_output.csv");
    csv << "t,h,v\n";
    csv << std::setprecision(12);

    csv << t << "," << s.h << "," << s.v << "\n";

    while (t < tFinal) {
        State sNext = rk4Step(s, dt);
        double g0 = eventFunction(s);
        double g1 = eventFunction(sNext);

        // Check for zero crossing while descending
        if (g0 > 0.0 && g1 <= 0.0 && s.v < 0.0) {
            double tEvent;
            State sMinus;
            refineEventBisection(t, s, t + dt, sNext, eventTol, tEvent, sMinus);

            csv << tEvent << "," << sMinus.h << "," << sMinus.v << "\n";

            if (std::abs(sMinus.v) < minBounceSpeed) {
                State rest{0.0, 0.0};
                csv << tEvent << "," << rest.h << "," << rest.v << "\n";
                break;
            }

            State sPlus = jumpMap(sMinus);
            t = tEvent;
            s = sPlus;
            csv << t << "," << s.h << "," << s.v << "\n";
            continue;
        }

        t += dt;
        s = sNext;
        csv << t << "," << s.h << "," << s.v << "\n";
    }

    csv.close();
    std::cout << "Simulation complete. Results written to Chapter15_Lesson4_cpp_output.csv\n";
    return 0;
}

The code writes a CSV file (Chapter15_Lesson4_cpp_output.csv) so the trajectory can be plotted externally. This pattern is common in embedded or real-time simulation pipelines where visualization is handled separately.

The Java example models a thermostat as a hybrid system with two discrete modes. The state (temperature) evolves continuously, while the mode switches when thresholds are crossed. Hysteresis prevents chattering and is numerically robust.

Chapter15_Lesson4.java

// Chapter15_Lesson4.java
// Hybrid continuous-discrete simulation in Java
// Example: Thermostat with hysteresis (mode switching)

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

public class Chapter15_Lesson4 {

    enum Mode { HEATER_ON, HEATER_OFF }

    static class State {
        double t;
        double T;
        Mode mode;

        State(double t, double T, Mode mode) {
            this.t = t;
            this.T = T;
            this.mode = mode;
        }
    }

    // Parameters
    static final double T_ENV = 18.0;     // ambient temperature
    static final double TAU = 12.0;       // thermal time constant
    static final double Q_HEAT = 8.0;     // heating gain
    static final double T_LOW = 20.0;     // lower threshold
    static final double T_HIGH = 22.0;    // upper threshold

    static double f(double T, Mode mode) {
        double heater = (mode == Mode.HEATER_ON) ? Q_HEAT : 0.0;
        return -(T - T_ENV) / TAU + heater / TAU;
    }

    static double rk4Temperature(double T, Mode mode, double dt) {
        double k1 = f(T, mode);
        double k2 = f(T + 0.5 * dt * k1, mode);
        double k3 = f(T + 0.5 * dt * k2, mode);
        double k4 = f(T + dt * k3, mode);
        return T + (dt / 6.0) * (k1 + 2*k2 + 2*k3 + k4);
    }

    static class EventResult {
        boolean switched;
        double tEvent;
        double TEvent;
        Mode nextMode;

        EventResult(boolean switched, double tEvent, double TEvent, Mode nextMode) {
            this.switched = switched;
            this.tEvent = tEvent;
            this.TEvent = TEvent;
            this.nextMode = nextMode;
        }
    }

    static EventResult detectAndRefineEvent(double t, double T, Mode mode, double dt) {
        double Tnext = rk4Temperature(T, mode, dt);

        if (mode == Mode.HEATER_ON && T < T_HIGH && Tnext >= T_HIGH) {
            // Bisection on phi(T) = T - T_HIGH
            double a = t, b = t + dt;
            double Ta = T, Tb = Tnext;
            for (int i = 0; i < 60; i++) {
                double m = 0.5 * (a + b);
                double Tm = rk4Temperature(Ta, mode, m - a);
                if (Math.abs(Tm - T_HIGH) < 1e-10 || (b - a) < 1e-10) {
                    return new EventResult(true, m, T_HIGH, Mode.HEATER_OFF);
                }
                if ((Ta - T_HIGH) * (Tm - T_HIGH) <= 0.0) {
                    b = m;
                    Tb = Tm;
                } else {
                    a = m;
                    Ta = Tm;
                }
            }
            return new EventResult(true, 0.5 * (a + b), T_HIGH, Mode.HEATER_OFF);
        }

        if (mode == Mode.HEATER_OFF && T > T_LOW && Tnext <= T_LOW) {
            // Bisection on phi(T) = T - T_LOW
            double a = t, b = t + dt;
            double Ta = T, Tb = Tnext;
            for (int i = 0; i < 60; i++) {
                double m = 0.5 * (a + b);
                double Tm = rk4Temperature(Ta, mode, m - a);
                if (Math.abs(Tm - T_LOW) < 1e-10 || (b - a) < 1e-10) {
                    return new EventResult(true, m, T_LOW, Mode.HEATER_ON);
                }
                if ((Ta - T_LOW) * (Tm - T_LOW) <= 0.0) {
                    b = m;
                    Tb = Tm;
                } else {
                    a = m;
                    Ta = Tm;
                }
            }
            return new EventResult(true, 0.5 * (a + b), T_LOW, Mode.HEATER_ON);
        }

        return new EventResult(false, t + dt, Tnext, mode);
    }

    public static void main(String[] args) {
        double tFinal = 48.0;
        double dt = 0.1;

        double t = 0.0;
        double T = 19.0;
        Mode mode = Mode.HEATER_ON;

        List<State> history = new ArrayList<>();
        history.add(new State(t, T, mode));

        while (t < tFinal) {
            EventResult er = detectAndRefineEvent(t, T, mode, dt);
            t = er.tEvent;
            T = er.TEvent;
            mode = er.nextMode;
            history.add(new State(t, T, mode));

            if (!er.switched) {
                continue;
            }
            // After a switch, continue integration with new mode.
        }

        System.out.println("t,Temperature,Mode");
        for (State s : history) {
            System.out.printf("%.6f,%.6f,%s%n", s.t, s.T, s.mode);
        }
    }
}

This example is useful for control logic prototyping and for demonstrating how hybrid modeling appears in cyber-physical software, not only in mechanical impact systems.

MATLAB provides event handling through ode45, ode15s, and related ODE solvers using an Events function. Simulink provides zero-crossing detection at the block level and can implement the reset logic with Triggered Subsystems or Stateflow.

The script below performs an event-driven bouncing-ball simulation and includes a short Simulink note.

Chapter15_Lesson4.m

% Chapter15_Lesson4.m
% Event handling and hybrid dynamics in MATLAB / Simulink-oriented style
% Example: Bouncing ball with ode45 Events + reset map loop

function Chapter15_Lesson4
    params.g = 9.81;
    params.e = 0.82;
    params.minBounceSpeed = 0.05;

    t0 = 0.0;
    tf = 8.0;
    y0 = [1.5; 0.0]; % [height; velocity]

    T_all = [];
    Y_all = [];
    impactTimes = [];

    opts = odeset('Events', @(t,y) groundEvent(t,y), ...
                  'RelTol', 1e-8, 'AbsTol', 1e-10, 'MaxStep', 0.03);

    tStart = t0;
    yStart = y0;

    while tStart < tf
        [T, Y, TE, YE, IE] = ode45(@(t,y) flow(t,y,params), [tStart tf], yStart, opts); %#ok<ASGLU>
        T_all = [T_all; T]; %#ok<AGROW>
        Y_all = [Y_all; Y]; %#ok<AGROW>

        if isempty(TE)
            break;
        end

        impactTimes(end+1,1) = TE(1); %#ok<AGROW>
        yMinus = YE(1,:).';

        if abs(yMinus(2)) < params.minBounceSpeed
            T_all = [T_all; TE(1)]; %#ok<AGROW>
            Y_all = [Y_all; 0.0, 0.0]; %#ok<AGROW>
            break;
        end

        yPlus = [0.0; -params.e * yMinus(2)];
        tStart = TE(1);
        yStart = yPlus + [0.0; 1e-12]; % avoid immediate retrigger
    end

    figure('Color', 'w');
    plot(T_all, Y_all(:,1), 'LineWidth', 1.5);
    grid on;
    xlabel('Time (s)');
    ylabel('Height (m)');
    title('Hybrid Simulation with Events (Bouncing Ball)');

    disp('Impact times:');
    disp(impactTimes);

    % Simulink note:
    % Use Integrator blocks with zero-crossing enabled, a Hit Crossing block
    % on h(t)=0, and a Triggered Subsystem / Stateflow chart to apply the reset
    % v^+ = -e v^- at each detected event.
end

function dydt = flow(~, y, params)
    dydt = [y(2); -params.g];
end

function [value, isterminal, direction] = groundEvent(~, y)
    value = y(1);      % event when height reaches zero
    isterminal = 1;    % stop integration and return control to script
    direction = -1;    % detect descending crossing only
end

In Simulink, the corresponding model uses two Integrator blocks for \( h \) and \( v \), a zero-crossing detector (e.g., Hit Crossing) for \( h=0 \), and reset/update logic for the velocity state.

Wolfram Mathematica supports hybrid event handling through WhenEvent inside NDSolve. The notebook content below uses event detection to stop integration at impact, apply the reset, and continue piecewise. This mirrors the event-driven strategy used in other solvers.

Chapter15_Lesson4.nb

Notebook[{
  Cell["Chapter15_Lesson4.nb - Event Handling, Switching, and Hybrid Dynamics", "Section"],
  Cell[BoxData[
    ToBoxes[
      Module[{g = 9.81, e = 0.82, tf = 8.0, x, v, impacts = {}, sol, tNow = 0., x0 = 1.5, v0 = 0.0, pieces = {}},
        While[tNow < tf,
          sol = NDSolveValue[
            {
              x'[t] == v[t],
              v'[t] == -g,
              x[tNow] == x0,
              v[tNow] == v0,
              WhenEvent[x[t] == 0 && v[t] < 0,
                {
                  AppendTo[impacts, t],
                  x[t] -> 0,
                  v[t] -> -e v[t]
                },
                "StopIntegration"
              ]
            },
            {x, v},
            {t, tNow, tf},
            Method -> {"EventLocator", "EventAction" -> "StopIntegration"}
          ];
          AppendTo[pieces, sol];
          If[impacts === {} || Last[impacts] <= tNow, Break[]];
          tNow = Last[impacts];
          x0 = 0.0;
          v0 = (sol[[2]][tNow]);
          If[Abs[v0] < 0.05, Break[]];
          v0 = v0 + 10^-12;
        ];
        Print["Impact times: ", impacts];
        Plot[
          Evaluate@Table[pieces[[k, 1]][t], {k, Length[pieces]}],
          {t, 0, tf},
          PlotRange -> All,
          GridLines -> {impacts, None},
          AxesLabel -> {"t", "h(t)"},
          PlotLabel -> "Hybrid Simulation with WhenEvent"
        ]
      ]
    ]
  ], "Input"]
}]

Problem 1 (Impact map derivation): Consider a point mass under gravity with state \( \mathbf{x}=[h,v]^\top \). At the ground, impose Newton restitution with coefficient \( e \in (0,1) \). Derive the jump map and show that kinetic energy immediately after impact is reduced by a factor \( e^2 \).

Solution: The impact law reverses the normal velocity and scales its magnitude:

\[ h^+ = 0,\qquad v^+ = -e v^-. \]

The pre- and post-impact kinetic energies are \( K^- = \frac{1}{2}m(v^-)^2 \) and \( K^+ = \frac{1}{2}m(v^+)^2 = \frac{1}{2}m e^2 (v^-)^2 \). Therefore

\[ K^+ = e^2 K^-, \]

so the impact dissipates a fraction \( 1-e^2 \) of kinetic energy.

Problem 2 (Transverse event uniqueness): Let \( \phi(t)=g(\mathbf{x}(t),t) \) be differentiable, and suppose \( \phi(t_e)=0 \) with \( \phi'(t_e)\neq 0 \). Explain why the event time is locally unique.

Solution: Since \( \phi'(t_e)\neq 0 \), the function \( \phi \) is strictly monotone in a small neighborhood of \( t_e \). A strictly monotone continuous function can cross zero at most once on a sufficiently small interval. Hence the event time is locally unique.

Problem 3 (Linear interpolation estimate of event time): Suppose within a step \( [t_n,t_{n+1}] \) the event function values satisfy \( g_n = g(t_n) > 0 \) and \( g_{n+1}=g(t_{n+1}) < 0 \). Derive the linear-interpolation estimate \( \hat{t}_e \).

Solution: Interpolate linearly: \( g(t) \approx g_n + \frac{t-t_n}{h}(g_{n+1}-g_n) \), where \( h=t_{n+1}-t_n \). Set this to zero and solve:

\[ \hat{t}_e = t_n - h \frac{g_n}{g_{n+1}-g_n}. \]

This estimate is first-order accurate for the event time unless enhanced by dense-output interpolation.

Problem 4 (Hysteresis avoids chattering): For a thermostat with switching thresholds \( T_{\text{low}} < T_{\text{high}} \), explain why the discrete mode cannot switch infinitely many times in an arbitrarily small neighborhood of one threshold (assuming bounded \( \dot{T} \)).

Solution: A switch from ON to OFF requires reaching \( T_{\text{high}} \). The next reverse switch requires the state to travel down to \( T_{\text{low}} \), which is separated by the positive gap \( \Delta T = T_{\text{high}} - T_{\text{low}} > 0 \). If \( |\dot{T}| \le M \), then the time to traverse this gap is at least \( \Delta T / M \). Therefore there is a positive lower bound on inter-switch time, which prevents arbitrarily fast repeated switching (chattering) near a single threshold.

Problem 5 (Hybrid simulation protocol design): Design a robust simulation protocol for a system with multiple event functions \( g_1,\dots,g_m \) and possible simultaneous events. Include tie-breaking and reset application order.

Solution: A robust protocol is:

1. Integrate one step and detect all sign changes
2. Refine candidate event times for all triggered functions
3. Select the smallest refined time \( t_e \)
4. Evaluate all event functions at \( t_e \) and form the active event set
5. Apply a pre-defined priority or a combined reset rule
6. Update mode, apply state reset, and restart integration

Priority rules are model-dependent, but they must be fixed in advance to ensure reproducibility.

We formalized hybrid continuous–discrete dynamics using event surfaces, mode transitions, and reset maps. The lesson emphasized transversality, accurate event-time localization, and the role of hysteresis in preventing chattering. We also connected these concepts to practical implementations in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica. These ideas are essential for the next lesson on numerical stability, error control, and verification of simulation results.