[ANN] Attractors.jl

Hello people,

I am super proud to announce Attractors.jl and the publication that accompanies the bulk of the package, which just appeared on arXiv. As eluded to in the DynamicalSystems.jl v3 announcement, Attractors.jl is a submodule of the DynamicalSystems.jl library that offers:

  • finding attractors of arbitrary dynamical systems
  • finding their basins of attraction or the state space fractions of the basins
  • “continuing” the attractors and their basins over a parameter range
  • finding the basin boundaries and analyzing their fractal properties
  • tipping points related functionality for systems with known dynamic rule
  • and more!

Arguably the most novel part is the third bullet point, that offers a completely new approach to ``continuation’', which means trying to figure out how the stability of the attractors of a dynamical system changes as a system parameter changes. We (the authors @Datseris @KalelR @awage ) believe that our new approach to continuation will become an essential tool in the analysis of stability of dynamical systems.

As an example of the capabilities of the framework, I am attaching here the code snippet that we use in our publication to illustrate framework usage:

using DynamicalSystems # includes Attractors.jl
using OrdinaryDiffEq   # high-accuracy ODE solvers

# initialize Lorenz84 as a `DynamicalSystem`
function lorenz84_rule(u, p, t)
    F, G, a, b = p
    x, y, z = u
    dx = -y^2 -z^2 -a*x + a*F
    dy = x*y - y - b*x*z + G
    dz = b*x*y + x*z - z
    return SVector(dx, dy, dz)
p0 = [6.886, 1.347, 0.255, 4.0]
diffeq = (alg = Vern9(), reltol = 1e-9, abstol = 1e-9)
ds = CoupledODEs(lorenz84_rule, ones(3), p0; diffeq)

# Provide state space box to search in
xg = yg = zg = range(-3, 3; length = 600)
grid = (xg, yg, zg)
# initialize recurrences-based algorithm
# and choose its metaparameters
mapper = AttractorsViaRecurrences(ds, grid;
    mx_chk_fnd_att = 1000, mx_chk_loc_att = 2000,
    mx_chk_lost = 100, mx_chk_safety = 1e8,
    Dt = 0.05, force_non_adaptive = true,

# find and continue attractors across a given
# parameter range for the `pidx`-th parameter
prange = range(1.34, 1.37; length = 101)
pidx = 2 # index of parameter
sampler = statespace_sampler(grid)[1]
rsc = RecurrencesFindAndMatch(mapper)

fractions_curves, attractors_info = continuation(
    rsc, prange, pidx, sampler

# Estimate Lyapunov spectra for all attractors
# by looping over the parameter range
lyapunovs_curves = map(eachindex(prange)) do index
    set_parameter!(ds, pidx, prange[index])
    attractor_dict = attractors_info[index]
    exponents = Dict(
        k => lyapunovspectrum(ds, 10000; u0 = A[1])
        for (k, A) in attractor_dict

In these 50 lines of code an incredible amount of information and automation is achieved. A dynamical system is defined from scratch. Then, the system attractors are found with the mapper = AttractorsViaRecurences(...) using no prior knowledge about the system besides a state space box that may contain the attractors. This mapper is used to construct a special “continuation” method rsc = RecurrecesFindAndMatch(...), that, when given to continuation, it provides how the attractors of the system, and their relative fractions of attraction, behave versus a parameter range. At the very end of the code snippet we highlight the importance of such tools being part of a large and coherent framework (here DynamicalSystems.jl): the same input that is used to find attractors (ds = CoupledODEs(...)) is used to compute the Lyapunov exponents of the attractors.

We can animate the central output of continuation using some of the pre-defined plotting functions offered by the framework. For example,

    ds, attractors_info, fractions_curves, prange, pidx;
    savename = "lorenz84.mp4", access = [1,3],
    limits = (-1,3,-2,2),
    markersize = 10,

or to make a nice plot of the basin fractions (global stability measure) and Lyapunov exponents:

fig = basins_curves_plot(fractions_curves, prange;
    axislegend_kwargs = (position = :lb,)
axl = Axis(fig[0,1])
axl.ylabel = "λ₁ + λ₂"
hidexdecorations!(axl; grid = false)

lyap_to_real(L) = L[1]+L[2]
attractors_curves_plot!(axl, lyapunovs_curves, lyap_to_real, prange;
    axislegend_kwargs = (position = :cb,)