Attractors.jl

Attractors.jl

A Julia module for finding attractors of dynamical systems, their basins and their boundaries, fractal properties of the boundaries, as well as continuing attractors and their basins across parameters. It can be used as a standalone package, or as part of DynamicalSystems.jl.

To install it, run import Pkg; Pkg.add("Attractors").

All further information is provided in the documentation, which you can either find online or build locally by running the docs/make.jl file.

Previously, Attractors.jl was part of ChaosTools.jl

DynamicalSystem

DynamicalSystem is an abstract supertype encompassing all concrete implementations of what counts as a "dynamical system" in the DynamicalSystems.jl library.

All concrete implementations of DynamicalSystem can be iteratively evolved in time via the step! function. Hence, most library functions that evolve the system will mutate its current state and/or parameters. See the documentation online for implications this has on for parallelization.

DynamicalSystem is further separated into two abstract types: ContinuousTimeDynamicalSystem, DiscreteTimeDynamicalSystem. The simplest and most common concrete implementations of a DynamicalSystem are DeterministicIteratedMap or CoupledODEs.

Description

A ds::DynamicalSystem representes a flow Φ in a state space. It mainly encapsulates three things:

1. A state, typically referred to as u, with initial value u0. The space that u occupies is the state space of ds and the length of u is the dimension of ds (and of the state space).
2. A dynamic rule, typically referred to as f, that dictates how the state evolves/changes with time when calling the step! function. f is a standard Julia function, see below.
3. A parameter container p that parameterizes f. p can be anything, but in general it is recommended to be a type-stable mutable container.

In sort, any set of quantities that change in time can be considered a dynamical system, however the concrete subtypes of DynamicalSystem are much more specific in their scope. Concrete subtypes typically also contain more information than the above 3 items.

In this scope dynamical systems have a known dynamic rule f defined as a standard Julia function. Observed or measured data from a dynamical system are represented using StateSpaceSet and are finite. Such data are obtained from the trajectory function or from an experimental measurement of a dynamical system with an unknown dynamic rule.

Construction instructions on f and u

Most of the concrete implementations of DynamicalSystem, with the exception of ArbitrarySteppable, have two ways of implementing the dynamic rule f, and as a consequence the type of the state u. The distinction is done on whether f is defined as an in-place (iip) function or out-of-place (oop) function.

• oop : f must be in the form f(u, p, t) -> out which means that given a state u::SVector{<:Real} and some parameter container p it returns the output of f as an SVector{<:Real} (static vector).
• iip : f must be in the form f!(out, u, p, t) which means that given a state u::AbstractArray{<:Real} and some parameter container p, it writes in-place the output of f in out::AbstractArray{<:Real}. The function must return nothing as a final statement.

t stands for current time in both cases. iip is suggested for systems with high dimension and oop for small. The break-even point is between 10 to 100 dimensions but should be benchmarked on a case-by-case basis as it depends on the complexity of f.

API

The API that the interface of DynamicalSystem employs is the functions listed below. Once a concrete instance of a subtype of DynamicalSystem is obtained, it can quieried or altered with the following functions.

The main use of a concrete dynamical system instance is to provide it to downstream functions such as lyapunovspectrum from ChaosTools.jl or basins_of_attraction from Attractors.jl. A typical user will likely not utilize directly the following API, unless when developing new algorithm implementations that use dynamical systems.

API - information

• ds(t) with ds an instance of DynamicalSystem: return the state of ds at time t. For continuous time systems this interpolates and extrapolates, while for discrete time systems it only works if t is the current time.
• current_state
• initial_state
• current_parameters
• initial_parameters
• isdeterministic
• isdiscretetime
• dynamic_rule
• current_time
• initial_time
• isinplace
• succesful_step

API - alter status

• reinit!
• set_state!
• set_parameter!
• set_parameters!

DeterministicIteratedMap <: DynamicalSystem
DeterministicIteratedMap(f, u0, p = nothing; t0 = 0)

A deterministic discrete time dynamical system defined by an iterated map as follows:

\[\vec{u}_{n+1} = \vec{f}(\vec{u}_n, p, n)\]

An alias for DeterministicIteratedMap is DiscreteDynamicalSystem.

Optionally configure the parameter container p and initial time t0.

For construction instructions regarding f, u0 see DynamicalSystem.

CoupledODEs <: ContinuousTimeDynamicalSystem
CoupledODEs(f, u0 [, p]; diffeq, t0 = 0.0)

A deterministic continuous time dynamical system defined by a set of coupled ordinary differential equations as follows:

\[\frac{d\vec{u}}{dt} = \vec{f}(\vec{u}, p, t)\]

An alias for CoupledODE is ContinuousDynamicalSystem.

Optionally provide the parameter container p and initial time as keyword t0.

For construction instructions regarding f, u0 see DynamicalSystem.

DifferentialEquations.jl keyword arguments and interfacing

The ODEs are evolved via the solvers of DifferentialEquations.jl. When initializing a CoupledODEs, you can specify the solver that will integrate f in time, along with any other integration options, using the diffeq keyword. For example you could use diffeq = (abstol = 1e-9, reltol = 1e-9). If you want to specify a solver, do so by using the keyword alg, e.g.: diffeq = (alg = Tsit5(), reltol = 1e-6). This requires you to have been first using OrdinaryDiffEq to access the solvers. The default diffeq is:

(alg = Tsit5(stagelimiter! = triviallimiter!, steplimiter! = triviallimiter!, thread = static(false)), abstol = 1.0e-6, reltol = 1.0e-6)

diffeq keywords can also include callback for event handling , however the majority of downstream functions in DynamicalSystems.jl assume that f is differentiable.

The convenience constructor CoupledODEs(prob::ODEProblem, diffeq) and CoupledODEs(ds::CoupledODEs, diffeq) are also available.

Dev note: CoupledODEs is a light wrapper of ODEIntegrator from DifferentialEquations.jl. The integrator is available as the field integ, and the ODEProblem is integ.sol.prob. The convenience syntax ODEProblem(ds::CoupledODEs, tspan = (t0, Inf)) is available.

StroboscopicMap <: DiscreteTimeDynamicalSystem
StroboscopicMap(ds::CoupledODEs, period::Real) → smap
StroboscopicMap(period::Real, f, u0, p = nothing; kwargs...)

A discrete time dynamical system that produces iterations of a time-dependent (non-autonomous) CoupledODEs system exactly over a given period. The second signature first creates a CoupledODEs and then calls the first.

StroboscopicMap follows the DynamicalSystem interface. In addition, the function set_period!(smap, period) is provided, that sets the period of the system to a new value (as if it was a parameter). As this system is in discrete time, current_time and initial_time are integers. The initial time is always 0, because current_time counts elapsed periods. Call these functions on the parent of StroboscopicMap to obtain the corresponding continuous time. In contrast, reinit! expects t0 in continuous time.

The convenience constructor

StroboscopicMap(T::Real, f, u0, p = nothing; diffeq, t0 = 0) → smap

is also provided.

See also PoincareMap.

PoincareMap <: DiscreteTimeDynamicalSystem
PoincareMap(ds::CoupledODEs, plane; kwargs...) → pmap

A discrete time dynamical system that produces iterations over the Poincaré map^{[DatserisParlitz2022]} of the given continuous time ds. This map is defined as the sequence of points on the Poincaré surface of section, which is defined by the plane argument.

See also StroboscopicMap, poincaresos.

Keyword arguments

• direction = -1: Only crossings with sign(direction) are considered to belong to the surface of section. Positive direction means going from less than $b$ to greater than $b$.
• u0 = nothing: Specify an initial state.
• rootkw = (xrtol = 1e-6, atol = 1e-8): A NamedTuple of keyword arguments passed to find_zero from Roots.jl.
• Tmax = 1e3: The argument Tmax exists so that the integrator can terminate instead of being evolved for infinite time, to avoid cases where iteration would continue forever for ill-defined hyperplanes or for convergence to fixed points, where the trajectory would never cross again the hyperplane. If during one step! the system has been evolved for more than Tmax, then step!(pmap) will terminate and error.

Description

The Poincaré surface of section is defined as sequential transversal crossings a trajectory has with any arbitrary manifold, but here the manifold must be a hyperplane. PoincareMap iterates over the crossings of the section.

If the state of ds is $\mathbf{u} = (u_1, \ldots, u_D)$ then the equation defining a hyperplane is

\[a_1u_1 + \dots + a_Du_D = \mathbf{a}\cdot\mathbf{u}=b\]

where $\mathbf{a}, b$ are the parameters of the hyperplane.

In code, plane can be either:

• A Tuple{Int, <: Real}, like (j, r): the plane is defined as when the jth variable of the system equals the value r.
• A vector of length D+1. The first D elements of the vector correspond to $\mathbf{a}$ while the last element is $b$.

PoincareMap uses ds, higher order interpolation from DifferentialEquations.jl, and root finding from Roots.jl, to create a high accuracy estimate of the section.

PoincareMap follows the DynamicalSystem interface with the following adjustments:

1. dimension(pmap) == dimension(ds), even though the Poincaré map is effectively 1 dimension less.
2. Like StroboscopicMap time is discrete and counts the iterations on the surface of section. initial_time is always 0 and current_time is current iteration number.
3. A new function current_crossing_time returns the real time corresponding to the latest crossing of the hyperplane, which is what the current_state(ds) corresponds to as well.
4. For the special case of plane being a Tuple{Int, <:Real}, a special reinit! method is allowed with input state of length D-1 instead of D, i.e., a reduced state already on the hyperplane that is then converted into the D dimensional state.

Example

using DynamicalSystemsBase
ds = Systems.rikitake(zeros(3); μ = 0.47, α = 1.0)
pmap = poincaremap(ds, (3, 0.0))
step!(pmap)
next_state_on_psos = current_state(pmap)

ProjectedDynamicalSystem <: DynamicalSystem
ProjectedDynamicalSystem(ds::DynamicalSystem, projection, complete_state)

A dynamical system that represents a projection of an existing ds on a (projected) space.

The projection defines the projected space. If projection isa AbstractVector{Int}, then the projected space is simply the variable indices that projection contains. Otherwise, projection can be an arbitrary function that given the state of the original system ds, returns the state in the projected space. In this case the projected space can be equal, or even higher-dimensional, than the original.

complete_state produces the state for the original system from the projected state. complete_state can always be a function that given the projected state returns a state in the original space. However, if projection isa AbstractVector{Int}, then complete_state can also be a vector that contains the values of the remaining variables of the system, i.e., those not contained in the projected space. In this case the projected space needs to be lower-dimensional than the original.

Notice that ProjectedDynamicalSystem does not require an invertible projection, complete_state is only used during reinit!. ProjectedDynamicalSystem is in fact a rather trivial wrapper of ds which steps it as normal in the original state space and only projects as a last step, e.g., during current_state.

Examples

Case 1: project 5-dimensional system to its last two dimensions.

ds = Systems.lorenz96(5)
projection = [4, 5]
complete_state = [0.0, 0.0, 0.0] # completed state just in the plane of last two dimensions
pds = projected_integrator(ds, projection, complete_state)
reinit!(pds, [0.2, 0.4])
step!(pds)
get_state(pds)

Case 2: custom projection to general functions of state. julia ds = Systems.lorenz96(5) projection(u) = [sum(u), sqrt(u[1]^2 + u[2]^2)] complete_state(y) = repeat(y[1]/5, 5) pds = # same as in above example...`

ArbitrarySteppable <: DiscreteTimeDynamicalSystem
ArbitrarySteppable(
    model, step!, extract_state, extract_parameters, reset_model!;
    isdeterministic = true, set_state = reinit!,
)

A dynamical system generated by an arbitrary "model" that can be stepped in-place with some function step!(model) for 1 step. The state of the model is extracted by the extract_state(model) -> u function The parameters of the model are extracted by the extract_parameters(model) -> p function. The system may be re-initialized, via reinit!, with the reset_model! user-provided function that must have the call signature

reset_model!(model, u, p)

given a (potentially new) state u and parameter container p, both of which will default to the initial ones in the reinit! call.

ArbitrarySteppable exists to provide the DynamicalSystems.jl interface to models from other packages that could be used within the DynamicalSystems.jl library. ArbitrarySteppable follows the DynamicalSystem interface with the following adjustments:

• initial_time is always 0, as time counts the steps the model has taken since creation or last reinit! call.
• set_state! is the same as reinit! by default. If not, the keyword argument set_state is a function set_state(model, u) that sets the state of the model to u.
• The keyword isdeterministic should be set properly, as it decides whether downstream algorithms should error or not.

StateSpaceSet{D, T} <: AbstractStateSpaceSet{D,T}

A dedicated interface for sets in a state space. It is an ordered container of equally-sized points of length D. Each point is represented by SVector{D, T}. The data are a standard Julia Vector{SVector}, and can be obtained with vec(ssset::StateSpaceSet). Typically the order of points in the set is the time direction, but it doesn't have to be.

When indexed with 1 index, StateSpaceSet is like a vector of points. When indexed with 2 indices it behaves like a matrix that has each of the columns be the timeseries of each of the variables. When iterated over, it iterates over its contained points. See description of indexing below for more.

StateSpaceSet also supports almost all sensible vector operations like append!, push!, hcat, eachrow, among others.

Description of indexing

In the following let i, j be integers, typeof(X) <: AbstractStateSpaceSet and v1, v2 be <: AbstractVector{Int} (v1, v2 could also be ranges, and for performance benefits make v2 an SVector{Int}).

• X[i] == X[i, :] gives the ith point (returns an SVector)
• X[v1] == X[v1, :], returns a StateSpaceSet with the points in those indices.
• X[:, j] gives the jth variable timeseries (or collection), as Vector
• X[v1, v2], X[:, v2] returns a StateSpaceSet with the appropriate entries (first indices being "time"/point index, while second being variables)
• X[i, j] value of the jth variable, at the ith timepoint

Use Matrix(ssset) or StateSpaceSet(matrix) to convert. It is assumed that each column of the matrix is one variable. If you have various timeseries vectors x, y, z, ... pass them like StateSpaceSet(x, y, z, ...). You can use columns(dataset) to obtain the reverse, i.e. all columns of the dataset in a tuple.