Attempts at reproducing phase-tipping plots of https://arxiv.org/abs/2101.12107

hassan_alkhayuon_statespace.jl

using DynamicalSystems, PyPlot, LinearAlgebra, OrdinaryDiffEq
# ODE solver options:
diffeq = (alg = Vern9(), reltol = 1e-10, abstol = 1e-10)

# non-dimensionalized Rozenweig-MacArthur model
function rozenweig_macarthur(u, p, t)
    k, m, c = p
    x, y = u
    common = m*x*y/(1+x)
    xdot = x*(1-x/k) - common
    ydot = - c*y + common
    return SVector(xdot, ydot)
end

# From MATLAB of Hassan Alkhayuon, specifically:
# https://github.com/hassanalkhayuon/phaseSensitiveTipping/blob/main/RM_model/BI_strip_RM.m
function hassan_alkhayuon(u, p, t)
    N, P = u
    RR = p[1]
    delta      = 2.2;
    C          = 0.19;
    gamma      = 0.004;
    beta       = 1.5;
    alpha      = 800;
    mu         = 0.03;
    nu         = 0.003;

    dN = RR * N *(1-((C/RR)*N))*((N - mu)/(nu + N)) - ((alpha*P*N)/(beta + N));
    dP = gamma * ((alpha*P*N)/(beta + N)) - (delta*P);
    return SVector(dN, dP)
end

# Toggle either of the two definitions to get either model.
# predator_pray = (u, p, t) -> rozenweig_macarthur(u, p, t)
predator_pray = (u, p, t) -> hassan_alkhayuon(u, p, t)

# This function creates a state space plot 
function statespaceplot(p)
    r, c, μ, ν, α, β, χ, δ = p
    # gx, gy are the grid in state space to make the streamplot of
    gx = range(0, 10; length = 10)
    gy = range(0, 15*1e3; length = 10) # using 15 or 15*1e3 makes no difference.
    vx = zeros(length(gx), length(gy))
    vy = copy(vx)
    # Make dynamical system, to use in `trajectory` to solve the ODE.
    ds = ContinuousDynamicalSystem(predator_pray, SVector(1.0, 1.0), p)
    figure(figsize = (10, 10))
    for (i, x) in enumerate(gx)
        for (j, y) in enumerate(gy)
            u0 = SVector(x, y)
            vx[i,j], vy[i,j] = predator_pray(u0, p, 0)
            # Also evolve initial condition `u0`, solving ODE:
            tr = trajectory(ds, 100.0, u0; diffeq...)
            plot(columns(tr)...; lw = 1.0)
        end
    end
    # plot state space stream plot:
    streamplot(Array(gx), Array(gy), vx', vy'; linewidth = 1.5)
end

# r, c, μ, ν, α, β, χ, δ = p
p = [2.47, 0.19, 0.03, 0.003, 800, 1.5, 0.004, 2.2]
statespaceplot(p) # only the first entry of the parameter container matters.