DynamicalSystemsBase.jl

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.

DynamicalSystem

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

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

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

Description

A DynamicalSystemrepresents the time evolution of a state in a state space. It mainly encapsulates three things:

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

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

In this scope dynamical systems have a known dynamic rule f. 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 (MTK)

Dynamical systems that have been constructed from DEProblems 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.

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 current and previous time step and is very accurate; 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
• named_variables

API - alter status

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

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.

initial_state(ds::DynamicalSystem) → u0

Return the initial state of ds. This state is never mutated and is set when initializing ds.

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

Return the state u of dsobserved 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 original state space state, not the projected one (this makes the most sense for allowing MTK-based indexing).

Use state_name for an accompanying name.

state_name(index)::String

Return a name that matches the outcome of observe_state with index.

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

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.

parameter_name(index)::String

Return a name that matches the outcome of current_parameter with index.

initial_parameters(ds::DynamicalSystem) → p0

Return the initial parameter container of ds. This is never mutated and is set when initializing ds.

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.

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.

dynamic_rule(ds::DynamicalSystem) → f

Return the dynamic rule of ds. This is never mutated and is set when initializing ds.

current_time(ds::DynamicalSystem) → t

Return the current time that ds is at. This is mutated when ds is evolved.

initial_time(ds::DynamicalSystem) → t0

Return the initial time defined for ds. This is never mutated and is set when initializing ds.

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.

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.

referrenced_sciml_model(ds::DynamicalSystem)

Return the ModelingToolkit.jl structurally-simplified model referrenced by ds. Return nothing if there is no referrenced model.

named_variables(ds::DynamicalSystem)

If ds is constructed via MTK, return a vector of the variable names (as symbols). Otherwise return nothing. If X is a StateSpaceSet, you can always do

X = StateSpaceSet(X; names = named_variables(ds))

in downstream functions to name a set coming from ds (if possible).

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.

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

1. to partially set only some state variables,
2. 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.

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.

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

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.

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. The returned time vector is t = (t0+Ttr):Δt:(t0+Ttr+T).

Optionally provide a starting state u0 which is current_state(ds) by default.

If time evolution diverged or in general failed before T, the remaining of the trajectory is set to the last valid point.

The dimensions of X are automatically named if ds referrences an MTK model and if save_idxs remains at its default value.

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: Which variables to output in X. By default it is nothing (all variables). It can be a vector of any type of index that can be given to observe_state. 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.

Description

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

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:

1. By giving in each individual columns of the state space set as Vector{<:Real}: StateSpaceSet(x, y, z, ...).
2. By giving in a matrix whose rows are the state space points: StateSpaceSet(m).
3. By giving in directly a vector of vectors (state space points): StateSpaceSet(v_of_v).

All constructors allow for two keywords:

• 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.
• names which can be an iterable of length D whose elements are Symbols. This allows assigning a name to each dimension and accessing the dimension by name, see below. names is nothing if not given. Use StateSpaceSet(s; names) to add names to an existing set s.

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

In all examples above, j can also be a Symbol, provided that names has been given when creating the state space set. This allows accessing a dimension by name. This is provided as a convenience and it is not an optimized operation, hence recommended to be used primarily with X[:, j::Symbol].

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 <: 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 <: 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.triviallimiter!), typeof(OrdinaryDiffEqCore.triviallimiter!), Static.False}(OrdinaryDiffEqCore.triviallimiter!, OrdinaryDiffEqCore.triviallimiter!, 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 <: DiscreteTimeDynamicalSystem
StroboscopicMap(ds::CoupledODEs, period::Real) → smap
StroboscopicMap(period::Real, f, u0, p = nothing; kwargs...)

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

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

The convenience constructor

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

is also provided.

See also PoincareMap.

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

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

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. If nothing it is the current_state(ds).
• rootkw = (xrtol = 1e-6, atol = 1e-8): A NamedTuple of keyword arguments passed to find_zero from Roots.jl.
• Tmax = 1e3: The argument Tmax exists so that the integrator can terminate instead of being evolved for infinite time, to avoid cases where iteration would continue forever for ill-defined hyperplanes or for convergence to fixed points, where the trajectory would never cross again the hyperplane. If during one step! the system has been evolved for more than Tmax, then step!(pmap) will terminate and error.

Description

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

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

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

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

In code, plane can be either:

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

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

PoincareMap follows the DynamicalSystem interface with the following adjustments:

1. dimension(pmap) == dimension(ds), even though the Poincaré map is effectively 1 dimension less.
2. Like StroboscopicMap time is discrete and counts the iterations on the surface of section. initial_time is always 0 and current_time is current iteration number.
3. A new function current_crossing_time returns the real time corresponding to the latest crossing of the hyperplane. The corresponding state on the hyperplane is current_state(pmap) as expected.
4. For the special case of plane being a Tuple{Int, <:Real}, a special reinit! method is allowed with input state of length D-1 instead of D, i.e., a reduced state already on the hyperplane that is then converted into the D dimensional state.
5. The initial_state(pmap) returns the state initial state of the map. This is not u0 because u0 is evolved forwards until it resides on the Poincaré plane.
6. In the reinit! function, the t0 keyword denotes the starting time of the continuous time dynamical system, as the starting time of the PoincareMap is by definition always 0.

Example

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

current_crossing_time(pmap::PoincareMap) → tcross

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

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.

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.

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

current_deviations(tands::TangentDynamicalSystem)

Return the deviation vectors of tands as a matrix with each column a vector.

set_deviations!(tands::TangentDynamicalSystem, Q)

Set the deviation vectors of tands to be Q, a matrix with each column a vector.

jacobian(ds::CoreDynamicalSystem)

Construct the Jacobian rule for the dynamical system ds. If the system already has a Jacobian rule constructed via ModelingToolkit.jl it returns this, otherwise it constructs the Jacobian rule with automatic differentiation using module ForwardDiff.

Description

For out-of-place systems, jacobian returns the Jacobian rule as a function Jf(u, p, t = 0) -> 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 = 0). Calling Jf!(J0, u, p) will compute the Jacobian at the state u, parameters p and time t and save the result in J0.

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 <: 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 <: 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 ith 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!.

initial_states(pds::ParallelDynamicalSystem)

Return an iterator over the initial parallel states of pds.

current_states(pds::ParallelDynamicalSystem)

Return an iterator over the parallel states of pds.

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.