It seems that main issue was the initial guess.
The algorithm was not able to reach an equilibrium point from the provided initial guess.
I'm no expert, but with minor tweaks I managed to get the solutions you mentioned.
The code:
using Revise, Parameters, Plots
using BifurcationKit
const BK = BifurcationKit
# vector field of the problem
function COm(u, p)
@unpack r,K,a,h,eps,mu = p
x, y = u
out = similar(u)
out[1] = r*x*(1.0 - x/K) - a*x*y/(1 + a*h*x)
out[2] = eps*a*x*y/(1 + a*h*x) - mu*y
out
end
####
# integrate ODE for sanity check
# a=0.1 stable
# a=0.3 oscillatory
using DifferentialEquations
z0 = [0.05, 0.1]
alg_ode = Rodas5()
# stable
@reset par_com.a = 0.1
prob_de = ODEProblem(COm, z0, (0.0, 300.0), par_com)
sol_ode = solve(prob_de, alg_ode);
plot(sol_ode)
# keep equilibrium point
zEq = sol_ode.u[end]
# oscillatory
@reset par_com.a = 0.3
prob_de = ODEProblem(COm, z0, (0.0, 300.0), par_com)
sol_ode = solve(prob_de, alg_ode);
plot!(sol_ode)
# parameters used in the model
par_com = (r = 1.0, K = 10.0, a = 0.1, h = 0.5, eps = 0.5, mu = 0.2)
# record variables for plotting
recordCO(x, p) = (x = x[1], y = x[2])
# initial condition
# z0 = [0.05, 0.1] ## fail to find equilibrium
# z0 = [5.0, 6.3] ## guess, from numerical data
z0 = zEq # initiate from numerical data
# Bifurcation Problem
prob = BifurcationProblem(COm, z0, par_com, (@optic _.a); record_from_solution = recordCO)
# continuation parameters
opts_br = ContinuationPar(p_min = 0.01, p_max = 0.5, dsmin=0.01, ds = 0.1, dsmax = 1.0)
# compute the branch of solutions
br = continuation(prob, PALC(), opts_br; plot = true, verbosity = 0, normC = norminf,bothside = true)
# plot the branch
scene = plot(br)
## inspect transient before and after the special points
pt = 3 # point of interest
prob_de = ODEProblem(COm, br.specialpoint[pt].x, (0.0, 1000.0),
@set par_com.a = br.specialpoint[pt].param - 0.01)
sol_ode = solve(prob_de, alg_ode);
plot(sol_ode)
prob_de = ODEProblem(COm, br.specialpoint[pt].x, (0.0, 1000.0),
@set par_com.a = br.specialpoint[pt].param + 0.01)
sol_ode = solve(prob_de, alg_ode);
plot!(sol_ode)
The continuation result:
julia> br
┌─ Curve type: EquilibriumCont
├─ Number of points: 272
├─ Type of vectors: Vector{Float64}
├─ Parameter a starts at 0.01, ends at 0.5
├─ Algo: PALC
└─ Special points:
- # 1, endpoint at a ≈ +0.01000000, step = 0
- # 2, bp at a ≈ +0.04983159 ∈ (+0.04983159, +0.05006280), |δp|=2e-04, [converged], δ = (-1, 0), step = 252
- # 3, hopf at a ≈ +0.30509470 ∈ (+0.29935407, +0.30509470), |δp|=6e-03, [converged], δ = ( 2, 2), step = 269
- # 4, endpoint at a ≈ +0.50000000,
I was going to add the images, but 10 reputation points are needed. It is my first post in here...
(Column A contains the formula, column B was manual input)
