# DynamicalSystemsBase.jl

`DynamicalSystemsBase`

— Module**DynamicalSystemsBase.jl**

A Julia package that defines the `DynamicalSystem`

interface and many concrete implementations used in the DynamicalSystems.jl ecosystem.

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

. Typically, you do not want to use `DynamicalSystemsBase`

directly, as downstream analysis packages re-export it.

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

file.

!!! note "Tutorial and examples at DynamicalSystems.jl docs! Please visit the documentation of the main DynamicalSystems.jl docs for a tutorial and examples on using the interface.

## The `DynamicalSystem`

API

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

A `DynamicalSystem`

**represents the time evolution of a state 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 typically a standard Julia function, see the online documentation for examples. - 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`

. Finite *measured* or *sampled* data from a dynamical system are represented using `StateSpaceSet`

. Such data are obtained from the `trajectory`

function or from an experimental measurement of a dynamical system with an unknown dynamic rule.

See also the DynamicalSystems.jl tutorial online for examples making dynamical systems.

**Integration with ModelingToolkit.jl**

Dynamical systems that have been constructed from `DEProblem`

s that themselves have been constructed from ModelingToolkit.jl keep a reference to the symbolic model and all symbolic variables. Accessing a `DynamicalSystem`

using symbolic variables is possible via the functions `observe_state`

, `set_state!`

, `current_parameter`

and `set_parameter!`

. The referenced MTK model corresponding to the dynamical system can be obtained with `model = referrenced_sciml_model(ds::DynamicalSystem)`

.

See also the DynamicalSystems.jl tutorial online for an example.

In ModelingToolkit.jl v9 the default `split`

behavior of the parameter container is `true`

. This means that the parameter container is no longer a `Vector{Float64}`

by default, which means that you cannot use integers to access parameters. It is recommended to keep `split = true`

(default) and only access parameters via their symbolic parameter binding. Use `structural_simplify(sys; split = false)`

to allow accessing parameters with integers again.

**API**

The API that `DynamicalSystem`

employs is composed of the functions listed below. Once a concrete instance of a subtype of `DynamicalSystem`

is obtained, it can queried 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 - obtain 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`

`observe_state`

`current_parameters`

`current_parameter`

`initial_parameters`

`isdeterministic`

`isdiscretetime`

`dynamic_rule`

`current_time`

`initial_time`

`isinplace`

`successful_step`

`referrenced_sciml_model`

**API - alter status**

`DynamicalSystemsBase.current_state`

— Function`current_state(ds::DynamicalSystem) → u::AbstractArray`

Return the current state of `ds`

. This state is mutated when `ds`

is mutated. See also `initial_state`

, `observe_state`

.

`DynamicalSystemsBase.initial_state`

— Function`initial_state(ds::DynamicalSystem) → u0`

Return the initial state of `ds`

. This state is never mutated and is set when initializing `ds`

.

`DynamicalSystemsBase.observe_state`

— Function`observe_state(ds::DynamicalSystem, i, u = current_state(ds)) → x::Real`

Return the state `u`

of `ds`

*observed* at "index" `i`

. Possibilities are:

`i::Int`

returns the`i`

-th dynamic variable.`i::Function`

returns`f(current_state(ds))`

.`i::SymbolLike`

returns the value of the corresponding symbolic variable. This is valid only for dynamical systems referrencing a ModelingToolkit.jl model which also has`i`

as one of its listed variables (either uknowns or observed). Here`i`

can be anything can be anything that could index the solution object`sol = ModelingToolkit.solve(...)`

, such as a`Num`

or`Symbol`

instance with the name of the symbolic variable. In this case, a last fourth optional positional argument`t`

defaults to`current_time(ds)`

and is the time to observe the state at.- Any symbolic expression involving variables present in the symbolic variables tracked by the system, e.g.,
`i = x^2 - y`

with`x, y`

symbolic variables.

For `ProjectedDynamicalSystem`

, this function assumes that the state of the system is the full state space state, not the projected one (this makes the most sense for allowing MTK-based indexing).

Use `state_name`

for an accompanying name.

`DynamicalSystemsBase.state_name`

— Function`state_name(index)::String`

Return a name that matches the outcome of `observe_state`

with `index`

.

`DynamicalSystemsBase.current_parameters`

— Function`current_parameters(ds::DynamicalSystem) → p`

Return the current parameter container of `ds`

. This is mutated in functions that need to evolve `ds`

across a parameter range.

See also `initial_parameters`

, `current_parameter`

, `set_parameter!`

.

`DynamicalSystemsBase.current_parameter`

— Function`current_parameter(ds::DynamicalSystem, index [,p])`

Return the specific parameter of `ds`

corresponding to `index`

, which can be anything given to `set_parameter!`

. `p`

defaults to `current_parameters`

and is the parameter container to extract the parameter from, which must match layout with its default value.

Use `parameter_name`

for an accompanying name.

`DynamicalSystemsBase.parameter_name`

— Function`parameter_name(index)::String`

Return a name that matches the outcome of `current_parameter`

with `index`

.

`DynamicalSystemsBase.initial_parameters`

— Function`initial_parameters(ds::DynamicalSystem) → p0`

Return the initial parameter container of `ds`

. This is never mutated and is set when initializing `ds`

.

`DynamicalSystemsBase.isdeterministic`

— Function`isdeterministic(ds::DynamicalSystem) → true/false`

Return `true`

if `ds`

is deterministic, i.e., the dynamic rule contains no randomness. This is information deduced from the type of `ds`

.

`DynamicalSystemsBase.isdiscretetime`

— Function`isdiscretetime(ds::DynamicalSystem) → true/false`

Return `true`

if `ds`

operates in discrete time, or `false`

if it is in continuous time. This is information deduced from the type of `ds`

.

`DynamicalSystemsBase.dynamic_rule`

— Function`dynamic_rule(ds::DynamicalSystem) → f`

Return the dynamic rule of `ds`

. This is never mutated and is set when initializing `ds`

.

`DynamicalSystemsBase.current_time`

— Function`current_time(ds::DynamicalSystem) → t`

Return the current time that `ds`

is at. This is mutated when `ds`

is evolved.

`DynamicalSystemsBase.initial_time`

— Function`initial_time(ds::DynamicalSystem) → t0`

Return the initial time defined for `ds`

. This is never mutated and is set when initializing `ds`

.

`SciMLBase.isinplace`

— Method`isinplace(ds::DynamicalSystem) → true/false`

Return `true`

if the dynamic rule of `ds`

is in-place, i.e., a function mutating the state in place. If `true`

, the state is typically `Array`

, if `false`

, the state is typically `SVector`

. A front-end user will most likely not care about this information, but a developer may care.

`DynamicalSystemsBase.successful_step`

— Function`successful_step(ds::DynamicalSystem) -> true/false`

Return `true`

if the last `step!`

call to `ds`

was successful, `false`

otherwise. For continuous time systems this uses DifferentialEquations.jl error checking, for discrete time it checks if any variable is `Inf`

or `NaN`

.

`DynamicalSystemsBase.referrenced_sciml_model`

— Function`referrenced_sciml_model(ds::DynamicalSystem)`

Return the ModelingToolkit.jl structurally-simplified model referrenced by `ds`

. Return `nothing`

if there is no referrenced model.

`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 state `u`

, you can also configure the keywords `t0 = initial_time(ds)`

and `p = current_parameters(ds)`

.

`reinit!(ds::DynamicalSystem, u::AbstractDict; kwargs...) → ds`

If `u`

is a `AbstractDict`

(for partially setting specific state variables in `set_state!`

), then the alterations are done in the state given by the keyword `reference_state = copy(initial_state(ds))`

.

`reinit!(ds, ::Nothing; kwargs...)`

This method 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.

`DynamicalSystemsBase.set_state!`

— Function`set_state!(ds::DynamicalSystem, u::AbstractArray{<:Real})`

Set the state of `ds`

to `u`

, which must match dimensionality with that of `ds`

. Also ensure that the change is notified to whatever integration protocol is used.

`set_state!(ds::DynamicalSystem, value::Real, i) → u`

Set the `i`

th variable of `ds`

to `value`

. The index `i`

can be an integer or a symbolic-like index for systems that reference a ModelingToolkit.jl model. For example:

```
i = :x # or `1` or `only(@variables(x))`
set_state!(ds, 0.5, i)
```

**Warning:** this function should not be used with derivative dynamical systems such as Poincare/stroboscopic/projected dynamical systems. Use the method below to manipulate an array and give that to `set_state!`

.

`set_state!(u::AbstractArray, value, index, ds::DynamicalSystem)`

Modify the given state `u`

and leave `ds`

untouched.

`set_state!(ds::DynamicalSystem, mapping::AbstractDict)`

Convenience version of `set_state!`

that iteratively calls `set_state!(ds, val, i)`

for all index-value pairs `(i, val)`

in `mapping`

. This is useful primarily in two cases:

- to partially set only some state variables,
- to set variables by name (if the system is created via ModelingToolkit.jl)

so that you don't have to keep track of the order of the dynamic variables.

`DynamicalSystemsBase.set_parameter!`

— Function`set_parameter!(ds::DynamicalSystem, index, value [, p])`

Change a parameter of `ds`

given the `index`

it has in the parameter container and the `value`

to set it to. This function works for any type of parameter container (array/dictionary/composite types) provided the `index`

is appropriate type.

The `index`

can be a traditional Julia index (integer for arrays, key for dictionaries, or symbol for composite types). It can also be a symbolic variable or `Symbol`

instance. This is valid only for dynamical systems referring a ModelingToolkit.jl model which also has `index`

as one of its parameters.

The last optional argument `p`

defaults to `current_parameters`

and is the parameter container whose value is changed at the given index. It must match layout with its default value.

`DynamicalSystemsBase.set_parameters!`

— Function`set_parameters!(ds::DynamicalSystem, p = initial_parameters(ds))`

Set the parameter values in the `current_parameters`

`(ds)`

to match those in `p`

. This is done as an in-place overwrite by looping over the keys of `p`

hence `p`

can be an arbitrary container mapping parameter indices to values (such as a `Vector{Real}`

, `Vector{Pair}`

, or `AbstractDict`

).

The keys of `p`

must be valid keys that can be given to `set_parameter!`

.

## Time evolution

`CommonSolve.step!`

— Method`step!(ds::DiscreteTimeDynamicalSystem [, n::Integer]) → ds`

Evolve the discrete time dynamical system for 1 or `n`

steps.

`step!(ds::ContinuousTimeDynamicalSystem, [, dt::Real [, stop_at_tdt]]) → ds`

Evolve the continuous time dynamical system for one integration step.

Alternatively, if a `dt`

is given, then progress the integration until there is a temporal difference `≥ dt`

(so, step *at least* for `dt`

time).

When `true`

is passed to the optional third argument, the integration advances for exactly `dt`

time.

`DynamicalSystemsBase.trajectory`

— Function`trajectory(ds::DynamicalSystem, T [, u0]; kwargs...) → X, t`

Evolve `ds`

for a total time of `T`

and return its trajectory `X`

, sampled at equal time intervals, and corresponding time vector. `X`

is a `StateSpaceSet`

. Optionally provide a starting state `u0`

which is `current_state(ds)`

by default.

The returned time vector is `t = (t0+Ttr):Δt:(t0+Ttr+T)`

.

If time evolution diverged, or in general failed, before `T`

, the remaining of the trajectory is set to the last valid point.

`trajectory`

is a very simple function provided for convenience. For continuous time systems, it doesn't play well with callbacks, use `DifferentialEquations.solve`

if you want a trajectory/timeseries that works with callbacks, or in general you want more flexibility in the generated trajectory (but remember to convert the output of `solve`

to a `StateSpaceSet`

).

**Keyword arguments**

`Δt`

: Time step of value output. For discrete time systems it must be an integer. Defaults to`0.1`

for continuous and`1`

for discrete time systems. If you don't have access to unicode, the keyword`Dt`

can be used instead.`Ttr = 0`

: Transient time to evolve the initial state before starting saving states.`t0 = initial_time(ds)`

: Starting time.`container = SVector`

: Type of vector that will represent the state space points that will be included in the`StateSpaceSet`

output. See`StateSpaceSet`

for valid options.`save_idxs::AbstractVector`

: Which variables to output in`X`

. It can be any type of index that can be given to`observe_state`

. Defaults to`1:dimension(ds)`

(all dynamic variables). Note: if you mix integer and symbolic indexing be sure to initialize the array as`Any`

so that integers`1, 2, ...`

are not converted to symbolic expressions.

`StateSpaceSets.StateSpaceSet`

— Type`StateSpaceSet{D, T, V} <: AbstractVector{V}`

A dedicated interface for sets in a state space. It is an **ordered container of equally-sized points** of length `D`

, with element type `T`

, represented by a vector of type `V`

. Typically `V`

is `SVector{D,T}`

or `Vector{T}`

and the data are always stored internally as `Vector{V}`

. `SSSet`

is an alias for `StateSpaceSet`

.

The underlying `Vector{V}`

can be obtained by `vec(ssset)`

, although this is almost never necessary because `StateSpaceSet`

subtypes `AbstractVector`

and extends its interface. `StateSpaceSet`

also supports almost all sensible vector operations like `append!, push!, hcat, eachrow`

, among others. When iterated over, it iterates over its contained points.

**Construction**

Constructing a `StateSpaceSet`

is done in three ways:

- By giving in each individual
**columns**of the state space set as`Vector{<:Real}`

:`StateSpaceSet(x, y, z, ...)`

. - By giving in a matrix whose rows are the state space points:
`StateSpaceSet(m)`

. - By giving in directly a vector of vectors (state space points):
`StateSpaceSet(v_of_v)`

.

All constructors allow for the keyword `container`

which sets the type of `V`

(the type of inner vectors). At the moment options are only `SVector`

, `MVector`

, or `Vector`

, and by default `SVector`

is used.

**Description of indexing**

When indexed with 1 index, `StateSpaceSet`

behaves exactly like its encapsulated vector. i.e., a vector of vectors (state space points). When indexed with 2 indices it behaves like a matrix where each row is a point.

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

th point (returns an`SVector`

)`X[v1] == X[v1, :]`

, returns a`StateSpaceSet`

with the points in those indices.`X[:, j]`

gives the`j`

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

th variable, at the`i`

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

`DeterministicIteratedMap`

`DynamicalSystemsBase.DeterministicIteratedMap`

— Type```
DeterministicIteratedMap <: DiscreteTimeDynamicalSystem
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 the DynamicalSystems.jl tutorial.

`CoupledODEs`

`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 the DynamicalSystems.jl tutorial.

**DifferentialEquations.jl 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`

(or smaller library package such as `OrdinaryDiffEqVerner`

) to access the solvers. The default `diffeq`

is:

(alg = OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial*limiter!), typeof(OrdinaryDiffEqCore.trivial*limiter!), Static.False}(OrdinaryDiffEqCore.trivial*limiter!, OrdinaryDiffEqCore.trivial*limiter!, static(false)), abstol = 1.0e-6, reltol = 1.0e-6)

`diffeq`

keywords can also include `callback`

for event handling .

The convenience constructors `CoupledODEs(prob::ODEProblem [, diffeq])`

and `CoupledODEs(ds::CoupledODEs [, diffeq])`

are also available. Use `ODEProblem(ds::CoupledODEs, tspan = (t0, Inf))`

to obtain the problem.

To integrate with ModelingToolkit.jl, the dynamical system **must** be created via the `ODEProblem`

(which itself is created via ModelingToolkit.jl), see the Tutorial for an example.

Dev note: `CoupledODEs`

is a light wrapper of `ODEIntegrator`

from DifferentialEquations.jl.

`StroboscopicMap`

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

.

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

Iterating `pmap`

also mutates `ds`

which is referrenced in `pmap`

.

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. The corresponding state on the hyperplane is`current_state(pmap)`

as expected. - 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.current_crossing_time`

— Function`current_crossing_time(pmap::PoincareMap) → tcross`

Return the time of the latest crossing of the Poincare section.

`DynamicalSystemsBase.poincaresos`

— Function`poincaresos(A::AbstractStateSpaceSet, plane; kwargs...) → P::StateSpaceSet`

Calculate the Poincaré surface of section of the given dataset with the given `plane`

by performing linear interpolation betweeen points that sandwich the hyperplane.

Argument `plane`

and keywords `direction, warning, save_idxs`

are the same as in `PoincareMap`

.

`poincaresos(ds::CoupledODEs, plane, T = 1000.0; kwargs...) → P::StateSpaceSet`

Return the iterations of `ds`

on the Poincaré surface of section with the `plane`

, by evolving `ds`

up to a total of `T`

. Return a `StateSpaceSet`

of the points that are on the surface of section.

This function initializes a `PoincareMap`

and steps it until its `current_crossing_time`

exceeds `T`

. You can also use `trajectory`

with `PoincareMap`

to get a sequence of `N::Int`

points instead.

The keywords `Ttr, save_idxs`

act as in `trajectory`

. See `PoincareMap`

for `plane`

and all other keywords.

`TangentDynamicalSystem`

`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.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 the Jacobian matrix used for the evolution of 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.current_deviations`

— Function`current_deviations(tands::TangentDynamicalSystem)`

Return the deviation vectors of `tands`

as a matrix with each column a vector.

`DynamicalSystemsBase.set_deviations!`

— Function`set_deviations!(tands::TangentDynamicalSystem, Q)`

Set the deviation vectors of `tands`

to be `Q`

, a matrix with each column a vector.

`DynamicalSystemsBase.jacobian`

— Function`jacobian(ds::CoreDynamicalSystem)`

Construct the Jacobian rule for the dynamical system `ds`

. This is done via automatic differentiation using module `ForwardDiff`

.

**Description**

For out-of-place systems, `jacobian`

returns the Jacobian rule as a function `Jf(u, p, t) -> J0::SMatrix`

. Calling `Jf(u, p, t)`

will compute the Jacobian at the state `u`

, parameters `p`

and time `t`

and return the result as `J0`

. For in-place systems, `jacobian`

returns the Jacobian rule as a function `Jf!(J0, u, p, t)`

. Calling `Jf!(J0, u, p, t)`

will compute the Jacobian at the state `u`

, parameters `p`

and time `t`

and save the result in `J0`

.

`StateSpaceSets.orthonormal`

— Function`orthonormal([T,] D, k) -> ws`

Return a matrix `ws`

with `k`

columns, each being an `D`

-dimensional orthonormal vector.

`T`

is the return type and can be either `SMatrix`

or `Matrix`

. If not given, it is `SMatrix`

if `D*k < 100`

, otherwise `Matrix`

.

`ProjectedDynamicalSystem`

`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
prods = ProjectedDynamicalSystem(ds, projection, complete_state)
reinit!(prods, [0.2, 0.4])
step!(prods)
current_state(prods)
```

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)
prods = # same as in above example...
```

`ParallelDynamicalSystem`

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

.

`ParallelDynamicalSystem(ds::DynamicalSystem, states::Vector{<:Dict})`

For a dynamical system referring a MTK model, one can specify states as a vector of dictionaries to alter the current state of `ds`

as in `set_state!`

.

`DynamicalSystemsBase.initial_states`

— Function`initial_states(pds::ParallelDynamicalSystem)`

Return an iterator over the initial parallel states of `pds`

.

`DynamicalSystemsBase.current_states`

— Function`current_states(pds::ParallelDynamicalSystem)`

Return an iterator over the parallel states of `pds`

.

`ArbitrarySteppable`

`DynamicalSystemsBase.ArbitrarySteppable`

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

## Parallelization

Since `DynamicalSystem`

s are mutable, one needs to copy them before parallelizing, to avoid having to deal with complicated race conditions etc. The simplest way is with `deepcopy`

. Here is an example block that shows how to parallelize calling some expensive function (e.g., calculating the Lyapunov exponent) over a parameter range using `Threads`

:

```
ds = DynamicalSystem(f, u, p) # some concrete implementation
parameters = 0:0.01:1
outputs = zeros(length(parameters))
# Since `DynamicalSystem`s are mutable, we need to copy to parallelize
systems = [deepcopy(ds) for _ in 1:Threads.nthreads()-1]
pushfirst!(systems, ds) # we can save 1 copy
Threads.@threads for i in eachindex(parameters)
system = systems[Threads.threadid()]
set_parameter!(system, 1, parameters[i])
outputs[i] = expensive_function(system, args...)
end
```

## Advanced example

This is an advanced example of making an in-place implementation of coupled standard maps. It will utilize a handcoded Jacobian, a sparse matrix for the Jacobinan, a default initial Jacobian matrix, as well as function-like-objects as the dynamic rule.

Coupled standard maps is a deterministic iterated map that can have arbitrary number of equations of motion, since you can couple `N`

standard maps which are 2D maps, like so:

\[\theta_{i}' = \theta_i + p_{i}' \\ p_{i}' = p_i + k_i\sin(\theta_i) - \Gamma \left[\sin(\theta_{i+1} - \theta_{i}) + \sin(\theta_{i-1} - \theta_{i}) \right]\]

To model this, we will make a dedicated `struct`

, which is parameterized on the number of coupled maps:

```
using DynamicalSystemsBase
struct CoupledStandardMaps{N}
idxs::SVector{N, Int}
idxsm1::SVector{N, Int}
idxsp1::SVector{N, Int}
end
```

(what these fields are will become apparent later)

We initialize the struct with the amount of standard maps we want to couple, and we also define appropriate parameters:

```
M = 5 # couple number
u0 = 0.001rand(2M) #initial state
ks = 0.9ones(M) # nonlinearity parameters
Γ = 1.0 # coupling strength
p = (ks, Γ) # parameter container
# Create struct:
SV = SVector{M, Int}
idxs = SV(1:M...) # indexes of thetas
idxsm1 = SV(circshift(idxs, +1)...) #indexes of thetas - 1
idxsp1 = SV(circshift(idxs, -1)...) #indexes of thetas + 1
# So that:
# x[i] ≡ θᵢ
# x[[idxsp1[i]]] ≡ θᵢ+₁
# x[[idxsm1[i]]] ≡ θᵢ-₁
csm = CoupledStandardMaps{M}(idxs, idxsm1, idxsp1)
```

`Main.CoupledStandardMaps{5}([1, 2, 3, 4, 5], [5, 1, 2, 3, 4], [2, 3, 4, 5, 1])`

We will now use this struct to define a function-like-object, a Type that also acts as a function

```
function (f::CoupledStandardMaps{N})(xnew::AbstractVector, x, p, n) where {N}
ks, Γ = p
@inbounds for i in f.idxs
xnew[i+N] = mod2pi(
x[i+N] + ks[i]*sin(x[i]) -
Γ*(sin(x[f.idxsp1[i]] - x[i]) + sin(x[f.idxsm1[i]] - x[i]))
)
xnew[i] = mod2pi(x[i] + xnew[i+N])
end
return nothing
end
```

We will use *the same*`struct`

to create a function for the Jacobian:

```
function (f::CoupledStandardMaps{M})(
J::AbstractMatrix, x, p, n) where {M}
ks, Γ = p
# x[i] ≡ θᵢ
# x[[idxsp1[i]]] ≡ θᵢ+₁
# x[[idxsm1[i]]] ≡ θᵢ-₁
@inbounds for i in f.idxs
cosθ = cos(x[i])
cosθp= cos(x[f.idxsp1[i]] - x[i])
cosθm= cos(x[f.idxsm1[i]] - x[i])
J[i+M, i] = ks[i]*cosθ + Γ*(cosθp + cosθm)
J[i+M, f.idxsm1[i]] = - Γ*cosθm
J[i+M, f.idxsp1[i]] = - Γ*cosθp
J[i, i] = 1 + J[i+M, i]
J[i, f.idxsm1[i]] = J[i+M, f.idxsm1[i]]
J[i, f.idxsp1[i]] = J[i+M, f.idxsp1[i]]
end
return nothing
end
```

This is possible because the system state is a `Vector`

while the Jacobian is a `Matrix`

, so multiple dispatch can differentiate between the two.

Notice in addition, that the Jacobian function accesses *only half the elements of the matrix*. This is intentional, and takes advantage of the fact that the other half is constant. We can leverage this further, by making the Jacobian a sparse matrix. Because the `DynamicalSystem`

constructors allow us to give in a pre-initialized Jacobian matrix, we take advantage of that and create:

```
using SparseArrays
J = zeros(eltype(u0), 2M, 2M)
# Set ∂/∂p entries (they are eye(M,M))
# And they dont change they are constants
for i in idxs
J[i, i+M] = 1
J[i+M, i+M] = 1
end
sparseJ = sparse(J)
csm(sparseJ, u0, p, 0) # apply Jacobian to initial state
sparseJ
```

```
10×10 SparseArrays.SparseMatrixCSC{Float64, Int64} with 40 stored entries:
3.9 -1.0 ⋅ ⋅ -1.0 1.0 ⋅ ⋅ ⋅ ⋅
-1.0 3.9 -1.0 ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ ⋅
⋅ -1.0 3.9 -1.0 ⋅ ⋅ ⋅ 1.0 ⋅ ⋅
⋅ ⋅ -1.0 3.9 -1.0 ⋅ ⋅ ⋅ 1.0 ⋅
-1.0 ⋅ ⋅ -1.0 3.9 ⋅ ⋅ ⋅ ⋅ 1.0
2.9 -1.0 ⋅ ⋅ -1.0 1.0 ⋅ ⋅ ⋅ ⋅
-1.0 2.9 -1.0 ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ ⋅
⋅ -1.0 2.9 -1.0 ⋅ ⋅ ⋅ 1.0 ⋅ ⋅
⋅ ⋅ -1.0 2.9 -1.0 ⋅ ⋅ ⋅ 1.0 ⋅
-1.0 ⋅ ⋅ -1.0 2.9 ⋅ ⋅ ⋅ ⋅ 1.0
```

Now we are ready to create our dynamical system

`ds = DeterministicIteratedMap(csm, u0, p)`

```
10-dimensional DeterministicIteratedMap
deterministic: true
discrete time: true
in-place: true
dynamic rule: CoupledStandardMaps
parameters: ([0.9, 0.9, 0.9, 0.9, 0.9], 1.0)
time: 0
state: [0.00041524111829434654, 0.0007229374394926134, 0.0006848502727239867, 0.0005907107844825832, 0.0006803903839332631, 0.0008921850182833862, 0.00012179045162710301, 0.0009690027609471248, 0.0004665835723713002, 0.00037448505102109466]
```

Of course, the reason we went through all this trouble was to make a `TangentDynamicalSystem`

, that can actually use the Jacobian function.

`tands = TangentDynamicalSystem(ds; J = csm, J0 = sparseJ, k = M)`

```
10-dimensional TangentDynamicalSystem
deterministic: true
discrete time: true
in-place: true
dynamic rule: CoupledStandardMaps
jacobian: CoupledStandardMaps
deviation vectors: 5
parameters: ([0.9, 0.9, 0.9, 0.9, 0.9], 1.0)
time: 0
state: [0.00041524111829434654, 0.0007229374394926134, 0.0006848502727239867, 0.0005907107844825832, 0.0006803903839332631, 0.0008921850182833862, 0.00012179045162710301, 0.0009690027609471248, 0.0004665835723713002, 0.00037448505102109466]
```

```
step!(tands, 5)
current_deviations(tands)
```

```
10×5 view(::Matrix{Float64}, :, 2:6) with eltype Float64:
3919.65 -2770.14 845.081 835.566 -2760.47
-2782.26 3943.87 -2784.3 847.221 845.205
846.645 -2773.97 3924.05 -2763.79 836.618
834.299 837.076 -2752.64 3900.81 -2749.87
-2747.17 834.944 833.672 -2745.91 3893.95
3263.48 -2344.26 733.337 723.854 -2334.64
-2356.34 3287.62 -2358.38 735.471 733.462
734.897 -2348.09 3267.89 -2337.95 724.907
722.584 725.36 -2326.83 3244.7 -2324.05
-2321.38 723.235 721.961 -2320.12 3237.88
```

(the deviation vectors will increase in magnitude rapidly because the dynamical system is chaotic)

- DatserisParlitz2022Datseris & Parlitz 2022,
*Nonlinear Dynamics: A Concise Introduction Interlaced with Code*, Springer Nature, Undergrad. Lect. Notes In Physics