ChaosTools.jl
ChaosTools
— ModuleChaosTools.jl
A Julia module that offers various tools for analysing nonlinear dynamics and chaotic behaviour. It can be used as a standalone package, or as part of DynamicalSystems.jl.
To install it, run import Pkg; Pkg.add("ChaosTools")
.
All further information is provided in the documentation, which you can either find online or build locally by running the docs/make.jl
file.
ChaosTools.jl is the jack-of-all-trades package of the DynamicalSystems.jl library: methods that are not extensive enough to be a standalone package are added here. You should see the full DynamicalSystems.jl library for other packages that may contain functionality you are looking for but did not find in ChaosTools.jl.
A good background for understanding the methods of ChaosTools.jl is the following textbook: Nonlinear Dynamics, Datseris & Parlitz, Springer 2022.
DynamicalSystemsBase.jl reference
As many docstrings in ChaosTools.jl point to the different DynamicalSystem
types, they are also provided here for reference.
DynamicalSystem
DeterministicIteratedMap
CoupledODEs
CoreDynamicalSystem
StroboscopicMap
PoincareMap
TangentDynamicalSystem
ParallelDynamicalSystem
ProjectedDynamicalSystem
reinit!
DynamicalSystemsBase.DynamicalSystem
— TypeDynamicalSystem
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
The documentation of DynamicalSystem
follows chapter 1 of Nonlinear Dynamics, Datseris & Parlitz, Springer 2022.
A ds::DynamicalSystem
representes a flow Φ in a state space. It mainly encapsulates three things:
- A state, typically referred to as
u
, with initial valueu0
. The space thatu
occupies is the state space ofds
and the length ofu
is the dimension ofds
(and of the state space). - A dynamic rule, typically referred to as
f
, that dictates how the state evolves/changes with time when calling thestep!
function.f
is a standard Julia function, see below. - A parameter container
p
that parameterizesf
.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 formf(u, p, t) -> out
which means that given a stateu::SVector{<:Real}
and some parameter containerp
it returns the output off
as anSVector{<:Real}
(static vector). - iip :
f
must be in the formf!(out, u, p, t)
which means that given a stateu::AbstractArray{<:Real}
and some parameter containerp
, it writes in-place the output off
inout::AbstractArray{<:Real}
. The function must returnnothing
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
.
Whether the dynamical system is autonomous (f
doesn't depend on time) or not, it is still necessary to include t
as an argument to f
. Some algorithms utilize this information, some do not, but we prefer to keep a consistent interface either way. You can also convert any system to autonomous by making time an additional variable. If the system is non-autonomous, its effective dimensionality is dimension(ds)+1
.
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)
withds
an instance ofDynamicalSystem
: return the state ofds
at timet
. For continuous time systems this interpolates and extrapolates, while for discrete time systems it only works ift
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
DynamicalSystemsBase.DeterministicIteratedMap
— TypeDeterministicIteratedMap <: 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
.
DynamicalSystemsBase.CoupledODEs
— TypeCoupledODEs <: 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.
DynamicalSystemsBase.CoreDynamicalSystem
— TypeCoreDynamicalSystem
Union type meaning either DeterministicIteratedMap
or CoupledODEs
, which are the core systems whose dynamic rule f
is known analytically.
This type is used for deciding whether a creation of a TangentDynamicalSystem
is possible or not.
DynamicalSystemsBase.StroboscopicMap
— TypeStroboscopicMap <: 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
.
DynamicalSystemsBase.PoincareMap
— TypePoincareMap <: 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 withsign(direction)
are considered to belong to the surface of section. Negative direction means going from less than $b$ to greater than $b$.u0 = nothing
: Specify an initial state.rootkw = (xrtol = 1e-6, atol = 1e-8)
: ANamedTuple
of keyword arguments passed tofind_zero
from Roots.jl.Tmax = 1e3
: The argumentTmax
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 onestep!
the system has been evolved for more thanTmax
, thenstep!(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 thej
th variable of the system equals the valuer
. - A vector of length
D+1
. The firstD
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:
dimension(pmap) == dimension(ds)
, even though the Poincaré map is effectively 1 dimension less.- Like
StroboscopicMap
time is discrete and counts the iterations on the surface of section.initial_time
is always0
andcurrent_time
is current iteration number. - A new function
current_crossing_time
returns the real time corresponding to the latest crossing of the hyperplane, which is what thecurrent_state(ds)
corresponds to as well. - For the special case of
plane
being aTuple{Int, <:Real}
, a specialreinit!
method is allowed with input state of lengthD-1
instead ofD
, i.e., a reduced state already on the hyperplane that is then converted into theD
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)
DynamicalSystemsBase.TangentDynamicalSystem
— TypeTangentDynamicalSystem <: DynamicalSystem
TangentDynamicalSystem(ds::CoreDynamicalSystem; kwargs...)
A dynamical system that bundles the evolution of ds
(which must be an CoreDynamicalSystem
) and k
deviation vectors that are evolved according to the dynamics in the tangent space (also called linearized dynamics or the tangent dynamics).
The state of ds
must be an AbstractVector
for TangentDynamicalSystem
.
TangentDynamicalSystem
follows the DynamicalSystem
interface with the following adjustments:
reinit!
takes an additional keywordQ0
(with same default as below)- The additional functions
current_deviations
andset_deviations!
are provided for the deviation vectors.
Keyword arguments
k
orQ0
:Q0
represents the initial deviation vectors (each column = 1 vector). Ifk::Int
is given, a matrixQ0
is created with the firstk
columns of the identity matrix. OtherwiseQ0
can be given directly as a matrix. It must hold thatsize(Q, 1) == dimension(ds)
. You can useorthonormal
for random orthonormal vectors. By defaultk = dimension(ds)
is used.u0 = current_state(ds)
: Starting state.J
andJ0
: See section "Jacobian" below.
Description
Let $u$ be the state of ds
, and $y$ a deviation (or perturbation) vector. These two are evolved in parallel according to
\[\begin{array}{rcl} \frac{d\vec{x}}{dt} &=& f(\vec{x}) \\ \frac{dY}{dt} &=& J_f(\vec{x}) \cdot Y \end{array} \quad \mathrm{or}\quad \begin{array}{rcl} \vec{x}_{n+1} &=& f(\vec{x}_n) \\ Y_{n+1} &=& J_f(\vec{x}_n) \cdot Y_n. \end{array}\]
for continuous or discrete time respectively. Here $f$ is the dynamic_rule
(ds)
and $J_f$ is the Jacobian of $f$.
Jacobian
The keyword J
provides the Jacobian function. It must be a Julia function in the same form as f
, the dynamic_rule
. Specifically, J(u, p, n) -> M::SMatrix
for the out-of-place version or J(M, u, p, n)
for the in-place version acting in-place on M
. in both cases M
is a matrix whose columns are the deviation vectors.
By default J = nothing
. In this case J
is constructed automatically using the module ForwardDiff
, hence its limitations also apply here. Even though ForwardDiff
is very fast, depending on your exact system you might gain significant speed-up by providing a hand-coded Jacobian and so it is recommended. Additionally, automatic and in-place Jacobians cannot be time dependent.
The keyword J0
allows you to pass an initialized Jacobian matrix J0
. This is useful for large in-place systems where only a few components of the Jacobian change during the time evolution. J0
can be a sparse or any other matrix type. If not given, a matrix of zeros is used. J0
is ignored for out of place systems.
DynamicalSystemsBase.ParallelDynamicalSystem
— TypeParallelDynamicalSystem <: DynamicalSystem
ParallelDynamicalSystem(ds::DynamicalSystem, states::Vector{<:AbstractArray})
A struct that evolves several states
of a given dynamical system in parallel at exactly the same times. Useful when wanting to evolve several different trajectories of the same system while ensuring that they share parameters and time vector.
This struct follows the DynamicalSystem
interface with the following adjustments:
- The function
current_state
is called ascurrent_state(pds, i::Int = 1)
which returns thei
th state. Same forinitial_state
. - Similarly,
set_state!
obtains a third argumenti::Int = 1
to set thei
-th state. current_states
andinitial_states
can be used to get all parallel states.reinit!
takes in a vector of states (likestates
) foru
.
DynamicalSystemsBase.ProjectedDynamicalSystem
— TypeProjectedDynamicalSystem <: 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 = ProjectedDynamicalSystem(ds, projection, complete_state)
reinit!(pds, [0.2, 0.4])
step!(pds)
get_state(pds)
Case 2: custom projection to general functions of state.
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...
SciMLBase.reinit!
— Methodreinit!(ds::DynamicalSystem, u = initial_state(ds); kwargs...) → ds
Reset the status of ds
, so that it is as if it has be just initialized with initial state u
. Practically every function of the ecosystem that evolves ds
first calls this function on it. Besides the new initial state u
, you can also configure the keywords t0 = initial_time(ds)
and p = current_parameters(ds)
.
Note the default settings: the state and time are the initial, but the parameters are the current.
The special method reinit!(ds, ::Nothing; kwargs...)
is also available, which does nothing and leaves the system as is. This is so that downstream functions that call reinit!
can still be used without resetting the system but rather continuing from its exact current state.
- DatserisParlitz2022Datseris & Parlitz 2022, Nonlinear Dynamics: A Concise Introduction Interlaced with Code, Springer Nature, Undergrad. Lect. Notes In Physics