Detecting & Categorizing Chaos

Being able to detect and distinguish chaotic from regular behavior is crucial in the study of dynamical systems. Most of the time a positive maximum lyapunov exponent and a bounded system indicate chaos.

However, the convergence of the Lyapunov exponent can be slow, or even misleading, as the types of chaotic behavior vary with respect to their predictability. There are many alternatives, some more efficient and some more accurate in characterizing chaotic and regular motion. Some of these methods are included in DynamicalSystems.jl.

Performance depends on the solver

Notice that the performance of functions that use ContinuousDynamicalSystems depend crucially on the chosen solver. Please see the documentation page on Choosing a solver for an in-depth discussion.

Generalized Alignment Index

"GALI" for sort, is a method that relies on the fact that initially orthogonal deviation vectors tend to align towards the direction of the maximum Lyapunov exponent for chaotic motion. It is one of the most recent and cheapest methods for distinguishing chaotic and regular behavior, introduced first in 2007 by Skokos, Bountis & Antonopoulos.

gali(ds::DynamicalSystem, tmax, k::Int | Q0; kwargs...) -> GALI_k, t

Compute $\text{GALI}_k$[Skokos2007] for a given k up to time tmax. Return $\text{GALI}_k(t)$ and time vector $t$.

The third argument, which sets the order of gali, can be an integer k, or a matrix with its columns being the deviation vectors (then k = size(Q0)[2]). In the first case random orthonormal vectors are chosen.

Keyword Arguments

  • threshold = 1e-12 : If GALI_k falls below the threshold iteration is terminated.
  • Δt = 1 : Time-step between deviation vector normalizations. For continuous systems this is approximate.
  • u0 : Initial state for the system. Defaults to get_state(ds).
  • diffeq is a NamedTuple (or Dict) of keyword arguments propagated into init of DifferentialEquations.jl. See trajectory for examples. Only valid for continuous systems.


The Generalized Alignment Index, $\text{GALI}_k$, is an efficient (and very fast) indicator of chaotic or regular behavior type in $D$-dimensional Hamiltonian systems ($D$ is number of variables). The asymptotic behavior of $\text{GALI}_k(t)$ depends critically on the type of orbit resulting from the initial condition. If it is a chaotic orbit, then

\[\text{GALI}_k(t) \sim \exp\left[\sum_{j=1}^k (\lambda_1 - \lambda_j)t \right]\]

with $\lambda_j$ being the j-th Lyapunov exponent (see lyapunov, lyapunovspectrum). If on the other hand the orbit is regular, corresponding to movement in $d$-dimensional torus with $1 \le d \le D/2$ then it holds

\[\text{GALI}_k(t) \sim \begin{cases} \text{const.}, & \text{if} \;\; 2 \le k \le d \; \; \text{and} \; \;d > 1 \\ t^{-(k - d)}, & \text{if} \;\; d < k \le D - d \\ t^{-(2k - D)}, & \text{if} \;\; D - d < k \le D \end{cases}\]

Traditionally, if $\text{GALI}_k(t)$ does not become less than the threshold until tmax the given orbit is said to be chaotic, otherwise it is regular.

Our implementation is not based on the original paper, but rather in the method described in[Skokos2016b], which uses the product of the singular values of $A$, a matrix that has as columns the deviation vectors.

Performance Notes

This function uses a tangent_integrator. For loops over initial conditions and/or parameter values one should use the low level method that accepts an integrator, and reinit! it to new initial conditions. See the "advanced documentation" for info on the integrator object. The low level method is

ChaosTools.gali(tinteg, tmax, Δt, threshold)

(section 5.3.1 and ref. [85] therein), Lecture Notes in Physics 915, Springer (2016)

Discrete Example

We will use 3 coupled standard maps as an example for a discrete system:

using DynamicalSystems
using PyPlot
M = 3; ks = 3ones(M); Γ = 0.1;
stable = [π, π, π, 0.01, 0, 0] .+ 0.1
chaotic = rand(2M)

ds = Systems.coupledstandardmaps(M, stable; ks, Γ)
6-dimensional discrete dynamical system
 state:       [3.24159, 3.24159, 3.24159, 0.11, 0.1, 0.1]
 rule f:      CoupledStandardMaps
 in-place?    true
 jacobian:    CoupledStandardMaps
 parameters:  [3.0, 3.0, 3.0, 0.1]

For this example we will see the behavior of GALI for a stable orbit

fig = figure(figsize = (9,5))
tr = trajectory(ds, 100000)

plot(tr[:,1], tr[:,1+M], label="stable", marker="o", ms=1, lw=0)

for k in [4, 5, 6]
    g, t = gali(ds, 1e5, k; threshold=1e-12)
    lt = log10.(t); lg = log10.(g)
    plot(lt, lg, label="GALI_$(k)")
lt = 2:0.5:5.5
plot(lt, -2(lt .- 3), label="slope -2")
plot(lt, -4(lt .- 3), label="slope -4")
plot(lt, -6(lt .- 3), label="slope -6")

xlim(2, 5.5)
ylim(-12, 2)
fig.tight_layout(pad=0.3); fig

Continuous Example

As an example of a continuous system, let's see the Henon-Heiles:

using DynamicalSystems
using PyPlot, OrdinaryDiffEq
Δt = 1.0
diffeq = (abstol=1e-9, retol=1e-9, alg = Vern9(), maxiters = typemax(Int))
sp = [0, .295456, .407308431, 0] # stable periodic orbit: 1D torus
qp = [0, .483000, .278980390, 0] # quasiperiodic orbit: 2D torus
ch = [0, -0.25, 0.42081, 0]      # chaotic orbit
ds = Systems.henonheiles(sp)
4-dimensional continuous dynamical system
 state:       [0.0, 0.295456, 0.407308, 0.0]
 rule f:      hhrule!
 in-place?    true
 jacobian:    hhjacob!
 parameters:  nothing

Let's see what happens with a quasi-periodic orbit:

fig = figure(figsize = (9,5))
tr = trajectory(ds, 10000.0, qp; Δt, diffeq...)
plot(tr[:,1], tr[:,3], alpha = 0.5,
label="qp",marker="o",markersize=2, linewidth=0)

for k in [2,3,4]
    g, t = gali(ds, 10000.0, k; u0 = qp, Δt, diffeq...)
    loglog(t, g, label="GALI_$(k)")
    if k == 2
        loglog(t, 1 ./ t.^(2k-4), label="slope -$(2k-4)")
        loglog(t, 100 ./ t.^(2k-4), label="slope -$(2k-4)")
ylim(1e-12, 2)
fig.tight_layout(pad=0.3); fig