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)
end
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
)
end
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,
animate_attractors_continuation(
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,)
)
fig
