Too many bifurcation points marked

I’m trying to use BifurcationKit.jl to continue a periodic orbit in a 2d dynamical system. Currently, I’m seeing “too many” bifurcation points marked, i.e. blue points marked “bp” appear where the system does not change stability. Looking at the eigenvalues, it looks like these blue dots come from places the eigenvalue at (0,0), which is always there because it is a periodic orbit, cross zero due to numerical error.

The only parameter change which “fixes” the problem that I have found is setting “tol_stability = 10^-1”.

Code below:

using Revise, Plots
using BifurcationKit
const BK = BifurcationKit
using DifferentialEquations

function cycle_2!(du, X, p, t = 0)
(;wEE_self, wEE, wEI, wIE, wII, tau_E, tau_I, b) = p
x1, x2 = X
du[1] =(1/tau_E)* (-x1 + max(0, wEE_selfx1 + wEI * x2 + b))
du[2] = (1/tau_I)
(-x2 + max(0, wIE * x1 + wII * x2))
du
end

par_cycle = (wEE_self =1.5, wEE = .75 , wEI = -1.5, wIE = 2.0, wII = -1.0, tau_E = 1.0, tau_I = 1.0, b = 1.0)

recordFromSolution(x, p; k…) = (u1 = x[1], u2 = x[2])

prob = BifurcationProblem(cycle_2!, z0, par_cycle,

specify the continuation parameter

(@optic _.tau_I), record_from_solution = recordFromSolution)

opts_br = ContinuationPar(p_min = 0.5, p_max = 8.0, dsmax = 0.01,

number of eigenvalues

nev = 2,

maximum number of continuation steps

max_steps = 10000,)

br = continuation(prob, PALC(tangent=Bordered()), opts_br; normC = norminf)

par_cycle = (wEE_self =1.5, wEE = .75 , wEI = -1.5, wIE = 2.0, wII = -1.0, tau_E = 1.0, tau_I = 3.9, b = 1.0)

z0 = [1.5, 1.5]

prob = BifurcationProblem(cycle_2!, z0, par_cycle, (@optic _.tau_I);
record_from_solution = (x, p; k…) → (x = x[1], y = x[2]))

prob_de = ODEProblem(cycle_2!, z0, (0,9.8), par_cycle)
sol = DifferentialEquations.solve(prob_de, Rodas5())

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 = “ext”, arg…, k…)
plot!(xtt.t, xtt[2,:]; label = “inh”, arg…, k…)
plot!(br; subplot = 1, putspecialptlegend = false)
end,

we use the supremum norm

normC = norminf)

probtrap, ci = BK.generate_ci_problem(PeriodicOrbitTrapProblem(M = 300),
prob, sol, 9.8)

opts_br = ContinuationPar(p_min = 0.01, p_max = 8.0, ds =-0.002, dsmax = 0.01)
opts_po_cont = ContinuationPar(opts_br, max_steps = 1000, tol_stability = 1e-3)
brpo_fold = continuation(probtrap, ci, PALC(), opts_po_cont;
verbosity = 3, plot = true,
args_po…
)

fu

Hi,

Your code is not runnable. You should use triple ``` like

a=1

Now, your example seems close to 🟢 Neural mass equation - MTK · Bifurcation Analysis in Julia

Finally, your vector field is not differentiable because of the max function. Hence, do not expect to get meaningfull results.

I think I fixed whatever the typo with the quotes was–this is running code, so it must be something introduced while copy/pasting here.

Should I really not expect to get meaningful results for a non-smooth system? It’s having no problem actually continuing the limit cycle.

You need to linearise (differentiate) around the periodic orbit to get its stability so…

Thanks. I’m fully aware that my system is not smooth. However, it is smooth (and even linear!) off of the two lines where either

wEE_self*x1 + wEI * x2 + b == 0

or wIE * x1 + wII * x2 == 0.

It’s even possible to write down the poincare map explicitly in terms of the amount of time the limit cycle spends in each linear chamber. This implies to me that it is possible to linearize around the limit cycle, so I must be missing something.
Sorry for the naieve question, I’m very new to bifurcation theory and to this software.

If your periodic orbit is stable (attracting), you can consider the global_continuation from Attractors.jl that doesn’t care about smoothness of the flow. It doesn’t give you bifurcation points however, if that’s what you are after; just the limit cycle itself, along with any other attractors in the system.

Its not clear from your post what you care about to achieve: continuing a stable limit cycle or finding bifurcations. My comment only helps for the first part.

I tried with Collocation and Shooting and the Floquet coefficients arent that great. I bet your piecewise linear system introduces jumps. I know that Stephen Coombes has worked out your case but BK wont be able to handle that gracefully. You can easily follow the PO but the stability is not that great. You can put a softmax though and be as close as you want from your orginal problem

If you want to work with the non-smooth system directly, the best software available at the moment is COCO by Dankowicz and Schilder. Look for hybrid periodic orbits in the documentation - it was originally designed for this sort of problem. Sadly, it’s in MATLAB.

As rveltz mentions, smoothing out the nonlinearity might be sufficient for you, however there can be some differences that occur if there are non-smooth bifurcations (e.g. grazing).

Thank you all for the software recommendations!