R-tipping functionality

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

@@ -5,6 +5,7 @@ pages = [
     "Tutorial" => "examples/tutorial.md",
     "Examples" => Any[
         "Defining stochastic systems" => "examples/stochastic-dynamics.md",
+        "Setting up a RateSystem" => "examples/RateSystem.md",
         "Large deviations: Maier-Stein system" => "examples/gMAM_Maierstein.md",
         "Simple gMAM: Kerr Parametric Oscillator" => "examples/sgMAM_KPO.md",
         "Transition Path Theory: Finite element method" => "examples/transition_path_theory_double_well.md",

docs/src/examples/RateSystem.md

@@ -0,0 +1,100 @@
+# Setting up a `RateSystem`: a prototypical model for R-tipping
+
+Let us explore an example of how to construct a `RateSystem` from an autonomous deterministic dynamical system (i.e. a `CoupledODEs`) and a time-dependent forcing protocol called `RateProtocol`.
+
+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+\lambda)^2 - 1
+\end{aligned}
+```
+The parameter ``\lambda`` 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]
+    λ = p[1]
+    dx = (x+λ)^2 - 1
+    return SVector{1}(dx)
+end
+
+lambda = 0.0 
+p = [lambda]
+x0 = [-1.]
+auto_sys = CoupledODEs(f,x0,p)
+```
+
+## Non-autonomous case
+
+Now, we want to explore a non-autonomous version of the system. 
+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_end]`` with a time-dependent parameter ramping given by the function ``\lambda(rt)``. Choosing different values of the parameter ``r`` allows us to vary the speed of the parameter ramping.
+
+We start by defining the function ``\lambda(t)``:
+```@example RateSystem
+function λ(p,t)
+    λ_max = p[1]
+    lambda = (λ_max/2)*(tanh(λ_max*t/2)+1)
+    return SVector{1}(lambda)
+end
+
+λ_max = 3.
+p_lambda = [λ_max] # parameter of the function lambda
+```
+
+
+We plot the function ``λ(p_lambda,t)``
+
+```@example RateSystem
+lambda_plotvals = [λ(p_lambda, t)[1] for t in -10.0:0.1:10.0]
+figλ = Figure(); axsλ = Axis(figλ[1,1],xlabel="t",ylabel=L"\lambda")
+lines!(axsλ,-10.0:0.1:10.0,lambda_plotvals,linewidth=2,label=L"\lambda(p_{\lambda}, t)")
+axislegend(axsλ,position=:rc,labelsize=10)
+figλ
+```
+
+
+Now, we define the RateProtocol that describes how and when the function ``λ(p_lambda,t)`` is used to
+make autonomous system non-autonomous:
+
+```@example RateSystem
+r = 4/3-0.02   # r just below critical rate 4/3
+t_start = -Inf # start time of non-autonomous part
+t_end = Inf    # end time of non-autonomous part
+
+rp = CriticalTransitions.RateProtocol(λ,p_lambda,r,t_start,t_end)
+```
+
+
+
+Now, we set up the combined system with autonomous past and future and non-autonomous ramping in between:
+
+```@example RateSystem
+t0 = -10.      # initial time of the system
+nonauto_sys = apply_ramping(auto_sys,rp,t0)
+```
+
+We can compute trajectories of this new system in the familiar way:
+```@example RateSystem
+T = 20.        # 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

@@ -57,6 +57,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 +83,8 @@ export MinimumActionPath
 export deterministic_orbit
 export transition, transitions
 
+export RateProtocol, apply_ramping
+
 # Error hint for extensions stubs
 function __init__()
     Base.Experimental.register_error_hint(

src/r_tipping/RateSystem.jl

@@ -0,0 +1,93 @@
+# 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
+
+Time-dependent forcing protocol specified by the following fields:
+- `λ::Function`: forcing function of the form `λ(p, t_start; kwargs...)``
+- `p_lambda::Vector`: parameters of forcing function
+- `r::Float64`: rate parameter
+- `t_start::Float64`: start time of protocol
+- `t_end::Float64`: end time of protocol
+
+If `t_start` and `t_end` are not provided, they are set to `t_start=-Inf`, and `t_end=Inf`.
+If `p_lambda` is not provided, it is set to an empty array `[]`.
+
+The forcing protocol contains all the information that allows to take an ODE of the form 
+`dxₜ/dt = f(xₜ,λ)` 
+with `λ ∈ Rᵐ` containing all the system parameters, and make it an ODE of the form 
+`dxₜ/dt = f(xₜ,λ(rt))`
+with `λ(t) ∈ Rᵐ`. 
+"""
+mutable struct RateProtocol
+    λ::Function
+    p_lambda::Vector
+    r::Float64
+    t_start::Float64
+    t_end::Float64
+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,
+    t_start::Float64,
+    t_end::Float64,
+    r::Float64;
+    kwargs...,
+)
+    if t_start > t_end
+        error("Please ensure that t_start ≤ t_end.")
+    end
+
+    p̃ = if r*t ≤ t_start
+        λ(p, t_start; kwargs...)
+    elseif t_start < r*t < t_end
+        λ(p, r*t; kwargs...) # the value(s) of λ(rt)
+    else
+        λ(p, t_end; kwargs...) # the value(s) of λ(rt)
+    end; # the value(s) of λ(rt)
+    return ds.integ.f(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`, 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_end` with the parameter shift given by the [`RateProtocol`](@def),
+and again autonomous from `t_end` to the end of the simulation:
+
+`t_0`  autonomous    `t_start`  non-autonomous   `t_end`  autonomous   `∞`
+"""
+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.t_start, rp.t_end, rp.r; kwargs...
+    )
+    prob = remake(auto_sys.integ.sol.prob; f, p=rp.p_lambda, tspan=(t0, Inf))
+    nonauto_sys = CoupledODEs(prob, auto_sys.diffeq)
+    return nonauto_sys
+end

src/r_tipping/critical_rate.jl

@@ -0,0 +1,147 @@
+using LinearAlgebra: norm
+
+## the following function integrates the nonautonomous system created by the ContinuousTimeDynamicalSystem ds and RateProtocol rp, with rate-parameter r 
+## the initial condition is e_start at time t_start
+## the system is integrated for T = t_end-t_start time units and the trajectory position at time T is then compared with the location of e_end 
+## if the trajectory position at time T is within a neighbourhood of radius rad of e_end, the function returns true, and false otherwise  
+
+function end_point_tracking(
+    ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, r::Float64, e_start::Vector, t_start::Float64, e_end::Vector, t_end::Float64, rad::Float64
+)
+
+    # note that should ensure that Δt << dλ(rt)/dt
+    # more generally, the suitability of the numerical integrator diffeq has to be ensured
+
+    T = t_end - t_start
+    rp.r = r # set the rate-parameter to the minimum value 
+    nonauto_sys = apply_ramping(ds, rp, t_start) # create the nonautonomous system
+    traj = trajectory(nonauto_sys, T, e_start) # simulate the nonautonomous system
+
+    return norm(traj[1][end,:]-e_end) < rad 
+
+end
+
+# the following function finds via brute force an interval [a,b] satisfying b-a ≤ tol satisfying f(a)*f(b) == false for a function f returning boolean values 
+# i.e. a small interval [a,b] for which the value of the function f evaluated at the points a and b returns both true and false  
+
+function binary_bisection(f::Function, a::Float64, b::Float64;
+    tol::Float64=1e-12, maxiter::Int64=10000)
+    a < b || error("Please enter a, b such that a < b.") 
+    fa = f(a)
+    fa*f(b) == false || error("The function does not satisfy the initial requirement that f(a)*f(b) == false")
+    ii = 0
+    local c
+    while b-a > tol
+        ii += 1
+        ii ≤ maxiter || error("The maximum number of iterations was exceeded.")
+        c = (a+b)/2
+        fc = f(c)
+        if fa*fc == true # reduce to the shrunken interval [c,b], since we must have f(c)*f(b) == false
+            a = c  
+            fa = fc 
+        else # reduce to the shrunken interval [a,c], since we have f(a)*f(c) == false
+            b = c  
+        end
+    end
+    return [a,b]
+end;
+
+## the following function aims to find a critical value for the rate-parameter between successful and non-successful "end-point tracking" 
+
+## the end-point tracking is currently defined with respect to a stable equilibria of the frozen-in system attained at time t = t_start...
+##...and a stable equilibria of the frozen-in system attained at time t = t_end 
+
+## there is currently no check that e_start and e_end are connected via a continuous branch of moving sinks attained across [t_start,t_end],...
+##... which would be a desirable check to perform although requires more work to implement (for a first draft I skipped this)
+## before doing this it probably makes sense to modify the moving_sinks function so as to return continuous branches of grouped equilibria...
+##...rather than just a vector of vectors containing all the (non-grouped) equilibria found for a range of frozen-in times  
+
+## there are check for ensuring that e_start and e_end are indeed stable equilibria of the frozen-in systems attained at times t_start and t_end respectively
+## there are also checks that the system succesfully end-point tracks for r = rmin and fails to end-point track for r = rmax (or vice versa, for exotic cases)
+
+## the function ultimately passes the end_point_tracking function to the binary_bisection function to find a critical rate for end-point tracking
+
+
+function critical_rate(
+    ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, e_start::Vector, t_start::Float64, e_end::Vector, t_end::Float64;  
+    rmin::Float64 = 1.e-2, 
+    rmax::Float64 = 1.e2,
+    rad::Float64 = 1.e-1,
+    tol::Float64 = 1.e-3,
+    maxiter::Int64 = Int64(1.e4)
+)
+
+    rmin < rmax || error("Please enter rmin, rmax such that rmin < rmax.") 
+
+    ##-----check no.1-----## 
+
+    ## that e_start is a stable equilibrium of the frozen-in system attained at t = t_start
+
+    λ_mod  = rp.λ
+
+    rp.λ = (p,τ) -> λ_mod(p,τ+t_start)
+    rate_sys_start = apply_ramping(ds, rp, rp.t_start)
+    box_start = [interval(x-0.1,x+0.1) for x ∈ e_start]
+    fp_start, ~, stab_start = fixedpoints(rate_sys_start, box_start)
+    if any(e -> isapprox(e,e_start;atol=1e-4),fp_start)
+        idx = findmin([norm(e-e_start) for e ∈ fp_start])[2] # the index in fp_start of the fixed point which is approximately equal to e_end
+        stab_start[idx] || @warn("The vector e_start = $(e_start) does not describe a stable equilibrium of the frozen-in system attained at t = t_start = $(t_start).")
+    else
+        @warn("The vector e_start = $(e_start) does not describe a stable equilibrium of the frozen-in system attained at t = t_start = $(t_start).")
+    end
+
+    ##-----check no.2-----##
+
+    ## that e_end is a stable equilibrium of the frozen-in system attained at t = t_end
+
+    rp.λ = (p,τ) -> λ_mod(p,τ+t_end)
+    rate_sys_end = apply_ramping(ds, rp, rp.t_start)
+    box_end = [interval(x-0.1,x+0.1) for x ∈ e_end]
+    fp_end, ~, stab_end = fixedpoints(rate_sys_end, box_end)
+    if any(e -> isapprox(e,e_end;atol=1e-4),fp_end)
+        idx = findmin([norm(e-e_end) for e ∈ fp_end])[2] # the index in fp_end of the fixed point which is approximately equal to e_end
+        stab_end[idx] || @warn("The vector e_end = $(e_end) does not describe a stable equilibrium of the frozen-in system attained at t = t_end = $(t_end).")
+    else
+        @warn("The vector e_end = $(e_end) does not describe a stable equilibrium of the frozen-in system attained at t = t_end = $(t_end).")
+    end
+
+    ##-----check no.3-----##
+
+    ## that
+    ## the system successfully end-point tracks for r = rmin and fails to end-point track for r = rmax
+    ## or
+    ## the system fails to end-point track for r = rmin and successfully end point tracks for r = rmax 
+
+    rp.λ = (p,τ) -> λ_mod(p,τ)
+
+    min_rate_track = end_point_tracking(ds,rp,rmin,e_start,t_start,e_end,t_end,rad) # true if the system end-point tracks for that rate
+
+    max_rate_track = end_point_tracking(ds,rp,rmax,e_start,t_start,e_end,t_end,rad) # true if the system end-point tracks for that rate
+
+    if min_rate_track && max_rate_track 
+        error("The system successfully end-point tracks for both r = rmin = $(rmin) and r = rmax = $(rmax).")
+    elseif !min_rate_track && !max_rate_track
+        error("The system fails to end-point track for both r = rmin = $(rmin) and r = rmax = $(rmax).")
+    elseif !min_rate_track && max_rate_track
+        @warn("The system fails to end-point track for r = rmin = $(rmin) but succesfully end-point tracks for r = rmax = $(rmax).")
+    end
+
+    ## all checks complete ##
+
+    ## performing the bisection to find a critical rate within the given tolerance ##
+
+    func(r) = end_point_tracking(ds,rp,r,e_start,t_start,e_end,t_end,rad)
+    rcrit_interval = binary_bisection(func,rmin,rmax;tol,maxiter)
+
+    return (rcrit_interval[1]+rcrit_interval[2])/2
+
+end
+
+# beginning notes
+
+# you are trying to locate the critical rate r for which the system end-point tracks the moving branch of sinks 
+
+# the first case is where after some travel time you are within a defined neighbourhood of a sink belonging to the frozen-in system at that time 
+# for this you should provide the tuple of tuples ((e₋,t₋),(e₊,t₊))
+# then you start at the position e⁻ at the time t⁻ and simulate t⁺-t⁻ time-units and compute the distance to e⁺
+# we assume that there exists some pairing (r₁,r₂) whereby the system lies within the neighbourhood of e⁺ under the rate r₁ and outside of the neighbourhood of e⁺ under the rate r₂