This post is related to my previous one.
I am new to bifurcation analyis and I am working with BifurcationKit to reproduce the bifurcation analysis performed by Grimbert Faugueras on the Jansen-Rit model.
Specifically, I want to reproduce these 2 figures:
I am not able to detect the fold bifurcation of limit cycles at p = 137.5.
I have tried continuation with detect_fold = true:
and computing branches of periodic orbits starting from the 1st Hopf point.
Since the branch I am trying to find is disconnected from the main branch, do I need to use deflated continuation? I donβt understand how to use deflated continuation.
My code so far is below, any help would be appreciated.
using Revise, Plots
using BifurcationKit
const BK = BifurcationKit
# Define the function Ο
function Ο(v, e0 = 2.5, v0 = 6, r=0.56)
2 * e0 / (1 + exp(r * (v0 - v)))
end
# vector field
function jansenrit(z, param)
(;A, a, B, b, C1, C2, C3, C4, p) = param
y0, y3, y1, y4, y2, y5 = z
[
y3
A * a * Ο(y1 - y2) - 2 * a * y3 - a^2 * y0
y4
A * a * (p + C2 * Ο(C1 * y0)) - 2 * a * y4 - a^2 * y1
y5
B * b * C4 * Ο(C3 * y0) - 2 * b * y5 - b^2 * y2
]
end
# parameter values
C= 135
par = (A = 3.25, a = 100, B = 22, b = 50, C1 = C, C2 = 0.8 * C, C3 = 0.25 * C, C4 = 0.25 * C, p = 50.0)
# initial condition
z0 = [0.0, 0.0, 15.0, 0.0, 10.0, 0.0]
# Bifurcation Problem
prob = BifurcationProblem(jansenrit, z0, par, (@optic _.p);
record_from_solution = (x, p; k...) -> (y0 = x[1],y3 = x[2],y1 = x[3],y4 = x[4],y2 = x[5],y5 = x[6],),)
optnewton = NewtonPar(tol = 1e-11, verbose = true)
# continuation options, we limit the parameter range for p
opts_br = ContinuationPar(p_min = -50.0, p_max = 400.0, max_steps = 8000, newton_options = optnewton, detect_fold = true)
# continuation of equilibria
br = continuation(prob, PALC(), opts_br;
# we want to compute both sides of the branch of the initial
# value of p = 100
bothside = true, normC = norminf)
#propertynames(br)
scene = plot(br, legend=:topleft)
print("PERIODIC ORBITS")
# arguments for periodic orbits
# one function to record information and one
# function for plotting
args_po = ( record_from_solution = (x, p; k...) -> begin
xtt = get_periodic_orbit(p.prob, x, p.p)
return (max = maximum(xtt[1,:]),
min = minimum(xtt[1,:]),
period = getperiod(p.prob, x, p.p))
end,
plot_solution = (x, p; k...) -> begin
xtt = get_periodic_orbit(p.prob, x, p.p)
arg = (marker = :d, markersize = 1)
plot!(xtt.t, xtt[1,:]; label = "y0", arg..., k...)
plot!(xtt.t, xtt[2,:]; label = "y3", arg..., k...)
plot!(xtt.t, xtt[3,:]; label = "y1", arg..., k...)
plot!(xtt.t, xtt[4,:]; label = "y4", arg..., k...)
plot!(xtt.t, xtt[5,:]; label = "y2", arg..., k...)
plot!(xtt.t, xtt[6,:]; label = "y5", arg..., k...)
plot!(br; subplot = 1, putspecialptlegend = false)
end,
# we use the supremum norm
normC = norminf)
# continuation parameters
opts_po_cont = ContinuationPar(opts_br, ds= 0.001, dsmin = 1e-4, dsmax = 0.1,
max_steps = 1000,
tol_stability = 1e-5)
br_pocoll = @time continuation(
# we want to branch from the 5th bif. point
br, 4, opts_po_cont,
# we want to use the Collocation method to locate PO, with polynomial degree 4
PeriodicOrbitOCollProblem(50, 4; meshadapt = true);
# regular continuation options
normC = norminf,
args_po...)
colour_stability = [stable ? :green : :red for stable in br_pocoll.stable]
plot(br, branchlabel = "equilibria")
plot!(br_pocoll.param, br_pocoll.max, label = "periodic orbits max", color = colour_stability)
plot!(br_pocoll.param, br_pocoll.min, color = colour_stability)
