# ChaosTools.jl

`ChaosTools`

— Module**ChaosTools.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`

— Type`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 ** 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`

can be iteratively evolved in time via the `step!`

function.`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 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). - 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. - 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`

.

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)`

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**

`DynamicalSystemsBase.DeterministicIteratedMap`

— Type```
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`

.

`DynamicalSystemsBase.CoupledODEs`

— Type```
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(; stage*limiter! = trivial*limiter!, step*limiter! = trivial*limiter!, 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`

— Type`CoreDynamicalSystem`

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`

— Type```
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`

.

`DynamicalSystemsBase.PoincareMap`

— Type```
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. 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)`

: 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`j`

th 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:

`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 always`0`

and`current_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 the`current_state(ds)`

corresponds to as well. - 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)
```

`DynamicalSystemsBase.TangentDynamicalSystem`

— Type```
TangentDynamicalSystem <: 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 keyword`Q0`

(with same default as below)- The additional functions
`current_deviations`

and`set_deviations!`

are provided for the deviation vectors.

**Keyword arguments**

`k`

or`Q0`

:`Q0`

represents the initial deviation vectors (each column = 1 vector). If`k::Int`

is given, a matrix`Q0`

is created with the first`k`

columns of the identity matrix. Otherwise`Q0`

can be given directly as a matrix. It must hold that`size(Q, 1) == dimension(ds)`

. You can use`orthonormal`

for random orthonormal vectors. By default`k = dimension(ds)`

is used.`u0 = current_state(ds)`

: Starting state.`J`

and`J0`

: 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`

— Type```
ParallelDynamicalSystem <: 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 as`current_state(pds, i::Int = 1)`

which returns the`i`

th state. Same for`initial_state`

. - Similarly,
`set_state!`

obtains a third argument`i::Int = 1`

to set the`i`

-th state. `current_states`

and`initial_states`

can be used to get all parallel states.`reinit!`

takes in a vector of states (like`states`

) for`u`

.

`DynamicalSystemsBase.ProjectedDynamicalSystem`

— Type```
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 = 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!`

— Method`reinit!(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