I am trying to construct a bifurcation diagram in Julia using Bifurcation.Kit package. I cannot get the graph I am expect to produce. For some context, I am modeling a system of ODEs (non linear, 6 equations, and 6 unknowns). I know that there are two alternative stable states, and therefore, I should expect a fold bifurcation. However, that is not showing up in my graph. I could really use any code debugging, insight, directions, or contacts of people who might be able to help.
I will be happy to clarify any aspect.
#Packages
using Revise, Parameters, Setfield, Plots, LinearAlgebra
using BifurcationKit
#Supremum Norm
norminf = x -> norm(x, Inf)
#Constants
const m = 0.250 #Mortality P and J
s = 0.500 #Stocking P
f = 1.000 #Fecundity
c = 0.100 #cannibalism
kₓ = 0.250 #Attack rate F on J in littoral
k_y = 0.500 #Attack rate F on J in pelagic
r = 0.500 #Forage fish growth rate
b = 0.050 #Density dependent constant forage fish
hₓ = 0.250 #Attack rate P on F in littoral
h_y = 0.500 #Attack rate P on F in pelagic
i = 1.000 #Stocking F
ϵ = 0.250 #Trophic Efficiency
γ = 0.500 #Area Ratio
σ₁ = 1.000 #Dispersal
σ₂ = 1.000 #Dispersal
σ₃ = 1.000 #Dispersal
σ₄ = 1.000 #Dispersal
σ₅ = 1.000 #Dispersal
σ₆ = 1.000 #Dispersal
tmat = 5.000 #Maturation
#Defining a vector field
function VF!(dz, z, p, t)
p₁, p₂ = p #Parameters
u₁, u₂, u₃, u₄, u₅, u₆ = z #State Variables
#Predator Equations
dz[1] = u₅/tmat - m*u₁ - p₁*u₁ + s + ϵ*hₓ*u₁*u₃ + σ₁*γ*p₂*u₂ - σ₂*(1/γ)*p₂*u₁ #Littoral
dz[2] = u₆/tmat - m*u₂ - p₁*u₂ + s + ϵ*h_y*u₂*u₄ - σ₁*γ*p₂*u₂ + σ₂*(1/γ)*p₂*u₁ #Pelagic
#Foragers Equations
dz[3] = r*u₃ - b*u₃*u₃ - hₓ*u₁*u₃ + i + kₓ*u₃*u₅ + σ₃*γ*p₂*u₄ - σ₄*(1/γ)*p₂*u₃ #Littoral
dz[4] = r*u₄ - b*u₄*u₄ - h_y*u₂*u₄ + i + k_y*u₄*u₆ - σ₃*γ*p₂*u₄ + σ₄*(1/γ)*p₂*u₃ #Pelagic
#Juveniles Equations
dz[5] = f*(u₁+γ*u₂) - u₅/tmat - m*u₅ - c*u₁*u₅ - kₓ*u₃*u₅ + σ₅*γ*p₂*u₆ - σ₆*(1/γ)*p₂*u₅ #Littoral
dz[6] = -u₆/tmat - m*u₆ - c*u₂*u₆ - k_y*u₄*u₆ - σ₅*γ*p₂*u₆ + σ₆*(1/γ)*p₂*u₅ #Pelagic
dz
end
#Initiliazing State Variables
u₀ = zeros(Float64,6)
u₀[1] = 10 #Predator; Littoral
u₀[2] = 10 #Predator; Palagic
u₀[3] = 10 #Forager; Littoral
u₀[4] = 10 #Forager; Palagic
u₀[5] = 10 #Juveniles; Littoral
u₀[6] = 10 #Juveniles; Palagic
#Initiliazing Parameters
p = zeros(Float64,2)
p[1] = 0.0 #Fishing effort
p[2] = 0.0 #Dispersal Rate
#Unclear what this step does
VF(z, p) = VF!(similar(z), z, p, 0) #Insert explanation here?????????
using OrdinaryDiffEq
#Defining Time Interval
tᵢ = 0.0 #Initial Time
t_f = 1000 #Final Time
tspan = (tᵢ,t_f)
#Defining ODE problem
prob = ODEProblem(
VF!, #Vector Field
u₀, #Initial Conditions State Variables
tspan, #Time Interval
p) #Initial Conditions Parameters
sol = solve(prob,Tsit5()) #Tsit5 --> Sitff ODE Solver
#Plotting Solution
solution = solve(prob, Rodas5())
plot(
solution,
tspan=(0,1000),
layout=(3,2),
label = "", title = ["Predator-Littoral" "Predator Littoral" "Forager-Littoral" "Forager-Pelagic" "Juvenile-littoral" "Juvenile-Pelagic"])
using BifurcationKit
opt_newton = NewtonPar(tol = 1e-11, maxIter = 10)
#Continuation Options
opts_br = ContinuationPar(
pMin = 0.0, #minimum parameter range for p
pMax = 1.0, #maximum paramter range for p
newtonOptions = opt_newton,
ds = 0.0002,
dsmax = 0.01, #maximum arclength allowed value.
dsmin = 0.0002,
detectBifurcation = 3, #If set to 3, a bisection algorithm is used to locate the bifurcations points (slower).
nInversion = 8, #number of sign inversions in bisection algorithm
maxBisectionSteps = 25, #Maximum number of bisection steps
nev = 6, #Number of eigenvalues to be computed.
maxSteps = 1000, #Maximum number of continuation steps
θ = 0.3) # θ is a parameter in the arclength constraint.
#Initilizaing Bifurcation Problem
#Calculating Guess
z₀ = last(sol)
#Parameters for Bifurcation
par_pp2 = (
p₁ = p[1], #fishing effort
p₂ = p[2]) #dispersal rate
return z₀, par_pp2
#function to record information from a solution
recordFromSolution(x, p) = (u1 = x[1], u2 = x[2],u3 = x[3],u4 = x[4], u5 = x[5], u6 = x[6]) #Insert Explanation here
prob = BifurcationProblem(
VF, #Vector field
z₀, #Initial guess
par_pp2, #parameters
(@lens _.p₁), #It specifies which parameter axis among params is used for continuation
recordFromSolution = recordFromSolution #Insert explanation here
)
#Compute the bifurcation diagram associated with the problem F(x, p) = 0 recursively.
diagram = continuation(
prob, #bifurcation problem
PALC(), # continuation algorithm
opts_br;
plot = true,
verbosity = 2,
normC = norminf
)
return nothing
