Chaos Detection

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 is often very slow and the computation costly. There are many alternatives that are both more efficient and more accurate in characterizing chaotic and regular motion, some of which are included in DynamicalSystems.jl.

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 chaos, introduced first in 2007 by Skokos, Bountis & Antonopoulos.

#ChaosTools.gali — Function.

gali(ds::DynamicalSystem, k::Int, tmax [, ws]; kwargs...) -> GALI_k, t

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

ws is an optional argument containing the deviation vectors $w_i$ for $i \in [1,k]$ , expected either as a matrix with each column a deviation vector, or as a vector of vectors. If not given, random orthonormal vectors are chosen.

Keyword Arguments

threshold : If GALI_k falls below the threshold iteration is terminated. Default value is 1e-12.
dt=1.0 : Time step between variational vector normalizations for continuous systems.
diff_eq_kwargs : See trajectory.

Description

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 assymptotic behavior of $\text{GALI}_k(t)$ depends critically of the type of orbit resulting from the initial condition ds.state. 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_1$ being the maximum lyapunov exponent. 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.

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

Performance Notes

If you want to do repeated evaluations of gali for many initial conditions and for continuous systems, you can take advantage of the function:

gali(integrator, k, W, tmax, dt, threshold)

in conjuction with reinit!(integrator, W) for various W=cat(2, state, ws). See the source code on how to set-up the integrator and W for the first time.

References

[1] : Skokos, C. H. et al., Physica D 231, pp 30–54 (2007)

[2] : Skokos, C. H. et al., Chaos Detection and Predictability - Chapter 5 (section 5.3.1 and ref. [85] therein), Lecture Notes in Physics 915, Springer (2016)

source

Discrete Example

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

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

ds = Systems.coupledstandardmaps(M, stable; ks=ks, Γ = Γ)
tr = trajectory(ds, 100000)

subplot(2,2,1)
plot(tr[:,1], tr[:,1+M], alpha = 0.5,
label="stable",marker="o", ms=1, linewidth=0)
legend()
#
subplot(2,2,2)
for k in [2,3,4, 5, 6]
    g, t = gali(ds, k, 1e5; threshold=1e-12)
    lt = log10.(t); lg = log10.(g)

    plot(lt, lg, label="GALI_$(k)")
end
lt = 2:0.5:5.5
plot(lt, zeros(lt), label="const")
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, 1)
legend(fontsize=12)
tight_layout()


ds = Systems.coupledstandardmaps(M, chaotic; ks=ks, Γ = Γ)
tr = trajectory(ds, 100000)
subplot(2,2,3)
plot(tr[:,1], tr[:,1+M], alpha = 0.5,
label="chaotic",marker="o", ms=1, linewidth=0)
legend()


subplot(2,2,4)
ls = lyapunovs(ds, 100000)
for k in [2,3,4,5 ,6]
    ex = sum(ls[1] - ls[j] for j in 2:k)
    g, t = gali(ds, k, 1000)
    semilogy(t, exp.(-ex.*t), label="exp. k=$k")
    semilogy(t, g, label="GALI_$(k)")
end
legend(fontsize=12)
xlim(0,50)
ylim(1e-12, 1)

GALI Discrete

Continuous Example

As an example of a continuous system, let's see the henonhelies:

using DynamicalSystems
using PyPlot
figure(figsize=(10, 12))
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
dt = 1.0

subplot(3,2,1)
ds = Systems.henonhelies(sp)
diffeq = Dict(:abstol=>1e-9, :reltol=>1e-9, :solver => Tsit5())
tr = trajectory(ds, 10000.0, dt=dt, diff_eq_kwargs = diffeq)
plot(tr[:,1], tr[:,3], alpha = 0.5,
label="sp",marker="o",markersize=2, linewidth=0)
legend()

subplot(3,2,2)
for k in [2,3,4]
    g, t = gali(ds, k, 50000.0; dt = dt, diff_eq_kwargs = diffeq)
    if k < 4
        loglog(t, 1./t.^(k-1), label="slope -$(k-1)")
    else
        loglog(t, 1./t.^(2k-4), label="slope -$(2k-4)")
    end
    loglog(t, g, label="GALI_$(k)")
end
legend(fontsize=12)

subplot(3,2,3)
ds = Systems.henonhelies(qp)
diffeq = Dict(:abstol=>1e-9, :reltol=>1e-9, :solver => Tsit5())
tr = trajectory(ds, 10000.0, dt=dt, diff_eq_kwargs = diffeq)
plot(tr[:,1], tr[:,3], alpha = 0.5,
label="qp",marker="o",markersize=2, linewidth=0)
legend()

subplot(3,2,4)
for k in [2,3,4]
    g, t = gali(ds, k, 10000.0; dt = dt, diff_eq_kwargs = diffeq)
    loglog(t, 1./t.^(2k-4), label="slope -$(2k-4)")
    loglog(t, g, label="GALI_$(k)")
end
legend(fontsize=12)
tight_layout()

ds = Systems.henonhelies(ch)
diffeq = Dict(:abstol=>1e-6, :reltol=>1e-6, :solver => Tsit5())
tr = trajectory(ds, 50000.0, dt=dt, diff_eq_kwargs = diffeq)
subplot(3,2,5)
plot(tr[:,1], tr[:,3], alpha = 0.5,
label="ch",marker="o",markersize=2, linewidth=0)
legend()

subplot(3,2,6)
ls = lyapunovs(ds, 5000.0, dt=dt)
for k in [2,3,4]
    ex = sum(ls[1] - ls[j] for j in 2:k)
    g, t = gali(ds, k, 1000; dt = dt)
    semilogy(t, exp.(-ex.*t), label="exp. k=$k")
    semilogy(t, g, label="GALI_$(k)")
end
legend(fontsize=12)
tight_layout()

GALI Continuous As you can see, the results of both discrete and continuous match very well the theory described in gali.

Using GALI

No-one in their right mind would try to fit power-laws in order to distinguish between chaotic and regular behavior, like the above examples. These were just demonstrations.

The most common usage of $\text{GALI}_k$ is to define a (sufficiently) small amount of time and a (sufficiently) small threshold and see whether $\text{GALI}_k$ stays below it, for a (sufficiently) big $k$ .

For example one could do

using DynamicalSystems
using PyPlot

chaoticness(ds) = gali(ds, 2, 500.0)[2][end]

function main(k)

    dens = 201
    chaoticity = zeros(dens,dens)
    θs = ps = linspace(0, 2π, dens+1)
    ds = Systems.standardmap(k = k)

    for (i, θ) ∈ enumerate(θs[1:dens])
        for (j, p) ∈ enumerate(ps[1:dens])
            ds.state = SVector{2}(θ, p)
            chaoticity[i, j] = chaoticness(ds)
        end
    end

    pcolormesh(θs .- (θs[2] - θs[1])/2, ps .- (ps[2] - ps[1])/2,
    chaoticity')
    colorbar()

end

main(0.9)

and after about two minutes you will get: Chaos detection