By the end of this lesson, you will be able to:

• Translate multi-domain physics into nonlinear state-space models \( \dot{\mathbf{x}} = \mathbf{f}(t,\mathbf{x}) \) with consistent units and parameters.
• Recognize and prove dissipativity (existence of an absorbing set) using energy-like Lyapunov arguments.
• Define the Poincaré map for periodically forced dynamics and use it as a discrete-time reduction for chaos analysis.
• Implement simulation + Poincaré sampling + crude Lyapunov exponent estimation (two-trajectory renormalization) in Python/C++/Java/MATLAB/Mathematica.

flowchart TD
  A["Define physics + assumptions"] --> B["Choose states x, inputs u, parameters p"]
  B --> C["Derive nonlinear model: xdot = f(t,x; p)"]
  C --> D["Check well-posedness: Lipschitz / smoothness"]
  D --> E["Prove boundedness or dissipativity (energy/Lyapunov)"]
  E --> F["Select forcing or operating regime (periodic input, saturation, etc.)"]
  F --> G["Simulate ODE (adaptive or fixed-step)"]
  G --> H["Build Poincare samples x(kT) for periodic forcing"]
  H --> I["Estimate invariants: fixed points, cycles, Lyapunov exponent"]
  I --> J["Validate numerics: step size, sensitivity, repeatability"]
        

Consider a \( T \)-periodic nonautonomous system \( \dot{\mathbf{x}} = \mathbf{f}(t,\mathbf{x}) \) with \( \mathbf{f}(t+T,\mathbf{x}) = \mathbf{f}(t,\mathbf{x}) \). Let \( \varphi(t;t_0,\mathbf{x}_0) \) denote its flow (solution map). The Poincaré map (stroboscopic map) is

\[ \mathbf{P}(\mathbf{x}_0) \;=\; \varphi(T;0,\mathbf{x}_0). \]

Iteration of the discrete-time system \( \mathbf{x}_{k+1}=\mathbf{P}(\mathbf{x}_k) \) captures fixed points (period-\( T \) solutions), periodic orbits (period-\( nT \) solutions), and chaotic invariant sets of the continuous-time dynamics.

Proposition (Well-defined Poincaré map). If \( \mathbf{f}(t,\mathbf{x}) \) is continuous in \( t \) and locally Lipschitz in \( \mathbf{x} \) uniformly on compact sets, then for any \( \mathbf{x}_0 \) in a domain where solutions exist on \( [0,T] \), the map \( \mathbf{P} \) is well-defined and continuous.

Sketch of proof. Local Lipschitz in \( \mathbf{x} \) implies existence and uniqueness of solutions; continuous dependence on initial conditions implies continuity of \( \mathbf{P}(\mathbf{x}_0)=\varphi(T;0,\mathbf{x}_0) \).

Lyapunov exponent connection. For a trajectory \( \mathbf{x}(t) \), the largest Lyapunov exponent is

\[ \lambda_{\max} \;=\; \limsup_{t\to\infty}\frac{1}{t}\log\frac{\|\delta\mathbf{x}(t)\|}{\|\delta\mathbf{x}(0)\|}, \]

where \( \delta\mathbf{x} \) satisfies the variational dynamics \( \delta\dot{\mathbf{x}} = \mathbf{J}(t)\delta\mathbf{x} \) with \( \mathbf{J}(t)=\frac{\partial \mathbf{f}}{\partial \mathbf{x}}(t,\mathbf{x}(t)) \). Numerically, a robust approach uses QR re-orthonormalization; a simpler (coarser) estimate uses two nearby trajectories with periodic renormalization (implemented in the provided code).

We model a DC motor with electrical inductance and back-EMF, driving a rigid load. Choose state \( \mathbf{x}=[\theta,\omega,i]^\top \) where \( \theta \) is shaft angle, \( \omega \) is angular speed, and \( i \) is armature current. With input voltage \( u(t) \) (saturated), the dynamics are:

\[ \begin{aligned} \dot{\theta} &= \omega, \\ J\dot{\omega} &= K_t i - \tau_f(\omega), \\ L\dot{i} &= u(t) - R i - K_e \omega, \end{aligned} \]

where \( J \) is inertia, \( R \) resistance, \( L \) inductance, \( K_t \) torque constant, \( K_e \) back-EMF constant. A smooth friction model combining viscous and Stribeck-like effects is:

\[ \tau_f(\omega) \;=\; b\omega \;+\; \Big(\tau_c + (\tau_s-\tau_c)\exp\!\big(-(\tfrac{|\omega|}{\omega_s})^2\big)\Big)\tanh(k\omega), \]

with viscous \( b \), Coulomb \( \tau_c \), static friction level \( \tau_s \), Stribeck speed \( \omega_s \), and smoothing gain \( k \). Input saturation is \( u(t)=\mathrm{sat}(u_0(t),u_{\max}) \).

Dissipativity (boundedness idea). Consider the storage function \( V(\omega,i) = \tfrac{1}{2}J\omega^2 + \tfrac{1}{2}L i^2 \). Using the dynamics (ignoring the integrator state \( \theta \) which is neutrally stable),

\[ \dot{V} = J\omega\dot{\omega} + Li\dot{i} = \omega(K_t i - \tau_f(\omega)) + i(u - Ri - K_e\omega). \]

If \( K_t=K_e \) (common in SI-consistent models), cross terms cancel:

\[ \dot{V} = -\omega\,\tau_f(\omega) - R i^2 + i\,u(t). \]

Since \( \omega\,\tau_f(\omega) \ge b\omega^2 \) for the viscous part and the smooth Coulomb term is nonnegative (it always opposes motion), we obtain a quadratic dissipation minus a bounded supply \( i\,u(t) \) with \( |u(t)|\le u_{\max} \). Applying Young's inequality, \( i\,u \le \tfrac{R}{2}i^2 + \tfrac{1}{2R}u^2 \), yields \( \dot{V} \le -b\omega^2 - \tfrac{R}{2}i^2 + \tfrac{1}{2R}u_{\max}^2 \), so solutions enter and remain in an absorbing set in \( (\omega,i) \).

We use the planar bicycle model with constant longitudinal speed \( U_x \) and states \( \mathbf{x}=[v_y,r]^\top \) (lateral velocity and yaw rate). For small angles, front and rear slip angles are

\[ \alpha_f \approx \frac{v_y + a r}{U_x} - \delta(t),\qquad \alpha_r \approx \frac{v_y - b r}{U_x}, \]

where \( a \), \( b \) are distances from center of gravity to front/rear axles, and \( \delta(t) \) is steering input. A smooth saturated tire model is

\[ F_y(\alpha) = C\,\alpha_s\,\tanh\!\Big(\frac{\alpha}{\alpha_s}\Big), \]

where \( C \) is (linear) cornering stiffness and \( \alpha_s \) sets saturation onset. Then the bicycle dynamics are:

\[ \begin{aligned} \dot{v}_y &= \frac{F_{yf}+F_{yr}}{m} - U_x r,\\ \dot{r} &= \frac{aF_{yf} - bF_{yr}}{I_z}, \end{aligned} \]

with \( F_{yf}=-F_y(\alpha_f) \), \( F_{yr}=-F_y(\alpha_r) \) (negative sign because lateral force opposes slip angle in this convention). A periodic steering can be taken as \( \delta(t) = \delta_0 \sin(2\pi t/T) \), defining a Poincaré map with period \( T \).

Linearized stability (local). Near \( \delta=0 \) and small slip (so \( \tanh(z)\approx z \)), the model becomes linear with effective stiffnesses \( C_f, C_r \). Stability depends on \( U_x \) (speed) and the understeer gradient; as \( U_x \) increases, eigenvalues can approach the imaginary axis, and nonlinear saturation + forcing can create complex responses.

Consider an irreversible exothermic reaction \( A \to B \) in a CSTR with constant volume \( V \). States: concentration \( C_A \) and temperature \( T \). Flow rate is \( q \), inlet concentration \( C_{Af} \). Reaction rate is Arrhenius:

\[ r_A(C_A,T) = k(T)C_A,\qquad k(T)=k_0\exp\!\Big(-\frac{E}{R_g T}\Big). \]

Mass balance:

\[ \dot{C}_A = \frac{q}{V}(C_{Af}-C_A) - r_A(C_A,T). \]

Energy balance (with coolant temperature \( T_c \) and heat-transfer \( UA \)):

\[ \dot{T} = \frac{q}{V}(T_f(t)-T) + \frac{-\Delta H}{\rho C_p}\,r_A(C_A,T) - \frac{UA}{\rho C_p V}(T-T_c). \]

We impose periodic inlet-temperature forcing \( T_f(t)=T_{f0}+A_f\sin(2\pi t/T) \), making the system \( T \)-periodic in time and enabling Poincaré analysis.

Dissipativity (absorbing set idea). From the mass balance, for \( C_A\ge 0 \) we have \( \dot{C}_A \le \frac{q}{V}C_{Af} - \frac{q}{V}C_A \), so \( C_A(t) \) is ultimately bounded by \( C_{Af} \). For temperature, the cooling term \( -\frac{UA}{\rho C_p V}(T-T_c) \) provides linear dissipation, while the reaction term is bounded when \( C_A \) is bounded and \( T \) is not too large; a standard comparison argument yields an ultimate bound for \( T(t) \) provided heat removal dominates for large \( T \).

For all three systems, a practical workflow is: (i) integrate one period to evaluate \( \mathbf{P}(\mathbf{x}) \), (ii) discard transients, (iii) retain a window of Poincaré points, and (iv) estimate a coarse largest Lyapunov exponent using two nearby trajectories with repeated renormalization.

flowchart TD
  S["Pick forcing period T and initial state x0"] --> I1["Integrate xdot=f(t,x) from t to t+T"]
  I1 --> P["Record x(t+T) as one Poincare sample"]
  P --> K{"Enough samples?"}
  K -->|no| I1
  K -->|yes| L0["Initialize nearby state x0 + d0 (small)"]
  L0 --> L1["Evolve both trajectories for dt steps"]
  L1 --> L2["Compute separation norm(d)"]
  L2 --> L3["Accumulate log(norm(d)/norm(d0))"]
  L3 --> L4["Renormalize d to norm(d0) and continue"]
        

Numerical caution. Poincaré sampling is sensitive to: (a) step size (fixed-step) or tolerances (adaptive), (b) event timing (must sample at exact multiples of \( T \)), and (c) transient length. Always verify convergence of qualitative conclusions under smaller steps/tighter tolerances.

The following script implements: models, Poincaré sampling, and a crude largest Lyapunov estimator.

# Chapter20_Lesson5.py
"""
System Dynamics — Chapter 20, Lesson 5
Integrated Case Studies and Projects:
(Mechatronic System, Vehicle Dynamics, Thermal–Fluid System)

This script provides:
1) Mechatronic DC motor + load with saturation + nonlinear friction (smooth Coulomb/Stribeck).
2) Vehicle planar bicycle model with nonlinear tire (tanh saturation) under periodic steering.
3) Thermal–fluid CSTR (mass + energy balance) with periodic inlet-temperature forcing.

For each case, we:
- define ODE f(t, x)
- simulate with solve_ivp
- compute a simple Poincaré map for periodically forced systems
- estimate a crude largest Lyapunov exponent via two-trajectory renormalization

Dependencies: numpy, scipy
"""
from __future__ import annotations
import numpy as np
from dataclasses import dataclass
from typing import Callable
from scipy.integrate import solve_ivp

def poincare_section(f: Callable[[float, np.ndarray], np.ndarray],
                     x0: np.ndarray,
                     T: float,
                     n_transient: int,
                     n_points: int,
                     t_per_period: int = 200) -> np.ndarray:
    t0 = 0.0
    x = x0.copy()
    samples = []
    for k in range(n_transient + n_points):
        sol = solve_ivp(f, (t0, t0 + T), x, max_step=T / t_per_period, rtol=1e-8, atol=1e-10)
        x = sol.y[:, -1]
        t0 += T
        if k >= n_transient:
            samples.append(x.copy())
    return np.array(samples)

def lyapunov_largest(f: Callable[[float, np.ndarray], np.ndarray],
                     x0: np.ndarray,
                     T: float,
                     n_steps: int,
                     dt: float,
                     delta0: float = 1e-8) -> float:
    rng = np.random.default_rng(0)
    d0 = rng.normal(size=x0.shape)
    d0 = delta0 * d0 / np.linalg.norm(d0)

    x = x0.copy()
    y = x0.copy() + d0
    s = 0.0
    t = 0.0

    def step(z, t, dt):
        k1 = f(t, z)
        k2 = f(t + 0.5*dt, z + 0.5*dt*k1)
        k3 = f(t + 0.5*dt, z + 0.5*dt*k2)
        k4 = f(t + dt, z + dt*k3)
        return z + (dt/6.0)*(k1 + 2*k2 + 2*k3 + k4)

    for _ in range(n_steps):
        x = step(x, t, dt)
        y = step(y, t, dt)
        d = y - x
        nd = np.linalg.norm(d)
        if nd > 0:
            s += np.log(nd / delta0)
            y = x + (delta0/nd)*d
        t += dt
        if T is not None and T > 0:
            t = t % T
    return s / (n_steps * dt)

# -----------------------------
# Case Study 1: DC motor
# -----------------------------
@dataclass
class MotorParams:
    J: float = 2e-3
    b: float = 1e-3
    Kt: float = 0.05
    Ke: float = 0.05
    R: float = 1.0
    L: float = 5e-3
    tau_c: float = 0.02
    tau_s: float = 0.03
    w_s: float = 2.0
    u_max: float = 12.0

def sat(u, umax):
    return np.clip(u, -umax, umax)

def stribeck_friction(w, p: MotorParams):
    sign = np.tanh(50.0 * w)
    tau = (p.tau_c + (p.tau_s - p.tau_c) * np.exp(-(np.abs(w)/p.w_s)**2)) * sign
    return tau

def motor_ode(p: MotorParams, u_fun: Callable[[float], float]):
    def f(t, x):
        th, w, i = x
        u = sat(u_fun(t), p.u_max)
        tau_f = p.b * w + stribeck_friction(w, p)
        dth = w
        dw = (p.Kt * i - tau_f) / p.J
        di = (u - p.R * i - p.Ke * w) / p.L
        return np.array([dth, dw, di], dtype=float)
    return f

# -----------------------------
# Case Study 2: Bicycle
# -----------------------------
@dataclass
class VehicleParams:
    m: float = 1500.0
    Iz: float = 2500.0
    a: float = 1.2
    b: float = 1.6
    Ux: float = 20.0
    Cf: float = 80000.0
    Cr: float = 90000.0
    alpha_sat: float = 0.15
    delta0: float = 0.06
    T: float = 1.0

def tire_force(alpha, C, alpha_sat):
    return C * alpha_sat * np.tanh(alpha / alpha_sat)

def bicycle_ode(p: VehicleParams):
    def delta(t):
        return p.delta0 * np.sin(2*np.pi*t/p.T)
    def f(t, x):
        vy, r = x
        alpha_f = (vy + p.a*r)/p.Ux - delta(t)
        alpha_r = (vy - p.b*r)/p.Ux
        Fyf = -tire_force(alpha_f, p.Cf, p.alpha_sat)
        Fyr = -tire_force(alpha_r, p.Cr, p.alpha_sat)
        dvy = (Fyf + Fyr)/p.m - p.Ux*r
        dr  = (p.a*Fyf - p.b*Fyr)/p.Iz
        return np.array([dvy, dr], dtype=float)
    return f

# -----------------------------
# Case Study 3: CSTR
# -----------------------------
@dataclass
class CSTRParams:
    V: float = 1.0
    rho: float = 1000.0
    Cp: float = 4180.0
    q: float = 1e-3
    CAf: float = 1.0
    Tf0: float = 300.0
    dTf: float = 5.0
    T: float = 2.0
    k0: float = 7.2e10
    E: float = 8.314e4
    Rg: float = 8.314
    dH: float = -5e7
    UA: float = 5e4
    Tc: float = 295.0

def cstr_ode(p: CSTRParams):
    def Tf(t):
        return p.Tf0 + p.dTf*np.sin(2*np.pi*t/p.T)
    def rate(CA, T):
        k = p.k0*np.exp(-p.E/(p.Rg*T))
        return k*CA
    def f(t, x):
        CA, T = x
        rA = rate(CA, T)
        dCA = (p.q/p.V)*(p.CAf - CA) - rA
        dT = (p.q/p.V)*(Tf(t) - T) \
             + (-p.dH/(p.rho*p.Cp))*rA \
             - (p.UA/(p.rho*p.Cp*p.V))*(T - p.Tc)
        return np.array([dCA, dT], dtype=float)
    return f

def main():
    # Motor
    mp = MotorParams()
    u_fun = lambda t: 8.0*np.sin(2*np.pi*0.8*t)
    f1 = motor_ode(mp, u_fun)
    x0_1 = np.array([0.0, 0.0, 0.0])
    T1 = 1.25
    P1 = poincare_section(f1, x0_1, T1, n_transient=50, n_points=20)
    lam1 = lyapunov_largest(f1, P1[-1], T1, n_steps=8000, dt=1e-3, delta0=1e-8)
    print("[Motor] Poincare last sample:", P1[-1])
    print("[Motor] crude largest Lyapunov ~", lam1, "1/s")

    # Vehicle
    vp = VehicleParams()
    f2 = bicycle_ode(vp)
    x0_2 = np.array([0.0, 0.0])
    P2 = poincare_section(f2, x0_2, vp.T, n_transient=200, n_points=50)
    lam2 = lyapunov_largest(f2, P2[-1], vp.T, n_steps=10000, dt=1e-3, delta0=1e-8)
    print("[Vehicle] Poincare last sample:", P2[-1])
    print("[Vehicle] crude largest Lyapunov ~", lam2, "1/s")

    # CSTR
    cp = CSTRParams()
    f3 = cstr_ode(cp)
    x0_3 = np.array([0.9, 305.0])
    P3 = poincare_section(f3, x0_3, cp.T, n_transient=400, n_points=60)
    lam3 = lyapunov_largest(f3, P3[-1], cp.T, n_steps=12000, dt=2e-3, delta0=1e-8)
    print("[CSTR] Poincare last sample:", P3[-1])
    print("[CSTR] crude largest Lyapunov ~", lam3, "1/s")

if __name__ == "__main__":
    main()

// Chapter20_Lesson5.cpp
/*
System Dynamics — Chapter 20, Lesson 5
Integrated Case Studies and Projects:
(Mechatronic DC motor, Vehicle bicycle dynamics, Thermal–fluid CSTR)

Self-contained RK4 + Poincaré sampling.
*/
#include <cmath>
#include <iostream>
#include <vector>
#include <array>
#include <iomanip>

static inline double sat(double u, double umax){
  if(u > umax) return umax;
  if(u < -umax) return -umax;
  return u;
}

template <size_t N, class F>
std::array<double,N> rk4_step(F f, double t, const std::array<double,N>& x, double dt){
  auto k1 = f(t, x);
  std::array<double,N> x2;
  for(size_t i=0;i<N;i++) x2[i] = x[i] + 0.5*dt*k1[i];
  auto k2 = f(t+0.5*dt, x2);
  for(size_t i=0;i<N;i++) x2[i] = x[i] + 0.5*dt*k2[i];
  auto k3 = f(t+0.5*dt, x2);
  for(size_t i=0;i<N;i++) x2[i] = x[i] + dt*k3[i];
  auto k4 = f(t+dt, x2);

  std::array<double,N> xn;
  for(size_t i=0;i<N;i++){
    xn[i] = x[i] + (dt/6.0)*(k1[i] + 2.0*k2[i] + 2.0*k3[i] + k4[i]);
  }
  return xn;
}

template <size_t N, class F>
std::array<double,N> integrate(F f, std::array<double,N> x0, double t0, double tf, double dt){
  double t = t0;
  auto x = x0;
  long steps = (long)std::ceil((tf - t0)/dt);
  for(long k=0;k<steps;k++){
    x = rk4_step<N>(f, t, x, dt);
    t += dt;
  }
  return x;
}

// Motor
struct MotorParams {
  double J=2e-3, b=1e-3, Kt=0.05, Ke=0.05, R=1.0, L=5e-3;
  double tau_c=0.02, tau_s=0.03, w_s=2.0, u_max=12.0;
};
static inline double stribeck(double w, const MotorParams& p){
  double s = std::tanh(50.0*w);
  double tau = (p.tau_c + (p.tau_s - p.tau_c)*std::exp(-std::pow(std::fabs(w)/p.w_s,2.0))) * s;
  return tau;
}
struct MotorODE {
  MotorParams p;
  double amp=8.0, freq=0.8;
  std::array<double,3> operator()(double t, const std::array<double,3>& x) const{
    double th=x[0], w=x[1], i=x[2];
    double u = sat(amp*std::sin(2.0*M_PI*freq*t), p.u_max);
    double tau_f = p.b*w + stribeck(w,p);
    return {w, (p.Kt*i - tau_f)/p.J, (u - p.R*i - p.Ke*w)/p.L};
  }
};

// Bicycle
struct VehicleParams {
  double m=1500.0, Iz=2500.0, a=1.2, b=1.6, Ux=20.0;
  double Cf=80000.0, Cr=90000.0, alpha_sat=0.15;
  double delta0=0.06, T=1.0;
};
static inline double tire(double alpha, double C, double alpha_sat){
  return C*alpha_sat*std::tanh(alpha/alpha_sat);
}
struct BicycleODE {
  VehicleParams p;
  std::array<double,2> operator()(double t, const std::array<double,2>& x) const{
    double vy=x[0], r=x[1];
    double delta = p.delta0*std::sin(2.0*M_PI*t/p.T);
    double alpha_f = (vy + p.a*r)/p.Ux - delta;
    double alpha_r = (vy - p.b*r)/p.Ux;
    double Fyf = -tire(alpha_f, p.Cf, p.alpha_sat);
    double Fyr = -tire(alpha_r, p.Cr, p.alpha_sat);
    return {(Fyf + Fyr)/p.m - p.Ux*r, (p.a*Fyf - p.b*Fyr)/p.Iz};
  }
};

// CSTR
struct CSTRParams {
  double V=1.0, rho=1000.0, Cp=4180.0, q=1e-3;
  double CAf=1.0, Tf0=300.0, dTf=5.0, T=2.0;
  double k0=7.2e10, E=8.314e4, Rg=8.314;
  double dH=-5e7, UA=5e4, Tc=295.0;
};
struct CSTRODE {
  CSTRParams p;
  std::array<double,2> operator()(double t, const std::array<double,2>& x) const{
    double CA=x[0], T=x[1];
    double Tf = p.Tf0 + p.dTf*std::sin(2.0*M_PI*t/p.T);
    double k = p.k0*std::exp(-p.E/(p.Rg*T));
    double rA = k*CA;
    double dCA = (p.q/p.V)*(p.CAf - CA) - rA;
    double dT  = (p.q/p.V)*(Tf - T)
               + (-p.dH/(p.rho*p.Cp))*rA
               - (p.UA/(p.rho*p.Cp*p.V))*(T - p.Tc);
    return {dCA,dT};
  }
};

template <size_t N, class F>
std::vector<std::array<double,N>> poincare(F f, std::array<double,N> x0, double period,
                                          int n_transient, int n_points, double dt){
  std::vector<std::array<double,N>> out;
  double t=0.0;
  auto x=x0;
  for(int k=0;k<n_transient+n_points;k++){
    x = integrate<N>(f, x, t, t+period, dt);
    t += period;
    if(k>=n_transient) out.push_back(x);
  }
  return out;
}

int main(){
  std::cout << std::setprecision(6) << std::fixed;

  MotorODE f1;
  double T1 = 1.25;
  auto P1 = poincare<3>(f1, {0.0,0.0,0.0}, T1, 50, 5, 1e-4);
  std::cout << "[Motor] Poincare samples (theta,w,i):\n";
  for(auto& s: P1) std::cout << "  " << s[0] << "  " << s[1] << "  " << s[2] << "\n";

  BicycleODE f2;
  auto P2 = poincare<2>(f2, {0.0,0.0}, f2.p.T, 200, 5, 1e-4);
  std::cout << "[Vehicle] Poincare samples (vy,r):\n";
  for(auto& s: P2) std::cout << "  " << s[0] << "  " << s[1] << "\n";

  CSTRODE f3;
  auto P3 = poincare<2>(f3, {0.9,305.0}, f3.p.T, 400, 5, 2e-4);
  std::cout << "[CSTR] Poincare samples (CA,T):\n";
  for(auto& s: P3) std::cout << "  " << s[0] << "  " << s[1] << "\n";

  return 0;
}

// Chapter20_Lesson5.java
/*
System Dynamics — Chapter 20, Lesson 5
Integrated Case Studies and Projects (Mechatronic / Vehicle / Thermal–Fluid)

Self-contained Java 17 code: RK4 + Poincaré sampling.
*/
import java.util.*;
import static java.lang.Math.*;

public class Chapter20_Lesson5 {

  interface ODE { double[] f(double t, double[] x); }

  static double[] rk4Step(ODE ode, double t, double[] x, double dt){
    double[] k1 = ode.f(t, x);
    double[] x2 = new double[x.length];
    for(int i=0;i<x.length;i++) x2[i] = x[i] + 0.5*dt*k1[i];
    double[] k2 = ode.f(t+0.5*dt, x2);
    for(int i=0;i<x.length;i++) x2[i] = x[i] + 0.5*dt*k2[i];
    double[] k3 = ode.f(t+0.5*dt, x2);
    for(int i=0;i<x.length;i++) x2[i] = x[i] + dt*k3[i];
    double[] k4 = ode.f(t+dt, x2);

    double[] xn = new double[x.length];
    for(int i=0;i<x.length;i++){
      xn[i] = x[i] + (dt/6.0)*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i]);
    }
    return xn;
  }

  static double[] integrate(ODE ode, double[] x0, double t0, double tf, double dt){
    double t=t0;
    double[] x = Arrays.copyOf(x0, x0.length);
    long steps = (long)ceil((tf - t0)/dt);
    for(long k=0;k<steps;k++){
      x = rk4Step(ode, t, x, dt);
      t += dt;
    }
    return x;
  }

  static List<double[]> poincare(ODE ode, double[] x0, double period,
                                int nTransient, int nPoints, double dt){
    List<double[]> out = new ArrayList<>();
    double t=0.0;
    double[] x = Arrays.copyOf(x0, x0.length);
    for(int k=0;k<nTransient+nPoints;k++){
      x = integrate(ode, x, t, t+period, dt);
      t += period;
      if(k>=nTransient) out.add(Arrays.copyOf(x, x.length));
    }
    return out;
  }

  static double sat(double u, double umax){ return max(-umax, min(umax, u)); }

  // Motor
  static class MotorParams{
    double J=2e-3, b=1e-3, Kt=0.05, Ke=0.05, R=1.0, L=5e-3;
    double tauC=0.02, tauS=0.03, wS=2.0, uMax=12.0;
  }
  static double stribeck(double w, MotorParams p){
    double s = tanh(50.0*w);
    return (p.tauC + (p.tauS - p.tauC)*exp(-pow(abs(w)/p.wS,2.0))) * s;
  }
  static ODE motorODE(MotorParams p, double amp, double freq){
    return (t, x) -> {
      double th=x[0], w=x[1], i=x[2];
      double u = sat(amp*sin(2.0*PI*freq*t), p.uMax);
      double tauF = p.b*w + stribeck(w,p);
      return new double[]{ w, (p.Kt*i - tauF)/p.J, (u - p.R*i - p.Ke*w)/p.L };
    };
  }

  // Bicycle
  static class VehicleParams{
    double m=1500.0, Iz=2500.0, a=1.2, b=1.6, Ux=20.0;
    double Cf=80000.0, Cr=90000.0, alphaSat=0.15;
    double delta0=0.06, T=1.0;
  }
  static double tire(double alpha, double C, double alphaSat){
    return C*alphaSat*tanh(alpha/alphaSat);
  }
  static ODE bicycleODE(VehicleParams p){
    return (t, x) -> {
      double vy=x[0], r=x[1];
      double delta = p.delta0*sin(2.0*PI*t/p.T);
      double alphaF = (vy + p.a*r)/p.Ux - delta;
      double alphaR = (vy - p.b*r)/p.Ux;
      double Fyf = -tire(alphaF, p.Cf, p.alphaSat);
      double Fyr = -tire(alphaR, p.Cr, p.alphaSat);
      double dvy = (Fyf + Fyr)/p.m - p.Ux*r;
      double dr  = (p.a*Fyf - p.b*Fyr)/p.Iz;
      return new double[]{dvy, dr};
    };
  }

  // CSTR
  static class CSTRParams{
    double V=1.0, rho=1000.0, Cp=4180.0, q=1e-3;
    double CAf=1.0, Tf0=300.0, dTf=5.0, T=2.0;
    double k0=7.2e10, E=8.314e4, Rg=8.314;
    double dH=-5e7, UA=5e4, Tc=295.0;
  }
  static ODE cstrODE(CSTRParams p){
    return (t, x) -> {
      double CA=x[0], Temp=x[1];
      double Tf = p.Tf0 + p.dTf*sin(2.0*PI*t/p.T);
      double k = p.k0*exp(-p.E/(p.Rg*Temp));
      double rA = k*CA;
      double dCA = (p.q/p.V)*(p.CAf - CA) - rA;
      double dT  = (p.q/p.V)*(Tf - Temp)
                 + (-p.dH/(p.rho*p.Cp))*rA
                 - (p.UA/(p.rho*p.Cp*p.V))*(Temp - p.Tc);
      return new double[]{dCA, dT};
    };
  }

  public static void main(String[] args){
    System.out.printf(Locale.US, "Chapter 20, Lesson 5 — Poincare samples%n");

    MotorParams mp = new MotorParams();
    ODE f1 = motorODE(mp, 8.0, 0.8);
    double T1 = 1.25;
    var P1 = poincare(f1, new double[]{0.0,0.0,0.0}, T1, 50, 5, 1e-4);
    System.out.println("[Motor] theta, w, i:");
    for(double[] s: P1) System.out.printf(Locale.US, "  %.6f  %.6f  %.6f%n", s[0], s[1], s[2]);

    VehicleParams vp = new VehicleParams();
    ODE f2 = bicycleODE(vp);
    var P2 = poincare(f2, new double[]{0.0,0.0}, vp.T, 200, 5, 1e-4);
    System.out.println("[Vehicle] vy, r:");
    for(double[] s: P2) System.out.printf(Locale.US, "  %.6f  %.6f%n", s[0], s[1]);

    CSTRParams cp = new CSTRParams();
    ODE f3 = cstrODE(cp);
    var P3 = poincare(f3, new double[]{0.9,305.0}, cp.T, 400, 5, 2e-4);
    System.out.println("[CSTR] CA, T:");
    for(double[] s: P3) System.out.printf(Locale.US, "  %.6f  %.6f%n", s[0], s[1]);
  }
}

% Chapter20_Lesson5.m
% System Dynamics — Chapter 20, Lesson 5
% Integrated Case Studies and Projects (Mechatronic / Vehicle / Thermal–Fluid)

clear; clc;

poincare = @(t, x, T, nTransient, nPoints) x(:, arrayfun(@(k) find(t >= k*T, 1, 'first'), (nTransient+1):(nTransient+nPoints)));

%% Case 1: DC motor
mp.J=2e-3; mp.b=1e-3; mp.Kt=0.05; mp.Ke=0.05; mp.R=1.0; mp.L=5e-3;
mp.tau_c=0.02; mp.tau_s=0.03; mp.w_s=2.0; mp.u_max=12.0;

u_amp=8.0; u_f=0.8;  T1=1.25;
motor_ode = @(t,x) motor_f(t,x,mp,u_amp,u_f);

tspan = [0, (50+20)*T1];
x0 = [0; 0; 0];
opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t1,x1] = ode45(motor_ode, tspan, x0, opts);
P1 = poincare(t1, x1', T1, 50, 20);
disp('[Motor] last Poincare sample [theta; w; i]:'), disp(P1(:,end));

%% Case 2: Bicycle
vp.m=1500; vp.Iz=2500; vp.a=1.2; vp.b=1.6; vp.Ux=20;
vp.Cf=80000; vp.Cr=90000; vp.alpha_sat=0.15; vp.delta0=0.06; vp.T=1.0;

bicycle_ode = @(t,x) bicycle_f(t,x,vp);

tspan = [0, (200+50)*vp.T];
x0 = [0;0];
[t2,x2] = ode45(bicycle_ode, tspan, x0, opts);
P2 = poincare(t2, x2', vp.T, 200, 50);
disp('[Vehicle] last Poincare sample [vy; r]:'), disp(P2(:,end));

%% Case 3: CSTR
cp.V=1.0; cp.rho=1000; cp.Cp=4180; cp.q=1e-3;
cp.CAf=1.0; cp.Tf0=300; cp.dTf=5; cp.T=2.0;
cp.k0=7.2e10; cp.E=8.314e4; cp.Rg=8.314;
cp.dH=-5e7; cp.UA=5e4; cp.Tc=295;

cstr_ode = @(t,x) cstr_f(t,x,cp);

tspan = [0, (400+60)*cp.T];
x0 = [0.9; 305];
[t3,x3] = ode45(cstr_ode, tspan, x0, opts);
P3 = poincare(t3, x3', cp.T, 400, 60);
disp('[CSTR] last Poincare sample [CA; T]:'), disp(P3(:,end));

%% Local functions
function dx = motor_f(t,x,p,amp,freq)
  w=x(2); i=x(3);
  u = sat(amp*sin(2*pi*freq*t), p.u_max);
  tau_f = p.b*w + stribeck(w,p);
  dx = [w; (p.Kt*i - tau_f)/p.J; (u - p.R*i - p.Ke*w)/p.L];
end

function tau = stribeck(w,p)
  s = tanh(50*w);
  tau = (p.tau_c + (p.tau_s - p.tau_c)*exp(-(abs(w)/p.w_s)^2)) * s;
end

function y = sat(u, umax)
  y = min(max(u, -umax), umax);
end

function F = tire(alpha, C, alpha_sat)
  F = C*alpha_sat*tanh(alpha/alpha_sat);
end

function dx = bicycle_f(t,x,p)
  vy=x(1); r=x(2);
  delta = p.delta0*sin(2*pi*t/p.T);
  alpha_f = (vy + p.a*r)/p.Ux - delta;
  alpha_r = (vy - p.b*r)/p.Ux;
  Fyf = -tire(alpha_f, p.Cf, p.alpha_sat);
  Fyr = -tire(alpha_r, p.Cr, p.alpha_sat);
  dx = [(Fyf + Fyr)/p.m - p.Ux*r; (p.a*Fyf - p.b*Fyr)/p.Iz];
end

function dx = cstr_f(t,x,p)
  CA=x(1); T=x(2);
  Tf = p.Tf0 + p.dTf*sin(2*pi*t/p.T);
  k = p.k0*exp(-p.E/(p.Rg*T));
  rA = k*CA;
  dCA = (p.q/p.V)*(p.CAf - CA) - rA;
  dT = (p.q/p.V)*(Tf - T) ...
     + (-p.dH/(p.rho*p.Cp))*rA ...
     - (p.UA/(p.rho*p.Cp*p.V))*(T - p.Tc);
  dx = [dCA; dT];
end

Simulink note. Each case maps naturally to an integrator chain with a nonlinear block for \( \tau_f(\omega) \), \( F_y(\alpha) \), or \( r_A(C_A,T) \), plus a periodic source for the forcing. Use a triggered subsystem (trigger at multiples of \( T \)) to log Poincaré samples.

(* Chapter20_Lesson5.nb (Mathematica package-style content) *)
ClearAll["Global`*"];

PoincareSamples[sol_, T_, nTransient_, nPoints_, vars_List] := Module[
  {ts = Table[(nTransient + k) T, {k, 1, nPoints}]},
  (vars /. sol /. t -> #) & /@ ts
];

(* Motor *)
mp = <|"J"->2.*^-3,"b"->1.*^-3,"Kt"->0.05,"Ke"->0.05,"R"->1.0,"L"->5.*^-3,
       "tauC"->0.02,"tauS"->0.03,"wS"->2.0,"uMax"->12.0|>;
sat[u_, umax_] := Max[-umax, Min[umax, u]];
stribeck[w_] := Module[{s = Tanh[50. w]},
  (mp["tauC"] + (mp["tauS"] - mp["tauC"]) Exp[-(Abs[w]/mp["wS"])^2]) s
];
uAmp=8.0; uF=0.8; T1=1.25;
u[t_] := sat[uAmp Sin[2 Pi uF t], mp["uMax"]];

eqMotor = {
  th'[t]==w[t],
  w'[t]==(mp["Kt"] i[t] - (mp["b"] w[t] + stribeck[w[t]]))/mp["J"],
  i'[t]==(u[t] - mp["R"] i[t] - mp["Ke"] w[t])/mp["L"],
  th[0]==0,w[0]==0,i[0]==0
};
solMotor = NDSolve[eqMotor, {th,w,i}, {t,0,(50+20) T1},
  Method->{"StiffnessSwitching"}, AccuracyGoal->12, PrecisionGoal->12][[1]];
P1 = PoincareSamples[solMotor, T1, 50, 10, {th[t], w[t], i[t]}];
Print["[Motor] last Poincare sample:", Last[P1]];

(* Bicycle *)
vp = <|"m"->1500.,"Iz"->2500.,"a"->1.2,"b"->1.6,"Ux"->20.,
       "Cf"->80000.,"Cr"->90000.,"alphaSat"->0.15,"delta0"->0.06,"T"->1.0|>;
tire[alpha_, C_] := C vp["alphaSat"] Tanh[alpha/vp["alphaSat"]];
delta[t_] := vp["delta0"] Sin[2 Pi t/vp["T"]];

eqVeh = {
  vy'[t]==(-tire[(vy[t]+vp["a"] r[t])/vp["Ux"]-delta[t], vp["Cf"]] -
           tire[(vy[t]-vp["b"] r[t])/vp["Ux"], vp["Cr"]])/vp["m"] - vp["Ux"] r[t],
  r'[t]==(vp["a"]*(-tire[(vy[t]+vp["a"] r[t])/vp["Ux"]-delta[t], vp["Cf"]]) -
          vp["b"]*(-tire[(vy[t]-vp["b"] r[t])/vp["Ux"], vp["Cr"]]))/vp["Iz"],
  vy[0]==0,r[0]==0
};
solVeh = NDSolve[eqVeh, {vy,r}, {t,0,(200+50) vp["T"]},
  Method->{"StiffnessSwitching"}, AccuracyGoal->12, PrecisionGoal->12][[1]];
P2 = PoincareSamples[solVeh, vp["T"], 200, 10, {vy[t], r[t]}];
Print["[Vehicle] last Poincare sample:", Last[P2]];

(* CSTR *)
cp = <|"V"->1.0,"rho"->1000.,"Cp"->4180.,"q"->1.*^-3,
       "CAf"->1.0,"Tf0"->300.,"dTf"->5.,"T"->2.0,
       "k0"->7.2*^10,"E"->8.314*^4,"Rg"->8.314,
       "dH"->-5.*^7,"UA"->5.*^4,"Tc"->295.|>;
Tf[t_] := cp["Tf0"] + cp["dTf"] Sin[2 Pi t/cp["T"]];
kArr[T_] := cp["k0"] Exp[-cp["E"]/(cp["Rg"] T)];
rA[CA_, T_] := kArr[T] CA;

eqCSTR = {
  CA'[t]==(cp["q"]/cp["V"])(cp["CAf"]-CA[t]) - rA[CA[t], Temp[t]],
  Temp'[t]==(cp["q"]/cp["V"])(Tf[t]-Temp[t]) +
            (-cp["dH"]/(cp["rho"] cp["Cp"])) rA[CA[t], Temp[t]] -
            (cp["UA"]/(cp["rho"] cp["Cp"] cp["V"])) (Temp[t]-cp["Tc"]),
  CA[0]==0.9, Temp[0]==305.
};
solCSTR = NDSolve[eqCSTR, {CA,Temp}, {t,0,(400+60) cp["T"]},
  Method->{"StiffnessSwitching"}, AccuracyGoal->12, PrecisionGoal->12][[1]];
P3 = PoincareSamples[solCSTR, cp["T"], 400, 10, {CA[t], Temp[t]}];
Print["[CSTR] last Poincare sample:", Last[P3]];

Problem 1 (Poincaré fixed point ↔ periodic orbit): Let \( \dot{\mathbf{x}}=\mathbf{f}(t,\mathbf{x}) \) be \( T \)-periodic in time. Prove that \( \mathbf{x}^\star \) is a fixed point of the Poincaré map \( \mathbf{P}(\mathbf{x})=\varphi(T;0,\mathbf{x}) \) if and only if there exists a \( T \)-periodic solution of the ODE with \( \mathbf{x}(0)=\mathbf{x}^\star \).

Solution: If \( \mathbf{P}(\mathbf{x}^\star)=\mathbf{x}^\star \), then by definition \( \varphi(T;0,\mathbf{x}^\star)=\mathbf{x}^\star \), so the solution returns to the same state after one period. Because the vector field is \( T \)-periodic, the continuation from time \( T \) repeats the same evolution as from time \( 0 \), yielding \( \mathbf{x}(t+T)=\mathbf{x}(t) \). Conversely, if a solution is \( T \)-periodic with \( \mathbf{x}(0)=\mathbf{x}^\star \), then \( \mathbf{x}(T)=\mathbf{x}(0)=\mathbf{x}^\star \) so \( \mathbf{P}(\mathbf{x}^\star)=\mathbf{x}^\star \).

Problem 2 (Motor dissipativity estimate): Assume \( K_t=K_e \) and the friction satisfies \( \omega\tau_f(\omega)\ge b\omega^2 \) for all \( \omega \). Show that with bounded input \( |u(t)|\le u_{\max} \), the energy-like function \( V(\omega,i)=\tfrac{1}{2}J\omega^2+\tfrac{1}{2}Li^2 \) satisfies \( \dot{V} \le -c_1 V + c_2 \) for some constants \( c_1>0 \), \( c_2>0 \).

Solution: From the derivation, \( \dot{V} = -\omega\tau_f(\omega) - Ri^2 + iu \le -b\omega^2 - Ri^2 + |i|u_{\max} \). Apply Young: \( |i|u_{\max} \le \tfrac{R}{2}i^2 + \tfrac{1}{2R}u_{\max}^2 \). Then \( \dot{V} \le -b\omega^2 - \tfrac{R}{2}i^2 + \tfrac{1}{2R}u_{\max}^2 \). Since \( V \ge \tfrac{1}{2}\min(J,L)(\omega^2+i^2) \) and \( b\omega^2 + \tfrac{R}{2}i^2 \ge c(\omega^2+i^2) \) for \( c=\min(b,R/2) \), we obtain \( \dot{V} \le -\frac{2c}{\min(J,L)}V + \tfrac{1}{2R}u_{\max}^2 \). Hence choose \( c_1=\frac{2c}{\min(J,L)} \), \( c_2=\tfrac{1}{2R}u_{\max}^2 \).

Problem 3 (Linearization of saturated tire model): Consider \( F_y(\alpha)=C\alpha_s\tanh(\alpha/\alpha_s) \). Show that for small \( |\alpha| \), \( F_y(\alpha)=C\alpha+O(\alpha^3) \), and compute \( \frac{dF_y}{d\alpha}(0) \).

Solution: Using the Taylor expansion \( \tanh(z)=z-\tfrac{1}{3}z^3+O(z^5) \) with \( z=\alpha/\alpha_s \), we get \( F_y=C\alpha_s\left(\alpha/\alpha_s - \tfrac{1}{3}(\alpha/\alpha_s)^3 + O(\alpha^5)\right)=C\alpha - \tfrac{C}{3\alpha_s^2}\alpha^3 + O(\alpha^5) \). Also \( \frac{dF_y}{d\alpha} = C\,\mathrm{sech}^2(\alpha/\alpha_s) \), so \( \frac{dF_y}{d\alpha}(0)=C \).

Problem 4 (CSTR equilibrium and local stability): For the unforced case \( T_f(t)=T_{f0} \), derive the equilibrium equations and explain why multiple equilibria can occur.

Solution: Equilibrium satisfies \( 0=\frac{q}{V}(C_{Af}-C_A^\star)-k(T^\star)C_A^\star \), giving \( C_A^\star = \frac{\frac{q}{V}C_{Af}}{\frac{q}{V}+k(T^\star)} \). Substitute into energy balance equilibrium \( 0=\frac{q}{V}(T_{f0}-T^\star)+\frac{-\Delta H}{\rho C_p}k(T^\star)C_A^\star-\frac{UA}{\rho C_p V}(T^\star-T_c) \). Because \( k(T)=k_0\exp(-E/(R_g T)) \) increases rapidly with \( T \), the heat generation term can intersect the heat removal line multiple times, yielding multiple equilibria and potential saddle-node bifurcations under parameter variation.

Problem 5 (Discrete-time reduction via stroboscopic sampling): Suppose you observe a periodically forced system at times \( t_k=kT \). Explain why chaos in the continuous system can be detected from the discrete sequence \( \{\mathbf{x}(t_k)\} \).

Solution: The sequence satisfies \( \mathbf{x}_{k+1}=\mathbf{P}(\mathbf{x}_k) \), where \( \mathbf{P} \) is the Poincaré map. If the continuous-time system has a strange attractor under periodic forcing, its intersection with the stroboscopic section is an invariant set for \( \mathbf{P} \). Sensitive dependence, dense periodic points, and topological mixing (properties associated with chaos) can manifest in the discrete map; thus analyzing \( \mathbf{P} \) (e.g., via Lyapunov exponent estimates or return maps) reveals chaos without continuous-time plotting.

We integrated the full nonlinear modeling and analysis loop across three domains: mechatronics (motor with saturation and nonlinear friction), vehicle dynamics (bicycle model with saturated tire forces), and thermal–fluid dynamics (forced CSTR with Arrhenius kinetics). We formalized periodic forcing via Poincaré maps, linked dissipativity to boundedness and invariant sets, and provided multi-language implementations for simulation, stroboscopic sampling, and coarse Lyapunov estimation.