RateSystem

Project.toml

@@ -66,4 +66,4 @@ julia = "1.10"
 
 [extras]
 Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
-ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
+ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"

docs/Project.toml

@@ -23,4 +23,4 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 
 [compat]
-Documenter = "^1.4.1"
+Documenter = "^1.4.1"

docs/pages.jl

@@ -9,6 +9,7 @@ pages = [
         "Simple gMAM: Kerr Parametric Oscillator" => "examples/sgMAM_KPO.md",
         "Transition Path Theory: Finite element method" => "examples/transition_path_theory_double_well.md",
         "Minimal action method: Optimal Control problem" => "examples/OC_mam.md",
+        "Studying R-Tipping" => "examples/RateSystem.md",
     ],
     "Manual" => Any[
         "Define your system" => "man/CoupledSDEs.md",

docs/src/examples/RateSystem.md

@@ -0,0 +1,101 @@
+# Studying R-Tipping
+
+Let us explore a simple prototypical example of how to use the R-tipping functionality of this package.
+We start with defining an autonomous deterministic dynamical system (i.e. a `CoupledODEs`). Then, we define a time-dependent forcing protocol called `RateConfig`, describing how the parameters of the previously defined autonomous `CoupledODEs` will be changed over time. For times `t` smaller than the time `t_start`, the system is autonomous, then non-autonomous for `t_start < t < t_start + ramp_t_length` with the paramter ramping given by the `RateConfig` and for `t > t_start + ramp_t_length` the system is autonomous again. This setting is a widely used and convenient for studying R-tipping.
+
+We first consider the following simple one-dimensional autonomous system with one attractor, given by the ordinary differential equation:
+```math
+\begin{aligned}
+    \dot{x} &= (x+p)^2 - 1
+\end{aligned}
+```
+The parameter ``p`` shifts the location of the extrema of the drift field. 
+We implement this system as follows:
+
+```@example RateSystem
+using CriticalTransitions
+using CairoMakie
+
+function f(u,p,t) # out-of-place
+    x = u[1]
+    dx = (x+p[1])^2 - 1
+    return SVector{1}(dx)
+end
+
+x0 = [-1.]
+auto_sys = CoupledODEs(f,x0,[0.0]);
+```
+
+## Applying the parameter ramping
+
+Now, we want to explore a non-autonomous version of the system by applying a parameter ramping. 
+As discussed, we consider a setting where in the past and in the future the system is autnonomous and in between there is a non-autonomous period `[t_start, t_start+ramp_t_length]` with a time-dependent parameter ramping given by the function ``p(t)``. Choosing different values of the parameter `ramp_t_length` allows us to vary the speed of the parameter ramping.
+
+We start by defining the function `p(p_parameters, t)`:
+```@example RateSystem
+function p(p_parameters, t)
+    p_max = p_parameters[1]
+    p_ = (p_max/2)*(tanh(p_max*t/2)+1)
+    return SVector{1}(p_)
+end
+
+p_max = 1.0
+p_parameters = [p_max] # parameter of the function p
+```
+
+
+We plot the function `p(p_parameters, t)`
+```@example RateSystem
+p_plotvals = [p(p_parameters, t)[1] for t in -10.0:0.1:10.0]
+figp = Figure(); axsp = Axis(figp[1,1],xlabel="t",ylabel=L"p")
+lines!(axsp,-10.0:0.1:10.0,p_plotvals,linewidth=2,label=L"p(p_{parameters}, t)")
+axislegend(axsp,position=:rc,labelsize=10)
+figp
+```
+Note that this function fulfills 
+```math
+    \begin{aligned}
+        p(t=-10)& =0 \ and \ p(t=10)=1.
+    \end{aligned}
+```
+We recommend this to be fulfilled for any ramping function to use the functionality of this package.
+
+
+Now, we define a `RateConfig`, which contains all the information to apply the parameter ramping given by 
+`p(p_parameters,t)` to the `auto_sys` during `[t_start, t_start+ramp_t_length]`:
+
+```@example RateSystem
+t_start = -10       # start time of non-autonomous part
+ramp_t_length = 5   # duration of non-autonomous part
+dp=3                # strength of the paramter ramping
+
+rc = CriticalTransitions.RateConfig(p, p_parameters, t_start,ramp_t_length,dp)
+```
+Note that `dp` is defined as a prefactor of the function `p`. Thus, changing `dp` will change the amount of the parameter ramping. If the user provides a function fulfilling ``p(t=-10) =0 and p(t=10)=1``, setting `dp=10` would result in a parameter ramping from ``p(t=-10) = 0 and p(t=10) = 10``.
+
+
+We set up the system with autonomous past and future and non-autonomous ramping in between:
+
+```@example RateSystem
+t0 = -10.0      # initial time of the system
+nonauto_sys = apply_ramping(auto_sys, rc, t0);
+```
+
+We can compute trajectories of this new system in the familiar way:
+```@example RateSystem
+T = 50.0        # final simulation time
+auto_traj = trajectory(auto_sys, T, x0)
+nonauto_traj = trajectory(nonauto_sys, T, x0);
+```
+
+We plot the two trajectories
+```@example RateSystem
+fig = Figure(); axs = Axis(fig[1,1],xlabel="t",ylabel="x")
+lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}")
+lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}")
+axislegend(axs,position=:rc,labelsize=10)
+fig
+```
+
+-----
+Author: Raphael Roemer; Date: 30 Jun 2025

src/CriticalTransitions.jl

@@ -19,7 +19,8 @@ using DynamicalSystemsBase:
     set_state!,
     trajectory,
     jacobian,
-    ContinuousTimeDynamicalSystem
+    ContinuousTimeDynamicalSystem,
+    referrenced_sciml_prob
 using ConstructionBase: ConstructionBase
 using StateSpaceSets: StateSpaceSets, dimension, StateSpaceSet
 using StochasticDiffEq: StochasticDiffEq
@@ -57,6 +58,8 @@ include("largedeviations/geometric_min_action_method.jl")
 include("largedeviations/sgMAM.jl")
 include("largedeviations/string_method.jl")
 
+include("r_tipping/RateSystem.jl")
+
 # Experimental features
 include("experimental/transition_path_theory/TransitionPathMesh.jl")
 include("experimental/transition_path_theory/langevin.jl")
@@ -81,6 +84,8 @@ export MinimumActionPath
 export deterministic_orbit
 export transition, transitions
 
+export RateConfig, apply_ramping
+
 # Error hint for extensions stubs
 function __init__()
     Base.Experimental.register_error_hint(

src/r_tipping/RateSystem.jl

@@ -0,0 +1,122 @@
+# we consider the ODE dxₜ/dt = f(xₜ,p(rt))
+# here p = p(t) ∈ Rᵐ is a function containing all the system parameters 
+
+# We ask the user to define: 
+#  1) a ContinuousTimeDynamicalSystem that should be investigated and
+#  2) a protocol for the time-dependent forcing with the struct RateConfig
+
+# Then we give back the ContinuousTimeDynamicalSystem with the parameter 
+# changing according to the rate protocol
+"""
+    RateConfig
+
+Time-dependent forcing protocol containing the information to apply a parameter shift to an autonomous system.
+
+Fields
+==============
+
+- p: monotonic function which describes the time-dependent parametric forcing 
+- p_parameters: the vector of parameters which are associated with p
+- t_pstart: the parameter values of the past limit system are given by p(p_parameteters,t_pstart)
+- t_pend: the parameter values of the future limit system are given by p(p_parameters,t_pend) 
+- t_start: the explicit value of time at which the nonautonomous dynamics starts 
+- t_ramp_length: the duration of time of nonautonomous dynamics [i.e. within (t_start,t_start+t_ramp_length)] 
+- dp: the difference in parameter values attained across the ramping
+
+Default values 
+==============
+
+- t_pstart = -100
+- t_pend = 100
+- p_parameters = []
+- t_start = -t_ramp_length/2
+- dp = 1
+"""
+mutable struct RateConfig
+    p::Function  
+    p_parameters::Vector
+    t_pstart::Float64
+    t_pend::Float64  
+    t_start::Float64
+    t_ramp_length::Float64
+    dp::Float64
+end
+
+## convenience functions
+
+RateConfig(p::Function, t_ramp_length::Float64) = RateConfig(p, [], -100.0, 100.0, -t_ramp_length/2, t_ramp_length, 1.0)
+RateConfig(p::Function, t_ramp_length::Float64, dp::Float64) = RateConfig(p, [], -100.0, 100.0, -t_ramp_length/2, t_ramp_length, dp)
+
+RateConfig(p::Function, p_parameters::Vector, t_ramp_length::Float64) = RateConfig(p, p_parameters, -100.0, 100.0, -t_ramp_length/2, t_ramp_length, 1.0)
+RateConfig(p::Function, p_parameters::Vector, t_ramp_length::Float64, dp::Float64) = RateConfig(p, p_parameters, -100.0, 100.0, -t_ramp_length/2, t_ramp_length, dp)
+
+## the following function creates a piecewise constant function with respect to t_pstart and t_pend...
+## ...which is written such that p̃(t_pend)-p̃(t_pstart) = dp, for the time-dependent entries (otherwise p̃(t_pend)-p̃(t_pstart) = 0 as p̃(t)≡c)
+
+function p_piecewise_scaled(p::Function,p_parameters::Vector,t_pstart::Float64,t_pend::Float64,dp::Float64,t::Float64)
+    if t ≤ t_pstart
+        return p(p_parameters,t_pstart)
+    else
+        np = length(p(p_parameters,t_pstart)) # the number of system parameters 
+        differences = p(p_parameters,t_pend) .- p(p_parameters,t_pstart)
+        nonautonomous_inds = [ii for ii ∈ 1:np if differences[ii] != 0]
+        inds = [ii ∈ nonautonomous_inds ? 1/abs(differences[ii]) : 0 for ii ∈ 1:np] # takes value 1/|p(...,t_pend)-p(...,t_start)| when p(...,t_pend)=/=p(...,t_start); zero otherwise
+        if t < t_pend
+            return p(p_parameters,t_pstart) .+ dp .* inds .* (p(p_parameters,t) .- p(p_parameters,t_pstart))
+        else 
+            return p(p_parameters,t_pstart) .+ dp .* inds .* (p(p_parameters,t_pend) .- p(p_parameters,t_pstart))
+        end
+    end
+end
+
+function modified_drift(
+    u,
+    p_parameters,
+    t,
+    ds::ContinuousTimeDynamicalSystem,
+    p::Function,
+    t_pstart::Float64,
+    t_pend::Float64,
+    t_start::Float64,
+    t_ramp_length::Float64,
+    dp::Float64;
+    kwargs...,
+)
+
+    q(t) = p_piecewise_scaled(p,p_parameters,t_pstart,t_pend,dp,t)
+
+    time_shift = ((t_pend-t_pstart)/t_ramp_length)*(t-t_start)+t_pstart # such that [t_start,t_start+t_ramp_length] is shifted into [t_pstart,t_pend]
+
+    q̃ = q(time_shift)
+
+    return dynamic_rule(ds)(u, q̃, t)
+end;
+
+"""
+    apply_ramping(sys::CoupledODEs, rc::RateConfig, t0=0.0; kwargs...)
+
+Applies a time-dependent [`RateConfig`](@def) to a given autonomous deterministic dynamical system
+`sys`, turning it into a non-autonomous dynamical system. The returned [`CoupledODEs`](@ref)
+has the explicit parameter time-dependence incorporated.
+
+The returned [`CoupledODEs`](@ref) is autonomous from `t_0` to `t_start`, 
+then non-autnonmous from `t_start` to `t_start+t_ramp_length` with the parameter shift given by the [`RateConfig`](@def),
+and again autonomous from `t_start+t_ramp_length` to the end of the simulation:
+
+`t_0`  autonomous    `t_start`  non-autonomous   `t_start+t_ramp_length`  autonomous   `∞`
+
+Computing trajectories of the returned [`CoupledODEs`](@ref) can then be done in the same way as for any other [`CoupledODEs`](@ref).
+"""
+function apply_ramping(auto_sys::CoupledODEs, rc::RateConfig, t0=0.0; kwargs...)
+    # we wish to return a continuous time dynamical system with modified drift field
+
+    f(u, p_parameters, t) = modified_drift(
+        u, p_parameters, t, auto_sys, rc.p, rc.t_pstart, rc.t_pend, rc.t_start, rc.t_ramp_length, rc.dp; kwargs...
+    )
+    prob = remake(referrenced_sciml_prob(auto_sys); f, p=rc.p_parameters, tspan=(t0, Inf))
+    nonauto_sys = CoupledODEs(prob, auto_sys.diffeq)
+    return nonauto_sys
+end
+
+
+

src/r_tipping/RateSystemTimeInterval.jl

@@ -0,0 +1,95 @@
+# we consider the ODE dxₜ/dt = f(xₜ,λ(rt))
+# here λ = λ(t) ∈ Rᵐ is a function containing all the system parameters 
+
+# We ask the user to define: 
+#  1) a ContinuousTimeDynamicalSystem that should be investigated and
+#  2) a protocol for the time-dependent forcing with the struct RateProtocol
+
+# Then we give back the ContinuousTimeDynamicalSystem with the parameter 
+# changing according to the rate protocol
+"""
+    RateProtocol
+
+The RateProtocol contains all the fields (information) that allows to take an ODE of the form 
+    `dxₜ/dt = f(xₜ, λ)` 
+with `λ ∈ Rᵐ` containing all system parameters, and make it an ODE of the form 
+    `dxₜ/dt = f(xₜ, λ(rt))`
+with `λ(t) ∈ Rᵐ` constant before time `time_interval[1]` and also after time `time_interval[2]`. 
+
+RateProtocol contains the following fields:
+- `λ::Function`: forcing function of the form `λ(p, t_start; kwargs...)``
+- `p_lambda::Vector`: parameters of the forcing function
+- `r::Float64`: rate parameter
+- `time_interval::Tuple`: start and end time of protocol
+
+Default values 
+==============
+
+time_interval = (-Inf, Inf)
+p_lambda = []
+"""
+mutable struct RateProtocol
+    λ::Function
+    p_lambda::Vector
+    r::Float64
+    time_interval::Tuple
+end
+
+# convenience functions
+
+function RateProtocol(λ::Function, p_lambda::Vector, r::Float64)
+    RateProtocol(λ, p_lambda, r, (-Inf, Inf))
+end
+RateProtocol(λ::Function, r::Float64)=RateProtocol(λ, [], r, (-Inf, Inf))
+#RateProtocol(λ::Function,p_lambda::Vector,r::Float64,t_start::Float64)=RateProtocol(λ,p_lambda,r,t_start,Inf)
+#RateProtocol(λ::Function,r::Float64,t_start::Float64)=RateProtocol(λ,[],r,t_start,Inf)	
+
+function modified_drift(
+    u,
+    p,
+    t,
+    ds::ContinuousTimeDynamicalSystem,
+    λ::Function,
+    time_interval::Tuple,
+    r::Float64;
+    kwargs...,
+)
+    if time_interval[1] > time_interval[2]
+        error("Please ensure that time_interval[1] ≤ time_interval[2].")
+    end
+
+    p̃ = if r*t ≤ time_interval[1]
+        λ(p, time_interval[1]; kwargs...)
+    elseif time_interval[1] < r*t < time_interval[2]
+        λ(p, r*t; kwargs...) # the value(s) of λ(rt)
+    else
+        λ(p, time_interval[2]; kwargs...) # the value(s) of λ(rt)
+    end; # the value(s) of λ(rt)
+    return dynamic_rule(ds)(u, p̃, t)
+end;
+
+"""
+    apply_ramping(sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs...)
+
+Applies a time-dependent [`RateProtocol`](@def) to a given autonomous deterministic dynamical system `sys`, 
+returning a non-autonomous dynamical system of type [`CoupledODEs`](@ref).
+
+The [`RateProtocol`](@def) replaces the parameters of `sys` by the function `λ(rt)` within the 
+time interval `time_interval`. Thus, the returned [`CoupledODEs`](@ref) has the explicit parameter time-dependence incorporated and is 
+autonomous from `t_0` to `time_interval[1]`, non-autnonmous from `time_interval[1]` to `time_interval[2]` with the parameter shift given by the [`RateProtocol`](@def),
+and autonomous from `time_interval[2]` to the end of the simulation:
+
+`t_0`  autonomous    `time_interval[1]`  non-autonomous   `time_interval[2]`  autonomous   `∞`
+
+Computing trajectories of the returned [`CoupledODEs`](@ref) can then be done in the same way as for any other [`CoupledODEs`](@ref).
+"""
+function apply_ramping(auto_sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs...)
+    # we wish to return a continuous time dynamical system with modified drift field
+
+    f(u, p, t) = modified_drift(
+        u, p, t, auto_sys, rp.λ, rp.time_interval, rp.r; kwargs...
+    )
+    prob = remake(referrenced_sciml_prob(auto_sys); f, p=rp.p_lambda, tspan=(t0, Inf))
+    nonauto_sys = CoupledODEs(prob, auto_sys.diffeq)
+    return nonauto_sys
+end

systems/CTLibrary.jl

@@ -13,8 +13,8 @@ include("stommel.jl")
 include("rivals.jl")
 
 export fitzhugh_nagumo!, fitzhugh_nagumo, stommel, rivals!, rivals, cessi, rooth_smooth
-modifiedtruscottbrindleywithdimensions!, modifiedtruscottbrindleywithdimensions
-originaltruscottbrindley!, originaltruscottbrindley
-rampedoriginaltruscottbrindley!, rampedoriginaltruscottbrindley
+#modifiedtruscottbrindleywithdimensions!, modifiedtruscottbrindleywithdimensions
+#originaltruscottbrindley!, originaltruscottbrindley
+#rampedoriginaltruscottbrindley!, rampedoriginaltruscottbrindley
 
 end