Dynamical System Definition

In DynamicalSystems.jl a Dynamical System can be either in continuous time

\[\frac{d\vec{u}}{dt} = \vec{f}(\vec{u}, p, t),\]

or discrete time

\[\vec{u}_{n+1} = \vec{f}(\vec{u}_n, p, n)\]

where $u$ is the state of the system and $p$ contains the parameters of the system. The function $f$ is called the dynamic rule of the system, also known as equations of motion.

In addition to $f$, information about the Jacobian of the system $J_f$ is also used throughout the library.


Keep in mind that almost all functions of DynamicalSystems.jl that use a DynamicalSystem assume that $f$ is differentiable!

Creating a Dynamical System


A central structure of DynamicalSystems.jl. All functions of the suite that can use known dynamic rule f (also called equations of motion or vector field) expect an instance of this type.

Constructing a DynamicalSystem

DiscreteDynamicalSystem(f, state, p [, jacobian [, J0]]; t0::Int = 0)
ContinuousDynamicalSystem(f, state, p [, jacobian [, J0]]; t0 = 0.0)

with f a Julia function (see below). p is a parameter container, which we highly suggest to be a mutable, concretely typed container. Pass nothing as p if your system does not have parameters. Use get_parameters, get_parameter and set_parameter! to interact with p.

t0, J0 allow you to choose the initial time and provide an initialized Jacobian matrix. Continuous systems are evolved via the solvers of DifferentialEquations.jl, see CDS_KWARGS for the default options and the discussion in trajectory.

Dynamic rule f

The are two "versions" for DynamicalSystem, depending on whether f is in-place (iip) or out-of-place (oop). Here is how to define them (1D systems are treated differently, see below):

  • oop : f must be in the form f(x, p, t) -> SVector which means that given a state x::SVector and some parameter container p it returns an SVector (from the StaticArrays module) containing the next state/rate-of-change.
  • iip : f must be in the form f!(xnew, x, p, t) which means that given a state x::Vector and some parameter container p, it writes in-place the new state/rate-of-change in xnew.

t stands for time (integer for discrete systems). iip is suggested for big systems, whereas oop is suggested for small systems. The break-even point at around 10 dimensions.

The constructor deduces automatically whether f is iip or oop. It is not possible however to deduce whether the system is continuous or discrete just from f, hence the 2 constructors.

Autonomous vs non-autonomous systems

Whether the dynamical system is autonomous (f doesn't depend on time) or not, it is still necessary to include t as an argument of f. While for some methods of DynamicalSystems.jl time-dependence is okay, the theoretical foundation of many functions of the library only makes sense with autonomous systems. If you use a non-autonomous system, it is your duty to know for which functions this is okay.


The optional argument jacobian for the constructors is a function and (if given) must also be of the same form as f, jacobian(x, p, n) -> SMatrix for the out-of-place version and jacobian!(J, x, p, n) for the in-place version.

The constructors also allows you to pass an initialized Jacobian matrix J0. This is useful for large iip systems where only a few components of the Jacobian change during the time evolution. J0 can have sparse structure for iip.

If jacobian is not given, it is constructed automatically using the module ForwardDiff. 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 we recommend it.

Comment on 1-D

One dimensional discrete systems expect the state always as a pure number, 0.8 instead of SVector(0.8). For continuous systems, the state can be in-place/out-of-place as in higher dimensions, however the derivative function must be always explicitly given.

Interface to DifferentialEquations.jl

Continuous systems are solved using DifferentialEquations.jl, by default using the keyword arguments contained in the constant CDS_KWARGS.

The following two interfaces are provided:

ContinuousDynamicalSystem(prob::ODEProblem [, jacobian [, J0]])
ODEProblem(continuous_dynamical_system, tspan, args...)

where in the second case args stands for the standard extra arguments of ODEProblem: callback, mass_matrix.

If you want to use callbacks with tangent_integrator or parallel_integrator, then invoke them with extra arguments as shown in the Advanced Documentation.

Relevant Functions

trajectory, set_parameter!.

Definition Table

Here is a handy table that summarizes in what form should be the functions required for the equations of motion and the Jacobian, for each system type:

System Typeequations of motionJacobian
in-place (big systems)eom!(du, u, p, t)jacobian!(J, u, p, t)
out-of-place (small systems)eom(u, p, t) -> SVectorjacobian(u, p, t) -> SMatrix
Use mutable containers for the parameters

It is highly suggested to use a subtype of Array, LMArray or a dictionary for the container of the model's parameters. Some functions offered by DynamicalSystems.jl, like e.g. orbitdiagram, assume that the parameters can be first accessed by p[x] with x some qualifier as well as that this value can be set by p[x] = newvalue.

The Labelled Arrays package offers Array implementations that can be accessed both by index as well as by some name.

Convenience functions

The following functions are defined for convenience for any dynamical system:

dimension(thing) -> D

Return the dimension of the thing, in the sense of state-space dimensionality.

jacobian(ds::DynamicalSystem, u = ds.u0, p = ds.p, t = ds.t0)

Return the jacobian of the system, always as a new matrix.

set_parameter!(ds_or_integ, index, value)

Change a parameter of the system given the index it has in the parameter container p and the value to set it to. This function works for both array/dictionary containers as well as composite types. In the latter case index needs to be a Symbol.

The same function also works for any integrator.

Example: continuous, out-of-place

Let's see an example for a small system, which is a case where out-of-place equations of motion are preferred.

using DynamicalSystems # also exports relevant StaticArrays names
# Lorenz system
# Equations of motion:
@inline @inbounds function loop(u, p, t)
    σ = p[1]; ρ = p[2]; β = p[3]
    du1 = σ*(u[2]-u[1])
    du2 = u[1]*(ρ-u[3]) - u[2]
    du3 = u[1]*u[2] - β*u[3]
    return SVector{3}(du1, du2, du3)
# Jacobian:
@inline @inbounds function loop_jac(u, p, t)
    σ, ρ, β = p
    J = @SMatrix [-σ  σ  0;
    ρ - u[3]  (-1)  (-u[1]);
    u[2]   u[1]  -β]
    return J

ds = ContinuousDynamicalSystem(loop, rand(3), [10.0, 28.0, 8/3], loop_jac)
3-dimensional continuous dynamical system
 state:       [0.495678, 0.17685, 0.972679]
 rule f:      loop
 in-place?    false
 jacobian:    loop_jac
 parameters:  [10.0, 28.0, 2.66667]

Example: discrete, in-place

The following example is only 2-dimensional, and thus once again it is "correct" to use out-of-place version with SVector. For the sake of example though, we use the in-place version.

# Henon map.
# equations of motion:
function hiip(dx, x, p, n)
    dx[1] = 1.0 - p[1]*x[1]^2 + x[2]
    dx[2] = p[2]*x[1]
# Jacobian:
function hiip_jac(J, x, p, n)
    J[1,1] = -2*p[1]*x[1]
    J[1,2] = 1.0
    J[2,1] = p[2]
    J[2,2] = 0.0
ds = DiscreteDynamicalSystem(hiip, zeros(2), [1.4, 0.3], hiip_jac)
2-dimensional discrete dynamical system
 state:       [0.0, 0.0]
 rule f:      hiip
 in-place?    true
 jacobian:    hiip_jac
 parameters:  [1.4, 0.3]

Or, if you don't want to write a Jacobian and want to use the auto-differentiation capabilities of DynamicalSystems.jl, which use the module ForwardDiff:

ds = DiscreteDynamicalSystem(hiip, zeros(2), [1.4, 0.3])
2-dimensional discrete dynamical system
 state:       [0.0, 0.0]
 rule f:      hiip
 in-place?    true
 jacobian:    ForwardDiff
 parameters:  [1.4, 0.3]

Complex Example

In this example we will go through the implementation of the coupled standard maps from our Predefined Dynamical Systems. It is the most complex implementation and takes full advantage of the flexibility of the constructors. The example will use a function-like-object as equations of motion, as well as a sparse matrix for the Jacobian.

Coupled standard maps is a big mapping that can have arbitrary number of equations of motion, since you can couple N standard maps which are 2D maps, like:

\[\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:

struct CoupledStandardMaps{N}
    idxs::SVector{N, Int}
    idxsm1::SVector{N, Int}
    idxsp1::SVector{N, Int}

(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])
    return nothing

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]]
    return nothing

The only reason that this is possible, is because the eom always takes a AbstractVector as first argument, while the Jacobian always takes an AbstractMatrix. Therefore we can take advantage of multiple dispatch!

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
sparseJ = sparse(J)

csm(sparseJ, u0, p, 0) # apply Jacobian to initial state

And finally, we are ready to create our dynamical system:

ds = DiscreteDynamicalSystem(csm, u0, p, csm, sparseJ)
10-dimensional discrete dynamical system
 state:       [0.000975375, 0.000235633, …, 0.00016233, 0.000931795]
 rule f:      CoupledStandardMaps
 in-place?    true
 jacobian:    CoupledStandardMaps
 parameters:  ([0.9, 0.9, 0.9, 0.9, 0.9], 1.0)

Automatic Jacobians

Notice that if you are using automatic differentiation for the Jacobian, you should take care to NOT define your equations of motion so that they explicitly use, or return, Float64 numbers. This is because ForwardDiff uses DualNumbers for differentiation. For example, if you did

function lorenz(u,p,t)
    σ, ρ, β = p
    dx = zeros(3)
    du1 = σ*(u[2] - u[1]) +
    du2 = u[1]*(ρ - u[3]) - u[2]
    du3 = u[1]*u[2] - β*u[3]
    return SVector{Float64, 3}(du1, du2, du3)

this function could not be used to auto-differentiate, as you would get an error when adding dual numbers to SVector{Float64}. Instead, leave the number type untyped, or use eltype(u) as the number type.

Time Evolution of Systems

DynamicalSystems.jl provides a convenient function for getting a trajectory of a system at equally spaced time points:

trajectory(ds::GeneralizedDynamicalSystem, T [, u]; kwargs...) → dataset

Return a dataset that will contain the trajectory of the system ds, after evolving it for total time T, optionally starting from state u. See Dataset for info on how to use this object.

The time vector is t = (t0+Ttr):Δt:(t0+Ttr+T) and is not returned (t0 is the starting time of ds which is by default 0).

trajectory also works with any type of dynamical system where it makes sense, see GeneralizedDynamicalSystem.

Keyword Arguments

  • Δt : Time step of value output. For discrete systems it must be an integer. Defaults to 0.01 for continuous and 1 for discrete.
  • Ttr=0 : Transient time to evolve the initial state before starting saving states.
  • save_idxs::AbstractVector{Int} : Which variables to output in the dataset (by default all).
  • diffeq is a NamedTuple (or Dict) of keyword arguments propagated into init of DifferentialEquations.jl. Only valid for continuous systems, see below.

DifferentialEquations.jl keyword arguments

Continuous dynamical systems are evolved via the solvers of DifferentialEquations.jl. Functions in DynamicalSystems.jl allow providing options to these solvers via 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(), maxiters = 100000). This requires you to have been first using OrdinaryDiffEq to access the solvers. See the CDS_KWARGS variable for the default values we use. Notice that diffeq keywords can also include callback for event handling.

Keep in mind that the default solver is SimpleATsit5, which only supports adaptive time-stepping. Use (alg = SimpleTsit5(), dt = your_step_size) as keywords for a non-adaptive time stepping solver.

Notice that if you want to do repeated evolutions of different states of a continuous system, you should use the integrator interface instead.


using DynamicalSystems
ds = Systems.towel()
3-dimensional discrete dynamical system
 state:       [0.085, -0.121, 0.075]
 rule f:      rule_towel
 in-place?    false
 jacobian:    jacob_towel
 parameters:  nothing
tr = trajectory(ds, 100)
3-dimensional Dataset{Float64} with 101 points
 0.085     -0.121       0.075
 0.285813  -0.0675286   0.238038
 0.76827   -0.038933    0.672094
 0.681871   0.0508933   0.825263
 0.837347   0.0372633   0.555269
 0.51969    0.0616256   0.940906
 0.966676  -0.00171595  0.2225
 0.112748   0.0674955   0.653573
 0.386547  -0.0886542   0.869349
 0.910741  -0.0316828   0.411607
 0.729743   0.0566836   0.25032
 0.739275   0.0308018   0.720691
 0.740845   0.047263    0.767058
 0.740185   0.0494093   0.684863
 0.738165   0.0466359   0.825703
 0.747373   0.0506512   0.553337
 0.719603   0.0437958   0.944377
 0.784243   0.0495777   0.20732
 0.631288   0.0375437   0.631114

Identical syntax is used for continuous systems

ds = Systems.lorenz()
3-dimensional continuous dynamical system
 state:       [0.0, 10.0, 0.0]
 rule f:      loop
 in-place?    false
 jacobian:    loop_jac
 parameters:  [10.0, 28.0, 2.66667]
tr = trajectory(ds, 10.0; Δt = 0.01)
3-dimensional Dataset{Float64} with 1001 points
  0.0       10.0       0.0
  0.951226  10.0352    0.0479007
  1.82768   10.3224    0.18683
  2.6593    10.839     0.416592
  3.47133   11.5678    0.744804
  4.28518   12.496     1.18585
  5.1191    13.6125    1.76047
  5.98858   14.9063    2.49569
  6.90648   16.3638    3.42487
  7.88291   17.9659    4.58749
 -8.3603    -9.01682  25.9478
 -8.42664   -9.09574  26.0152
 -8.49391   -9.16911  26.0931
 -8.56145   -9.23601  26.1809
 -8.62859   -9.29558  26.278
 -8.69462   -9.34701  26.3833
 -8.75883   -9.38954  26.4959
 -8.82052   -9.4225   26.6146
 -8.879     -9.44531  26.7382

And a final example controlling the integrator accuracy:

ds = Systems.lorenz()
tr = trajectory(ds, 10.0; Δt = 0.1, abstol = 1e-9, reltol = 1e-9)
3-dimensional Dataset{Float64} with 101 points
  0.0      10.0       0.0
  8.92478  19.6843    6.02839
 20.0553   24.8278   39.9785
  9.89835  -7.67212  42.312
 -2.67297  -9.86486  30.4531
 -6.73013  -8.71064  27.4841
 -7.91     -8.61271  26.6731
 -8.44189  -8.84913  26.7275
 -8.72492  -8.86503  27.1909
 -8.70356  -8.53165  27.5098
 -7.90028  -8.12207  25.9823
 -8.35252  -8.96846  25.9951
 -8.94418  -9.4035   26.9655
 -9.09673  -8.90402  27.9987
 -8.6238   -7.96535  28.007
 -8.01079  -7.54848  27.0382
 -7.83288  -7.96251  26.0288
 -8.23289  -8.84589  25.8466
 -8.87894  -9.44531  26.7381