In previous lessons of this chapter, you saw that nonlinear systems can generate trajectories that appear irregular and unpredictable despite being deterministic. The precise mathematical signature of chaos is that two trajectories starting extremely close can diverge rapidly. Let \( \mathbf{x}(t) \in \mathbb{R}^n \) be a system state and \( \delta\mathbf{x}(t) \) a small perturbation. If typical perturbations grow approximately like \( \|\delta\mathbf{x}(t)\| \approx \|\delta\mathbf{x}(0)\| e^{\lambda t} \), then \( \lambda \) measures predictability loss. The quantity \( \lambda \) is a (largest) Lyapunov exponent.

flowchart TD
  A["Choose system model (map or ODE)"] --> B["Pick two nearby initial states"]
  B --> C["Evolve both states forward in time"]
  C --> D["Measure separation d(t)=||x2(t)-x1(t)||"]
  D --> E["Fit growth rate: d(t) ~ d0 * exp(lambda * t)"]
  E --> F["Interpret: lambda > 0 implies sensitivity and finite predictability"]
        

Important: Lyapunov exponents are asymptotic quantities. Transient behaviors can mimic growth/decay. Proper estimation requires discarding transients and using sufficiently long time horizons.

Consider a discrete-time nonlinear map \( \mathbf{x}_{k+1} = \mathbf{f}(\mathbf{x}_k) \) or a continuous-time ODE \( \dot{\mathbf{x}} = \mathbf{F}(\mathbf{x}) \). Sensitivity to initial conditions means that arbitrarily close initial conditions may separate to a macroscopic distance after some time.

A common formal definition (topological form) says: a dynamical system on a metric space has sensitivity if there exists \( \varepsilon > 0 \) such that for any point and any neighborhood, some point in that neighborhood eventually separates by at least \( \varepsilon \).

In engineering and computation, we quantify sensitivity through rates of separation. Define the separation \( d(t) \): as

\[ d(t) \equiv \|\mathbf{x}(t;\mathbf{x}_0+\delta\mathbf{x}_0) - \mathbf{x}(t;\mathbf{x}_0)\|. \]

For many systems, for sufficiently small \( \|\delta\mathbf{x}_0\| \) and for times before saturation, one observes approximate exponential behavior \( d(t) \):

\[ d(t) \approx d(0)\, e^{\lambda t} \quad \text{(ODEs)}, \qquad d_k \approx d_0\, e^{\lambda k} \quad \text{(maps)}. \]

If \( \lambda > 0 \), errors double roughly every \( (\ln 2)/\lambda \) units (time for ODEs; iterations for maps), producing a short predictability horizon.

Let \( \mathbf{x}_{k+1}=\mathbf{f}(\mathbf{x}_k) \) with \( \mathbf{f}:\mathbb{R}^n \to \mathbb{R}^n \) differentiable. Linearizing around a trajectory gives the tangent dynamics \( \delta\mathbf{x}_k \):

\[ \delta\mathbf{x}_{k+1} = \mathbf{J}(\mathbf{x}_k)\,\delta\mathbf{x}_k, \qquad \mathbf{J}(\mathbf{x}_k) = \frac{\partial \mathbf{f}}{\partial \mathbf{x}}(\mathbf{x}_k). \]

Iterating, \( \delta\mathbf{x}_k \):

\[ \delta\mathbf{x}_{k} = \left(\prod_{i=0}^{k-1}\mathbf{J}(\mathbf{x}_i)\right)\delta\mathbf{x}_0. \]

The largest Lyapunov exponent can be defined (when the limit exists for typical initial conditions) as \( \lambda_{\max} \):

\[ \lambda_{\max} = \lim_{k\to\infty}\frac{1}{k}\ln\frac{\|\delta\mathbf{x}_k\|}{\|\delta\mathbf{x}_0\|}. \]

1D special case (important and computationally cheap). For \( x_{k+1}=f(x_k) \), the tangent dynamics is \( \delta x_{k+1}=f'(x_k)\delta x_k \), so

\[ \frac{\delta x_k}{\delta x_0} = \prod_{i=0}^{k-1} f'(x_i) \quad \Rightarrow \quad \lambda = \lim_{k\to\infty}\frac{1}{k}\sum_{i=0}^{k-1}\ln|f'(x_i)|. \]

This turns Lyapunov estimation into a time-average along the orbit. You already met an example in Lesson 1: the logistic map \( x_{k+1}=r x_k(1-x_k) \) has \( f'(x) = r(1-2x) \):.

Derivation (short proof for 1D formula). Starting from \( \delta x_{k+1}=f'(x_k)\delta x_k \), recursively:

\[ \delta x_k = \left(\prod_{i=0}^{k-1} f'(x_i)\right)\delta x_0 \Rightarrow \ln\left|\frac{\delta x_k}{\delta x_0}\right| = \sum_{i=0}^{k-1}\ln|f'(x_i)|. \]

Divide by \( k \) and take \( k\to\infty \) (when the average converges) to obtain the formula.

Consider an autonomous ODE \( \dot{\mathbf{x}} = \mathbf{F}(\mathbf{x}) \) with a trajectory \( \mathbf{x}(t) \). Perturb the initial condition by \( \delta\mathbf{x}(0) \). For small perturbations, the perturbation evolves according to the variational equation: \( \delta\mathbf{x}(t) \):

\[ \frac{d}{dt}\delta\mathbf{x}(t) = \mathbf{A}(t)\,\delta\mathbf{x}(t), \qquad \mathbf{A}(t) = \frac{\partial \mathbf{F}}{\partial \mathbf{x}}\big(\mathbf{x}(t)\big). \]

The state transition matrix \( \mathbf{\Phi}(t) \) solves \( \dot{\mathbf{\Phi}}(t)=\mathbf{A}(t)\mathbf{\Phi}(t),\ \mathbf{\Phi}(0)=\mathbf{I} \), giving \( \delta\mathbf{x}(t)=\mathbf{\Phi}(t)\delta\mathbf{x}(0) \). The largest exponent is defined as \( \lambda_{\max} \):

\[ \lambda_{\max} = \lim_{t\to\infty}\frac{1}{t}\ln\frac{\|\delta\mathbf{x}(t)\|}{\|\delta\mathbf{x}(0)\|}. \]

Key linear benchmark (proof). If \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \) with constant matrix \( \mathbf{A} \), then \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}(0) \). For generic initial conditions, the dominant growth rate is governed by the eigenvalue of \( \mathbf{A} \) with the largest real part:

\[ \lambda_{\max} = \max_i \operatorname{Re}(\mu_i), \quad \text{where } \mu_i \text{ are eigenvalues of } \mathbf{A}. \]

Sketch. Put \( \mathbf{A}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1} \) (diagonalizable case), then \( e^{\mathbf{A}t}=\mathbf{V}e^{\mathbf{\Lambda}t}\mathbf{V}^{-1} \). Components along eigenvectors grow like \( e^{\operatorname{Re}(\mu_i)t} \); the largest \( \operatorname{Re}(\mu_i) \) dominates asymptotically.

For maps and ODEs, the clean definition uses tangent dynamics. Numerically, a very common introductory estimator (robust and easy to implement) is the two-trajectory renormalization method:

flowchart TD
  S["Start with reference x and perturbed x+delta"] --> I["Integrate both for m steps"]
  I --> M["Measure d = ||(x+delta) - x||"]
  M --> A["Accumulate sum += log(d/d0)"]
  A --> R["Renormalize: delta = d0 * ((x+delta)-x)/d"]
  R --> L["Repeat many times; lambda ~= sum / total_time"]
        

The method enforces \( \|\delta\| \approx d_0 \) to stay in the linearized regime while accumulating the average exponential stretching. If the estimator converges to \( \lambda_{\max} > 0 \), the system exhibits strong sensitivity.

Numerical cautions (engineering view):

• Discard transients: compute exponents on the attractor, not during approach.
• Choose \( d_0 \) small (linear regime) but not so small that finite precision causes collapse.
• Use stable integrators and verify step-size convergence (especially for ODEs).
• Expect saturation: the separation cannot grow forever because trajectories remain in a bounded region; this is why renormalization is used.

Predictability horizon. If an initial uncertainty is \( d(0) \) and you can tolerate error up to \( d_{\text{tol}} \), then a crude horizon estimate is \( T_p \):

\[ T_p \approx \frac{1}{\lambda_{\max}}\ln\left(\frac{d_{\text{tol}}}{d(0)}\right), \qquad (\lambda_{\max} > 0). \]

We implement: (i) the 1D time-average formula for the logistic map, and (ii) the two-trajectory renormalization estimator for the Lorenz ODE. These provide an introductory but practical Lyapunov workflow.

6.1 Python (NumPy) Implementation

File: Chapter20_Lesson3.py

# Chapter20_Lesson3.py
# System Dynamics (Control Engineering) - Chapter 20, Lesson 3
# Sensitivity to Initial Conditions and Lyapunov Exponents (Intro)

import math
import numpy as np

# -----------------------------
# Part A) Discrete-time map: Logistic map
# x_{n+1} = r x_n (1 - x_n)
# Largest Lyapunov exponent (LLE) approximation:
#   lambda ≈ (1/N) * sum_{n=0}^{N-1} log | f'(x_n) |
# with f'(x) = r (1 - 2x)
# -----------------------------

def logistic_map_step(x: float, r: float) -> float:
    return r * x * (1.0 - x)

def lyapunov_logistic(r: float, x0: float = 0.2, n: int = 200_000, discard: int = 5_000) -> float:
    x = float(x0)
    # transient removal
    for _ in range(discard):
        x = logistic_map_step(x, r)

    s = 0.0
    for _ in range(n):
        # derivative along the orbit
        fp = r * (1.0 - 2.0 * x)
        s += math.log(abs(fp) + 1e-300)  # tiny offset avoids log(0) in degenerate cases
        x = logistic_map_step(x, r)
    return s / n

# -----------------------------
# Part B) Continuous-time system: Lorenz (classic chaotic ODE)
#   xdot = sigma (y - x)
#   ydot = x (rho - z) - y
#   zdot = x y - beta z
#
# Largest Lyapunov exponent approximation by "two-trajectory renormalization":
# 1) Integrate a reference state X(t)
# 2) Integrate a perturbed state Xp(t) = X(t) + d0 * u (small separation)
# 3) Every m steps: measure separation d = ||Xp - X||,
#    accumulate log(d/d0), then renormalize: Xp = X + d0 * (Xp - X)/d
# 4) lambda ≈ (1/T) * sum log(d/d0)
# -----------------------------

def lorenz_rhs(t: float, X: np.ndarray, sigma: float = 10.0, rho: float = 28.0, beta: float = 8.0/3.0) -> np.ndarray:
    x, y, z = X
    return np.array([
        sigma * (y - x),
        x * (rho - z) - y,
        x * y - beta * z
    ], dtype=float)

def rk4_step(rhs, t: float, y: np.ndarray, dt: float, *rhs_args, **rhs_kwargs) -> np.ndarray:
    k1 = rhs(t, y, *rhs_args, **rhs_kwargs)
    k2 = rhs(t + 0.5*dt, y + 0.5*dt*k1, *rhs_args, **rhs_kwargs)
    k3 = rhs(t + 0.5*dt, y + 0.5*dt*k2, *rhs_args, **rhs_kwargs)
    k4 = rhs(t + dt, y + dt*k3, *rhs_args, **rhs_kwargs)
    return y + (dt/6.0) * (k1 + 2*k2 + 2*k3 + k4)

def lyapunov_lorenz_lle(
    X0 = (1.0, 1.0, 1.0),
    dt: float = 0.01,
    T: float = 100.0,
    transient: float = 10.0,
    renorm_every: int = 10,
    d0: float = 1e-8,
    sigma: float = 10.0,
    rho: float = 28.0,
    beta: float = 8.0/3.0,
    seed: int = 0
) -> float:
    rng = np.random.default_rng(seed)
    X = np.array(X0, dtype=float)

    # random unit direction for initial perturbation
    u = rng.normal(size=3)
    u = u / np.linalg.norm(u)
    Xp = X + d0 * u

    # integrate out transient without measuring
    t = 0.0
    n_trans = int(transient / dt)
    for _ in range(n_trans):
        X = rk4_step(lorenz_rhs, t, X, dt, sigma=sigma, rho=rho, beta=beta)
        Xp = rk4_step(lorenz_rhs, t, Xp, dt, sigma=sigma, rho=rho, beta=beta)
        t += dt
        # keep perturbation small to avoid collapse or blow-up
        d = np.linalg.norm(Xp - X)
        if d == 0.0:
            Xp = X + d0 * u
        else:
            Xp = X + d0 * (Xp - X) / d

    # main accumulation
    s = 0.0
    steps = int(T / dt)
    count_renorm = 0
    for k in range(steps):
        X = rk4_step(lorenz_rhs, t, X, dt, sigma=sigma, rho=rho, beta=beta)
        Xp = rk4_step(lorenz_rhs, t, Xp, dt, sigma=sigma, rho=rho, beta=beta)
        t += dt

        if (k + 1) % renorm_every == 0:
            d = np.linalg.norm(Xp - X)
            if d == 0.0:
                # extremely rare; re-seed perturbation
                u = rng.normal(size=3)
                u = u / np.linalg.norm(u)
                Xp = X + d0 * u
                continue
            s += math.log(d / d0)
            count_renorm += 1
            Xp = X + d0 * (Xp - X) / d

    total_time = count_renorm * renorm_every * dt
    return s / total_time

def demo():
    print("Logistic map LLE examples:")
    for r in [3.2, 3.5, 3.9, 4.0]:
        lam = lyapunov_logistic(r, x0=0.234, n=100_000, discard=5_000)
        print(f"  r={r:>3}: lambda ≈ {lam:.6f}")

    print("\nLorenz LLE example (classic sigma=10, rho=28, beta=8/3):")
    lam_lor = lyapunov_lorenz_lle(T=120.0, transient=20.0, dt=0.01, renorm_every=10, d0=1e-8)
    print(f"  lambda_max ≈ {lam_lor:.4f} 1/time-unit")

if __name__ == "__main__":
    demo()

6.2 C++ Implementation (RK4 + Renormalization)

File: Chapter20_Lesson3.cpp

// Chapter20_Lesson3.cpp
// System Dynamics (Control Engineering) - Chapter 20, Lesson 3
// Sensitivity to Initial Conditions and Lyapunov Exponents (Intro)

#include <iostream>
#include <iomanip>
#include <cmath>
#include <array>
#include <random>

// -----------------------------
// Part A) Logistic map LLE
// x_{n+1} = r x_n (1 - x_n)
// lambda ≈ (1/N) sum log |r (1 - 2 x_n)|
// -----------------------------
static inline double logistic_step(double x, double r) {
    return r * x * (1.0 - x);
}

double lyapunov_logistic(double r, double x0=0.2, int n=200000, int discard=5000) {
    double x = x0;
    for (int i=0; i<discard; ++i) x = logistic_step(x, r);

    long double s = 0.0L;
    for (int i=0; i<n; ++i) {
        double fp = r * (1.0 - 2.0*x);
        s += std::log(std::fabs(fp) + 1e-300);
        x = logistic_step(x, r);
    }
    return (double)(s / (long double)n);
}

// -----------------------------
// Part B) Lorenz LLE (two-trajectory renormalization) with RK4
// -----------------------------
using Vec3 = std::array<double,3>;

static inline Vec3 lorenz_rhs(const Vec3& X, double sigma=10.0, double rho=28.0, double beta=8.0/3.0) {
    const double x = X[0], y = X[1], z = X[2];
    return Vec3{
        sigma * (y - x),
        x * (rho - z) - y,
        x * y - beta * z
    };
}

static inline Vec3 add(const Vec3& a, const Vec3& b) {
    return Vec3{a[0]+b[0], a[1]+b[1], a[2]+b[2]};
}
static inline Vec3 scale(const Vec3& a, double s) {
    return Vec3{a[0]*s, a[1]*s, a[2]*s};
}
static inline Vec3 sub(const Vec3& a, const Vec3& b) {
    return Vec3{a[0]-b[0], a[1]-b[1], a[2]-b[2]};
}
static inline double norm2(const Vec3& a) {
    return std::sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
}

Vec3 rk4_step_lorenz(const Vec3& X, double dt, double sigma=10.0, double rho=28.0, double beta=8.0/3.0) {
    const Vec3 k1 = lorenz_rhs(X, sigma, rho, beta);
    const Vec3 k2 = lorenz_rhs(add(X, scale(k1, 0.5*dt)), sigma, rho, beta);
    const Vec3 k3 = lorenz_rhs(add(X, scale(k2, 0.5*dt)), sigma, rho, beta);
    const Vec3 k4 = lorenz_rhs(add(X, scale(k3, dt)), sigma, rho, beta);
    return add(X, scale(add(add(k1, scale(k2,2.0)), add(scale(k3,2.0), k4)), dt/6.0));
}

double lyapunov_lorenz_lle(
    Vec3 X0 = {1.0,1.0,1.0},
    double dt = 0.01,
    double T = 100.0,
    double transient = 10.0,
    int renorm_every = 10,
    double d0 = 1e-8,
    double sigma = 10.0,
    double rho = 28.0,
    double beta = 8.0/3.0,
    unsigned seed = 0
) {
    std::mt19937 gen(seed);
    std::normal_distribution<double> N01(0.0, 1.0);

    Vec3 X = X0;
    Vec3 u{N01(gen), N01(gen), N01(gen)};
    double un = norm2(u);
    u = scale(u, 1.0/un);
    Vec3 Xp = add(X, scale(u, d0));

    // transient
    int n_trans = (int)std::round(transient/dt);
    for (int i=0; i<n_trans; ++i) {
        X = rk4_step_lorenz(X, dt, sigma, rho, beta);
        Xp = rk4_step_lorenz(Xp, dt, sigma, rho, beta);
        Vec3 dvec = sub(Xp, X);
        double d = norm2(dvec);
        if (d == 0.0) {
            Xp = add(X, scale(u, d0));
        } else {
            Xp = add(X, scale(dvec, d0/d));
        }
    }

    // accumulate
    long double s = 0.0L;
    int steps = (int)std::round(T/dt);
    int count = 0;
    for (int k=0; k<steps; ++k) {
        X  = rk4_step_lorenz(X, dt, sigma, rho, beta);
        Xp = rk4_step_lorenz(Xp, dt, sigma, rho, beta);

        if ((k+1) % renorm_every == 0) {
            Vec3 dvec = sub(Xp, X);
            double d = norm2(dvec);
            if (d == 0.0) continue;
            s += std::log(d/d0);
            count += 1;
            Xp = add(X, scale(dvec, d0/d));
        }
    }

    double total_time = (double)count * (double)renorm_every * dt;
    return (double)(s / (long double)total_time);
}

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

    std::cout << "Logistic map LLE examples:\n";
    for (double r : {3.2, 3.5, 3.9, 4.0}) {
        double lam = lyapunov_logistic(r, 0.234, 100000, 5000);
        std::cout << "  r=" << r << " : lambda ~= " << lam << "\n";
    }

    std::cout << "\nLorenz LLE example (sigma=10, rho=28, beta=8/3):\n";
    double lamL = lyapunov_lorenz_lle({1.0,1.0,1.0}, 0.01, 120.0, 20.0, 10, 1e-8, 10.0, 28.0, 8.0/3.0, 0);
    std::cout << "  lambda_max ~= " << lamL << " 1/time-unit\n";

    return 0;
}

6.3 Java Implementation

File: Chapter20_Lesson3.java

// Chapter20_Lesson3.java
// System Dynamics (Control Engineering) - Chapter 20, Lesson 3
// Sensitivity to Initial Conditions and Lyapunov Exponents (Intro)

import java.util.Random;

public class Chapter20_Lesson3 {

    // -----------------------------
    // Part A) Logistic map LLE
    // -----------------------------
    static double logisticStep(double x, double r) {
        return r * x * (1.0 - x);
    }

    static double lyapunovLogistic(double r, double x0, int n, int discard) {
        double x = x0;
        for (int i = 0; i < discard; i++) x = logisticStep(x, r);

        double s = 0.0;
        for (int i = 0; i < n; i++) {
            double fp = r * (1.0 - 2.0 * x);
            s += Math.log(Math.abs(fp) + 1e-300);
            x = logisticStep(x, r);
        }
        return s / n;
    }

    // -----------------------------
    // Part B) Lorenz LLE via two-trajectory renormalization with RK4
    // -----------------------------
    static double[] lorenzRhs(double[] X, double sigma, double rho, double beta) {
        double x = X[0], y = X[1], z = X[2];
        return new double[] {
            sigma * (y - x),
            x * (rho - z) - y,
            x * y - beta * z
        };
    }

    static double[] rk4StepLorenz(double[] X, double dt, double sigma, double rho, double beta) {
        double[] k1 = lorenzRhs(X, sigma, rho, beta);

        double[] X2 = new double[] { X[0] + 0.5*dt*k1[0], X[1] + 0.5*dt*k1[1], X[2] + 0.5*dt*k1[2] };
        double[] k2 = lorenzRhs(X2, sigma, rho, beta);

        double[] X3 = new double[] { X[0] + 0.5*dt*k2[0], X[1] + 0.5*dt*k2[1], X[2] + 0.5*dt*k2[2] };
        double[] k3 = lorenzRhs(X3, sigma, rho, beta);

        double[] X4 = new double[] { X[0] + dt*k3[0], X[1] + dt*k3[1], X[2] + dt*k3[2] };
        double[] k4 = lorenzRhs(X4, sigma, rho, beta);

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

    static double norm(double[] v) {
        return Math.sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
    }

    static double lyapunovLorenzLLE(
        double[] X0,
        double dt,
        double T,
        double transient,
        int renormEvery,
        double d0,
        double sigma,
        double rho,
        double beta,
        long seed
    ) {
        Random rng = new Random(seed);
        double[] X = new double[] { X0[0], X0[1], X0[2] };

        // random unit direction
        double[] u = new double[] { rng.nextGaussian(), rng.nextGaussian(), rng.nextGaussian() };
        double un = norm(u);
        u[0] /= un; u[1] /= un; u[2] /= un;

        double[] Xp = new double[] { X[0] + d0*u[0], X[1] + d0*u[1], X[2] + d0*u[2] };

        // transient integration (keep separation ~ d0 but don't accumulate)
        int nTrans = (int)Math.round(transient / dt);
        for (int i = 0; i < nTrans; i++) {
            X  = rk4StepLorenz(X,  dt, sigma, rho, beta);
            Xp = rk4StepLorenz(Xp, dt, sigma, rho, beta);

            double[] dvec = new double[] { Xp[0]-X[0], Xp[1]-X[1], Xp[2]-X[2] };
            double d = norm(dvec);
            if (d == 0.0) {
                Xp[0] = X[0] + d0*u[0]; Xp[1] = X[1] + d0*u[1]; Xp[2] = X[2] + d0*u[2];
            } else {
                double s = d0 / d;
                Xp[0] = X[0] + s*dvec[0]; Xp[1] = X[1] + s*dvec[1]; Xp[2] = X[2] + s*dvec[2];
            }
        }

        // main accumulation
        double sum = 0.0;
        int steps = (int)Math.round(T / dt);
        int count = 0;
        for (int k = 0; k < steps; k++) {
            X  = rk4StepLorenz(X,  dt, sigma, rho, beta);
            Xp = rk4StepLorenz(Xp, dt, sigma, rho, beta);

            if ((k+1) % renormEvery == 0) {
                double[] dvec = new double[] { Xp[0]-X[0], Xp[1]-X[1], Xp[2]-X[2] };
                double d = norm(dvec);
                if (d == 0.0) continue;
                sum += Math.log(d / d0);
                count += 1;

                double s = d0 / d;
                Xp[0] = X[0] + s*dvec[0]; Xp[1] = X[1] + s*dvec[1]; Xp[2] = X[2] + s*dvec[2];
            }
        }

        double totalTime = count * renormEvery * dt;
        return sum / totalTime;
    }

    public static void main(String[] args) {
        System.out.println("Logistic map LLE examples:");
        double[] rList = new double[] {3.2, 3.5, 3.9, 4.0};
        for (double r : rList) {
            double lam = lyapunovLogistic(r, 0.234, 100000, 5000);
            System.out.printf("  r=%.1f : lambda ~= %.6f%n", r, lam);
        }

        System.out.println("\nLorenz LLE example (sigma=10, rho=28, beta=8/3):");
        double lamL = lyapunovLorenzLLE(new double[]{1,1,1}, 0.01, 120.0, 20.0, 10, 1e-8, 10.0, 28.0, 8.0/3.0, 0);
        System.out.printf("  lambda_max ~= %.4f 1/time-unit%n", lamL);
    }
}

6.4 MATLAB (and Simulink-style scripting) Implementation

File: Chapter20_Lesson3.m

% Chapter20_Lesson3.m
% System Dynamics (Control Engineering) - Chapter 20, Lesson 3
% Sensitivity to Initial Conditions and Lyapunov Exponents (Intro)

clear; clc;

%% Part A) Logistic map LLE
% x_{n+1} = r x_n (1 - x_n)
% lambda ≈ (1/N) sum log | r (1 - 2 x_n) |

rList = [3.2, 3.5, 3.9, 4.0];
for r = rList
    lam = lyapunov_logistic(r, 0.234, 100000, 5000);
    fprintf('Logistic: r=%.1f, lambda ~= %.6f\n', r, lam);
end

%% Part B) Lorenz LLE (two-trajectory renormalization with RK4)
sigma = 10; rho = 28; beta = 8/3;
dt = 0.01;
T = 120;
transient = 20;
renormEvery = 10;
d0 = 1e-8;

lamL = lyapunov_lorenz_lle([1;1;1], dt, T, transient, renormEvery, d0, sigma, rho, beta, 0);
fprintf('Lorenz: lambda_max ~= %.4f 1/time-unit\n', lamL);

%% --- Functions ---
function xnext = logistic_step(x, r)
    xnext = r*x*(1-x);
end

function lam = lyapunov_logistic(r, x0, N, discard)
    x = x0;
    for i = 1:discard
        x = logistic_step(x, r);
    end

    s = 0.0;
    for i = 1:N
        fp = r*(1 - 2*x);
        s = s + log(abs(fp) + 1e-300);
        x = logistic_step(x, r);
    end
    lam = s / N;
end

function dX = lorenz_rhs(X, sigma, rho, beta)
    x = X(1); y = X(2); z = X(3);
    dX = [ sigma*(y-x);
           x*(rho - z) - y;
           x*y - beta*z ];
end

function Xn = rk4_step(rhs, X, dt, sigma, rho, beta)
    k1 = rhs(X, sigma, rho, beta);
    k2 = rhs(X + 0.5*dt*k1, sigma, rho, beta);
    k3 = rhs(X + 0.5*dt*k2, sigma, rho, beta);
    k4 = rhs(X + dt*k3, sigma, rho, beta);
    Xn = X + (dt/6)*(k1 + 2*k2 + 2*k3 + k4);
end

function lam = lyapunov_lorenz_lle(X0, dt, T, transient, renormEvery, d0, sigma, rho, beta, seed)
    rng(seed);

    X = X0(:);
    u = randn(3,1); u = u / norm(u);
    Xp = X + d0*u;

    nTrans = round(transient/dt);
    for i = 1:nTrans
        X  = rk4_step(@lorenz_rhs, X,  dt, sigma, rho, beta);
        Xp = rk4_step(@lorenz_rhs, Xp, dt, sigma, rho, beta);
        dvec = Xp - X;
        d = norm(dvec);
        if d == 0
            Xp = X + d0*u;
        else
            Xp = X + d0*(dvec/d);
        end
    end

    steps = round(T/dt);
    s = 0.0;
    count = 0;
    for k = 1:steps
        X  = rk4_step(@lorenz_rhs, X,  dt, sigma, rho, beta);
        Xp = rk4_step(@lorenz_rhs, Xp, dt, sigma, rho, beta);

        if mod(k, renormEvery) == 0
            dvec = Xp - X;
            d = norm(dvec);
            if d == 0
                continue;
            end
            s = s + log(d/d0);
            count = count + 1;
            Xp = X + d0*(dvec/d);
        end
    end

    totalTime = count * renormEvery * dt;
    lam = s / totalTime;
end

Simulink note (workflow-only, no extra diagram): You can implement the Lorenz model in Simulink using three Integrator blocks and algebra blocks for the right-hand side, then run two copies with slightly different initial conditions. Use a MATLAB script to compute \( d(t)=\|\mathbf{x}_2(t)-\mathbf{x}_1(t)\| \) and apply the renormalization loop externally (post-processing), which matches the estimator in this lesson.

6.5 Wolfram Mathematica Implementation

File: Chapter20_Lesson3.nb

(* Chapter20_Lesson3.nb
   System Dynamics (Control Engineering) - Chapter 20, Lesson 3
   Sensitivity to Initial Conditions and Lyapunov Exponents (Intro)

   This is a plain-text Mathematica notebook expression.
   Open it in Mathematica; it should render as a notebook.
*)

Notebook[{
  Cell["Chapter 20 - Lesson 3: Sensitivity to Initial Conditions and Lyapunov Exponents (Intro)", "Title"],

  Cell["Part A) Logistic map: LLE via time-average of log|f'(x)|", "Section"],
  Cell[
"(* Logistic map: x_{n+1} = r x_n (1 - x_n)
   LLE: lambda ~ (1/N) Sum_{n=0}^{N-1} Log[ Abs[ r (1 - 2 x_n ) ] ]
*)

ClearAll[logisticStep, lyapLogistic];
logisticStep[x_, r_] := r x (1 - x);

lyapLogistic[r_, x0_:0.234, N_:100000, discard_:5000] := Module[{x=x0, s=0.0, fp},
  Do[x = logisticStep[x, r], {discard}];
  Do[
    fp = r (1 - 2 x);
    s = s + Log[Abs[fp] + 10^-300];
    x = logisticStep[x, r];
  , {N}];
  s/N
];

rList = {3.2, 3.5, 3.9, 4.0};
TableForm[Table[{r, lyapLogistic[r]}, {r, rList}], TableHeadings -> {None, {\"r\", \"lambda\"}}]
", "Input"],

  Cell["Part B) Lorenz system: LLE via two-trajectory renormalization", "Section"],
  Cell[
"(* Lorenz ODE:
   x' = sigma (y - x)
   y' = x (rho - z) - y
   z' = x y - beta z
*)

ClearAll[lorenz, rk4Step, lleLorenz];
lorenz[{x_, y_, z_}, sigma_:10.0, rho_:28.0, beta_:8.0/3.0] := {
  sigma (y - x),
  x (rho - z) - y,
  x y - beta z
};

rk4Step[f_, y_, dt_] := Module[{k1, k2, k3, k4},
  k1 = f[y];
  k2 = f[y + 0.5 dt k1];
  k3 = f[y + 0.5 dt k2];
  k4 = f[y + dt k3];
  y + (dt/6) (k1 + 2 k2 + 2 k3 + k4)
];

lleLorenz[X0_:{1.0,1.0,1.0}, dt_:0.01, T_:120.0, transient_:20.0, renormEvery_:10, d0_:10^-8,
          sigma_:10.0, rho_:28.0, beta_:8.0/3.0, seed_:0] := Module[
  {X = N@X0, u, Xp, steps, nTrans, s=0.0, count=0, k, dvec, d, f},
  SeedRandom[seed];
  u = Normalize[RandomVariate[NormalDistribution[], 3]];
  Xp = X + d0 u;

  f[v_] := lorenz[v, sigma, rho, beta];

  nTrans = Round[transient/dt];
  Do[
    X  = rk4Step[f, X, dt];
    Xp = rk4Step[f, Xp, dt];
    dvec = Xp - X;
    d = Norm[dvec];
    If[d == 0, Xp = X + d0 u, Xp = X + d0 dvec/d];
  , {nTrans}];

  steps = Round[T/dt];
  Do[
    X  = rk4Step[f, X, dt];
    Xp = rk4Step[f, Xp, dt];

    If[Mod[k, renormEvery] == 0,
      dvec = Xp - X;
      d = Norm[dvec];
      If[d != 0,
        s = s + Log[d/d0];
        count = count + 1;
        Xp = X + d0 dvec/d;
      ];
    ];
  , {k, 1, steps}];

  s/(count * renormEvery * dt)
];

lleLorenz[]
", "Input"]
}]

Problem 1 (1D Map Lyapunov Formula): Consider \( x_{k+1} = f(x_k) \) where \( f \) is differentiable. Starting from the linearized perturbation \( \delta x_{k+1}=f'(x_k)\delta x_k \), derive \( \lambda = \lim_{k\to\infty}\frac{1}{k}\sum_{i=0}^{k-1}\ln|f'(x_i)| \).

Solution: Recursively expand: \( \delta x_1=f'(x_0)\delta x_0 \), \( \delta x_2=f'(x_1)\delta x_1=f'(x_1)f'(x_0)\delta x_0 \), and in general \( \delta x_k=\left(\prod_{i=0}^{k-1}f'(x_i)\right)\delta x_0 \). Taking logs:

\[ \ln\left|\frac{\delta x_k}{\delta x_0}\right| = \sum_{i=0}^{k-1}\ln|f'(x_i)|. \]

Divide by \( k \) and take the limit \( k\to\infty \) (assuming convergence) to obtain the formula.

Problem 2 (Logistic Map at r=4): For the logistic map \( x_{k+1}=4x_k(1-x_k) \), show that the Lyapunov exponent is \( \lambda = \ln 2 \) for typical initial conditions (a classical result).

Solution (guided derivation): The derivative is \( f'(x)=4(1-2x) \), hence \( \lambda=\lim_{k\to\infty}\frac{1}{k}\sum \ln|4(1-2x_k)| \). For \( r=4 \), the map has an invariant density \( \rho(x)=\frac{1}{\pi\sqrt{x(1-x)}} \) on \( (0,1) \) (ergodic case), so the time average equals the space average:

\[ \lambda = \int_0^1 \ln|4(1-2x)|\,\rho(x)\,dx = \ln 2. \]

This is a deep fact relying on ergodicity/mixing at \( r=4 \). In practice, your code in this lesson should numerically approach \( \ln 2 \approx 0.6931 \).

Problem 3 (Linear ODE Benchmark): For \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \) with constant \( \mathbf{A} \), prove that \( \lambda_{\max}=\max_i \operatorname{Re}(\mu_i) \), where \( \mu_i \) are eigenvalues of \( \mathbf{A} \).

Solution: If \( \mathbf{A} \) is diagonalizable, \( \mathbf{A}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1} \), then \( e^{\mathbf{A}t}=\mathbf{V}e^{\mathbf{\Lambda}t}\mathbf{V}^{-1} \). Each eigen-direction scales like \( e^{\mu_i t} \); thus the dominant asymptotic growth rate in norm is governed by the largest real part:

\[ \lim_{t\to\infty}\frac{1}{t}\ln\|e^{\mathbf{A}t}\| = \max_i \operatorname{Re}(\mu_i). \]

This matches Lyapunov exponents for linear systems and provides a correctness check for numerical estimators.

Problem 4 (Renormalization Estimator): Suppose you renormalize every \( \Delta t \) (for ODEs) so that after each interval the separation is reset from \( d \) to \( d_0 \). Show that the estimator

\[ \hat{\lambda} = \frac{1}{N\Delta t}\sum_{i=1}^{N}\ln\left(\frac{d_i}{d_0}\right) \]

is consistent with \( \lambda_{\max} \) under typical assumptions (ergodicity and sufficiently small \( d_0 \)).

Solution (explanation): Over each interval, the local stretching factor is approximately \( d_i/d_0 \approx \exp(\lambda_i \Delta t) \), where \( \lambda_i \) is the local finite-time rate. Summing logs:

\[ \sum_{i=1}^{N}\ln\left(\frac{d_i}{d_0}\right) \approx \sum_{i=1}^{N}\lambda_i \Delta t. \]

Dividing by total time \( N\Delta t \) yields the time-average of finite-time rates, which converges to the asymptotic Lyapunov exponent when averages along typical trajectories converge.

• Sensitivity to initial conditions is quantified by exponential growth of perturbations: \( d(t)\approx d(0)e^{\lambda t} \) (ODEs) or \( d_k\approx d_0 e^{\lambda k} \) (maps).
• For 1D maps, \( \lambda=\lim_{k\to\infty}\frac{1}{k}\sum \ln|f'(x_k)| \) gives a direct numerical estimator.
• For ODEs, Lyapunov exponents arise from variational dynamics \( \dot{\delta\mathbf{x}}=\mathbf{A}(t)\delta\mathbf{x} \) and can be estimated via two-trajectory renormalization.
• Positive largest Lyapunov exponent (\( \lambda_{\max} > 0 \)) indicates strong sensitivity and finite predictability.