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) — the average exponential rate at which two initially close trajectories pull apart — is the gold-standard definition of chaos, converting the intuition of "butterfly effects" into a precise, computable number.

TL;DR: Lyapunov exponents measure how fast nearby trajectories in a system diverge; a positive value signals chaos. This article covers the full spectrum, the Kaplan–Yorke dimension formula for fractal attractors, three computation methods (direct formula, Wolf–Swinney–Vastano, Gram–Schmidt), Takens' embedding for real data, and the link to predictability horizons like weather forecasting limits.

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 →