Hello!
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.
I have managed to get the continuation graph.
I also managed to get the periodic orbit starting from the 1st Hopf bifurcation.
Now, I am trying to get the periodic orbits produced by the 2nd Hopf Bifurcation at p = 89. I am using PeriodicOrbitOCollProblem, same as I used to get the graph above. The calculation fails to converge.
Error: Failure to converge with given tolerances.
│ We reached the smallest value [dsmin] valid for ds, namely 1.0e-8.
│ Stopping continuation at continuation step 1.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/RaJtn/src/continuation/Contbase.jl:67
24.719435 seconds (25.34 M allocations: 20.385 GiB, 20.35% gc time, 45.28% compilation time)
I also get logs saying the precision on the FloquetMultiplier is not enough. I have been decreasing the tol_stability, but these warnings keep appearing.
┌ Warning: The precision on the Floquet multipliers is 0.01921577015310253.
│ It may be not enough to allow for precise bifurcation detection.
│ Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetColl` for computing Floquet coefficients.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/RaJtn/src/periodicorbit/Floquet.jl:532
Should I be using another type of PeriodOrbit problem? I am confused about when I should be using which.
Do I need to use a different method than FloquetColl? Which one? I have looked at the docs but can’t figure out how to change it and what the other options are.
My code is below. Any pointers would be really 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)
# 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)
scene = plot(br, legend=:topleft, ylabel = "y")
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.00001, dsmin = 1e-8, dsmax = 0.1,
max_steps = 100,
tol_stability = 1e-10)
br_pocoll = @time continuation(
# we want to branch from the 5th bif. point
br, 5, opts_po_cont,
# we want to use the Collocation method to locate PO, with polynomial degree 4
PeriodicOrbitOCollProblem(60, 4; meshadapt = true);
# regular continuation options
normC = norminf,
args_po...)
br_pocoll