# Integrators

When writing algorithms for nonlinear dynamics in code, it is almost always more convenient to work with integrators. DynamicalSystems.jl provides several different integrators of dynamical systems and a unified API to treat all of them. The API is based on the same API that DifferentialEquations.jl uses.

## Available integrators

`DynamicalSystemsBase.integrator`

— Function`integrator(ds::DynamicalSystem [, u0]; diffeq) -> integ`

Return an integrator object that can be used to evolve a system interactively using `step!(integ [, Δt])`

. Optionally specify an initial state `u0`

.

The state of this integrator is a vector.

`diffeq`

is a`NamedTuple`

(or`Dict`

) of keyword arguments propagated into`init`

of DifferentialEquations.jl. See`trajectory`

for examples. Only valid for continuous systems.

`DynamicalSystemsBase.parallel_integrator`

— Function`parallel_integrator(ds::DynamicalSystem, states; kwargs...)`

Return an integrator object that can be used to evolve many `states`

of a system in parallel at the *exact same times*, using `step!(integ [, Δt])`

.

`states`

are expected as vectors of vectors.

**Keyword Arguments**

`diffeq`

is a`NamedTuple`

(or`Dict`

) of keyword arguments propagated into`init`

of DifferentialEquations.jl. See`trajectory`

for examples. Only valid for continuous systems.

It is *heavily* advised to use the functions `get_state`

and `set_state!`

to manipulate the integrator. Provide `i`

as a second argument to change the `i`

-th state.

`DynamicalSystemsBase.tangent_integrator`

— Function`tangent_integrator(ds::DynamicalSystem, Q0 | k::Int; kwargs...)`

Return an integrator object that evolves in parallel both the system as well as deviation vectors living on the tangent space, also called linearized space.

`Q0`

is a *matrix* whose columns are initial values for deviation vectors. If instead of a matrix `Q0`

an integer `k`

is given, then `k`

random orthonormal vectors are choosen as initial conditions.

It is *heavily* advised to use the functions `get_state`

, `get_deviations`

, `set_state!`

, `set_deviations!`

to manipulate the integrator.

**Keyword Arguments**

`u0`

: Optional different initial state.`diffeq`

is a`NamedTuple`

(or`Dict`

) of keyword arguments propagated into`init`

of DifferentialEquations.jl. See`trajectory`

for examples. Only valid for continuous systems.

**Description**

If $J$ is the jacobian of the system then the *tangent dynamics* are the equations that evolve in parallel the system as well as a deviation vector (or matrix) $w$:

\[\begin{aligned} \dot{u} &= f(u, p, t) \\ \dot{w} &= J(u, p, t) \times w \end{aligned}\]

with $f$ being the equations of motion and $u$ the system state. Similar equations hold for the discrete case.

`DynamicalSystemsBase.projected_integrator`

— Function`projected_integrator(ds::DynamicalSystem, projection, complete_state; kwargs...) → integ`

Return an integrator that produces iterations of the dynamical system `ds`

on a projected state space. See Integrator API for handling integrators.

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, 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 the full state. 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. Obviously in this case the projected space needs to be lower-dimensional than the original.

Notice that `projected_integrator`

does not require an invertible projection, `complete_state`

is only used during `reinit!`

. Of course the internal integrator operates in the full state space, where the dynamic rule has been defined. The projection always happens as a last step, e.g., during `get_state`

.

**Keyword Arguments**

`u0`

: initial state`diffeq`

is a`NamedTuple`

(or`Dict`

) of keyword arguments propagated into`init`

of DifferentialEquations.jl.

**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
pinteg = projected_integrator(ds, projection, complete_state)
reinit!(pinteg, [0.2, 0.4])
step!(pinteg)
get_state(pinteg)
```

Case 2: custom projection to general functions of state. `julia ds = Systems.lorenz96(5) projection(u) = [sum(u), sqrt(u[1]^2 + u[2]^2)] complete_state(y) = repeat(y[1]/5, 5) pinteg = # same as in above example...`

`

`DynamicalSystemsBase.stroboscopicmap`

— Function`stroboscopicmap(ds::ContinuousDynamicalSystem, T; kwargs...) → smap`

Return a map (integrator) that produces iterations over a period `T`

of the `ds`

, known as a stroboscopic map. See Integrator API for handling integrators.

See also `poincaremap`

.

**Keyword Arguments**

`u0`

: initial state`diffeq`

is a`NamedTuple`

(or`Dict`

) of keyword arguments propagated into`init`

of DifferentialEquations.jl.

**Example**

```
f = 0.27; ω = 0.1
ds = Systems.duffing(zeros(2); ω, f, d = 0.15, β = -1)
smap = stroboscopicmap(ds, 2π/ω; diffeq = (;reltol = 1e-8))
reinit!(smap, [1.0, 1.0])
u = step!(smap)
u = step!(smap, 4) # step 4 iterations forward
```

`ChaosTools.poincaremap`

— Function`poincaremap(ds::ContinuousDynamicalSystem, plane, Tmax=1e3; kwargs...) → pmap`

Return a map (integrator) that produces iterations over the Poincaré map of `ds`

. This map is defined as the sequence of points on the Poincaré surface of section. See `poincaresos`

for details on `plane`

and all other `kwargs`

. Keyword `idxs`

does not apply to `poincaremap`

, as it doesn't save any states.

Notice that while in theory the Poincaré map has one less dimension than `ds`

, in code the map operates on the full `D`

-dimensional state of `ds`

because that is the only way to accommodate planes with generic orientation.

The output `pmap`

follows the Integrator API, i.e., `step!`

and `reinit!`

. `current_time(pmap)`

returns the time of the last crossing. For the special case of `plane`

being a `Tuple{Int, <:Real}`

, a special `reinit!`

method is allowed with input state with length `D-1`

instead of `D`

, i.e., a reduced state already on the hyperplane that is then converted into the `D`

dimensional state.

**Notice**: 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. If during one `step!`

the system has been evolved for more than `Tmax`

, then `step!(pmap)`

will terminate and return `nothing`

.

**Example**

```
ds = Systems.rikitake([0.,0.,0.], μ = 0.47, α = 1.0)
pmap = poincaremap(ds, (3, 0.0))
next_state_on_psos = step!(pmap)
# Change initial condition
reinit!(pmap, [1.0, 0]) # 3rd variable gets value 0 from the plane
next_state_on_psos = step!(pmap)
```

`DynamicalSystemsBase.GeneralizedDynamicalSystem`

— Type`GeneralizedDynamicalSystem`

An umbrella term that includes any instance of `DynamicalSystem`

but also `stroboscopicmap`

, `poincaremap`

and `projeted_integrator`

.

It is used as type declaration in functions that work for any conceivable system type, such as `trajectory`

or `AttractorMapper`

.

## Integrator API

After you have initialized any integrator, you can use the same functions to handle them. The most important function is `step!`

, that will simply progress the integrator. `get_state`

will return its current state. `reinit!`

can be used to re-start the integrator at a possibly different new state. `current_time`

returns the time the integrator is currently at.

Especially for the `tangent_integrator`

, there are also two more functions: `get_deviations`

, `set_deviations`

.

`SciMLBase.step!`

— Method`step!(discrete_integ [, dt::Int])`

Progress the discrete-time integrator for 1 or `dt`

steps.

`step!(continuous_integ, [, dt [, stop_at_tdt]])`

Progress the continuous-time integrator for one step.

Alternative, if a `dt`

is given, then `step!`

the integrator until there is a temporal difference `≥ dt`

(so, step *at least* for `dt`

time).

When `true`

is passed to the optional third argument, the integrator advances for exactly `dt`

time.

`DynamicalSystemsBase.get_state`

— Function`get_state(ds::DynamicalSystem)`

Return the state of `ds`

.

`get_state(integ [, i::Int = 1])`

Return the state of the integrator, in the sense of the state of the dynamical system. This function works for any integrator of the library.

If the integrator is a `parallel_integrator`

, passing `i`

will return the `i`

-th state. The function also correctly returns the true state of the system for tangent integrators.

`SciMLBase.reinit!`

— Method`reinit!(integ [, state])`

Re-initialize the integrator to the new given `state`

. The `state`

type should match the integrator type. Specifically:

`integrator, stroboscopicmap, poincaremap`

: it needs to be a vector of size the same as the dimension of the dynamical system that produced the integrator.`projected_integrator`

: it needs to be a vector of same size as the projected space.`parallel_integrator`

: it needs to be a vector of vectors of size the same as the dimension of the dynamical system that produced the integrator.`tangent_integrator`

: there the special signature`reinit!(integ, state, Q0::AbstractMatrix)`

should be used, with`Q0`

containing the deviation vectors. Alternatively, use`set_deviations!`

if you only want to change the deviation vectors.

`DynamicalSystemsBase.current_time`

— Function`current_time(integ) → t`

Return the current time of the integrator.

`DynamicalSystemsBase.get_deviations`

— Function`get_deviations(tang_integ)`

Return the deviation vectors of the `tangent_integrator`

in a form of a matrix with columns the vectors.

`DynamicalSystemsBase.set_deviations!`

— Function`set_deviations!(tang_integ, Q)`

Set the deviation vectors of the `tangent_integrator`

to `Q`

, which must be a matrix with each column being a deviation vector.

Automatically does `u_modified!(tang_integ, true)`

.

### Re-initializing an integrator

It is more efficient to re-initialize an integrator using `reinit!`

than to create a new one. This can be very helpful when looping over initial conditions and/or parameter values.

All high-level functions from `ChaosTools`

have a set-up part that creates an integrator, and a low-level part that does the computation. The low level part is your friend! Use it! See the Using GALI page for an example as well as the section below.

#### Example: Re-init of continuous tangent integrator

Here we compute the `lyapunovspectrum`

for many different initial conditions.

```
ds = Systems.lorenz()
tinteg = tangent_integrator(ds, 2)
ics = [rand(3) for i in 1:100]
for ic in ics
reinit!(tinteg, ic, orthonormal(3, 2))
λ = lyapunovspectrum(tinteg, 1000, 0.1, 10.0)
# reminder: lyapunovspectrum(tinteg, N, Δt::Real, Ttr::Real = 0.0)
end
```

## Using callbacks with integrators

For the case of continuous systems you can add callbacks from the event handling of **DifferentialEquations.jl**. This is done simply as a keyword argument to the initializers.

In this example we use a simple `SavingCallback`

to save the distance between the two states of a `parallel_integrator`

.

```
using DynamicalSystems, DiffEqCallbacks
using LinearAlgebra: norm
ds = Systems.lorenz()
d0 = 1e-9
T = 100.0
save_func(u, t, integrator) = norm(u[1] - u[2])
saved_values = SavedValues(eltype(ds.t0), eltype(get_state(ds)))
cb = SavingCallback(save_func, saved_values)
diffeq = (abstol=1e-14, reltol=1e-14, maxiters=1e9, callback = cb)
u0 = get_state(ds)
pinteg = parallel_integrator(ds, [u0, u0 + rand(SVector{3})*d0*√3]; diffeq)
step!(pinteg, T)
t = saved_values.t
n = saved_values.saveval
```

`Float64[]`

As expected you can see that the recorded distance between two states is increasing.

## Choosing a solver

`ContinuousDynamicalSystem`

s are evolved using solvers from DifferentialEquations.jl. In this page we discuss the importance of which solver to choose.

### Default Solver

The default solver is:

```
using DynamicalSystems
DynamicalSystemsBase.DEFAULT_SOLVER
```

`SimpleATsit5()`

which is a Runge-Kutta-like solver. The number in the solver's name is the "order" of the solver.

### Speed of a solver

Estimating a given solver's performance for a particular problem is not trivial. The following are general rules of thumb:

- Higher order solvers call the equations of motion function more times per step.
- Higher order solvers can cover larger timespans per step.
- Higher order solvers do better at small tolerances.

This means that there is a delicate balance between how expensive is your function and how large of a step a solver can take while it is still efficient. In general you want to strike a point of taking large steps but also not calling the function exceedingly often.

### How do I pick?

The answer to this question is easy: **benchmarks!**

Here is a simple case: let's compute the Lyapunov spectrum of the Lorenz system using `lyapunovspectrum`

:

```
ds = Systems.lorenz()
tols = (abstol = 1e-6, reltol = 1e-6)
lyapunovspectrum(ds, 2000; Ttr = 100.0, diffeq = tols)
```

```
3-element Vector{Float64}:
0.9138290308099412
0.00020584718966332785
-14.580623290765311
```

The above uses the default solver. Let's now benchmark using two different solvers, `SimpleATsit5`

and `Vern9`

. Since the `SimpleATsit5`

case is of lower order, naively one might think it is faster because it makes less function calls. This argument is not necessarily true though.

It is important to understand that when calling `lyapunovspectrum(ds, 2000)`

you want the system (and the tangent space) to be evolved so that it reaches a total time of `2000*Δt`

, which by default is `2000.0`

units of time. Even though `SimpleATsit5`

requires less function calls per step, `Vern9`

can cover larger timespans per step.

Here are the numbers:

```
using BenchmarkTools, OrdinaryDiffEq, SimpleDiffEq, Statistics
b1 = @benchmark lyapunovspectrum(ds, 2000; diffeq = (alg = SimpleATsit5(), tols...), Ttr = 100.0);
b2 = @benchmark lyapunovspectrum(ds, 2000; diffeq = (alg = Vern9(), tols...), Ttr = 100.0);
println("Timing for SimpleATsit5:")
println(mean(b1))
println("Timing for Vern9:")
println(mean(b2))
```

```
Timing for SimpleATsit5:
TrialEstimate(53.517 ms)
Timing for Vern9:
TrialEstimate(27.511 ms)
```

As you can see `Vern9`

is faster in doing the *entire* computation! Of course this does not have to be universally true. It is true for the Lorenz system, but for your specific system you should do dedicated benchmarks!

### DifferentialEquations.jl

For more info about the possible solvers be sure to head over to the documentation of DifferentialEquations.jl!