Defining a forced dynamical system
This section explains how to specify your dynamical system and forcing of interest.
To specify a system, CriticalTransitions
imports two core system types from DynamicalSystemsBase.jl
:
- CoupledODEs - used to define deterministic systems of ordinary differential equations of the form $\frac{\text{d}\mathbf{u}}{\text{d}t} = \mathbf{f}(\mathbf{u}, p, t)$.
- CoupledSDEs - used to define systems of stochastic differential equations of the form $\text{d}\mathbf{u} = \mathbf{f}(\mathbf{u}, p, t) \text{d}t + \mathbf{g}(\mathbf{u}, p, t) \text{d}\mathcal{N}_t$.
Both CoupledODEs
and CoupledSDEs
support nonautonomous (i.e. explicitly time-dependent) dynamics. A time-dependent parameter change can easily be specified:
RateProtocol
- allows to specify the time evolution of a system parameter.apply_ramping
- used to apply aRateProtocol
to a given dynamical system.
CoupledODEs
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
(or smaller library package such as OrdinaryDiffEqVerner
) to access the solvers. The default diffeq
is:
(alg = OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.triviallimiter!), typeof(OrdinaryDiffEqCore.triviallimiter!), Static.False}(OrdinaryDiffEqCore.triviallimiter!, OrdinaryDiffEqCore.triviallimiter!, 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. Use ODEProblem(ds::CoupledODEs, tspan = (t0, Inf))
to obtain the problem.
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.
CoupledSDEs
DynamicalSystemsBase.CoupledSDEs
— TypeCoupledSDEs <: ContinuousTimeDynamicalSystem
CoupledSDEs(f, u0 [, p]; kwargs...)
A stochastic continuous time dynamical system defined by a set of coupled stochastic differential equations (SDEs) as follows:
\[\text{d}\mathbf{u} = \mathbf{f}(\mathbf{u}, p, t) \text{d}t + \mathbf{g}(\mathbf{u}, p, t) \text{d}\mathcal{N}_t\]
where $\mathbf{u}(t)$ is the state vector at time $t$, $\mathbf{f}$ describes the deterministic dynamics, and the noise term $\mathbf{g}(\mathbf{u}, p, t) \text{d}\mathcal{N}_t$ describes the stochastic forcing in terms of a noise function (or diffusion function) $\mathbf{g}$ and a noise process $\mathcal{N}_t$. The parameters of the functions $\mathcal{f}$ and $\mathcal{g}$ are contained in the vector $p$.
There are multiple ways to construct a CoupledSDEs
depending on the type of stochastic forcing.
The only required positional arguments are the deterministic dynamic rule f(u, p, t)
, the initial state u0
, and optinally the parameter container p
(by default p = nothing
). For construction instructions regarding f, u0
see the DynamicalSystems.jl tutorial .
By default, the noise term is standard Brownian motion, i.e. additive Gaussian white noise with identity covariance matrix. To construct different noise structures, see below.
Noise term
The noise term can be specified via several keyword arguments. Based on these keyword arguments, the noise function g
is constructed behind the scenes unless explicitly given.
- The noise strength (i.e. the magnitude of the stochastic forcing) can be scaled with
noise_strength
(defaults to1.0
). This factor is multiplied with the whole noise term. - For non-diagonal and correlated noise, a covariance matrix can be provided via
covariance
(defaults to identity matrix of sizelength(u0)
.) - For more complicated noise structures, including state- and time-dependent noise, the noise function
g
can be provided explicitly as a keyword argument (defaults tonothing
). For construction instructions, continue reading.
The function g
interfaces to the diffusion function specified in an SDEProblem
of DynamicalSystems.jl. g
must follow the same syntax as f
, i.e., g(u, p, t)
for out-of-place (oop) and g!(du, u, p, t)
for in-place (iip).
Unless g
is of vector form and describes diagonal noise, a prototype type instance for the output of g
must be specified via the keyword argument noise_prototype
. It can be of any type A
that has the method LinearAlgebra.mul!(Y, A, B) -> Y
defined. Commonly, this is a matrix or sparse matrix. If this is not given, it defaults to nothing
, which means the g
should be interpreted as being diagonal.
The noise process can be specified via noise_process
. It defaults to a standard Wiener process (Gaussian white noise). For details on defining noise processes, see the docs of DiffEqNoiseProcess.jl . A complete list of the pre-defined processes can be found here. Note that DiffEqNoiseProcess.jl
also has an interface for defining custom noise processes.
By combining g
and noise_process
, you can define different types of stochastic systems. Examples of different types of stochastic systems are listed on the StochasticDiffEq.jl tutorial page. A quick overview of common types of stochastic systems can also be found in the online docs for CoupledSDEs
.
Keyword arguments
g
: noise function (defaultnothing
)noise_strength
: scaling factor for noise strength (default1.0
)covariance
: noise covariance matrix (defaultnothing
)noise_prototype
: prototype instance for the output ofg
(defaultnothing
)noise_process
: stochastic process as provided by DiffEqNoiseProcess.jl (defaultnothing
, i.e. standard Wiener process)t0
: initial time (default0.0
)diffeq
: DiffEq solver settings (see below)seed
: random number seed (defaultUInt64(0)
)
DifferentialEquations.jl interfacing
The CoupledSDEs
is evolved using the solvers of DifferentialEquations.jl. To specify a solver via the diffeq
keyword argument, use the flag alg
, which can be accessed after loading StochasticDiffEq.jl (using StochasticDiffEq
). The default diffeq
is:
(alg = SOSRA(), abstol = 1.0e-2, reltol = 1.0e-2)
diffeq
keywords can also include a callback
for event handling .
Dev note: CoupledSDEs
is a light wrapper of SDEIntegrator
from StochasticDiffEq.jl. The integrator is available as the field integ
, and the SDEProblem
is integ.sol.prob
. The convenience syntax SDEProblem(ds::CoupledSDEs, tspan = (t0, Inf))
is available to extract the problem.
Converting between CoupledSDEs
and CoupledODEs
You can convert a CoupledSDEs
system to CoupledODEs
to analyze its deterministic part using the function CoupledODEs(ds::CoupledSDEs; diffeq, t0)
. Similarly, use CoupledSDEs(ds::CoupledODEs, p; kwargs...)
to convert a CoupledODEs
into a CoupledSDEs
.
Check out some examples of how to construct different types of stochastic dynamics here.
Note that nonlinear mixings of the Noise Process $\mathcal{W}$ fall into the class of random ordinary differential equations (RODEs) which have a separate set of solvers. See this example of DifferentialEquations.jl.
Accessing CoupledSDEs
properties
CriticalTransitions.solver
— Functionsolver(ds::ContinuousTimeDynamicalSystem) -> Any
Returns the SDE solver specified in the diffeq
settings of the CoupledSDEs
.
CriticalTransitions.drift
— Functiondrift(
sys::Union{CoupledODEs{IIP}, CoupledSDEs{IIP}},
x;
t
) -> Any
Returns the deterministic drift $f(x)$ of the CoupledSDEs sys
at state x
. For time-dependent systems, the time can be specified as a keyword argument t
(by default t=0
).
CriticalTransitions.div_drift
— Functiondiv_drift(sys::ContinuousTimeDynamicalSystem, x) -> Any
div_drift(sys::ContinuousTimeDynamicalSystem, x, t) -> Any
Computes the divergence of the drift field $f(x)$ at state x
. For time- dependent systems, the time can be specified as a keyword argument t
(by default t=0
).
StochasticSystemsBase.covariance_matrix
— Functioncovariance_matrix(ds::CoupledSDEs)
Returns the covariance matrix of the stochastic term of the CoupledSDEs
ds
, provided that the diffusion function g
can be expressed as a constant invertible matrix. If this is not the case, returns nothing
.
See also diffusion_matrix
.
StochasticSystemsBase.diffusion_matrix
— Functiondiffusion_matrix(ds::CoupledSDEs)
Returns the diffusion matrix of the stochastic term of the CoupledSDEs
ds
, provided that the diffusion function g
can be expressed as a constant invertible matrix. If this is not the case, returns nothing
.
Note: The diffusion matrix $Σ$ is the square root of the noise covariance matrix $Q$ (see covariance_matrix
), defined via the Cholesky decomposition $Q = Σ Σ^\top$.
CriticalTransitions.noise_process
— Functionnoise_process(sys::CoupledSDEs) -> Any
Fetches the stochastic process $\mathcal{N}$ specified in the intergrator of sys
. Returns the type DiffEqNoiseProcess.NoiseProcess
.
Non-autonomous systems
RateProtocol apply_ramping
Converting between system types
The deterministic part of a CoupledSDEs
can be extracted as a CoupledODEs
, making it compatible with functionality of DynamicalSystems.jl
. In turn, a CoupledODEs
can easily be extended into a CoupledSDEs
.
DynamicalSystemsBase.CoupledODEs
— MethodCoupledODEs(ds::CoupledSDEs; kwargs...)
Converts a CoupledSDEs
into a CoupledODEs
by extracting the deterministic part of ds
.
DynamicalSystemsBase.CoupledSDEs
— MethodCoupledSDEs(ds::CoupledODEs, p; kwargs...)
Converts a CoupledODEs
system into a CoupledSDEs
.
For example, the Lyapunov spectrum of a CoupledSDEs
in the absence of noise, here exemplified by the FitzHugh-Nagumo model, can be computed by typing:
using CriticalTransitions
using ChaosTools: lyapunovspectrum
function fitzhugh_nagumo(u, p, t)
x, y = u
ϵ, β = p
dx = (-x^3 + x - y)/ϵ
dy = -β*y + x
return SVector{2}([dx, dy])
end
p = [1.,3.]
# Define the CoupledSDEs
sys = CoupledSDEs(fitzhugh_nagumo, ones(2), p; noise_strength=0.1)
# Compute Lyapunov spectrum of deterministic flow
ls = lyapunovspectrum(CoupledODEs(sys), 10000)
2-element Vector{Float64}:
0.0037846226653203657
-0.028280268988271094