GitHub

Defining a forced dynamical system

This section explains how to specify your dynamical system and forcing of interest.

To specify a system, CriticalTransitions provides three core system types:

  • CoupledODEs - used to define a deterministic system 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 a system 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$.
  • RateSystem - used to define a non-autonomous system with parametric forcing of the form $\frac{\text{d}\mathbf{u}}{\text{d}t} = \mathbf{f}(\mathbf{u}(t), p(t))$.

The CoupledODEs and CoupledSDEs system types are inherited from DynamicalSystemsBase.jl. The RateSystem type is added in CriticalTransitions.jl to enable easy construction of non-autonomous dynamical systems in which a parameter changes over time.

Deterministic: CoupledODEs

DynamicalSystemsBase.CoupledODEs — Type
CoupledODEs <: 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.

source

Stochastic: CoupledSDEs

DynamicalSystemsBase.CoupledSDEs — Type
CoupledSDEs <: 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 to 1.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 size length(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 to nothing). 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 (default nothing)
  • noise_strength: scaling factor for noise strength (default 1.0)
  • covariance: noise covariance matrix (default nothing)
  • noise_prototype: prototype instance for the output of g (default nothing)
  • noise_process: stochastic process as provided by DiffEqNoiseProcess.jl (default nothing, i.e. standard Wiener process)
  • t0: initial time (default 0.0)
  • diffeq: DiffEq solver settings (see below)
  • seed: random number seed (default UInt64(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.

source
Tip

Check out some examples of how to construct different types of stochastic dynamics here.

Info

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 — Function
solver(ds::ContinuousTimeDynamicalSystem) -> Any

Returns the SDE solver specified in the diffeq settings of the CoupledSDEs.

source
CriticalTransitions.drift — Function
drift(
    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).

source
CriticalTransitions.div_drift — Function
div_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).

source
StochasticSystemsBase.diffusion_matrix — Function
diffusion_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$.

source
CriticalTransitions.noise_process — Function
noise_process(sys::CoupledSDEs) -> Any

Fetches the stochastic process $\mathcal{N}$ specified in the intergrator of sys. Returns the type DiffEqNoiseProcess.NoiseProcess.

source

Non-autonomous: RateSystem

CriticalTransitions.RateSystem — Type
RateSystem <: ContinuousTimeDynamicalSystem
RateSystem(ds::ContinuousTimeDynamicalSystem, fp::ForcingProfile, pidx; kwargs...)

A non-autonomous dynamical system constructed by applying a time-dependent parametric forcing protocol (ForcingProfile) to an underlying autonomous continuous-time dynamical system ds.

The ForcingProfile, applied to the parameter with index pidx, describes the functional form of the parameter evolution over a given interval. This interval is rescaled in system time units by defining its start time (forcing_start_time) and duration (forcing_duration). The magnitude of the forcing can be adjusted by scaling the forcing profile by a factor forcing_scale.

The integer pidx refers to the relevant index of the parameter container p of ds.

Before t = forcing_start_time, the system parameters correspond to that of the underlying autonomous system ds. At the end of the forcing interval (t = forcing_start_time + forcing_duration), the parameters are frozen at their current value and thereafter the system is again autonomous.

Keyword arguments

  • forcing_start_time = fp.interval[1]: Time when parameter change starts
  • forcing_duration = fp.interval[2] - fp.interval[1]: Duration of the parameter change (in system time units)
  • forcing_scale = 1.0: Amplitude multiplication factor of the forcing protocol
  • t0 = 0.0: Initial time of the system
source
CriticalTransitions.ForcingProfile — Type
ForcingProfile(profile::Function, interval)

Time-dependent forcing profile $p(t)$ describing the evolution of a parameter over a domain interval interval = (start, end). Used to define a parametric forcing when constructing a non-autonomous RateSystem.

Keyword arguments

  • profile: function $p(t)$ from $ℝ → ℝ$ describing the parameter change
  • interval: domain interval over which the profile is considered

Note: The units of t are arbitrary; the forcing profile is rescaled in system time units and magnitude when using it to construct a RateSystem.

Convenience functions

  • ForcingProfile(:linear): Creates a linear ramp from 0 to 1.
  • ForcingProfile(:tanh): Creates a hyperbolic tangent ramping from 0 to 1.
source

Converting between systems

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.

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