Detecting & Categorizing Chaos

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... : Keyword arguments propagated into init of DifferentialEquations.jl. See trajectory for examples. Only valid for continuous systems.

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 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)

predictability(ds::DynamicalSystem; kwargs...) -> chaos_type, ν, C

Determine whether ds displays strongly chaotic, partially-predictable chaotic or regular behaviour, using the method by Wernecke et al. described in^{[Wernecke2017]}.

Return the type of the behavior, the cross-distance scaling coefficient ν and the correlation coefficient C. Typical values for ν, C and chaos_type are given in Table 2 of^{[Wernecke2017]}:

Keyword Arguments

• Ttr = 200 : Extra "transient" time to evolve the system before sampling from the trajectory. Should be Int for discrete systems.
• T_sample = 1e4 : Time to evolve the system for taking samples. Should be Int for discrete systems.
• n_samples = 500 : Number of samples to take for use in calculating statistics.
• λ_max = lyapunov(ds, 5000) : Value to use for largest Lyapunov exponent for finding the Lyapunov prediction time. If it is less than zero a regular result is returned immediatelly.
• d_tol = 1e-3 : tolerance distance to use for calculating Lyapunov prediction time.
• T_multiplier = 10 : Multiplier from the Lyapunov prediction time to the evaluation time.
• T_max = Inf : Maximum time at which to evaluate trajectory distance. If the internally computed evaluation time is larger than T_max, stop at T_max instead.
• δ_range = 10.0 .^ (-9:-6) : Range of initial condition perturbation distances to use to determine scaling ν.
• diffeq... : Keyword arguments propagated into init of DifferentialEquations.jl. See trajectory for examples. Only valid for continuous systems.

Description

Samples points from a trajectory of the system to be used as initial conditions. Each of these initial conditions is randomly perturbed by a distance δ, and the trajectories for both the original and perturbed initial conditions are computed to the 'evaluation time' T.

The average (over the samples) distance and cross-correlation coefficient of the state at time T is computed. This is repeated for a range of δ (defined by δ_range), and linear regression is used to determine how the distance and cross-correlation scale with δ, allowing for identification of chaos type.

The evaluation time T is calculated as T = T_multiplier*Tλ, where the Lyapunov prediction time Tλ = log(d_tol/δ)/λ_max. This may be very large if the λ_max is small, e.g. when the system is regular, so this internally computed time T can be overridden by a smaller T_max set by the user.

Performance Notes

For continuous systems, it is likely that the maxiters used by the integrators needs to be increased, e.g. to 1e9. This is part of the diffeq kwargs. In addition, be aware that this function does a lot of internal computations. It is operating in a different speed than e.g. lyapunov.

testchaos01(φ::Vector [, cs, N0]) -> chaotic?

Perform the so called "0-1" test for chaos introduced by Gottwald and Melbourne^{[Gottwald2016]} on the timeseries φ. Return true if φ is chaotic, false otherwise.

Description

This method tests if the given timeseries is chaotic or not by transforming it into a two-dimensional diffusive process. If the timeseries is chaotic, the mean square displacement of the process grows as sqrt(length(φ)), while it stays constant if the timeseries is regular. The implementation here computes K, the correlation coefficient (median of Kc for c ∈ cs), and simply checks if K > 0.5.

If you want to access the various Kc you should call the method testchaos01(φ, c::Real, N0) which returns Kc.

cs defaults to 3π/5*rand(10) + π/4 and N0, the length of the two-dimensional process, is N0 = length(φ)/10.

Notice that for data sampled from continous dynamical systems, some care must be taken regarding the values of cs, see^{[Gottwald2016]}.

expansionentropy(ds::DynamicalSystem, sampler, restraining; kwargs...)

Calculate the expansion entropy^[Hunt2015] of ds, in the restraining region $S$ defined by restraining, by estimating the slope of the biggest linear region of the curve $\log E_{t0+T, t0}(f, S)$ versus $T$ (using linear_region). This is an approximation of the expansion entropy $H_0$, according to^[Hunt2015].

sampler is a 0-argument function that generates a random initial condition (a sample) of ds. restraining is a 1-argument function restraining(u) that given the state u it returns true if the state is inside the restraining region $S$.

Use boxregion for an easy way to define sampler and restraining on a multidimension box.

Keyword Arguments

• N = 1000 : Number of samples taken at each batch (same as $N$ of [1]).
• steps = 40 : The maximal steps for which the system will be run.
• Ttr = 0 : Transient time to evolve each initial condition before starting to comute $E$. This is t0 of [1] and of the following notation.
• batches = 100 : Number of batches to run the calculation, see below.
• Δt = 1 : Integration step size.
• diffeq... : Other keywords are propagated to the solvers of DifferentialEquations.jl.

Description

N samples are initialized and propagated forwards in time (along with their tangent space). At every time $t$ in [t0+Δt, t0+2dt, ... t0+steps*Δt] we calculate $H$:

\[H[t] = \log E_{t0+T, t0}(f, S),\]

with

\[E_{t0+T, t0}(f, S) = \frac 1 N \sum_{i'} G(Df_{t0+t, t0}(x_i))\]

(using same notation as ^[Hunt2015]). In principle $E$ is the average largest possible growth ratio within the restraining region (sampled by the initial conditions). The summation is only over $x_i$ that stay inside the region $S$ defined by the boolean function restraining. This process is done by the expansionentropy_sample function.

Then, this is repeated for batches amount of times, as recommended in^[Hunt2015]. From all these batches, the mean and std of $H$ is computed at every time point. This is done by the expansionentropy_batch function. When plotted versus $t$, these create the curves and error bars of e.g. Figs 2, 3 of [1].

This function expansionentropy simply returns the slope of the biggest linear region of the curve $H$ versus $t$, which approximates the expansion entropy $H_0$. It is therefore recommended to use expansionentropy_batch directly and evaluate the result yourself, as this step is known to be inaccurate for non-chaotic systems (where $H$ fluctuates strongly around 0).

boxregion(as, bs) -> sampler, restraining

Define a box in $\mathbb{R}^d$ with edges the as and bs and then return two functions: sampler, which generates a random initial condition in that box and restraining that returns true if a given state is in the box.

expansionentropy_sample(ds, sampler, restraining; kwargs...)

Return times, H for one sample of ds (see expansionentropy). Accepts the same argumets as expansionentropy, besides batches.

expansionentropy_batch(ds, sampler, restraining; kwargs...)

Run expansionentropy_sample batch times, and return times, mean(H), std(H) for all resulting H, see expansionentropy.

Accepts the same arguments as expansionentropy.