Hi!
I am using BifurcationKit to compute the periodic orbits of the Jansen-Rit model from the first Hopf Point. I seem to be having issues computing them for y1 - y2. I am using record_from_solution to compute them for y1-y2 because it seems to me that if I don’t use record_from_solution only the values for y0 are contained in the result.
I am getting this (and the branch should not be flicking back up)
It should stop at p =113, as shown here:
Here is my full code:
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 #These are the parameters. p will be our control paramter
y0, y3, y1, y4, y2, y5 = z #These are our functions
[
#We define the derivative functions in the same order as above
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]
# 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], y = x[3]-x[5]),)
prob = BifurcationProblem(jansenrit, z0, par, (@optic _.p);
record_from_solution = (x, p; k...) -> (y = x[3] - x[5]),)
opts_br_po = ContinuationPar(p_min =-50.0, p_max = 400.0, max_steps = 3000, dsmin = 1e-7, ds = 1e-3, newton_options = NewtonPar(tol = 1e-10, verbose = true), detect_bifurcation = 3, tol_stability = 1e-4)
br_po = @time continuation(br, 4, opts_br_po,
PeriodicOrbitOCollProblem(50,5; jacobian = BK.DenseAnalytical(), meshadapt = true);
verbosity = 0,
alg = PALC(tangent = Bordered()),
linear_algo = COPBLS(),
normC = norminf,
callback_newton = BK.cbMaxNorm(1e2), #limit residual to avoid Inf or NaN
record_from_solution = (x, p; k...) -> (y = x[3] - x[5])
)
#Plotting
# Colour by stability of each periodic orbit bracn
stability_po = [stable ? :green : :red for stable in br_po.stable]
# Plot the equilibria branch
plot(br, branchlabel = "equilibria")
# Plot the max and min of the periodic orbit branch
plot!(br_po, c = colour_stability)
After some experimenting, I believe this is due to “record_from_solution”. If I set record_from_solution to y0, but don’t specify it in the Periodic Orbit Continuation, I get the result I expect for y0. However, if I use “record_from_solution” to record y0, I also get strange behaviour, with the branch flipping back up and joining with the bifurcation point.
— Without setting record_from_solution:
br_po = @time continuation(br, 4, opts_br_po,
PeriodicOrbitOCollProblem(50,5; jacobian = BK.DenseAnalytical(), meshadapt = true);
verbosity = 0,
alg = PALC(tangent = Bordered()),
linear_algo = COPBLS(),
normC = norminf,
callback_newton = BK.cbMaxNorm(1e2), #limit residual to avoid Inf or NaN
)
—With record_from_solution set to y0:
br_po = @time continuation(br, 4, opts_br_po,
PeriodicOrbitOCollProblem(50,5; jacobian = BK.DenseAnalytical(), meshadapt = true);
verbosity = 0,
alg = PALC(tangent = Bordered()),
linear_algo = COPBLS(),
normC = norminf,
callback_newton = BK.cbMaxNorm(1e2), #limit residual to avoid Inf or NaN
record_from_solution = (x, p; k...) -> (y = x[1])
)
Lastly, it seems like I can’t get both the min and max values of the periodic orbit with the way I set y0 or y1-y2 in “record_from_solution”. Is there a way to do so?
