# Producing Fold Bifurcation Diagram using BirfurcationKit

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 = u₅/tmat - m*u₁ - p₁*u₁ + s + ϵ*hₓ*u₁*u₃              + σ₁*γ*p₂*u₂ - σ₂*(1/γ)*p₂*u₁    #Littoral
dz = u₆/tmat - m*u₂ - p₁*u₂ + s + ϵ*h_y*u₂*u₄             - σ₁*γ*p₂*u₂ + σ₂*(1/γ)*p₂*u₁    #Pelagic

#Foragers Equations
dz = r*u₃ - b*u₃*u₃ - hₓ*u₁*u₃ +  i + kₓ*u₃*u₅ + σ₃*γ*p₂*u₄ - σ₄*(1/γ)*p₂*u₃             #Littoral
dz = r*u₄ - b*u₄*u₄ - h_y*u₂*u₄ + i + k_y*u₄*u₆ - σ₃*γ*p₂*u₄ + σ₄*(1/γ)*p₂*u₃            #Pelagic

#Juveniles Equations
dz = f*(u₁+γ*u₂) - u₅/tmat - m*u₅ - c*u₁*u₅ - kₓ*u₃*u₅ + σ₅*γ*p₂*u₆ - σ₆*(1/γ)*p₂*u₅       #Littoral
dz = -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₀ = 10 #Predator; Littoral
u₀ = 10 #Predator; Palagic
u₀ = 10 #Forager; Littoral
u₀ = 10 #Forager; Palagic
u₀ = 10 #Juveniles; Littoral
u₀ = 10 #Juveniles; Palagic

#Initiliazing Parameters
p = zeros(Float64,2)
p = 0.0 #Fishing effort
p = 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).
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, #fishing effort
p₂ = p) #dispersal rate

return z₀, par_pp2

#function to record information from a solution
recordFromSolution(x, p) = (u1 = x, u2 = x,u3 = x,u4 = x, u5 = x, u6 = x) #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
``````
``````#Unclear what this step does
VF(z, p) = VF!(similar(z), z, p, 0) #Insert explanation here?????????
``````

BifurcationKit does not use inplace formulation for now but DifferentialEquations.jl does. You use both in your code so we need both formulations.

Your code seems fine. I use deflated continuation but got a mess, I am wondering if your model specification is right. Also, why dont you use namedtuple for the parameters (or else) instead of const?

What makes you think you have 2 equilibria?

you have many branches, only the red is stable (and I cut negative u1) Hi @rveltz. Thank you for taking the time to respond.

I believe that I have 2 equilibria because when p_1 = 0 and p_2 = 0.8, it seems to stabilizing on a two different equilibria respectively. (Imgur: The magic of the Internet)

My understanding is that there should be a fold bifurcation.

no necessarily. It may happen that you have two disconnected branches of solutions (see picture above).

In any case, if you know the two equilibria, you can continue them and see if the branches merge