Tutorial

To give you an idea of how our package works, this tutorial provides some example code with explanations.

Example: FitzHugh-Nagumo model

Let's consider a simple 2-dimensional dynamical system - the FitzHugh-Nagumo model:

\[\begin{aligned} \frac{du}{dt} &= \frac{1}{\epsilon} \left( -\alpha u^3 + \gamma u - \kappa v + I \right) \\ \frac{dv}{dt} &= -\beta v + u \ , \end{aligned}\]

where $\epsilon$ is the parameter of time scale separation between the state variables $u$ and $v$. The parameters $\alpha >0$, $\beta >1$, $\gamma>0$, and $\kappa>0$ are real constants, and $I$ denotes a driving term.

Let's investigate this system under stochastic forcing.

System definition

First, we need to translate the system equations above into Julia code.

This works exactly as in DynamicalSystems.jl by defining a function f(u,p,t) which takes as input a vector u of state variables ($u$,$v$), a vector p of parameters, and time t. The function must return an array of flow increments ($\text{d}u$, $\text{d}v$). For performance reasons, it is advisable to return a StaticArray SA[du, dv] rather than just a Vector [du, dv].

using CriticalTransitions

function fitzhugh_nagumo(u,p,t)
    u, v = u
    ϵ, β, α, γ, κ, I = p

    du = (-α*u^3 + γ*u - κ*v + I)/ϵ
    dv = -β*v + u

    SA[du, dv]
end

fitzhugh_nagumo (generic function with 1 method)

CoupledSDE

Next, we construct a stochastic system with the fitzhugh_nagumo equation as the deterministic part. Suppose we would like to force both state variables $u$ and $v$ with additive, uncorrelated Gaussian noise of intensity $\sigma$. This is the default case. We simply write

p = [1., 3., 1., 1., 1., 0.] # Parameters (ϵ, β, α, γ, κ, I)
σ = 0.2 # noise strength

# CoupledSDE
sys = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=σ)

2-dimensional CoupledSDEs
 deterministic: false
 discrete time: false
 in-place:      false
 dynamic rule:  fitzhugh_nagumo
 SDE solver:    SOSRA
 SDE kwargs:    (abstol = 0.01, reltol = 0.01, dt = 0.1)
 Noise type:    (additive = true, autonomous = true, linear = true, invertible = true)
 parameters:    [1.0, 3.0, 1.0, 1.0, 1.0, 0.0]
 time:          0.0
 state:         [0.0, 0.0]

Here the first field fitzhugh_nagumo specifies the deterministic dynamics f (see Define a CoupledSDEs system). We have chosen zeros(2) as the initial state of the system, which is the second field. The length of this vector must match the system's dimensionality. In the (optional) third field, we specify the parameter vector p, which includes the parameters of f followed by the parameters of g (in this case, there are no parameters for g). Lastly, noise_strength sets the noise strength. Since we have not specified a noise process, the default case of an uncorrelated Wiener process is used.

That's it! Now we can apply the toolbox of CriticalTransitions to our stochastic FitzHugh-Nagumo system sys.

Find stable equilibria

For the parameters chosen above, the FitzHugh-Nagumo system is bistable. Let's compute the fixed points using the fixedpoints function. This function is borrowed from ChaosTools.jl and is loaded as an extension when we write using ChaosTools.

using ChaosTools
# Calculate fixed points and store the stable ones
eqs, eigs, stab = fixedpoints(sys, [-2,-2], [2,2])
fp1, fp2 = eqs[stab]

2-dimensional StateSpaceSet{Float64} with 2 points
 -0.816497  -0.272166
  0.816497   0.272166

Stochastic simulation

Using the trajectory function, we now run a simulation of our system for 100 time units starting out from the fixed point fp1:

sim = trajectory(sys, 100, fp1)

(2-dimensional StateSpaceSet{Float64} with 1001 points, 0.0:0.1:100.0)

In the keyword arguments, we have specified at which interval the solution is saved. Further keyword arguments can be used to change the solver (the default is SOSRA() for stochastic integration) and other settings.

The simulated trajectory is stored in sim in the usual output format of the solve method of DifferentialEquations.jl, including the solution sim.u and the vector of time points sim.t. The solution can also be accessed as a matrix sim[i, t], where i is the i-th component of u and t the time index.

Let's plot the result. Did the trajectory transition to the other attractor?

using Plots
plt = plot(sim[1][:,1], sim[1][:,2]; xlabel="u", ylabel="v", legend=false)
scatter!([fp1[1], fp2[1]], [fp1[2], fp2[2]], color=:red, markersize=4)
xlims!(-1.2, 1.2)
ylims!(-0.6, 0.6)
plt

Hopefully, this helped you to get started. For more info, check out the Manual section of these docs.