# 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 allows you to partially set only some state 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. 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 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.

**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.`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} <: 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`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`

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 .

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

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

are also available. 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. 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 to extract the problem.

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

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.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.0004910500876909757, 0.0005451362693489591, 0.00045952843209030284, 0.00038184661251345, 0.00030669932287036305, 0.0008870399208510343, 0.000372312112955555, 0.0004520998624870806, 0.00014542948644631715, 0.0003366323274470259]
```

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.0004910500876909757, 0.0005451362693489591, 0.00045952843209030284, 0.00038184661251345, 0.00030669932287036305, 0.0008870399208510343, 0.000372312112955555, 0.0004520998624870806, 0.00014542948644631715, 0.0003366323274470259]
```

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

```
10×5 view(::Matrix{Float64}, :, 2:6) with eltype Float64:
3915.25 -2765.67 842.737 836.352 -2759.15
-2780.17 3943.21 -2785.32 848.572 843.472
848.42 -2790.46 3953.81 -2791.01 848.971
846.477 849.07 -2789.31 3950.25 -2786.73
-2764.39 836.905 846.347 -2774.03 3924.85
3259.13 -2339.82 730.999 724.649 -2333.34
-2354.25 3286.94 -2359.38 736.813 731.731
736.662 -2364.5 3297.5 -2365.05 737.213
734.725 737.312 -2363.36 3293.95 -2360.78
-2338.55 725.198 734.595 -2348.13 3268.67
```

(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