Overarching tutorial for DynamicalSystems.jl
This page serves as a short but to-the-point introduction to the DynamicalSystems.jl library. It outlines the core components, and how they establish an interface that is used by the rest of the library. It also provides a couple of usage examples to connect the various packages of the library together.
Going through this tutorial should take you about 20 minutes.
Installation
To install DynamicalSystems.jl, simply do:
using Pkg; Pkg.add("DynamicalSystems")
This installs several packages for the Julia language. These are the sub-modules/packages that comprise DynamicalSystems.jl, see contents for more. All of the functionality is brought into scope when doing:
using DynamicalSystems
in your Julia session.
Package versions used
import Pkg
Pkg.status(["DynamicalSystems", "CairoMakie", "GLMakie", "OrdinaryDiffEq", "BenchmarkTools"]; mode = Pkg.PKGMODE_MANIFEST)
Status `~/work/DynamicalSystems.jl/DynamicalSystems.jl/docs/Manifest.toml`
[6e4b80f9] BenchmarkTools v1.5.0
[13f3f980] CairoMakie v0.12.7
[61744808] DynamicalSystems v3.3.20 `~/work/DynamicalSystems.jl/DynamicalSystems.jl`
[1dea7af3] OrdinaryDiffEq v6.87.0
Core components
The individual packages that compose DynamicalSystems
interact flawlessly with each other because of the following two components:
- The
StateSpaceSet
, which represents numerical data. They can be observed or measured from experiments, sampled trajectories of dynamical systems, or just unordered sets in a state space. AStateSpaceSet
is a container of equally-sized points, representing multivariate timeseries or multivariate datasets. Timeseries, which are univariate sets, are represented by theAbstractVector{<:Real}
Julia base type. - The
DynamicalSystem
, which is the abstract representation of a dynamical system with a known dynamic evolution rule.DynamicalSystem
defines an extendable interface, but typically one uses existing implementations such asDeterministicIteratedMap
orCoupledODEs
.
Making dynamical systems
In the majority of cases, to make a dynamical system one needs three things:
- The dynamic rule
f
: A Julia function that provides the instructions of how to evolve the dynamical system in time. - The state
u
: An array-like container that contains the variables of the dynamical system and also defines the starting state of the system. - The parameters
p
: An arbitrary container that parameterizesf
.
For most concrete implementations of DynamicalSystem
there are two ways of defining f, u
. The distinction is done on whether f
is defined as an in-place (iip) function or out-of-place (oop) function.
- oop :
f
must be in the formf(u, p, t) -> out
which means that given a stateu::SVector{<:Real}
and some parameter containerp
it returns the output off
as anSVector{<:Real}
(static vector). - iip :
f
must be in the formf!(out, u, p, t)
which means that given a stateu::AbstractArray{<:Real}
and some parameter containerp
, it writes in-place the output off
inout::AbstractArray{<:Real}
. The function must returnnothing
as a final statement.
t
stands for current time in both cases. iip is suggested for systems with high dimension and oop for small. The break-even point is between 10 to 100 dimensions but should be benchmarked on a case-by-case basis as it depends on the complexity of f
.
Whether the dynamical system is autonomous (f
doesn't depend on time) or not, it is still necessary to include t
as an argument to f
. Some algorithms utilize this information, some do not, but we prefer to keep a consistent interface either way.
Example: Henon map
Let's make the Henon map, defined as
\[\begin{aligned} x_{n+1} &= 1 - ax^2_n+y_n \\ y_{n+1} & = bx_n \end{aligned}\]
with parameters $a = 1.4, b = 0.3$.
First, we define the dynamic rule as a standard Julia function. Since the dynamical system is only two-dimensional, we should use the out-of-place form that returns an SVector
with the next state:
using DynamicalSystems
function henon_rule(u, p, n) # here `n` is "time", but we don't use it.
x, y = u # system state
a, b = p # system parameters
xn = 1.0 - a*x^2 + y
yn = b*x
return SVector(xn, yn)
end
henon_rule (generic function with 1 method)
Then, we define initial state and parameters
u0 = [0.2, 0.3]
p0 = [1.4, 0.3]
2-element Vector{Float64}:
1.4
0.3
Lastly, we give these three to the DeterministicIteratedMap
:
henon = DeterministicIteratedMap(henon_rule, u0, p0)
2-dimensional DeterministicIteratedMap
deterministic: true
discrete time: true
in-place: false
dynamic rule: henon_rule
parameters: [1.4, 0.3]
time: 0
state: [0.2, 0.3]
henon
is a DynamicalSystem
, one of the two core structures of the library. They can evolved interactively, and queried, using the interface defined by DynamicalSystem
. The simplest thing you can do with a DynamicalSystem
is to get its trajectory:
total_time = 10_000
X, t = trajectory(henon, total_time)
X
2-dimensional StateSpaceSet{Float64} with 10001 points
0.2 0.3
1.244 0.06
-1.10655 0.3732
-0.341035 -0.331965
0.505208 -0.102311
0.540361 0.151562
0.742777 0.162108
0.389703 0.222833
1.01022 0.116911
-0.311842 0.303065
⋮
-0.582534 0.328346
0.853262 -0.17476
-0.194038 0.255978
1.20327 -0.0582113
-1.08521 0.36098
-0.287758 -0.325562
0.558512 -0.0863275
0.476963 0.167554
0.849062 0.143089
X
is a StateSpaceSet
, the second of the core structures of the library. We'll see below how, and where, to use a StateSpaceset
, but for now let's just do a scatter plot
using CairoMakie
scatter(X[:, 1], X[:, 2])
Example: Lorenz96
Let's also make another dynamical system, the Lorenz96 model:
\[\frac{dx_i}{dt} = (x_{i+1}-x_{i-2})x_{i-1} - x_i + F\]
for $i \in \{1, \ldots, N\}$ and $N+j=j$.
Here, instead of a discrete time map we have $N$ coupled ordinary differential equations. However, creating the dynamical system works out just like above, but using CoupledODEs
instead of DeterministicIteratedMap
.
First, we make the dynamic rule function. Since this dynamical system can be arbitrarily high-dimensional, we prefer to use the in-place form for f
, overwriting in place the rate of change in a pre-allocated container. It is customary to append the name of functions that modify their arguments in-place with a bang (!
).
function lorenz96_rule!(du, u, p, t)
F = p[1]; N = length(u)
# 3 edge cases
du[1] = (u[2] - u[N - 1]) * u[N] - u[1] + F
du[2] = (u[3] - u[N]) * u[1] - u[2] + F
du[N] = (u[1] - u[N - 2]) * u[N - 1] - u[N] + F
# then the general case
for n in 3:(N - 1)
du[n] = (u[n + 1] - u[n - 2]) * u[n - 1] - u[n] + F
end
return nothing # always `return nothing` for in-place form!
end
lorenz96_rule! (generic function with 1 method)
then, like before, we define an initial state and parameters, and initialize the system
N = 6
u0 = range(0.1, 1; length = N)
p0 = [8.0]
lorenz96 = CoupledODEs(lorenz96_rule!, u0, p0)
6-dimensional CoupledODEs
deterministic: true
discrete time: false
in-place: true
dynamic rule: lorenz96_rule!
ODE solver: Tsit5
ODE kwargs: (abstol = 1.0e-6, reltol = 1.0e-6)
parameters: [8.0]
time: 0.0
state: [0.1, 0.28, 0.46, 0.64, 0.82, 1.0]
and, again like before, we may obtain a trajectory the same way
total_time = 12.5
sampling_time = 0.02
Y, t = trajectory(lorenz96, total_time; Ttr = 2.2, Δt = sampling_time)
Y
6-dimensional StateSpaceSet{Float64} with 626 points
3.15368 -4.40493 0.0311581 0.486735 1.89895 4.15167
2.71382 -4.39303 0.395019 0.66327 2.0652 4.32045
2.25088 -4.33682 0.693967 0.879701 2.2412 4.46619
1.7707 -4.24045 0.924523 1.12771 2.42882 4.58259
1.27983 -4.1073 1.08656 1.39809 2.62943 4.66318
0.785433 -3.94005 1.18319 1.6815 2.84384 4.70147
0.295361 -3.74095 1.2205 1.96908 3.07224 4.69114
-0.181932 -3.51222 1.20719 2.25296 3.3139 4.62628
-0.637491 -3.25665 1.154 2.5267 3.56698 4.50178
-1.06206 -2.9781 1.07303 2.7856 3.82827 4.31366
⋮ ⋮
3.17245 2.3759 3.01796 7.27415 7.26007 -0.116002
3.29671 2.71146 3.32758 7.5693 6.75971 -0.537853
3.44096 3.09855 3.66908 7.82351 6.13876 -0.922775
3.58387 3.53999 4.04452 8.01418 5.39898 -1.25074
3.70359 4.03513 4.45448 8.1137 4.55005 -1.5042
3.78135 4.57879 4.89677 8.09013 3.61125 -1.66943
3.80523 5.16112 5.36441 7.90891 2.61262 -1.73822
3.77305 5.7684 5.84318 7.53627 1.59529 -1.71018
3.6934 6.38507 6.30923 6.94454 0.61023 -1.59518
We can't scatterplot something 6-dimensional but we can visualize all timeseries
fig = Figure()
ax = Axis(fig[1, 1]; xlabel = "time", ylabel = "variable")
for var in columns(Y)
lines!(ax, t, var)
end
fig
ODE solving and choosing solver
Continuous time dynamical systems are evolved through DifferentialEquations.jl. In this sense, the above trajectory
function is a simplified version of DifferentialEquations.solve
. If you only care about evolving a dynamical system forwards in time, you are probably better off using DifferentialEquations.jl directly. DynamicalSystems.jl can be used to do many other things that either occur during the time evolution or after it, see the section below on using dynamical systems.
When initializing a CoupledODEs
you can tune the solver properties to your heart's content using any of the ODE solvers and any of the common solver options. For example:
using OrdinaryDiffEq # accessing the ODE solvers
diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)
lorenz96_vern = ContinuousDynamicalSystem(lorenz96_rule!, u0, p0; diffeq)
6-dimensional CoupledODEs
deterministic: true
discrete time: false
in-place: true
dynamic rule: lorenz96_rule!
ODE solver: Vern9
ODE kwargs: (abstol = 1.0e-9, reltol = 1.0e-9)
parameters: [8.0]
time: 0.0
state: [0.1, 0.28, 0.46, 0.64, 0.82, 1.0]
Y, t = trajectory(lorenz96_vern, total_time; Ttr = 2.2, Δt = sampling_time)
Y[end]
6-element SVector{6, Float64} with indices SOneTo(6):
3.8390248122407535
6.155709531170093
6.080625689047345
7.278588308970321
1.258215221210072
-1.5297062916990858
The choice of the solver algorithm can have huge impact on the performance and stability of the ODE integration! We will showcase this with two simple examples
Higher accuracy, higher order
The solver Tsit5
(the default solver) is most performant when medium-high error tolerances are requested. When we require very small errors, choosing a different solver can be more accurate. This can be especially impactful for chaotic dynamical systems. Let's first expliclty ask for a given accuracy when solving the ODE by passing the keywords abstol, reltol
(for absolute and relative tolerance respectively), and compare performance to a naive solver one would use:
using BenchmarkTools: @btime
using OrdinaryDiffEq: BS3 # equivalent of odeint23
for alg in (BS3(), Vern9())
diffeq = (; alg, abstol = 1e-12, reltol = 1e-12)
lorenz96 = CoupledODEs(lorenz96_rule!, u0, p0; diffeq)
@btime step!($lorenz96, 100.0) # evolve for 100 time units
end
┌ Warning: Assignment to `diffeq` in soft scope is ambiguous because a global variable by the same name exists: `diffeq` will be treated as a new local. Disambiguate by using `local diffeq` to suppress this warning or `global diffeq` to assign to the existing global variable.
└ @ tutorial.md:224
┌ Warning: Assignment to `lorenz96` in soft scope is ambiguous because a global variable by the same name exists: `lorenz96` will be treated as a new local. Disambiguate by using `local lorenz96` to suppress this warning or `global lorenz96` to assign to the existing global variable.
└ @ tutorial.md:225
626.921 ms (0 allocations: 0 bytes)
2.797 ms (0 allocations: 0 bytes)
The performance difference is dramatic!
Stiff problems
A "stiff" ODE problem is one that can be numerically unstable unless the step size (or equivalently, the step error tolerances) are extremely small. There are several situations where a problem may be come "stiff":
- The derivative values can get very large for some state values.
- There is a large timescale separation between the dynamics of the variables
- There is a large speed separation between different state space regions
One must be aware whether this is possible for their system and choose a solver that is better suited to tackle stiff problems. If not, a solution may diverge and the ODE integrator will throw an error or a warning.
Many of the problems in DifferentialEquations.jl are suitable for dealing with stiff problems. We can create a stiff problem by using the well known Van der Pol oscillator with a timescale separation:
\[\begin{aligned} \dot{x} & = y \\ \dot{y} / \mu &= (1-x^2)y - x \end{aligned}\]
with $\mu$ being the timescale of the $y$ variable in units of the timescale of the $x$ variable. For very large values of $\mu$ this problem becomes stiff.
Let's compare
function vanderpol_rule(u, μ, t)
x, y = u
dx = y
dy = μ*((1-x^2)*y - x)
return SVector(dx, dy)
end
μ = 1e6
for alg in (Tsit5(), Rodas5P()) # default vs specialized solver
diffeq = (; alg, abstol = 1e-12, reltol = 1e-12, maxiters = typemax(Int))
vdp = CoupledODEs(vanderpol_rule, SVector(1.0, 1.0), μ; diffeq)
@btime step!($vdp, 100.0)
end
┌ Warning: Assignment to `diffeq` in soft scope is ambiguous because a global variable by the same name exists: `diffeq` will be treated as a new local. Disambiguate by using `local diffeq` to suppress this warning or `global diffeq` to assign to the existing global variable.
└ @ tutorial.md:266
8.578 s (0 allocations: 0 bytes)
1.952 s (19170892 allocations: 562.55 MiB)
We see that the stiff solver Rodas5P
is much faster than the default Tsit5
when there is a large timescale separation. This happened because Rodas5P
required much less steps to integrated the same total amount of time. In fact, there are cases where regular solvers will fail to integrate the ODE if the problem is very stiff, e.g. in the ROBER example.
So using an appropriate solver really does matter! For more information on choosing solvers consult the DifferentialEquations.jl documentation.
Using dynamical systems
You may use the DynamicalSystem
interface to develop algorithms that utilize dynamical systems with a known evolution rule. The two main packages of the library that do this are ChaosTools
and Attractors
. For example, you may want to compute the Lyapunov spectrum of the Lorenz96 system from above. This is as easy as calling the lyapunovspectrum
function with lorenz96
steps = 10_000
lyapunovspectrum(lorenz96, steps)
6-element Vector{Float64}:
0.9480080466096642
0.00012614364764332677
-0.1568159270952503
-0.747794014329633
-1.4066157568705158
-4.636905163494868
As expected, there is at least one positive Lyapunov exponent, because the system is chaotic, and at least one zero Lyapunov exponent, because the system is continuous time.
Alternatively, you may want to estimate the basins of attraction of a multistable dynamical system. The Henon map is "multistable" in the sense that some initial conditions diverge to infinity, and some others converge to a chaotic attractor. Computing these basins of attraction is simple with Attractors
, and would work as follows:
# define a state space grid to compute the basins on:
xg = yg = range(-2, 2; length = 201)
# find attractors using recurrences in state space:
mapper = AttractorsViaRecurrences(henon, (xg, yg); sparse = false)
# compute the full basins of attraction:
basins, attractors = basins_of_attraction(mapper; show_progress = false)
([-1 -1 … -1 -1; -1 -1 … -1 -1; … ; -1 -1 … -1 -1; -1 -1 … -1 -1], Dict{Int64, StateSpaceSet{2, Float64}}(1 => 2-dimensional StateSpaceSet{Float64} with 416 points))
Let's visualize the result
heatmap_basins_attractors((xg, yg), basins, attractors)
One last thing to highlight in this short overview are the interactive GUI apps one can launch to examine a DynamicalSystem
. The simplest is interactive_trajectory_timeseries
. To actually make it interactive one needs to enable GLMakie.jl as a backend:
import GLMakie
GLMakie.activate!()
and then launch the app:
u0s = [10rand(5) for _ in 1:3]
parameter_sliders = Dict(1 => 0:0.01:32)
fig, dsobs = interactive_trajectory_timeseries(
lorenz96, [1, 2, 3, 4, 5], u0s;
Δt = 0.02, parameter_sliders
)
fig
Developing new algorithms
You could also be using a DynamicalSystem
instance directly to build your own algorithm if it isn't already implemented (and then later contribute it so it is implemented ;) ). A dynamical system can be evolved forwards in time using step!
:
henon
2-dimensional DeterministicIteratedMap
deterministic: true
discrete time: true
in-place: false
dynamic rule: henon_rule
parameters: [1.4, 0.3]
time: 5
state: [-1.5266434026801804e8, -3132.7519146699206]
Notice how the time is not 0, because henon
has already been stepped when we called the function basins_of_attraction
with it. We can step it more:
step!(henon)
2-dimensional DeterministicIteratedMap
deterministic: true
discrete time: true
in-place: false
dynamic rule: henon_rule
parameters: [1.4, 0.3]
time: 6
state: [-3.262896110526e16, -4.579930208040541e7]
step!(henon, 2)
2-dimensional DeterministicIteratedMap
deterministic: true
discrete time: true
in-place: false
dynamic rule: henon_rule
parameters: [1.4, 0.3]
time: 8
state: [-3.110262842032839e66, -4.4715262317959936e32]
For more information on how to directly use DynamicalSystem
instances, see the documentation of DynamicalSystemsBase
.
State space sets
Let's recall that the output of the trajectory
function is a StateSpaceSet
:
X
2-dimensional StateSpaceSet{Float64} with 10001 points
0.2 0.3
1.244 0.06
-1.10655 0.3732
-0.341035 -0.331965
0.505208 -0.102311
0.540361 0.151562
0.742777 0.162108
0.389703 0.222833
1.01022 0.116911
-0.311842 0.303065
⋮
-0.582534 0.328346
0.853262 -0.17476
-0.194038 0.255978
1.20327 -0.0582113
-1.08521 0.36098
-0.287758 -0.325562
0.558512 -0.0863275
0.476963 0.167554
0.849062 0.143089
It is printed like a matrix where each column is the timeseries of each dynamic variable. In reality, it is a vector of statically sized vectors (for performance reasons). When indexed with 1 index, it behaves like a vector of vectors
X[1]
2-element SVector{2, Float64} with indices SOneTo(2):
0.2
0.3
X[2:5]
2-dimensional StateSpaceSet{Float64} with 4 points
1.244 0.06
-1.10655 0.3732
-0.341035 -0.331965
0.505208 -0.102311
When indexed with two indices, it behaves like a matrix
X[7:13, 2] # 2nd column
7-element Vector{Float64}:
0.1621081681101694
0.22283309461548204
0.11691103950545975
0.30306503631282444
-0.09355263057973214
0.35007640234803744
-0.29998206408499634
When iterated, it iterates over the contained points
for (i, point) in enumerate(X)
@show point
i > 5 && break
end
point = [0.2, 0.3]
point = [1.244, 0.06]
point = [-1.1065503999999997, 0.3732]
point = [-0.34103530283622296, -0.3319651199999999]
point = [0.5052077711071681, -0.10231059085086688]
point = [0.5403605603672313, 0.1515623313321504]
map(point -> point[1] + 1/(point[2]+0.1), X)
10001-element Vector{Float64}:
2.7
7.494
1.006720944040575
-4.652028265916192
-432.28452634595817
4.515518505137784
4.557995741581754
3.4872793228808314
5.620401692451571
2.1691470931310732
⋮
1.752024323892008
-12.522817477226136
2.6151207744662046
25.133206775921803
1.0840850677362275
-4.721137673509619
73.69787164232238
4.214533120554577
4.9627832739756785
The columns of the set are obtained with the convenience columns
function
x, y = columns(X)
summary.((x, y))
("10001-element Vector{Float64}", "10001-element Vector{Float64}")
Using state space sets
Several packages of the library deal with StateSpaceSets
.
You could use ComplexityMeasures
to obtain the entropy, or other complexity measures, of a given set. Below, we obtain the entropy of the natural density of the chaotic attractor by partitioning into a histogram of approximately 50
bins per dimension:
prob_est = ValueHistogram(50)
entropy(prob_est, X)
7.825799208736613
Or, obtain the permutation and sample entropies of the two columns of X
:
pex = entropy_permutation(x; m = 4)
sey = entropy_sample(y; m = 2)
pex, sey
(3.15987571159201, 0.02579132263914716)
Alternatively, you could use FractalDimensions
to get the fractal dimensions of the chaotic attractor of the henon map using the Grassberger-Procaccia algorithm:
grassberger_proccacia_dim(X)
1.2232922815092426
Or, you could obtain a recurrence matrix of a state space set with RecurrenceAnalysis
R = RecurrenceMatrix(Y, 8.0)
Rg = grayscale(R)
rr = recurrencerate(R)
heatmap(Rg; colormap = :grays,
axis = (title = "recurrence rate = $(round(rr; digits = 3))", aspect = 1)
)
More nonlinear timeseries analysis
A trajectory
of a known dynamical system is one way to obtain a StateSpaceSet
. However, another common way is via a delay coordinates embedding of a measured/observed timeseries. For example, we could use optimal_separated_de
from DelayEmbeddings
to create an optimized delay coordinates embedding of a timeseries
w = Y[:, 1] # first variable of Lorenz96
𝒟, τ, e = optimal_separated_de(w)
𝒟
5-dimensional StateSpaceSet{Float64} with 558 points
3.15369 -2.40036 1.60497 2.90499 5.72572
2.71384 -2.24811 1.55832 3.04987 5.6022
2.2509 -2.02902 1.50499 3.20633 5.38629
1.77073 -1.75077 1.45921 3.37699 5.07029
1.27986 -1.42354 1.43338 3.56316 4.65003
0.785468 -1.05974 1.43672 3.76473 4.12617
0.295399 -0.673567 1.47423 3.98019 3.50532
-0.181891 -0.280351 1.54635 4.20677 2.80048
-0.637447 0.104361 1.64932 4.44054 2.03084
-1.06201 0.465767 1.77622 4.67654 1.22067
⋮
7.42111 9.27879 -1.23936 5.15945 3.25618
7.94615 9.22663 -1.64222 5.24344 3.34749
8.40503 9.13776 -1.81947 5.26339 3.46932
8.78703 8.99491 -1.77254 5.22631 3.60343
9.08701 8.77963 -1.51823 5.13887 3.72926
9.30562 8.47357 -1.08603 5.00759 3.82705
9.4488 8.06029 -0.514333 4.83928 3.88137
9.52679 7.52731 0.153637 4.6414 3.88458
9.55278 6.86845 0.873855 4.42248 3.83902
and compare
fig = Figure()
axs = [Axis3(fig[1, i]) for i in 1:2]
for (S, ax) in zip((Y, 𝒟), axs)
lines!(ax, S[:, 1], S[:, 2], S[:, 3])
end
fig
Since 𝒟
is just another state space set, we could be using any of the above analysis pipelines on it just as easily.
The last package to mention here is TimeseriesSurrogates
, which ties with all other observed/measured data analysis by providing a framework for confidence/hypothesis testing. For example, if we had a measured timeseries but we were not sure whether it represents a deterministic system with structure in the state space, or mostly noise, we could do a surrogate test. For this, we use surrogenerator
and RandomFourier
from TimeseriesSurrogates
, and the generalized_dim
from FractalDimensions
(because it performs better in noisy sets)
x # Henon map timeseries
# contaminate with noise
using Random: Xoshiro
rng = Xoshiro(1234)
x .+= randn(rng, length(x))/100
# compute noise-contaminated fractal dim.
Δ_orig = generalized_dim(embed(x, 2, 1))
1.3801073957979793
And we do the surrogate test
surrogate_method = RandomFourier()
sgen = surrogenerator(x, surrogate_method, rng)
Δ_surr = map(1:1000) do i
s = sgen()
generalized_dim(embed(s, 2, 1))
end
1000-element Vector{Float64}:
1.8297218640747606
1.8449246422083758
1.827998413768852
1.8301078704767055
1.8107042515928138
1.83249783593657
1.8300954188512213
1.8416570396197014
1.8575804541517287
1.821435647618282
⋮
1.8768262063118628
1.8287938131103985
1.8435474417451545
1.8129026565561648
1.8551700924676993
1.8283378225272842
1.821134647531232
1.8369834750866052
1.8368699339310084
and visualize the test result
fig, ax = hist(Δ_surr)
vlines!(ax, Δ_orig)
fig
since the real value is outside the distribution we have confidence the data are not pure noise.
Integration with ModelingToolkit.jl
DynamicalSystems.jl understands when a model has been generated via ModelingToolkit.jl. The symbolic variables used in ModelingToolkit.jl can be used to access the state or parameters of the dynamical system.
To access this functionality, the DynamicalSystem
must be created from a DEProblem
of the SciML ecosystem, and the DEProblem
itself must be created from a ModelingToolkit.jl model.
ProcessBasedModelling.jl is an extension to ModelingToolkit.jl for creating models from a set of equations. It has been designed to be useful for scenarios applicable to a typical nonlinear dynamics analysis workflow, and provides better error messages during system construction than MTK. Have a look at its docs!
Let's create a the Roessler system as an MTK model:
using ModelingToolkit
@variables t # use unitless time
D = Differential(t)
@mtkmodel Roessler begin
@parameters begin
a = 0.2
b = 0.2
c = 5.7
end
@variables begin
x(t) = 1.0
y(t) = 0.0
z(t) = 0.0
nlt(t) # nonlinear term
end
@equations begin
D(x) ~ -y -z
D(y) ~ x + a*y
D(z) ~ b + nlt
nlt ~ z*(x - c)
end
end
@mtkbuild model = Roessler()
\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} &= - y\left( t \right) - z\left( t \right) \\ \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} &= x\left( t \right) + a y\left( t \right) \\ \frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} &= b + \mathrm{nlt}\left( t \right) \end{align} \]
this model can then be made into an ODEProblem
:
prob = ODEProblem(model)
ODEProblem with uType Vector{Float64} and tType Nothing. In-place: true
timespan: (nothing, nothing)
u0: 3-element Vector{Float64}:
1.0
0.0
0.0
(notice that because we specified initial values for all parameters and variables during the model creation we do need to provide additional initial values)
Now, this problem can be made into a CoupledODEs
:
roessler = CoupledODEs(prob)
3-dimensional CoupledODEs
deterministic: true
discrete time: false
in-place: true
dynamic rule: f
ODE solver: Tsit5
ODE kwargs: (abstol = 1.0e-6, reltol = 1.0e-6)
parameters: ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Any, Vector{Any}}, Tuple{}, Tuple{}}([0.2, 0.2, 5.7], Any[], (), ())
time: 0.0
state: [1.0, 0.0, 0.0]
This dynamical system instance can be used in the rest of the library like anything else. Additionally, you can "observe" referenced symbolic variables:
observe_state(roessler, model.x)
1.0
observe_state(roessler, :nlt) # can use `Symbol`s as well
-0.0
These observables can also be used in the GUI visualization interactive_trajectory_timeseries
.
You can also symbolically alter parameters
current_parameter(roessler, :c)
5.7
set_parameter!(roessler, :c, 5.0)
current_parameter(roessler, :c)
5.0
This symbolic indexing can be given anywhere in the ecosystem where you would be altering the parameters.
Core components reference
StateSpaceSets.StateSpaceSet
— TypeStateSpaceSet{D, T} <: AbstractStateSpaceSet{D,T}
A dedicated interface for sets in a state space. It is an ordered container of equally-sized points of length D
. Each point is represented by SVector{D, T}
. The data are a standard Julia Vector{SVector}
, and can be obtained with vec(ssset::StateSpaceSet)
. Typically the order of points in the set is the time direction, but it doesn't have to be.
When indexed with 1 index, StateSpaceSet
is like a vector of points. When indexed with 2 indices it behaves like a matrix that has each of the columns be the timeseries of each of the variables. When iterated over, it iterates over its contained points. See description of indexing below for more.
StateSpaceSet
also supports almost all sensible vector operations like append!, push!, hcat, eachrow
, among others.
Description of indexing
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 thei
th point (returns anSVector
)X[v1] == X[v1, :]
, returns aStateSpaceSet
with the points in those indices.X[:, j]
gives thej
th variable timeseries (or collection), asVector
X[v1, v2], X[:, v2]
returns aStateSpaceSet
with the appropriate entries (first indices being "time"/point index, while second being variables)X[i, j]
value of thej
th variable, at thei
th timepoint
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.
DynamicalSystemsBase.DynamicalSystem
— TypeDynamicalSystem
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 DynamicalSystem
represents the time evolution of a state in a state space. It mainly encapsulates three things:
- A state, typically referred to as
u
, with initial valueu0
. The space thatu
occupies is the state space ofds
and the length ofu
is the dimension ofds
(and of the state space). - A dynamic rule, typically referred to as
f
, that dictates how the state evolves/changes with time when calling thestep!
function.f
is typically a standard Julia function, see the online documentation for examples. - A parameter container
p
that parameterizesf
.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
Dynamical systems that have been constructed from DEProblem
s 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.
In ModelingToolkit.jl v9 the default split
behavior of the parameter container is true
. This means that the parameter container is no longer a Vector{Float64}
by default, which means that you cannot use integers to access parameters. It is recommended to keep split = true
(default) and only access parameters via their symbolic parameter binding. Use structural_simplify(sys; split = false)
to allow accessing parameters with integers again.
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)
withds
an instance ofDynamicalSystem
: return the state ofds
at timet
. For continuous time systems this interpolates and extrapolates, while for discrete time systems it only works ift
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
API - alter status
Dynamical system implementations
DynamicalSystemsBase.DeterministicIteratedMap
— TypeDeterministicIteratedMap <: 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.
DynamicalSystemsBase.CoupledODEs
— TypeCoupledODEs <: 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
to access the solvers. The default diffeq
is:
(alg = Tsit5(; stagelimiter! = triviallimiter!, steplimiter! = triviallimiter!, thread = 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. 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. The integrator is available as the field integ
, and the ODEProblem
is integ.sol.prob
. The convenience syntax ODEProblem(ds::CoupledODEs, tspan = (t0, Inf))
is available to extract the problem.
DynamicalSystemsBase.StroboscopicMap
— TypeStroboscopicMap <: 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
.
DynamicalSystemsBase.PoincareMap
— TypePoincareMap <: 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.
See also StroboscopicMap
, poincaresos
.
Keyword arguments
direction = -1
: Only crossings withsign(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.rootkw = (xrtol = 1e-6, atol = 1e-8)
: ANamedTuple
of keyword arguments passed tofind_zero
from Roots.jl.Tmax = 1e3
: The argumentTmax
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 onestep!
the system has been evolved for more thanTmax
, thenstep!(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 thej
th variable of the system equals the valuer
. - A vector of length
D+1
. The firstD
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:
dimension(pmap) == dimension(ds)
, even though the Poincaré map is effectively 1 dimension less.- Like
StroboscopicMap
time is discrete and counts the iterations on the surface of section.initial_time
is always0
andcurrent_time
is current iteration number. - A new function
current_crossing_time
returns the real time corresponding to the latest crossing of the hyperplane, which is what thecurrent_state(ds)
corresponds to as well. - For the special case of
plane
being aTuple{Int, <:Real}
, a specialreinit!
method is allowed with input state of lengthD-1
instead ofD
, i.e., a reduced state already on the hyperplane that is then converted into theD
dimensional state.
Example
using DynamicalSystemsBase
ds = Systems.rikitake(zeros(3); μ = 0.47, α = 1.0)
pmap = poincaremap(ds, (3, 0.0))
step!(pmap)
next_state_on_psos = current_state(pmap)
DynamicalSystemsBase.ProjectedDynamicalSystem
— TypeProjectedDynamicalSystem <: 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...
DynamicalSystemsBase.ArbitrarySteppable
— TypeArbitrarySteppable <: 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 lastreinit!
call.set_state!
is the same asreinit!
by default. If not, the keyword argumentset_state
is a functionset_state(model, u)
that sets the state of the model tou
.- The keyword
isdeterministic
should be set properly, as it decides whether downstream algorithms should error or not.
Learn more
To learn more, you need to visit the documentation pages of the modules that compose DynamicalSystems.jl. See the contents page for more!
- DatserisParlitz2022Datseris & Parlitz 2022, Nonlinear Dynamics: A Concise Introduction Interlaced with Code, Springer Nature, Undergrad. Lect. Notes In Physics