Lyapunov Exponents

Lyapunov exponents measure rates of separation of nearby trajectories in the flow of a dynamical system. The Wikipedia and the Scholarpedia entries have a lot of valuable information about the history and usage of these quantities.

The naming comes after Aleksandr M. Lyapunov, a Russian mathematician/physicist that had major impact on the analysis of the stability of systems.

This page treats systems where the equations of motion are known. If instead you have numerical data, see the nonlinear timeseries analysis page.

Lyapunov Spectrum

The function lyapunovs calculates the entire spectrum of the Lyapunov exponents of a system:

#ChaosTools.lyapunovs — Function.

lyapunovs(ds::DynamicalSystem, N; kwargs...) -> [λ1, λ2, ..., λD]

Calculate the spectrum of Lyapunov exponents [1] of ds by applying a QR-decomposition on the parallelepiped matrix N times. Return the spectrum sorted from maximum to minimum.

Keyword Arguments

  • Ttr = 0 : Extra "transient" time to evolve the system before application of the algorithm. Should be Int for discrete systems.
  • dt = 1.0 : (only for continuous) Time of individual evolutions between successive orthonormalization steps.
  • diff_eq_kwargs = Dict() : (only for continuous) Keyword arguments passed into the solvers of the DifferentialEquations package (see trajectory for more info).

Description

The method we employ is "H2" of [2], originally stated in [3]. The vectors defining a D-dimensional parallepiped are evolved using the tangent dynamics (known also as variational equations) of the system. A QR-decomposition at each step yields the local growth rate for each dimension of the parallepiped. The growth rates are then averaged over N successive steps, yielding the lyapunov exponent spectrum.

For discrete systems the QR-decomposition is performed at every step i ∈ 1:N.

References

[1] : A. M. Lyapunov, The General Problem of the Stability of Motion, Taylor & Francis (1992)

[2] : K. Geist et al., Progr. Theor. Phys. 83, pp 875 (1990)

[3] : G. Benettin et al., Meccanica 15, pp 9-20 & 21-30 (1980)

source

As you can see, the documentation string is detailed and self-contained. For example, the lyapunov spectrum of the folded towel map is calculated as:

using DynamicalSystems

ds = Systems.towel()
λλ = lyapunovs(ds, 10000)

3-element SVector{3,Float64}:
  0.432253
  0.371617
 -3.29632

Similarly, for a continuous system, e.g. the Lorenz system, you would do:

lor = Systems.lorenz(ρ = 32.0) #this is not the original parameter!
λλ = lyapunovs(lor, 10000,
dt = 0.1, diff_eq_kwargs = Dict(:abstol => 1e-9, :reltol => 1e-9))

3-element Array{Float64,1}:
   0.985688
   0.00271333
 -14.6551

Maximum Lyapunov Exponent

The function lyapunov calculates the maximum lyapunov exponent of a system, much more efficiently than getting the first result of lyapunovs:

#ChaosTools.lyapunov — Function.

lyapunov(ds::DynamicalSystem, Τ; kwargs...)

Calculate the maximum Lyapunov exponent λ using a method due to Benettin [1], which simply evolves two neighboring trajectories (one called "given" and one called "test") while constantly rescaling the test one. T denotes the total time of evolution (should be Int for discrete systems).

Keyword Arguments

  • Ttr = 0 : Extra "transient" time to evolve the system before application of the algorithm. Should be Int for discrete systems.
  • d0 = 1e-9 : Initial & rescaling distance between the two neighboring trajectories.
  • threshold = 1e-5 : Distance threshold for rescaling.
  • diff_eq_kwargs = Dict(:abstol=>d0, :reltol=>d0) : (only for continuous) Keyword arguments passed into the solvers of the DifferentialEquations package (see trajectory for more info).
  • dt = 0.1 : (only for continuous) Time of evolution between each check of distance exceeding the threshold.
  • inittest = (st1, d0) -> st1 .+ d0/sqrt(D) : A function that given (st1, d0) initializes the test state with distance d0 from the given state st1 (D is the dimension of the system). This function can be used when you want to avoid the test state appearing in a region of the phase-space where it would have e.g. different energy or escape to infinity.

Description

Two neighboring trajectories with initial distance d0 are evolved in time. At time d(t_i) their distance exceeds the threshold, which initializes a rescaling of the test trajectory back to having distance d0 from the given one, while the rescaling keeps the distance vector along the maximal expansion direction.

The maximum Lyapunov exponent is the average of the time-local Lyapunov exponents

\lambda = \frac{1}{t_{n}}\sum_{i=1}^{n} \ln\left( a_i \right),\quad a_i = \frac{d(t_{i})}{d_0}.

Performance Notes

For the continuous case, the algorithm becomes faster with increasing dt, since integration is interrupted less frequenty. For the fastest performance you want to fine-tune dt, d0, threshold such that you have the minimum amount of rescalings while still being well within the linearized dynamics region.

You can easily modify the source code to return the convergence timeseries of the exponent, if need be.

References

[1] : G. Benettin et al., Phys. Rev. A 14, pp 2338 (1976)

source

For example:

using DynamicalSystems

henon = Systems.henon()
λ = lyapunov(henon, 10000, d0 = 1e-7, threshold = 1e-4, Ttr = 100)

0.42007471604734054

The same is done for continuous systems:

using DynamicalSystems, OrdinaryDiffEq

ross = Systems.roessler(a = 0.1, b = 0.1, c = 14.0) #not original parameters
λ = lyapunov(ross, 100000, dt = 10.0,
diff_eq_kwargs = Dict(:solver => Vern8(), :abstol=>1e-9, :reltol=>1e-9),
Ttr = 100.0)

0.07076483526799626