Chaos and ensembles¶

Lyapunov exponents with error bars, covariant Lyapunov vectors, Poincaré maps, and the fate of many trajectories integrated as one batch.

Lyapunov spectrum¶

The exponents are computed from the tangent dynamics of the same discrete map you integrate, over several independent trajectories, with a standard error from block averages.

import numpy as np

from conatus.bifurcation import ode_problem
from conatus.chaos import lyapunov_spectrum
from conatus.discretize import discretize
from conatus.timestep import IVPSolver

lorenz = ode_problem(["x", "y", "z"], {"sigma": 10.0, "r": 28.0, "b": 8 / 3},
                     lambda X, P: [P["sigma"] * (X["y"] - X["x"]),
                                   P["r"] * X["x"] - X["y"] - X["x"] * X["z"],
                                   X["x"] * X["y"] - P["b"] * X["z"]], name="Lorenz")
s = IVPSolver(discretize(lorenz), "RK443", 0.01)      # scheme and step: your choice

rng = np.random.default_rng(0)
X0 = np.column_stack([rng.uniform(-10, 10, 16), rng.uniform(-10, 10, 16),
                      rng.uniform(10, 40, 16)])          # 16 independent trajectories
res = lyapunov_spectrum(s, X0, n_vectors=3, transient_steps=1000,
                        tangent_transient_steps=1000, steps=5000, reortho_every=10,
                        inner="euclidean", n_blocks=5, seed=1)
print(res.explain())                                     # exponents with standard errors
D, dD = res.kaplan_yorke()

inner= declares the inner product of the Gram–Schmidt steps; res.kaplan_yorke() returns the Kaplan–Yorke dimension with its error. With clv_window=... the result also carries the covariant Lyapunov vectors. poincare_section records the crossings of a surface \(g(x) = 0\), located to a declared tolerance, and their return map.

From a native simulation¶

sim.lyapunov computes the same spectrum, and the covariant vectors, from a native model; lengths are durations of the problem (each a whole number of steps or of intervals reortho_every * dt):

import numpy as np

import conatus as cn

L = 22.0
m = cn.Model(cn.Box(x=cn.periodic(L)), name="Kuramoto-Sivashinsky L=22")
m.scalar("u")
m.equations("dt(u) + lap(u) + lap(lap(u)) = -u*dx(u)")
sim = m.simulate(N=32, scheme="RK443", dt=0.05)
sim.initial(u=lambda x: 0.3 * np.cos(2 * np.pi * x / L) + 0.2 * np.sin(4 * np.pi * x / L))
res = sim.lyapunov(12, transient=100.0, align=50.0, T=200.0, reortho_every=10,
                   inner="euclidean", n_blocks=4, seed=1, zero_mean="u", covariant=True,
                   clv_window=40.0, clv_backward=40.0, clv_checkpoints=8)
print(res.explain())
dom = res.domination(5.0)
print(dom.explain())

zero_mean="u" keeps the tangent vectors in the invariant zero-mean subspace (all axes periodic); project=f declares another invariant subspace. starts=[{...}, ...] runs several trajectories in one batch (each dict as for sim.initial); the simulation’s own state is not advanced.
covariant=True computes the covariant Lyapunov vectors (Ginelli et al. 2007) over a window of length clv_window that ends clv_backward before the end of the run. clv_checkpoints=s stores s checkpoints instead of every interval and recomputes them once (memory and recomputation check in explain()).
explain() also gives the backward convergence (the covariant vectors from two independent random starts of the backward iteration, compared: the error left by clv_backward), the covariance residual (each vector advanced by the exact tangent over one interval is parallel to the vector there, to its rounding floor) and the smallest principal angle at every split.
res.angles(j), res.pair_angle(i, j): angles at every kept time; res.clv_grid("u", i): a covariant vector as grid values.
res.domination(tau): the fraction of windows of length \(\tau\) where the local growth rates of \(\gamma_i\) and \(\gamma_j\) (\(i < j\)) are in the wrong order; modes joined by such violations form clusters, the leading one is the set of physical modes (dom.n_physical; eight for KS \(L = 22\), tutorial 15. Covariant Lyapunov vectors of Kuramoto–Sivashinsky: tangencies and the physical modes).

Ensembles¶

An ensemble integrates many members of one model as a single batched state, on CPU or GPU. ensemble_fate classifies each member with your criterion and stops when every member is decided or the horizon is reached; members still undecided are reported as such, never assigned.

from conatus.ensembles import Ensemble, ensemble_fate

s20 = IVPSolver(discretize(ode_problem(
    ["x", "y", "z"], {"sigma": 10.0, "r": 20.0, "b": 8 / 3},
    lambda X, P: [P["sigma"] * (X["y"] - X["x"]), P["r"] * X["x"] - X["y"] - X["x"] * X["z"],
                  X["x"] * X["y"] - P["b"] * X["z"]], name="Lorenz r=20")), "RK443", 0.01)
c = np.sqrt(8 / 3 * 19.0)
Cp, Cm = np.array([c, c, 19.0]), np.array([-c, -c, 19.0])   # the two stable equilibria
starts = np.column_stack([rng.uniform(-20, 20, 256), rng.uniform(-20, 20, 256),
                          np.full(256, 19.0)])
ens = Ensemble(s20, len(starts), rhs="batched")             # 256 members in one batch
ens.set_states(starts)

def attractor(e):                                           # declared classifier
    X = e.states()
    return np.where(np.linalg.norm(X - Cp, axis=1) < 1, 1,
                    np.where(np.linalg.norm(X - Cm, axis=1) < 1, -1, 0))

fate = ensemble_fate(ens, attractor, check_every=100, min_steps=100, max_steps=10000,
                     classifier="|X - C+-| < 1")
print(fate.explain())                                       # fractions with intervals

The proportions come with confidence intervals (proportion_ci).

See also: tutorial 12. Chaos and ensembles: Lyapunov spectrum, covariant vectors, Poincaré map, basins (Lorenz attractor, Sprott’s exponents); conatus.chaos.lyapunov_spectrum, conatus.chaos.poincare_section, Ensemble, ensemble_fate.