Discrete Systems

Discrete systems are of the form:

DynamicalSystems.jl categorizes discrete systems in three cases, due to the extremely performant handling that StaticArrays offers for small dimensionalities.

When defining a system using Functors, be sure to use EomVector and EomMatrix (see the continuous system page for more info).

High-Dimensional

At around D=10 dimensions, Static Arrays start to become less efficient than Julia's base Arrays, provided that the latter use in-place operations. For cases of discrete systems with much high dimensionality, we offer a type called BigDiscreteDS:

#DynamicalSystemsBase.BigDiscreteDS — Type.

BigDiscreteDS(state, eom! [, jacob! [, J]]) <: DynamicalSystem

D-dimensional discrete dynamical system (used for big D). The equations for this system perform all operations in-place.

Fields:

• state::Vector{T} : Current state-vector of the system. Do ds.state .= u to change the state.
• eom! (function) : The function that represents the system's equations of motion (also called vector field). The function is of the format: eom!(xnew, x) which means that given a state-vector x and another similar one xnew, it writes in-place the new state in xnew.
• jacob! (function) : A function that calculates the system's jacobian matrix, based on the format: jacob!(J, x) which means that given a state-vector x it writes in-place the Jacobian in J.
• J::Matrix{T} : Initialized Jacobian matrix (optional).
• dummystate::Vector{T} : Dummy vector, which most of the time fills the role of the previous state in e.g. evolve.

Only the first two fields of this type are displayed during print.

As mentioned in our official documentation, it is preferred to use Functors for both the equations of motion and the Jacobian.

If the jacob is not provided by the user, it is created automatically using the module ForwardDiff.

source

This system is identical to ContinuousDS as far as definition is concerned. All operations are done in place, and the type is immutable. The same suggestions about using Functors and initialized Jacobians also apply here.

See the source code of the pre-defined coupled standard maps for an example of a BigDiscreteDS definition.

Low-dimensional

The definition of low-dimensional discrete systems differs fundamentally from high dimensional ones, because everything is much more efficiently done with statically sized vectors. The struct representing such systems is called DiscreteDS:

#DynamicalSystemsBase.DiscreteDS — Type.

DiscreteDS(state, eom [, jacob]) <: DynamicalSystem

D-dimensional discrete dynamical system.

Fields:

• state::SVector{D} : Current state-vector of the system, stored in the data format of StaticArray's SVector. Use ds.state = newstate to set a new state.
• eom (function) : The function that represents the system's equations of motion (also called vector field). The function is of the format: eom(u) -> SVector which means that given a state-vector u it returns an SVector containing the next state.
• jacob (function) : A function that calculates the system's jacobian matrix, based on the format: jacob(u) -> SMatrix which means that given a state-vector u it returns an SMatrix containing the Jacobian at that state.

Only the first two fields of this type are displayed during print.

If the jacob is not provided by the user, it is created automatically using the module ForwardDiff.

source

The documentation string of the constructor is perfectly self-contained, but for the sake of clarity we will go through all the steps in the following.

state is simply the state the system starts (a.k.a. initial conditions) and is always of type SVector from StaticArrays.jl. eom is a function that takes a state as an input and returns the next state as an output. The jacob is also a function that takes a state as an input and returns the Jacobian matrix of the system (at this state). So far this is actually different than BigDiscreteDS where the functions where in-place.

Something important to note when defining a DiscreteDS using Functors: since the function calls take only one argument (always a state), it is impossible to use multiple dispatch to differentiate between a call to the e.o.m. or the Jacobian functions.

However, it is very easy to still define both function calls using a single struct, by using a 2 argument function given to the constructor. For example:

using DynamicalSystems
using StaticArrays # only necessary when defining a system

mutable struct HénonMap
    a::Float64
    b::Float64
end
(f::HénonMap)(x::EomVector) = SVector{2}(1.0 - f.a*x[1]^2 + x[2], f.b*x[1])
(f::HénonMap)(x::EomVector, no::Void) = @SMatrix [-2*f.a*x[1] 1.0; f.b 0.0]

function henon(u0=zeros(2); a = 1.4, b = 0.3)
    he = HénonMap(a,b)
    # The jacobian function: (still uses the HénonMap)
    @inline jacob_henon(x) = he(x, nothing)
    return DiscreteDS(u0, he, jacob_henon)
end

ds = henon()

2-dimensional discrete system
 state: [0.0, 0.0]
 eom: ex-henon.HénonMap

Here the example uses the type Void for dispatch, but you could use any other bittype like e.g. ::Float64 and pass in zero(T).

In this example case, doing ds.eom.a = 2.5 would still affect both the equations of motion as well as the Jacobian, making everything work perfectly!

Defining a DynamicalSystem without Functors

As an example of defining a system without using Functors, let's create another one of the Predefined Systems offered by this package, the Folded Towel map.

Because this map doesn't have any parameters, it is unnecessary to associate a Functor with it.

using DynamicalSystems
using StaticArrays # only necessary when defining a system

@inline @inbounds function eom_towel(x)
x1, x2, x3 = x[1], x[2], x[3]
SVector( 3.8*x1*(1-x1) - 0.05*(x2+0.35)*(1-2*x3),
0.1*( (x2+0.35)*(1-2*x3) - 1 )*(1 - 1.9*x1),
3.78*x3*(1-x3)+0.2*x2 )
end

@inline @inbounds function jacob_towel(x)
    @SMatrix [3.8*(1 - 2x[1]) -0.05*(1-2x[3]) 0.1*(x[2] + 0.35);
    -0.19((x[2] + 0.35)*(1-2x[3]) - 1)  0.1*(1-2x[3])*(1-1.9x[1])  -0.2*(x[2] + 0.35)*(1-1.9x[1]);
    0.0  0.2  3.78(1-2x[3]) ]
end

u0=[0.085, -0.121, 0.075]
towel =  DiscreteDS(u0, eom_towel, jacob_towel)

3-dimensional discrete system
 state: [0.085, -0.121, 0.075]
 eom: ex-2.#eom_towel

If we did not want to write a Jacobian for it, we could do

towl_nojac = DiscreteDS(rand(3), eom_towel)

and the Jacobian function is created automatically.

One-Dimensional

In the case of maps, there a special structure for one-dimensional systems. The syntax is DiscreteDS1D(state, eom [, deriv]). In this one-dimensional case, you don't need to worry about StaticArrays.jl because everything is in plain numbers. For example:

using DynamicalSystems

@inline eom_logistic(r) = (x) -> r*x*(1-x)  # this is a closure
@inline deriv_logistic(r) = (x) -> r*(1-2x) # this is a closure
r = 3.7
logistic = DiscreteDS1D(rand(), eom_logistic(r), deriv_logistic(r))

1-dimensional discrete dynamical system:
state: 0.4505133861785231
eom: ex-3.##1#2{Float64}

Once again, if you skip the derivative functions it will be calculated automatically using ForwardDiff.jl.