# 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

`DynamicalSystemsBase.DynamicalSystem`

— Type`DynamicalSystem`

The central structure of **DynamicalSystems.jl**. All functions of the suite that can use known equations of motion expect an instance of this type.

**Constructing a DynamicalSystem**

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

with `eom`

the equations of motion function (see below). `p`

is a parameter container, which we highly suggest to use a mutable object like `Array`

, `LMArray`

or a dictionary. Pass `nothing`

in the place of `p`

if your system does not have parameters.

`t0`

, `J0`

allow you to choose the initial time and provide an initialized Jacobian matrix. See `CDS_KWARGS`

for the default options used to evolve continuous systems (through `OrdinaryDiffEq`

).

**Equations of motion**

The are two "versions" for `DynamicalSystem`

, depending on whether the equations of motion (`eom`

) are in-place (iip) or out-of-place (oop). Here is how to define them (1-D systems are treated differently, see below):

**oop**: The`eom`

**must**be in the form`eom(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.**iip**: The`eom`

**must**be in the form`eom!(xnew, x, p, t)`

which means that given a state`x::Vector`

and some parameter container`p`

, it writes in-place the new state 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 100 dimensions, and for using functions that use the tangent space (like e.g. `lyapunovs`

or `gali`

), the break-even point is at around 10 dimensions.

The constructor deduces automatically whether `eom`

is iip or oop. It is not possible however to deduce whether the system is continuous or discrete just from the equations of motion, hence the 2 constructors.

**Jacobian**

The optional argument `jacobian`

for the constructors is a *function* and (if given) must also be of the same form as the `eom`

, `jacobian(x, p, n) -> SMatrix`

for the out-of-place version and `jacobian!(Jnew, x, p, n)`

for the in-place version.

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

### 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 Type | equations of motion | Jacobian |
---|---|---|

in-place (big systems) | `eom!(du, u, p, t)` | `jacobian!(J, u, p, t)` |

out-of-place (small systems) | `eom(u, p, t) -> SVector` | `jacobian(u, p, t) -> SMatrix` |

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:

`DelayEmbeddings.dimension`

— Function`dimension(thing) -> D`

Return the dimension of the `thing`

, in the sense of state-space dimensionality.

`DynamicalSystemsBase.jacobian`

— Function`jacobian(ds::DynamicalSystem, u = ds.u0, t = ds.t0)`

Return the jacobian of the system at `u`

, at `t`

.

`DynamicalSystemsBase.set_parameter!`

— Function```
set_parameter!(ds::DynamicalSystem, index, value)
set_parameter!(ds::DynamicalSystem, values)
```

Change one or many parameters of the system by setting `p[index] = value`

in the first case and `p .= values`

in the second.

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)
end
# Jacobian:
@inline @inbounds function loop_jac(u, p, t)
σ, ρ, β = p
J = @SMatrix [-σ σ 0;
ρ - u[3] (-1) (-u[1]);
u[2] u[1] -β]
return J
end
ds = ContinuousDynamicalSystem(loop, rand(3), [10.0, 28.0, 8/3], loop_jac)
```

```
3-dimensional continuous dynamical system
state: [0.068248, 0.828095, 0.0743729]
e.o.m.: loop
in-place? false
jacobian: loop_jac
```

## 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]
return
end
# 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
return
end
ds = DiscreteDynamicalSystem(hiip, zeros(2), [1.4, 0.3], hiip_jac)
```

```
2-dimensional discrete dynamical system
state: [0.0, 0.0]
e.o.m.: hiip
in-place? true
jacobian: hiip_jac
```

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]
e.o.m.: hiip
in-place? true
jacobian: ForwardDiff
```

## 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}
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);
```

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

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:

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

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

`ds = DiscreteDynamicalSystem(csm, u0, p, csm, sparseJ)`

```
10-dimensional discrete dynamical system
state: [0.000803001, 0.00092095, 0.000313022, …, 3.07769e-5, 0.000670152]
e.o.m.: CoupledStandardMaps
in-place? true
jacobian: CoupledStandardMaps
parameters: Tuple
```

## 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)
end
```

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:

`DynamicalSystemsBase.trajectory`

— Function`trajectory(ds::DynamicalSystem, T [, u]; kwargs...) -> dataset`

Return a dataset that will contain the trajectory of the system, after evolving it for total time `T`

, optionally starting from state `u`

. See `Dataset`

for info on how to use this object.

A `W×D`

dataset is returned, with `W = length(t0:dt:T)`

with `t0:dt:T`

representing the time vector (*not* returned) and `D`

the system dimension. For discrete systems both `T`

and `dt`

must be integers.

**Keyword Arguments**

`dt`

: Time step of value output during the solving of the continuous system. For discrete systems it must be an integer. Defaults to`0.01`

for continuous and`1`

for discrete.`Ttr`

: Transient time to evolve the initial state before starting saving states.`diffeq...`

: Keyword arguments propagated into`init`

of DifferentialEquations.jl. For example`abstol = 1e-9`

. Only valid for continuous systems. If you want to specify a solver, do so by using the name`alg`

, e.g.:`alg = Tsit5(), maxiters = 1000`

. This requires you to have been first`using OrdinaryDiffEq`

to access the solvers. See`DynamicalSystemsBase.CDS_KWARGS`

for default values. These keywords can also include`callback`

for event handling. Using a`SavingCallback`

with`trajectory`

will lead to unexpected behavior!

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

interface instead.

### Example

```
using DynamicalSystems
ds = Systems.towel()
```

3-dimensional discrete dynamical system state: [0.085, -0.121, 0.075] e.o.m.: 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

To get every 3-rd point of the trajectory, do

`tr = trajectory(ds, 100; dt = 3)`

3-dimensional Dataset{Float64} with 34 points 0.085 -0.121 0.075 0.681871 0.0508933 0.825263 0.966676 -0.00171595 0.2225 0.910741 -0.0316828 0.411607 0.545332 0.0508239 0.819404 0.544182 -0.0940112 0.270849 0.690701 -0.0679162 0.727605 0.939111 0.010943 0.736301 0.841034 0.0322099 0.689976 0.140788 0.0876455 0.918065 ⋮ 0.223785 0.0902301 0.75948 0.510326 0.0726233 0.682717 0.529268 -0.069801 0.962927 0.435893 -0.0686173 0.92088 0.583619 -0.0769684 0.76526 0.729705 -0.0680752 0.508933 0.739275 0.0308018 0.720691 0.738165 0.0466359 0.825703 0.784243 0.0495777 0.20732

Identical syntax is used for continuous systems

`ds = Systems.lorenz()`

3-dimensional continuous dynamical system state: [0.0, 10.0, 0.0] e.o.m.: loop in-place? false jacobian: loop_jac parameters: [10.0, 28.0, 2.6666666666666665]

`tr = trajectory(ds, 10.0; dt = 0.01)`

3-dimensional Dataset{Float64} with 1001 points 0.0 10.0 0.0 0.951226 10.0352 0.0479002 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.82053 -9.4225 26.6146 -8.879 -9.44532 26.7382

And a final example controlling the integrator accuracy:

```
ds = Systems.lorenz()
tr = trajectory(ds, 10.0; dt = 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

### Solution precision for continuous systems

A numerical solution of an ODE is not the "true" solution, uniquely defined by a (well-defined) ODE and an initial condition. Especially for chaotic systems, where deviations are amplified exponentially, one is left worried if the numerical solutions truly are part of the system and can truly give insight in understanding the system.

DifferentialEquations.jl offers a tool, called Uncertainty Quantification, which allows users to asses up to what time-scales the numerical solution is close to the "true" solution. For example, using the default solving parameters of **DynamicalSystems.jl**, the Lorenz system is accurate up to time `t = 50.0`

.

However, fortunately for us, there is not too much worry about the numerical solution diverging from the true solution. That is because of the shadowing theorem (or shadowing lemma):

Although a numerically computed chaotic trajectory diverges exponentially from the true trajectory with the same initial coordinates, there exists an errorless trajectory with a slightly different initial condition that stays near ("shadows") the numerically computed one.

This simply means that one can always numerically study chaos not only qualitatively but also quantitatively. For more information, see the book *Chaos in Dynamical Systems* by E. Ott, or the scholarpedia entry.