Measuring Chaos: Lyapunov Exponents and Fractal Dimensions#
“Chaotic” is a quantitative claim: nearby trajectories separate
exponentially, at a rate given by the largest Lyapunov exponent
\(\lambda_1 > 0\), and the attractor they settle on has a
non-integer dimension. This tutorial measures both with
physicskit.chaos, on three systems of increasing complexity: the
logistic map, the Hénon map and the Lorenz flow. At each step the
numerical answer is checked against a result known in closed form, so
you can tell a real measurement from an artefact of the method.
Everything here is dimensionless. The maps count time in iterations and the Lorenz equations in their own natural time unit.
The logistic map: from order to chaos#
The logistic map \(x_{n+1} = r x_n (1 - x_n)\)
(LogisticMap) is the simplest system with a
period-doubling route to chaos (May, Nature 261, 459, 1976). Iterate
it past its transient and count the distinct points it keeps visiting,
then estimate \(\lambda\) with
map_lyapunov_spectrum(), which evolves a
tangent vector along the orbit (the Benettin QR method):
import numpy as np
import physicskit as pk
from physicskit.chaos.utils import map_lyapunov_spectrum
for r in [2.8, 3.2, 3.5, 3.56, 4.0]:
logistic = pk.chaos.LogisticMap(r=r)
orbit = logistic.trajectory(np.array([0.2]), n_iter=2000)
period = len(np.unique(np.round(orbit[-64:, 0], 6)))
lam = map_lyapunov_spectrum(logistic, np.array([0.2]), n_iter=20000)[0]
print(f"r={r:<5} distinct points={period:<3} lambda={lam:+.4f}")
# r=2.8 distinct points=1 lambda=-0.2231
# r=3.2 distinct points=2 lambda=-0.9163
# r=3.5 distinct points=4 lambda=-0.8725
# r=3.56 distinct points=8 lambda=-0.0771
# r=4.0 distinct points=64 lambda=+0.6933
The period doubles 1 → 2 → 4 → 8, and \(\lambda\) stays negative (stable cycles attract nearby orbits) until the cascade ends and \(\lambda\) turns positive. Three of these rows can be checked exactly:
At \(r = 2.8\) the orbit sits on the fixed point \(x^* = 1 - 1/r\), where \(f'(x^*) = 2 - r\), so \(\lambda = \ln|2 - r| = \ln 0.8 = -0.2231\).
At \(r = 3.2\) the 2-cycle has multiplier \(f'(x_1) f'(x_2) = 4 + 2r - r^2\), so \(\lambda = \tfrac12 \ln|4 + 2r - r^2| = \tfrac12\ln 0.16 = -0.9163\).
At \(r = 4\) the map is conjugate to the tent map, giving \(\lambda = \ln 2 = 0.6931\) (Ulam and von Neumann, 1947).
Near \(r = 3.56\), \(\lambda\) is close to zero because the cascade is about to reach the Feigenbaum point \(r_\infty = 3.5699...\), where the period becomes infinite.
The Hénon map: a strange attractor in two dimensions#
The Hénon map (HenonMap, Hénon, Commun. Math.
Phys. 50, 69, 1976) has two Lyapunov exponents. Its Jacobian has the
constant determinant \(-b\), so every area shrinks by \(|b|\)
per step. That fixes the sum of the exponents exactly,
\(\lambda_1 + \lambda_2 = \ln b\), which gives a sharp test of the
numerical spectrum:
from physicskit.chaos.utils import box_counting_dimension
henon = pk.chaos.HenonMap(a=1.4, b=0.3)
spectrum = map_lyapunov_spectrum(henon, np.array([0.1, 0.1]), n_iter=50000)
print(spectrum.round(4), spectrum.sum().round(4), np.log(0.3).round(4))
# [ 0.4187 -1.6226] -1.204 -1.204
One exponent is positive, so the map stretches, and the sum is negative, so it contracts area overall. An attractor with both properties is a fractal. The Kaplan-Yorke conjecture predicts its dimension from the spectrum alone, \(D_{KY} = 1 + \lambda_1/|\lambda_2|\) (Kaplan and Yorke, 1979). Compare that with a direct box count of the attractor:
kaplan_yorke = 1 + spectrum[0] / abs(spectrum[1])
orbit = henon.trajectory(np.array([0.1, 0.1]), n_iter=200000)[1000:]
D0, eps, counts = box_counting_dimension(orbit)
print(round(kaplan_yorke, 3), round(D0, 3))
# 1.258 1.303
Both land near the accepted \(D \approx 1.26\). The box count comes
out slightly high because a finite orbit covers the fine boxes unevenly.
box_counting_dimension() returns
eps and counts so you can plot \(\log N(\epsilon)\) against
\(\log(1/\epsilon)\) and check that the fit really covers a straight
section.
The Lorenz attractor: flows and the zero exponent#
For a continuous flow, benettin_lyapunov_spectrum()
integrates the trajectory and its tangent space together. The Lorenz
system (Lorenz) has divergence
\(-(\sigma + 1 + \beta)\) everywhere, so again the sum is known
exactly, and any bounded, non-fixed-point trajectory of an autonomous
flow must have one exponent equal to zero (the direction along the flow
itself):
from physicskit.chaos.utils import benettin_lyapunov_spectrum, correlation_dimension
lorenz = pk.chaos.Lorenz(sigma=10.0, rho=28.0, beta=8.0 / 3.0)
spectrum = benettin_lyapunov_spectrum(
lorenz, np.array([1.0, 1.0, 1.0]), dt=0.01, n_steps=100000, n_transient=2000
)
print(spectrum.round(3), spectrum.sum().round(3))
# [ 9.040e-01 -1.000e-03 -1.457e+01] -13.667
The spectrum \((0.904, 0, -14.57)\) matches the literature values (e.g. Sprott, Chaos and Time-Series Analysis, 2003), and the sum is \(-(10 + 1 + 8/3) = -13.667\) to the digits shown.
The attractor’s correlation dimension \(D_2\) comes from the Grassberger-Procaccia correlation sum (Phys. Rev. Lett. 50, 346, 1983), as the slope of \(\log C(\epsilon)\) against \(\log\epsilon\). That slope is only meaningful over the scaling region. At large \(\epsilon\) the sum saturates at the size of the attractor, and at small \(\epsilon\) it runs out of point pairs:
t, states = lorenz.trajectory(np.array([1.0, 1.0, 1.0]), dt=0.01, n_steps=200000)
points = states[5000::50][:3000]
D2_default, _, _ = correlation_dimension(points)
D2, eps, C = correlation_dimension(points, eps_min=0.5, eps_max=5.0)
print(round(D2_default, 3), round(D2, 3))
# 1.681 2.047
With the default range (the 5th to 50th percentile of pair distances)
the fit reaches into the saturated region and returns 1.68. Restricting
it to \(0.5 \le \epsilon \le 5\) gives 2.05, Grassberger and
Procaccia’s own value (\(2.05 \pm 0.01\)). Before trusting a
dimension estimate, look at the returned eps and C on a log-log
plot.
Lyapunov exponents from a single time series#
Experiments rarely give you the equations of motion, only a signal.
Takens’ theorem (1981) says that delay coordinates
\((x(t), x(t+\tau), x(t+2\tau), \dots)\) reconstruct an attractor
equivalent to the original. Rosenstein’s method (Physica D 65, 117,
1993) then follows nearest neighbours in that reconstruction and
averages their log separation,
average_log_divergence(). Its slope is
\(\lambda_1\). Here we use only the \(x\) component of the
Lorenz run:
from physicskit.chaos.utils import average_log_divergence
lorenz = pk.chaos.Lorenz()
t, states = lorenz.trajectory(np.array([1.0, 1.0, 1.0]), dt=0.01, n_steps=65000)
x = states[5000:, 0] # observe only x(t), sampled every dt = 0.01
steps, divergence = average_log_divergence(x, dim=3, tau=10, min_tsep=100, max_iter=300)
for lo, hi in [(0, 50), (50, 200)]:
slope = np.polyfit(steps[lo:hi] * 0.01, divergence[lo:hi], 1)[0]
print(f"fit over steps {lo:>3}-{hi:<3}: lambda = {slope:.3f}")
# fit over steps 0-50 : lambda = 2.266
# fit over steps 50-200: lambda = 0.915
The first half time unit is a transient. The neighbour pair starts
separated along every direction, and it takes a while to rotate onto the
most unstable one, so the curve is steeper there. Fitting that part
overestimates \(\lambda_1\) by more than a factor of two. Fitting
the straight section after it gives 0.915, within 1% of the Benettin
value computed from the equations. The convenience wrapper
rosenstein_lyapunov() fits from step 0, so
for a signal like this one plot the curve and choose the linear region
yourself.