In Chapters 15–16 you learned numerical simulation basics (time stepping, stability, error control) and in Lessons 1–3 of this chapter you met canonical chaotic systems (Logistic map and Lorenz system), bifurcations, and Lyapunov exponents at an introductory level. Here we focus on computational tooling: how to encode models, select solvers, validate simulations, and run systematic analyses.

The three tool families differ mainly in model representation:

• MATLAB (ODE suite): user provides \( \dot{\mathbf{x}}=\mathbf{f}(t,\mathbf{x}) \) as a function; solver chooses adaptive steps.
• Simulink: block-diagram realization of integrators, algebraic blocks, and event logic; solver integrates the implied ODE/DAE.
• Python (SciPy): user provides \( \mathbf{f} \) to solve_ivp (adaptive explicit or implicit methods); ecosystems for reproducibility (NumPy, pandas, Jupyter).
• Modelica: user writes equations (often DAEs), the compiler performs symbolic processing (index reduction, tearing), then a numeric solver runs.

flowchart TD
  A["Start: choose model type (ODE / DAE / map)"] --> B["Write model equations and parameters"]
  B --> C["Choose tool: MATLAB, Simulink, Python, Modelica"]
  C --> D["Select solver class: explicit RK / implicit stiff / DAE"]
  D --> E["Set tolerances: rtol, atol, max step"]
  E --> F["Run simulation; log states and events to files"]
  F --> G["Verification: step refinement, invariants, sensitivity checks"]
  G --> H["Analysis: Poincare section, bifurcation sweep, Lyapunov estimate"]
  H --> I["Report: plots, tables, reproducible code + seeds"]
        

The scientific goal is not “getting a trajectory” but producing a reliable computational experiment: controlled error, reproducible settings, and validation checks that separate numerical artifacts from dynamical behavior.

Most simulation tools act on one of the following model classes.

ODE (flow) model. A parameterized initial value problem (IVP) is

\[ \dot{\mathbf{x}}(t)=\mathbf{f}\!\left(t,\mathbf{x}(t);\boldsymbol{\theta}\right), \quad \mathbf{x}(t_0)=\mathbf{x}_0, \quad \mathbf{x}\in\mathbb{R}^n. \]

The exact solution defines a flow map \( \Phi_{t,t_0} \) such that \( \mathbf{x}(t)=\Phi_{t,t_0}(\mathbf{x}_0) \). In chaotic systems, the flow exhibits sensitive dependence; therefore numerical error control is essential.

DAE (equation-based) model. Many physical models lead to algebraic constraints plus dynamics:

\[ \mathbf{F}\!\left(t,\mathbf{x}(t),\dot{\mathbf{x}}(t)\right)=\mathbf{0}. \]

Here \( \mathbf{x} \) may include algebraic variables, and initialization requires a consistent pair \( (\mathbf{x}(t_0),\dot{\mathbf{x}}(t_0)) \). Modelica is built around this class; toolchains often apply symbolic transformations before numeric solution.

Discrete-time map. A nonlinear map (e.g., Logistic map from Lesson 1) is

\[ \mathbf{x}_{k+1}=\mathbf{G}\!\left(\mathbf{x}_k;\boldsymbol{\mu}\right), \quad k\in\mathbb{Z}_{\ge 0}. \]

A discrete map can be produced by sampling a flow at fixed period \( T_s \): \( \mathbf{x}_{k+1}=\Phi_{t_0+(k+1)T_s,t_0+kT_s}(\mathbf{x}_k) \), or by Poincaré return maps. Tooling must support event detection to compute such maps reliably.

Most general-purpose solvers can be described as one-step methods producing \( \mathbf{x}_{n+1}=\Psi_h(t_n,\mathbf{x}_n) \) with step size \( h \). For an order-\( p \) method, the local truncation error scales like \( \mathcal{O}(h^{p+1}) \).

Embedded Runge–Kutta pairs. Popular solvers (MATLAB ode45, SciPy RK45) use an embedded pair producing two approximations of different orders, e.g., order \( p \) and \( p-1 \): \( \mathbf{x}_{n+1}^{(p)} \) and \( \mathbf{x}_{n+1}^{(p-1)} \). The difference serves as an error estimator: \( \mathbf{e}_n \approx \mathbf{x}_{n+1}^{(p)}-\mathbf{x}_{n+1}^{(p-1)} \).

A typical adaptive controller chooses the next step size \( h_{\text{new}} \) using (componentwise) scaling by absolute and relative tolerances.

\[ h_{\text{new}} = s\, h \left(\frac{\text{tol}}{\|\mathbf{e}_n\|_{\text{scaled}}}\right)^{1/p}, \quad \|\mathbf{e}_n\|_{\text{scaled}} = \left\|\frac{\mathbf{e}_n}{\mathbf{a}+\mathbf{r}\odot \max(|\mathbf{x}_n|,|\mathbf{x}_{n+1}|)}\right\| . \]

Convergence theorem (sketch). Suppose \( \mathbf{f} \) is globally Lipschitz in \( \mathbf{x} \) with constant \( L \). Assume a one-step method has local error \( \|\boldsymbol{\tau}_{n+1}\|\le C h^{p+1} \). Then the global error satisfies \( \|\mathbf{x}(t_n)-\mathbf{x}_n\|=\mathcal{O}(h^p) \).

Proof idea: Let \( \mathbf{e}_n=\mathbf{x}(t_n)-\mathbf{x}_n \). One derives a recursion \( \|\mathbf{e}_{n+1}\|\le (1+Lh)\|\mathbf{e}_n\| + C h^{p+1} \). Applying the discrete Grönwall inequality yields \( \|\mathbf{e}_n\|\le \exp(L(t_n-t_0))\,\mathcal{O}(h^p) \). In chaotic flows, the factor \( \exp(Lt) \) reflects exponential separation; thus long-time accuracy requires careful diagnostics.

Toolboxes expose multiple solver families because nonlinear systems may show stiff and nonstiff regimes, sometimes within the same model (hybrid events, constraints, fast/slow modes). Stability of time integration is understood via the scalar test equation \( y'(t)=\lambda y(t) \), whose exact step is \( y_{n+1}=e^{\lambda h}y_n \).

For a one-step method, the numerical update is \( y_{n+1}=R(z)\,y_n \), where \( z=\lambda h \) and \( R(z) \) is the stability function. The absolute stability region is \( \{z\in\mathbb{C}:\;|R(z)|<1\} \). Stiff problems need methods stable for large negative \( z \).

A-stability. A method is A-stable if it is stable for all \( \operatorname{Re}(z)<0 \). Many implicit methods (e.g., BDF family and implicit Runge–Kutta) have strong stability properties, hence MATLAB offers ode15s (variable-order BDF) and SciPy offers Radau/BDF in solve_ivp.

Practical implication: use explicit RK (ode45/RK45) for smooth nonstiff dynamics (common for Lorenz at standard parameters), and switch to stiff solvers for DAEs or problems with widely separated time scales.

MATLAB solvers accept a function handle for \( \mathbf{f}(t,\mathbf{x}) \) and options via odeset. A key capability for chaos workflows is event detection, used to define Poincaré sections or terminate integration on constraints.

An event is defined by a scalar function \( g(t,\mathbf{x}) \) (e.g., \( g(t,\mathbf{x})=x_1 \)). The solver locates times \( t_e \) where \( g(t_e,\mathbf{x}(t_e))=0 \) using root-finding on an interpolant.

Key point: Event detection interacts with step-size control; tight tolerances improve event time accuracy, but cost increases.

File: Chapter20_Lesson4.m

% Chapter20_Lesson4.m
%{
Chapter 20 — Chaos, Complex Dynamics, and Computational Tools
Lesson 4 — Computational Tools

MATLAB demonstration:
1) Lorenz ODE with ode45 (nonstiff) and ode15s (stiff-capable)
2) Event detection for Poincaré section (x=0 crossings, increasing direction)
3) Logistic map bifurcation tail data + Lyapunov exponent curve
4) Optional: Programmatic creation of a simple Simulink model for Lorenz

Outputs: CSV files written to a local folder "outputs_ch20_l4_matlab".
%}

function Chapter20_Lesson4()
  outDir = "outputs_ch20_l4_matlab";
  if ~exist(outDir, 'dir'), mkdir(outDir); end

  % Parameters
  p.sigma = 10.0;
  p.rho   = 28.0;
  p.beta  = 8.0/3.0;

  x0 = [1;1;1];
  tspan = [0, 40];

  % --- ODE45 with event detection ---
  opts = odeset('RelTol',1e-9,'AbsTol',1e-12,'MaxStep',0.02,'Events',@(t,x) ev_x0(t,x,p));
  [t, x, te, xe] = ode45(@(t,x) lorenz_rhs(t,x,p), tspan, x0, opts);

  writematrix([t, x], fullfile(outDir, "lorenz_traj.csv"));

  if ~isempty(xe)
    writematrix(xe, fullfile(outDir, "lorenz_poincare_x0.csv"));
  end

  % --- Compare with ode15s (useful when stiffness appears) ---
  % Lorenz is not stiff in this parameter set, but we show usage:
  opts15 = odeset('RelTol',1e-9,'AbsTol',1e-12,'MaxStep',0.02);
  [t15, x15] = ode15s(@(t,x) lorenz_rhs(t,x,p), tspan, x0, opts15); %#ok<ASGLU>
  % (Not exported to keep files minimal.)

  % --- Logistic map bifurcation data ---
  rVals = linspace(2.5, 4.0, 400);
  nTrans = 800; nKeep = 120;
  bif = zeros(numel(rVals)*nKeep, 2);
  idx = 1;
  for i = 1:numel(rVals)
    r = rVals(i);
    x = 0.2;
    for k = 1:nTrans, x = r*x*(1-x); end
    for k = 1:nKeep
      x = r*x*(1-x);
      bif(idx,:) = [r, x];
      idx = idx + 1;
    end
  end
  writematrix(bif, fullfile(outDir, "logistic_bifurcation.csv"));

  % --- Logistic map Lyapunov curve ---
  rVals2 = linspace(2.8, 4.0, 200);
  ly = zeros(numel(rVals2), 2);
  for i = 1:numel(rVals2)
    r = rVals2(i);
    ly(i,:) = [r, logistic_lyapunov(r, 0.2, 1000, 4000)];
  end
  writematrix(ly, fullfile(outDir, "logistic_lyapunov.csv"));

  % --- Optional: build a Simulink model programmatically ---
  % If you prefer manual building, create 3 Integrator blocks and one MATLAB Function
  % that outputs dx,dy,dz and wire them as xdot -> Integrator -> x.
  if exist('simulink', 'file')
    try
      build_simulink_lorenz(p);
    catch ME
      fprintf("Simulink model build failed: %s\n", ME.message);
    end
  end

  disp("Done. Outputs written to: " + outDir);
end

function dx = lorenz_rhs(~, x, p)
  X = x(1); Y = x(2); Z = x(3);
  dx = [ p.sigma*(Y - X);
         X*(p.rho - Z) - Y;
         X*Y - p.beta*Z ];
end

function [value, isterminal, direction] = ev_x0(~, x, ~)
  value = x(1);        % x == 0
  isterminal = 0;      % do not stop
  direction = +1;      % negative -> positive only
end

function lam = logistic_lyapunov(r, x0, nTrans, n)
  x = x0;
  for k = 1:nTrans, x = r*x*(1-x); end
  s = 0.0;
  for k = 1:n
    x = r*x*(1-x);
    d = abs(r*(1 - 2*x));
    if d < 1e-300, d = 1e-300; end
    s = s + log(d);
  end
  lam = s / n;
end

function build_simulink_lorenz(p)
  mdl = "Chapter20_Lesson4_Lorenz";
  if bdIsLoaded(mdl), close_system(mdl, 0); end
  new_system(mdl); open_system(mdl);

  % Add 3 integrators for x,y,z
  add_block('simulink/Continuous/Integrator', mdl + "/IntX", 'Position',[200 60 230 90]);
  add_block('simulink/Continuous/Integrator', mdl + "/IntY", 'Position',[200 130 230 160]);
  add_block('simulink/Continuous/Integrator', mdl + "/IntZ", 'Position',[200 200 230 230]);

  % Set initial conditions
  set_param(mdl + "/IntX", 'InitialCondition', '1');
  set_param(mdl + "/IntY", 'InitialCondition', '1');
  set_param(mdl + "/IntZ", 'InitialCondition', '1');

  % MATLAB Function block for derivatives
  add_block('simulink/User-Defined Functions/MATLAB Function', mdl + "/LorenzF", 'Position',[60 110 140 180]);
  code = sprintf([ ...
    "function dx = f(x)\n" ...
    "%% x = [X;Y;Z]\n" ...
    "sigma = %.17g; rho = %.17g; beta = %.17g;\n" ...
    "X = x(1); Y = x(2); Z = x(3);\n" ...
    "dx = [ sigma*(Y-X); X*(rho-Z)-Y; X*Y - beta*Z ];\n" ...
    "end\n"], p.sigma, p.rho, p.beta);
  set_param(mdl + "/LorenzF", 'Script', code);

  % Mux and Demux
  add_block('simulink/Signal Routing/Mux', mdl + "/Mux", 'Inputs','3','Position',[160 110 180 180]);
  add_block('simulink/Signal Routing/Demux', mdl + "/Demux", 'Outputs','3','Position',[150 60 170 230]);

  % Connect integrator outputs -> mux -> function -> demux -> integrator inputs
  add_line(mdl, "IntX/1", "Mux/1");
  add_line(mdl, "IntY/1", "Mux/2");
  add_line(mdl, "IntZ/1", "Mux/3");
  add_line(mdl, "Mux/1", "LorenzF/1");
  add_line(mdl, "LorenzF/1", "Demux/1");
  add_line(mdl, "Demux/1", "IntX/1");
  add_line(mdl, "Demux/2", "IntY/1");
  add_line(mdl, "Demux/3", "IntZ/1");

  % Add scopes
  add_block('simulink/Sinks/Scope', mdl + "/Scope", 'Position',[270 120 310 160]);
  add_line(mdl, "Mux/1", "Scope/1");

  set_param(mdl, 'StopTime', '40');
  save_system(mdl);
end
      

In SciPy, solve_ivp implements adaptive explicit solvers (e.g., RK45) and implicit solvers for stiff systems (e.g., Radau, BDF). The workflow we use for chaos studies is:

• Integrate trajectories with tight tolerances and bounded maximum step.
• Use events to compute Poincaré sections / return maps.
• Export to CSV so downstream plotting (Python/MATLAB) is decoupled from simulation.
• Use structured code (functions, dataclasses) for parameter sweeps and reproducibility.

File: Chapter20_Lesson4.py

# Chapter20_Lesson4.py
"""
Chapter 20 — Chaos, Complex Dynamics, and Computational Tools
Lesson 4 — Computational Tools: MATLAB/Simulink, Python (SciPy), and Modelica for System Dynamics

This script demonstrates:
1) Continuous-time ODE simulation (Lorenz system) with SciPy solve_ivp (adaptive RK and stiff solvers)
2) Event detection for Poincaré sections
3) Largest Lyapunov exponent estimation (Benettin-style) for the Lorenz flow (intro level)
4) Discrete-time map simulation (logistic map), bifurcation data, and Lyapunov exponent of the map
5) Reproducible export to CSV for post-processing / plotting in any tool

Dependencies:
- numpy
- scipy
"""

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

import numpy as np
from scipy.integrate import solve_ivp


# ----------------------------
# 1) Lorenz system (ODE)
# ----------------------------
@dataclass
class LorenzParams:
    sigma: float = 10.0
    rho: float = 28.0
    beta: float = 8.0 / 3.0


def lorenz_rhs(t: float, x: np.ndarray, p: LorenzParams) -> np.ndarray:
    """dx/dt = f(x)."""
    X, Y, Z = x
    return np.array([
        p.sigma * (Y - X),
        X * (p.rho - Z) - Y,
        X * Y - p.beta * Z
    ], dtype=float)


def lorenz_jacobian(x: np.ndarray, p: LorenzParams) -> np.ndarray:
    """Jacobian J = df/dx."""
    X, Y, Z = x
    return np.array([
        [-p.sigma, p.sigma, 0.0],
        [p.rho - Z, -1.0, -X],
        [Y, X, -p.beta]
    ], dtype=float)


def event_poincare_x0(t: float, x: np.ndarray, p: LorenzParams) -> float:
    """Event: x crosses 0 with positive direction."""
    return x[0]


event_poincare_x0.terminal = False
event_poincare_x0.direction = 1.0


def simulate_lorenz(
    x0: np.ndarray,
    t_span: Tuple[float, float],
    p: LorenzParams,
    method: str = "RK45",
    rtol: float = 1e-9,
    atol: float = 1e-12,
    max_step: float = 0.05,
) -> solve_ivp:
    """Adaptive integration with optional event handling."""
    fun = lambda t, x: lorenz_rhs(t, x, p)
    ev = lambda t, x: event_poincare_x0(t, x, p)
    sol = solve_ivp(
        fun=fun,
        t_span=t_span,
        y0=x0,
        method=method,
        rtol=rtol,
        atol=atol,
        max_step=max_step,
        events=ev,
        dense_output=True
    )
    if not sol.success:
        raise RuntimeError(f"Integration failed: {sol.message}")
    return sol


def export_trajectory_csv(path: str, t: np.ndarray, y: np.ndarray) -> None:
    """Write trajectory to CSV: t,x,y,z."""
    with open(path, "w", newline="") as f:
        w = csv.writer(f)
        w.writerow(["t", "x", "y", "z"])
        for i in range(t.size):
            w.writerow([float(t[i]), float(y[0, i]), float(y[1, i]), float(y[2, i])])


def export_points_csv(path: str, pts: np.ndarray, header=("x", "y", "z")) -> None:
    with open(path, "w", newline="") as f:
        w = csv.writer(f)
        w.writerow(list(header))
        for row in pts:
            w.writerow([float(v) for v in row])


# ----------------------------
# 2) Largest Lyapunov exponent (Lorenz, intro level)
#    Benettin-style: integrate tangent dynamics and renormalize.
# ----------------------------
def lce_lorenz_benettin(
    x0: np.ndarray,
    p: LorenzParams,
    t_transient: float = 5.0,
    t_total: float = 50.0,
    dt_renorm: float = 0.05,
    method: str = "RK45",
) -> float:
    """
    Estimate largest Lyapunov exponent λ_max.

    We integrate the extended system:
      x' = f(x)
      v' = J(x) v
    and periodically renormalize v to avoid overflow.

    Notes:
    - This is an educational implementation; for research-grade estimation,
      use multiple tangent vectors and QR orthonormalization.
    """
    # Transient integration to approach attractor
    sol_tr = simulate_lorenz(x0, (0.0, t_transient), p, method=method, max_step=dt_renorm)
    x = sol_tr.y[:, -1].astype(float)

    # Initial tangent vector
    v = np.array([1.0, 0.0, 0.0], dtype=float)
    v /= np.linalg.norm(v)

    # Integrate in chunks of dt_renorm
    t = 0.0
    sum_log = 0.0
    n = 0

    def ext_rhs(t_local: float, z: np.ndarray) -> np.ndarray:
        x_local = z[0:3]
        v_local = z[3:6]
        dx = lorenz_rhs(t_local, x_local, p)
        J = lorenz_jacobian(x_local, p)
        dv = J @ v_local
        return np.concatenate([dx, dv])

    while t < t_total:
        z0 = np.concatenate([x, v])
        sol = solve_ivp(
            fun=ext_rhs,
            t_span=(0.0, dt_renorm),
            y0=z0,
            method=method,
            rtol=1e-9,
            atol=1e-12,
            max_step=dt_renorm
        )
        if not sol.success:
            raise RuntimeError(f"Extended integration failed: {sol.message}")
        z1 = sol.y[:, -1]
        x = z1[0:3]
        v = z1[3:6]
        norm_v = float(np.linalg.norm(v))
        if norm_v == 0.0:
            # extremely unlikely; re-seed
            v = np.array([1.0, 0.0, 0.0], dtype=float)
            norm_v = 1.0
        sum_log += math.log(norm_v)
        n += 1
        v = v / norm_v
        t += dt_renorm

    return sum_log / (n * dt_renorm)


# ----------------------------
# 3) Logistic map (discrete-time)
# ----------------------------
def logistic_map(r: float, x0: float, n: int) -> np.ndarray:
    """x_{k+1} = r x_k (1 - x_k), returns array length n+1."""
    xs = np.empty(n + 1, dtype=float)
    x = float(x0)
    xs[0] = x
    for k in range(n):
        x = r * x * (1.0 - x)
        xs[k + 1] = x
    return xs


def logistic_bifurcation_data(
    r_values: np.ndarray,
    x0: float = 0.2,
    n_transient: int = 1000,
    n_keep: int = 200,
) -> np.ndarray:
    """
    Return (r, x) pairs after discarding transient.
    Useful for bifurcation diagrams.
    """
    rows = []
    for r in r_values:
        xs = logistic_map(float(r), x0, n_transient + n_keep)
        tail = xs[-n_keep:]
        for x in tail:
            rows.append((float(r), float(x)))
    return np.array(rows, dtype=float)


def logistic_lyapunov(r: float, x0: float = 0.2, n_transient: int = 1000, n: int = 5000) -> float:
    """
    Largest Lyapunov exponent for logistic map:
      λ = lim (1/N) Σ ln | f'(x_k) |, with f'(x)=r(1-2x)
    """
    xs = logistic_map(r, x0, n_transient + n)
    xs = xs[n_transient:]  # drop transient
    deriv = np.abs(r * (1.0 - 2.0 * xs))
    # Avoid log(0)
    deriv = np.maximum(deriv, 1e-300)
    return float(np.mean(np.log(deriv)))


# ----------------------------
# 4) Main demo
# ----------------------------
def main() -> None:
    out_dir = "outputs_ch20_l4"
    import os
    os.makedirs(out_dir, exist_ok=True)

    p = LorenzParams(sigma=10.0, rho=28.0, beta=8.0/3.0)
    x0 = np.array([1.0, 1.0, 1.0], dtype=float)

    # Lorenz simulation with adaptive RK (Dormand–Prince-like)
    sol = simulate_lorenz(
        x0=x0,
        t_span=(0.0, 40.0),
        p=p,
        method="RK45",
        max_step=0.02
    )
    export_trajectory_csv(os.path.join(out_dir, "lorenz_traj.csv"), sol.t, sol.y)

    # Poincaré section points (events at x=0 crossing)
    if sol.t_events and len(sol.t_events[0]) > 0:
        t_ev = sol.t_events[0]
        y_ev = sol.sol(t_ev) if sol.sol is not None else None
        if y_ev is not None:
            pts = y_ev.T  # rows
            export_points_csv(os.path.join(out_dir, "lorenz_poincare_x0.csv"), pts, header=("x", "y", "z"))

    # Largest Lyapunov exponent (educational)
    lce = lce_lorenz_benettin(x0=x0, p=p, t_transient=5.0, t_total=60.0, dt_renorm=0.05)
    with open(os.path.join(out_dir, "lorenz_lce.txt"), "w") as f:
        f.write(f"lambda_max ~ {lce:.6f}\n")

    # Logistic map: bifurcation data and Lyapunov
    r_vals = np.linspace(2.5, 4.0, 400)
    bif = logistic_bifurcation_data(r_vals, x0=0.2, n_transient=800, n_keep=120)
    np.savetxt(os.path.join(out_dir, "logistic_bifurcation.csv"), bif, delimiter=",", header="r,x", comments="")

    lyap_rows = []
    for r in np.linspace(2.8, 4.0, 200):
        lyap_rows.append((float(r), logistic_lyapunov(float(r), x0=0.2, n_transient=1000, n=4000)))
    np.savetxt(os.path.join(out_dir, "logistic_lyapunov.csv"), np.array(lyap_rows), delimiter=",", header="r,lambda", comments="")

    print("Done. Outputs written to:", out_dir)
    print("Lorenz LCE estimate:", lce)


if __name__ == "__main__":
    main()
      

Exercise implementation (adaptive RKF45 from scratch): this mirrors the embedded-pair principle behind many production solvers.

File: Chapter20_Lesson4_Ex1.py

# Chapter20_Lesson4_Ex1.py
"""
Exercise 1: Implement an embedded Runge–Kutta–Fehlberg 4(5) method (RKF45)
with adaptive step size control for a scalar test problem.

Goal:
- Understand how error estimators drive step-size selection (a core idea behind
  MATLAB ODE45 / SciPy RK45 style solvers).

We solve: y' = -y, y(0)=1 on [0, T].
Exact solution: y(t)=exp(-t).
"""

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


def rkf45_step(f: Callable[[float, float], float], t: float, y: float, h: float) -> Tuple[float, float]:
    """One RKF45 step. Returns (y5, err_est) where err_est = y5 - y4."""
    k1 = h * f(t, y)
    k2 = h * f(t + h/4, y + k1/4)
    k3 = h * f(t + 3*h/8, y + 3*k1/32 + 9*k2/32)
    k4 = h * f(t + 12*h/13, y + 1932*k1/2197 - 7200*k2/2197 + 7296*k3/2197)
    k5 = h * f(t + h, y + 439*k1/216 - 8*k2 + 3680*k3/513 - 845*k4/4104)
    k6 = h * f(t + h/2, y - 8*k1/27 + 2*k2 - 3544*k3/2565 + 1859*k4/4104 - 11*k5/40)

    # 4th order
    y4 = y + (25*k1/216) + (1408*k3/2565) + (2197*k4/4104) - (k5/5)

    # 5th order
    y5 = y + (16*k1/135) + (6656*k3/12825) + (28561*k4/56430) - (9*k5/50) + (2*k6/55)

    return y5, (y5 - y4)


def integrate_rkf45(
    f: Callable[[float, float], float],
    t0: float,
    y0: float,
    T: float,
    h0: float = 0.1,
    rtol: float = 1e-6,
    atol: float = 1e-9,
) -> Tuple[List[float], List[float]]:
    t = t0
    y = y0
    h = h0

    ts = [t]
    ys = [y]

    safety = 0.9
    p = 5  # higher order solution used
    h_min = 1e-10

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

        y_new, err = rkf45_step(f, t, y, h)
        scale = atol + rtol * max(abs(y), abs(y_new))
        err_norm = abs(err) / scale

        if err_norm <= 1.0:
            # accept
            t += h
            y = y_new
            ts.append(t)
            ys.append(y)

        # adapt step
        if err_norm == 0.0:
            h *= 2.0
        else:
            h = safety * h * (1.0 / err_norm) ** (1.0 / p)

        if h < h_min:
            raise RuntimeError("Step size underflow; tolerance too strict?")

    return ts, ys


def main() -> None:
    f = lambda t, y: -y
    ts, ys = integrate_rkf45(f, 0.0, 1.0, 10.0, h0=0.2, rtol=1e-6, atol=1e-9)

    # simple report
    y_exact = math.exp(-ts[-1])
    print("t_final =", ts[-1])
    print("y_final =", ys[-1])
    print("y_exact =", y_exact)
    print("abs error =", abs(ys[-1] - y_exact))
    print("steps =", len(ts)-1)


if __name__ == "__main__":
    main()
      

For system dynamics engineers, C++ is often used for embedded simulation, real-time HIL, or custom solvers. Here we implement a minimal fixed-step RK4 integrator (for learning) and produce data compatible with MATLAB/Python plotting. Production systems typically replace this with tested libraries (e.g., SUNDIALS, Boost.Odeint), but implementing RK4 once is valuable.

File: Chapter20_Lesson4.cpp

// Chapter20_Lesson4.cpp
/*
Chapter 20 — Chaos, Complex Dynamics, and Computational Tools
Lesson 4 — Computational Tools

C++ (standard library only) demonstration:
- Lorenz ODE simulated with fixed-step RK4
- Simple Poincaré section detection (x crosses 0 with positive direction)
- Logistic map bifurcation tail data and Lyapunov exponent for the map
- CSV export for plotting in MATLAB/Python/etc.

Build:
  g++ -O2 -std=c++17 Chapter20_Lesson4.cpp -o ch20_l4

Run:
  ./ch20_l4
*/

#include <cmath>
#include <cstdint>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>

struct LorenzParams {
  double sigma = 10.0;
  double rho   = 28.0;
  double beta  = 8.0 / 3.0;
};

struct State3 {
  double x, y, z;
};

static inline State3 add(const State3& a, const State3& b) {
  return {a.x + b.x, a.y + b.y, a.z + b.z};
}
static inline State3 mul(double s, const State3& a) {
  return {s * a.x, s * a.y, s * a.z};
}

static inline State3 lorenz_rhs(const State3& s, const LorenzParams& p) {
  return {
    p.sigma * (s.y - s.x),
    s.x * (p.rho - s.z) - s.y,
    s.x * s.y - p.beta * s.z
  };
}

static inline State3 rk4_step(const State3& s, double h, const LorenzParams& p) {
  State3 k1 = lorenz_rhs(s, p);
  State3 k2 = lorenz_rhs(add(s, mul(h * 0.5, k1)), p);
  State3 k3 = lorenz_rhs(add(s, mul(h * 0.5, k2)), p);
  State3 k4 = lorenz_rhs(add(s, mul(h, k3)), p);

  State3 incr = mul(h / 6.0, add(add(k1, mul(2.0, k2)), add(mul(2.0, k3), k4)));
  return add(s, incr);
}

static inline double logistic_next(double r, double x) {
  return r * x * (1.0 - x);
}

static double logistic_lyapunov(double r, double x0, int n_transient, int n) {
  double x = x0;
  for (int i = 0; i < n_transient; ++i) x = logistic_next(r, x);
  double sum = 0.0;
  for (int i = 0; i < n; ++i) {
    x = logistic_next(r, x);
    double d = std::abs(r * (1.0 - 2.0 * x));
    if (d < 1e-300) d = 1e-300;
    sum += std::log(d);
  }
  return sum / static_cast<double>(n);
}

int main() {
  const std::string outDir = "outputs_ch20_l4_cpp";

  // --- Lorenz simulation ---
  LorenzParams p;
  State3 s{1.0, 1.0, 1.0};
  double t = 0.0;
  const double T = 40.0;
  const double h = 0.001; // fixed step for simplicity

  std::ofstream traj(outDir + "_lorenz_traj.csv");
  traj << "t,x,y,z\n";

  std::ofstream poinc(outDir + "_lorenz_poincare_x0.csv");
  poinc << "x,y,z\n";

  double prev_x = s.x;
  for (int k = 0; t <= T; ++k) {
    traj << t << "," << s.x << "," << s.y << "," << s.z << "\n";

    // detect crossing x: negative -> positive
    if (prev_x < 0.0 && s.x >= 0.0) {
      poinc << s.x << "," << s.y << "," << s.z << "\n";
    }
    prev_x = s.x;

    s = rk4_step(s, h, p);
    t += h;
  }

  std::cout << "Wrote Lorenz CSV files: " << outDir << "_lorenz_*.csv\n";

  // --- Logistic map: bifurcation tail data ---
  std::ofstream bif(outDir + "_logistic_bifurcation.csv");
  bif << "r,x\n";
  const double r_min = 2.5;
  const double r_max = 4.0;
  const int n_r = 400;
  const int n_trans = 800;
  const int n_keep = 120;

  for (int i = 0; i < n_r; ++i) {
    double r = r_min + (r_max - r_min) * (static_cast<double>(i) / (n_r - 1));
    double x = 0.2;
    for (int k = 0; k < n_trans; ++k) x = logistic_next(r, x);
    for (int k = 0; k < n_keep; ++k) {
      x = logistic_next(r, x);
      bif << r << "," << x << "\n";
    }
  }

  // --- Logistic map: Lyapunov curve ---
  std::ofstream ly(outDir + "_logistic_lyapunov.csv");
  ly << "r,lambda\n";
  const int n_r2 = 200;
  for (int i = 0; i < n_r2; ++i) {
    double r = 2.8 + (4.0 - 2.8) * (static_cast<double>(i) / (n_r2 - 1));
    double lam = logistic_lyapunov(r, 0.2, 1000, 4000);
    ly << r << "," << lam << "\n";
  }

  std::cout << "Wrote logistic CSV files: " << outDir << "_logistic_*.csv\n";
  return 0;
}
      

Java is common in cross-platform engineering tools and large codebases. The numerical structure is identical to C++ here: implement RK4, run, export CSV.

File: Chapter20_Lesson4.java

// Chapter20_Lesson4.java
/*
Chapter 20 — Chaos, Complex Dynamics, and Computational Tools
Lesson 4 — Computational Tools

Java demonstration:
- Lorenz ODE simulated with fixed-step RK4
- Simple Poincaré section detection (x crosses 0 with positive direction)
- Logistic map bifurcation tail data and Lyapunov exponent for the map
- CSV export

Compile:
  javac Chapter20_Lesson4.java

Run:
  java Chapter20_Lesson4
*/

import java.io.PrintWriter;
import java.io.IOException;

public class Chapter20_Lesson4 {

  static class LorenzParams {
    double sigma = 10.0;
    double rho   = 28.0;
    double beta  = 8.0 / 3.0;
  }

  static double[] lorenzRhs(double[] s, LorenzParams p) {
    double x = s[0], y = s[1], z = s[2];
    return new double[] {
      p.sigma * (y - x),
      x * (p.rho - z) - y,
      x * y - p.beta * z
    };
  }

  static double[] add(double[] a, double[] b) {
    return new double[] { a[0] + b[0], a[1] + b[1], a[2] + b[2] };
  }

  static double[] mul(double s, double[] a) {
    return new double[] { s * a[0], s * a[1], s * a[2] };
  }

  static double[] rk4Step(double[] s, double h, LorenzParams p) {
    double[] k1 = lorenzRhs(s, p);
    double[] k2 = lorenzRhs(add(s, mul(h * 0.5, k1)), p);
    double[] k3 = lorenzRhs(add(s, mul(h * 0.5, k2)), p);
    double[] k4 = lorenzRhs(add(s, mul(h, k3)), p);

    double[] incr = mul(h / 6.0, add(add(k1, mul(2.0, k2)), add(mul(2.0, k3), k4)));
    return add(s, incr);
  }

  static double logisticNext(double r, double x) {
    return r * x * (1.0 - x);
  }

  static double logisticLyapunov(double r, double x0, int nTransient, int n) {
    double x = x0;
    for (int i = 0; i < nTransient; i++) x = logisticNext(r, x);
    double sum = 0.0;
    for (int i = 0; i < n; i++) {
      x = logisticNext(r, x);
      double d = Math.abs(r * (1.0 - 2.0 * x));
      if (d < 1e-300) d = 1e-300;
      sum += Math.log(d);
    }
    return sum / (double)n;
  }

  public static void main(String[] args) throws IOException {
    String outPrefix = "outputs_ch20_l4_java";

    // --- Lorenz simulation ---
    LorenzParams p = new LorenzParams();
    double[] s = new double[] {1.0, 1.0, 1.0};
    double t = 0.0;
    double T = 40.0;
    double h = 0.001;

    try (PrintWriter traj = new PrintWriter(outPrefix + "_lorenz_traj.csv");
         PrintWriter poinc = new PrintWriter(outPrefix + "_lorenz_poincare_x0.csv")) {

      traj.println("t,x,y,z");
      poinc.println("x,y,z");

      double prevX = s[0];

      while (t <= T) {
        traj.printf("%.10f,%.10f,%.10f,%.10f%n", t, s[0], s[1], s[2]);

        if (prevX < 0.0 && s[0] >= 0.0) {
          poinc.printf("%.10f,%.10f,%.10f%n", s[0], s[1], s[2]);
        }
        prevX = s[0];

        s = rk4Step(s, h, p);
        t += h;
      }
    }

    // --- Logistic map: bifurcation tail data ---
    try (PrintWriter bif = new PrintWriter(outPrefix + "_logistic_bifurcation.csv")) {
      bif.println("r,x");
      double rMin = 2.5, rMax = 4.0;
      int nR = 400, nTrans = 800, nKeep = 120;

      for (int i = 0; i < nR; i++) {
        double r = rMin + (rMax - rMin) * ((double)i / (double)(nR - 1));
        double x = 0.2;
        for (int k = 0; k < nTrans; k++) x = logisticNext(r, x);
        for (int k = 0; k < nKeep; k++) {
          x = logisticNext(r, x);
          bif.printf("%.10f,%.10f%n", r, x);
        }
      }
    }

    // --- Logistic map: Lyapunov curve ---
    try (PrintWriter ly = new PrintWriter(outPrefix + "_logistic_lyapunov.csv")) {
      ly.println("r,lambda");
      int nR2 = 200;
      for (int i = 0; i < nR2; i++) {
        double r = 2.8 + (4.0 - 2.8) * ((double)i / (double)(nR2 - 1));
        double lam = logisticLyapunov(r, 0.2, 1000, 4000);
        ly.printf("%.10f,%.10f%n", r, lam);
      }
    }

    System.out.println("Done. Wrote CSV files with prefix: " + outPrefix + "_*.csv");
  }
}
      

Modelica is an equation-based, acausal modeling language designed for physical systems (multi-domain). Unlike MATLAB/Python where you typically write \( \dot{\mathbf{x}}=\mathbf{f}(t,\mathbf{x}) \), in Modelica you write equations and the tool derives a solvable form (often a DAE) automatically.

This changes the computational problem: before integrating, the compiler performs structural analysis, possibly index reduction and tearing, then generates code for a numeric DAE/ODE solver.

flowchart TD
  M["Write equations (acausal)"] --> S["Structural analysis (DAE form)"]
  S --> IR["Index reduction + consistent init"]
  IR --> T["Tearing / symbolic simplification"]
  T --> CG["Code generation (C/C++)"]
  CG --> NS["Numeric solver (ODE/DAE, events)"]
  NS --> OUT["Results (trajectories, events, logs)"]
        

Even though the Lorenz system is an ODE, writing it in Modelica illustrates the equation-based style.

File: Chapter20_Lesson4.mo

// Chapter20_Lesson4.mo
// Chapter 20 — Chaos, Complex Dynamics, and Computational Tools
// Lesson 4 — Modelica model (Lorenz system)
//
// This is a minimal Modelica model. Simulate it in OpenModelica or Dymola.
//
// Suggested experiment: StopTime=40, Interval=0.01

model Chapter20_Lesson4_Lorenz
  parameter Real sigma = 10.0;
  parameter Real rho   = 28.0;
  parameter Real beta  = 8.0/3.0;

  Real x(start=1.0);
  Real y(start=1.0);
  Real z(start=1.0);

equation
  der(x) = sigma*(y - x);
  der(y) = x*(rho - z) - y;
  der(z) = x*y - beta*z;

annotation(
  experiment(StopTime=40, Interval=0.01)
);
end Chapter20_Lesson4_Lorenz;
      

DAE note: In more realistic multi-domain models, constraints (e.g., kinematic loops, ideal electrical elements) yield DAEs. Tool robustness depends on consistent initialization and solver selection (DAE solvers such as DASSL-family).

Mathematica is especially useful for system dynamics because it combines symbolic manipulation (Jacobians, stability analysis, simplification) with robust numerical solvers (NDSolve) and a strong data/export pipeline. Below is a notebook-style file that solves Lorenz, computes a simple Poincaré section by sign-change detection, and generates logistic bifurcation/lyapunov CSVs.

File: Chapter20_Lesson4.nb

(* Chapter20_Lesson4.nb *)

Notebook[{
  Cell["Chapter 20 — Chaos, Complex Dynamics, and Computational Tools", "Title"],
  Cell["Lesson 4 — Computational Tools (Wolfram Language)", "Subtitle"],

  Cell["1) Lorenz system via NDSolve and CSV export", "Section"],
  Cell["lorenz = {x'[t] == 10 (y[t] - x[t]), y'[t] == x[t] (28 - z[t]) - y[t], z'[t] == x[t] y[t] - (8/3) z[t], x[0] == 1, y[0] == 1, z[0] == 1};\nsol = NDSolveValue[lorenz, {x, y, z}, {t, 0, 40}, MaxStepFraction -> 1/2000];\n(* Sample the solution *)\nts = Range[0, 40, 0.01];\ntraj = Transpose[{ts, sol[[1]] /@ ts, sol[[2]] /@ ts, sol[[3]] /@ ts}];\nExport[\"lorenz_traj.csv\", traj];", "Input"],

  Cell["2) Poincaré section: x(t)=0 crossings (simple sign-change detection)", "Section"],
  Cell["xs = sol[[1]] /@ ts; ys = sol[[2]] /@ ts; zs = sol[[3]] /@ ts;\nidx = Flatten@Position[Partition[Sign[xs], 2, 1], {-1, 1}];\npoinc = Transpose[{xs[[idx + 1]], ys[[idx + 1]], zs[[idx + 1]]}];\nExport[\"lorenz_poincare_x0.csv\", poinc];", "Input"],

  Cell["3) Logistic map: bifurcation tail data and Lyapunov exponent curve", "Section"],
  Cell["logistic[r_?NumericQ, x_?NumericQ] := r x (1 - x);\nlogisticIter[r_?NumericQ, x0_?NumericQ, n_Integer] := NestList[logistic[r, #] &, x0, n];\n\nrVals = N@Subdivide[2.5, 4.0, 399];\ntrans = 800; keep = 120;\nbif = Flatten[Table[\n  xs = logisticIter[r, 0.2, trans + keep];\n  Transpose[{ConstantArray[r, keep], Take[xs, -keep]}]\n, {r, rVals}], 1];\nExport[\"logistic_bifurcation.csv\", bif];\n\nlyap[r_?NumericQ] := Module[{xs, d},\n  xs = Take[logisticIter[r, 0.2, 1000 + 4000], -4000];\n  d = Abs[r (1 - 2 xs)];\n  Mean[Log[Max[d, 10^-300]]]\n];\n\nrVals2 = N@Subdivide[2.8, 4.0, 199];\nly = Table[{r, lyap[r]}, {r, rVals2}];\nExport[\"logistic_lyapunov.csv\", ly];", "Input"]
},
WindowTitle -> "Chapter20_Lesson4",
WindowSize -> {900, 700}
]
      

Problem 1 (Embedded RK step-size law): Assume an embedded pair produces an error estimate with scaling \( \|\mathbf{e}_n\| \approx K h^{p+1} \). Derive a step update that targets a tolerance \( \text{tol} \) for the next step.

Solution: Set \( \text{tol} \approx K h_{\text{new}}^{p+1} \). Using the current estimate \( \|\mathbf{e}_n\| \approx K h^{p+1} \), eliminate \( K \):

\[ \frac{\text{tol}}{\|\mathbf{e}_n\|} \approx \left(\frac{h_{\text{new}}}{h}\right)^{p+1} \;\Rightarrow\; h_{\text{new}} \approx h\left(\frac{\text{tol}}{\|\mathbf{e}_n\|}\right)^{1/(p+1)}. \]

Practical solvers use the exponent \( 1/p \) or \( 1/(p+1) \) depending on conventions (which solution is accepted), plus a safety factor and bounds.

Problem 2 (Absolute stability of RK4 on \(y'=\lambda y\)): For the classical RK4 method, show that the stability function is \( R(z)=1+z+z^2/2+z^3/6+z^4/24 \). Conclude why explicit RK4 can require tiny step sizes for large negative \( \lambda \).

Solution: On the test equation, each stage is a polynomial in \( z \), and RK4 matches the 4th-order Taylor polynomial of \( e^z \):

\[ y_{n+1} = R(z) y_n,\quad R(z)=\sum_{k=0}^4 \frac{z^k}{k!}. \]

Stability requires \( |R(z)|<1 \). For stiff modes with large negative \( z=\lambda h \), this fails unless \( h \) is very small, motivating stiff (often implicit) solvers.

Problem 3 (Event accuracy and step refinement): Consider an event \( g(t,\mathbf{x}(t))=0 \) located by root-finding on an interpolant. Explain why reducing the maximum step size improves the accuracy of event times, and propose a verification test.

Solution: Interpolants are constructed from step data, so coarser stepping produces larger interpolation error and a noisier root. A verification test is a step-halving study: run with \( h_{\max} \) and \( h_{\max}/2 \) (same tolerances), compare event times \( t_e \) and section points \( \mathbf{x}(t_e) \). Convergence under refinement suggests numerical robustness.

Problem 4 (Logistic map Lyapunov exponent): For \( x_{k+1}=r x_k(1-x_k) \), show that the largest Lyapunov exponent is \( \lambda = \lim_{N\to\infty}\frac{1}{N}\sum_{k=0}^{N-1}\ln|r(1-2x_k)| \).

Solution: Linearize the map: \( \delta x_{k+1}=f'(x_k)\delta x_k \) with \( f'(x)=r(1-2x) \). Then \( |\delta x_N| = |\delta x_0|\prod_{k=0}^{N-1}|f'(x_k)| \). Taking logs and dividing by \( N \) gives the stated formula.

Problem 5 (DAE initialization concept): For \( \mathbf{F}(t,\mathbf{x},\dot{\mathbf{x}})=\mathbf{0} \), explain what it means for initial conditions to be “consistent” and why equation-based tools care.

Solution: Consistency means that at \( t_0 \) there exist values \( \mathbf{x}(t_0) \) and \( \dot{\mathbf{x}}(t_0) \) satisfying \( \mathbf{F}(t_0,\mathbf{x}(t_0),\dot{\mathbf{x}}(t_0))=\mathbf{0} \). If constraints are violated at initialization, the solver starts off the constraint manifold and may fail immediately or produce nonphysical transients. Modelica tools often solve a nonlinear system at start-up to enforce consistency.

We formalized how system dynamics models appear as ODEs, DAEs, and discrete-time maps; connected solver choices to stability and stiffness; and showed how modern tools support event handling and reproducible analysis (Poincaré sections, bifurcation sweeps, Lyapunov estimation). The provided MATLAB/Simulink, Python/SciPy, Modelica, C++, Java, and Mathematica implementations are consistent, export common CSV artifacts, and illustrate a tool-agnostic scientific workflow.