# Too many bifurcation points marked

**URL:** <https://discourse.julialang.org/t/too-many-bifurcation-points-marked/123504>\
**Category:** General Usage\
**Tags:** question, package, dynamical-systems\
**Created:** [December 5, 2024, 4:17pm UTC](https://discourse.julialang.org/t/too-many-bifurcation-points-marked/123504 "2024-12-05T16:17:49Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Caitlin\_Lienkaemper](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caitlin_lienkaemper/32/50996_2.png) [@Caitlin\_Lienkaemper](https://discourse.julialang.org/u/Caitlin_Lienkaemper)\
**Post date:** [December 5, 2024, 4:17pm UTC](https://discourse.julialang.org/t/too-many-bifurcation-points-marked/123504/1 "2024-12-05T16:17:49Z")

</div>

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\_self_x1 + 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…  
> )

 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/6/b65362b44a8014837115c3056976c18957a70c79.png)fu

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 5, 2024, 4:26pm UTC](https://discourse.julialang.org/t/too-many-bifurcation-points-marked/123504/2 "2024-12-05T16:26:49Z")

</div>

Hi,

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

```julia
a=1

```

Now, your example seems close to [🟢 Neural mass equation - MTK · Bifurcation Analysis in Julia](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/ode/NME-MTK/#Neural-mass-equation-MTK)

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

---

<div class="post-metadata">

**Author:** ![Caitlin\_Lienkaemper](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caitlin_lienkaemper/32/50996_2.png) [@Caitlin\_Lienkaemper](https://discourse.julialang.org/u/Caitlin_Lienkaemper)\
**Post date:** [December 5, 2024, 4:41pm UTC](https://discourse.julialang.org/t/too-many-bifurcation-points-marked/123504/3 "2024-12-05T16:41:17Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 5, 2024, 4:44pm UTC](https://discourse.julialang.org/t/too-many-bifurcation-points-marked/123504/4 "2024-12-05T16:44:05Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Caitlin\_Lienkaemper](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caitlin_lienkaemper/32/50996_2.png) [@Caitlin\_Lienkaemper](https://discourse.julialang.org/u/Caitlin_Lienkaemper)\
**Post date:** [December 5, 2024, 4:55pm UTC](https://discourse.julialang.org/t/too-many-bifurcation-points-marked/123504/5 "2024-12-05T16:55:48Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [December 5, 2024, 5:59pm UTC](https://discourse.julialang.org/t/too-many-bifurcation-points-marked/123504/6 "2024-12-05T17:59:43Z")

</div>

> [@Caitlin\_Lienkaemper](#):
>
> continue a periodic orbit in a 2d dynamical system

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.

> **[Attractors.jl Tutorial · Attractors.jl](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/attractors/stable/tutorial/)**
>
> Documentation for Attractors.jl.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 5, 2024, 6:54pm UTC](https://discourse.julialang.org/t/too-many-bifurcation-points-marked/123504/7 "2024-12-05T18:54:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dawbarton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dawbarton/32/215461_2.png) [@dawbarton](https://discourse.julialang.org/u/dawbarton)\
**Post date:** [December 5, 2024, 7:09pm UTC](https://discourse.julialang.org/t/too-many-bifurcation-points-marked/123504/8 "2024-12-05T19:09:11Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![Caitlin\_Lienkaemper](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caitlin_lienkaemper/32/50996_2.png) [@Caitlin\_Lienkaemper](https://discourse.julialang.org/u/Caitlin_Lienkaemper)\
**Post date:** [December 6, 2024, 5:18am UTC](https://discourse.julialang.org/t/too-many-bifurcation-points-marked/123504/9 "2024-12-06T05:18:59Z")

</div>

Thank you all for the software recommendations!