The hallmark of a chaotic system is sensitive dependence on initial conditions: two trajectories starting arbitrarily close diverge exponentially. Lyapunov exponents quantify that divergence rate. A positive maximal Lyapunov exponent (MLE) is the gold-standard definition of chaos, converting the intuition of "butterfly effects" into a precise, computable number.

1. Definition and Geometric Meaning

Consider a smooth dynamical system ẋ = f(x) (flow) or xₙ₊₁ = F(xₙ) (map). Starting from x₀, perturb the initial condition by a tiny vector δ₀ with |δ₀| = ε. After time t (or n iterations) the perturbation evolves to δ(t) via the linearised (tangent) dynamics:

δ(t) = DFᵗ(x₀) · δ₀,    DFᵗ ≡ ∂Fᵗ/∂x|_{x₀}

The Lyapunov exponent for direction δ₀ is:

λ(x₀, δ₀) = lim_{t→∞} (1/t) · ln |δ(t)| / |δ₀|

By Oseledets' multiplicative ergodic theorem (1968), for almost every x₀ (with respect to an ergodic invariant measure), this limit converges to one of at most d values (where d is the state-space dimension), independent of the starting direction δ₀ for a full-measure set of directions. These values λ₁ ≥ λ₂ ≥ … ≥ λ_d form the Lyapunov spectrum.

Exponent sign Geometric meaning Typical system
λ₁ > 0 Exponential divergence — chaos Lorenz, logistic r > r∞
λ₁ = 0 Neutral / marginal — bifurcation or quasi-periodic Limit cycle, KAM tori
λ₁ < 0 Convergence to attractor Stable fixed point, stable limit cycle
Σλᵢ < 0 Phase-volume contraction (dissipative) All physical attractors
Σλᵢ = 0 Volume-preserving (Hamiltonian) Pendulum, N-body

2. Kaplan–Yorke Dimension

The fractal (information) dimension of a strange attractor is linked to the Lyapunov spectrum by the Kaplan–Yorke conjecture (1979), confirmed numerically for many systems:

D_KY = j + Σᵢ₌₁ʲ λᵢ / |λⱼ₊₁|

where j is the largest index such that Σᵢ₌₁ʲ λᵢ ≥ 0. For the Lorenz system with σ=10, ρ=28, β=8/3 the spectrum is approximately λ₁ ≈ 0.906, λ₂ = 0, λ₃ ≈ −14.57, giving D_KY ≈ 2 + 0.906/14.57 ≈ 2.062. The famous Lorenz attractor is indeed a non-integer-dimensional fractal set.

Constraints on the Lyapunov spectrum

• One exponent is always 0 along the flow direction for autonomous continuous systems.
• For Hamiltonian systems exponents pair symmetrically: λᵢ = −λ_{d+1−i}.
• For a dissipative system: Σλᵢ = ∫ ∇·f dμ (time-averaged divergence of the vector field).
• For the Lorenz system: Σλᵢ = −(σ + 1 + β) = −13.67.

3. Computing the MLE from a 1-D Map

For a one-dimensional map xₙ₊₁ = F(xₙ) the MLE reduces to:

λ = lim_{N→∞} (1/N) · Σₙ₌₀^{N−1} ln |F′(xₙ)|

For the logistic map F(x) = rx(1−x):

F′(x) = r(1 − 2x), λ(r) = lim_{N→∞} (1/N) Σ ln |r(1 − 2xₙ)|
Parameter r Regime MLE λ
r = 2.0 Stable fixed point ln(r−2) → −∞ for fixed pt
r = 3.2 Period-2 cycle ≈ −0.345
r = 3.5 Period-4 cycle ≈ −0.163
r ≈ 3.5699 Onset of chaos (r∞) = 0
r = 3.8 Chaos (interior crisis) ≈ +0.433
r = 4.0 Full chaos (tent map conjugacy) = ln 2 ≈ 0.693

The bifurcation points of the logistic map are exactly the zeroes of λ(r), confirming that λ = 0 is both necessary and sufficient for parameter-space bifurcations in this family.

4. Wolf–Swinney–Vastano Algorithm

For a continuous flow (Lorenz, Rössler, etc.) or an experimental time series, direct differentiation is unavailable. The WSV algorithm (Wolf, Swift, Swinney & Vastano, 1985) estimates the maximal Lyapunov exponent from a single trajectory:

WSV Algorithm — pseudocode

1. Integrate orbit to find fiducial trajectory x(t)
2. Choose nearby point x̃(0) with |δ₀| = |x̃(0) − x(0)| = ε (small)
3. Integrate both trajectories for evolution time T_ev
4. Measure new separation d₁ = |x̃(T_ev) − x(T_ev)|
5. Accumulate: L += ln(d₁ / ε)
6. Rescale: replace x̃ → x(T_ev) + ε · δ₁ / d₁  (renormalise, preserve direction)
7. Repeat from step 3; after N cycles:
   λ₁ = (1 / N·T_ev) · L

Key parameters: T_ev should be short enough to avoid orbit tangencies but long enough to reduce statistical noise; ε must be small relative to attractor size yet large relative to floating-point noise. Well-tested defaults: ε ~ 10⁻⁸, T_ev ~ 0.5–2.0 time units.

5. The Gram–Schmidt Full Spectrum Method

To compute all d exponents simultaneously, integrate the tangent (variational) equations alongside the trajectory. The tangent vectors must be periodically Gram–Schmidt orthonormalised to prevent collapse onto the maximal direction.

dδᵢ/dt = Df(x(t)) · δᵢ for i = 1, …, d

After each renormalisation step at time τ, record the log-stretching of each basis vector before orthonormalisation. The accumulated sum / total time → λᵢ. This is the foundation of the Benettin et al. (1980) algorithm and its numerical implementation in most scientific computing packages.

Gram–Schmidt QR interpretation

At each step the d tangent vectors form a matrix Φ. Factor Φ = QR with Q orthogonal; the diagonal entries rᵢᵢ of R give the stretching along each Oseledets direction. After N steps:
λᵢ = (1 / N·τ) · Σₖ ln |rᵢᵢ⁽ᵏ⁾|

6. Lyapunov Exponents from Experimental Time Series

When only a scalar measurement series {s₁, s₂, …, sₙ} is available (e.g. ECG, EEG, stock prices), the phase space must first be reconstructed by Takens' delay embedding:

x(t) = [s(t), s(t + τ), s(t + 2τ), …, s(t + (m−1)τ)] ∈ ℝᵐ

By Takens' theorem the reconstructed attractor is diffeomorphic to the original if m ≥ 2d_A + 1 (where d_A is the attractor dimension). Rosenstein et al. (1993) and Kantz (1994) provide numerically robust algorithms for the MLE on embedded data — the Rosenstein method tracks nearest-neighbour divergence across all initial conditions simultaneously:

d_j(i) = 〈ln Δⱼ(i)〉_j,   slope of d_j(i) vs. i·Δt ≈ λ₁

Needed parameters: embedding dimension m (from false nearest neighbours, FNN), delay τ (from autocorrelation zero or first minimum of mutual information), theiler window W (to exclude temporal neighbours).

7. Connection to Entropy and Predictability

The Kolmogorov–Sinai (KS) entropy h_KS quantifies the rate of information production of a dynamical system. Pesin's formula connects it directly to positive Lyapunov exponents:

h_KS = Σ_{λᵢ > 0} λᵢ (Pesin's identity)

Predictability horizon: an initial uncertainty of ε grows to order-1 uncertainty after time t_pred ≈ −ln(ε) / λ₁. For the Lorenz system with λ₁ ≈ 0.906 and ε = 10⁻⁵ the predictability horizon is ≈ 13 time units — explaining why numerical weather prediction becomes unreliable beyond ~2 weeks.

System MLE λ₁ Predictability (ε=10⁻⁵)
Lorenz (σ=10, ρ=28, β=8/3) ≈ 0.906 ≈ 12.7 time units
Rössler (a=0.2, b=0.2, c=5.7) ≈ 0.071 ≈ 162 time units
Logistic r = 4.0 = ln 2 ≈ 0.693 ≈ 16.6 iterations
Double pendulum (energy ≈ 10E_min) ≈ 3.5 rad⁻¹ ≈ 3.2 radian-time
Solar system (outer planets) ≈ 1/(5 Myr) ≈ 58 Myr (1/λ₁)

8. Interactive: Real-Time MLE for the Logistic Map

Drag the r slider to watch the Lyapunov exponent update in real time. The upper panel shows the orbit diagram (last 200 iterates from a warm-up of 400); the lower panel plots λ(r) computed on-the-fly via the sum-of-log-derivatives formula. Positive λ is red (chaos), zero crossings mark bifurcations, and negative λ is green (periodicity).

λ ≈ ?

The full-scan mode sweeps r from 2.5 to 4.0 and renders λ(r) as a curve — bifurcation points appear as sharp zero-crossings, and the period-3 window at r≈3.828 shows a dip to strongly negative values.

← Bifurcation Diagrams Takens' Theorem →