# BifurcationKit: detecting fold bifurcation
**URL:** https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096
**Category:** Modelling & Simulations
**Tags:** question, package
**Created:** [November 26, 2024, 12:18pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096 "2024-11-26T12:18:56Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@mmplac14](https://discourse.julialang.org/u/mmplac14)
#### Post date: [November 26, 2024, 12:18pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/1 "2024-11-26T12:18:56Z")
This post is related to my previous [one](https://discourse.julialang.org/t/bifurcationkit-periodicorbitocollproblem-failure-to-converge-with-given-tolerances/123048).
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](https://www.researchgate.net/publication/252117635_Analysis_of_Jansen's_model_of_a_single_cortical_column).
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](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/DeflatedContinuation/)? I don’t understand how to use deflated continuation.
My code so far is below, any help would be appreciated.
```julia
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)
```
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [November 26, 2024, 1:12pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/3 "2024-11-26T13:12:50Z")
The problem is again the tolerance as for the fixed points. I dont know why but it can’t get to 1e-11 in sup norm. You can see that by looking at the newton iterations:
```julia
opts_br_po = ContinuationPar(p_min =-100.0, p_max = 400.0, max_steps = 3000, newton_options = NewtonPar(tol = 1e-12, verbose = true), plot_every_step = 150, detect_bifurcation = 0)
br_po2 = @time continuation(br, 3, opts_br_po,
PeriodicOrbitOCollProblem(50,5; jacobian = BK.DenseAnalyticalInplace(), meshadapt = true);
verbosity = 1,
plot = true,
alg = PALC(tangent = Bordered()),
linear_algo = COPBLS(),
normC = norminf,
```
which gives
```julia
Continuation step 310
Step size = 1.0000e-04
Parameter p = 1.4571e+01 ⟶ 1.4571e+01 [guess]
┌─────────────────────────────────────────────────────┐
│ Newton step residual linear iterations │
├─────────────┬──────────────────────┬────────────────┤
│ 0 │ 3.3916e-05 │ 0 │
│ 1 │ 1.2987e-11 │ 1 │
│ 2 │ 2.9181e-12 │ 1 │
│ 3 │ 3.4021e-12 │ 1 │
│ 4 │ 1.1385e-11 │ 1 │
│ 5 │ 7.6389e-12 │ 1 │
│ 6 │ 2.2868e-12 │ 1 │
│ 7 │ 1.7537e-11 │ 1 │
│ 8 │ 1.4163e-11 │ 1 │
│ 9 │ 1.0234e-11 │ 1 │
│ 10 │ 8.0037e-12 │ 1 │
│ 11 │ 1.0922e-11 │ 1 │
│ 12 │ 1.2774e-11 │ 1 │
│ 13 │ 1.9409e-11 │ 1 │
│ 14 │ 1.2409e-11 │ 1 │
│ 15 │ 7.0216e-12 │ 1 │
│ 16 │ 7.3513e-12 │ 1 │
│ 17 │ 1.6162e-11 │ 1 │
│ 18 │ 3.1005e-12 │ 1 │
│ 19 │ 3.1916e-12 │ 1 │
│ 20 │ 9.8766e-12 │ 1 │
│ 21 │ 2.0967e-11 │ 1 │
│ 22 │ 6.6990e-12 │ 1 │
│ 23 │ 2.3850e-12 │ 1 │
│ 24 │ 1.5516e-11 │ 1 │
│ 25 │ 2.2145e-11 │ 1 │
└─────────────┴──────────────────────┴────────────────┘
```
Instead, if you use smaller tolerances:
```julia
br_po2 = @time continuation(br, 3, opts_br_po,
PeriodicOrbitOCollProblem(50,5; jacobian = BK.DenseAnalyticalInplace(), meshadapt = true);
verbosity = 1,
alg = PALC(tangent = Bordered()),
normC = norminf,
callback_newton = BK.cbMaxNorm(1e2), #limit residual to avoid Inf or NaN
)
plot(br, br_po, br_po2)
plot!(br_po2, vars=(:param, :min))
```
---
### Author:  [@mmplac14](https://discourse.julialang.org/u/mmplac14)
#### Post date: [November 26, 2024, 2:13pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/4 "2024-11-26T14:13:46Z")
Thanks for looking into this @rveltz
Is it meant to be BK.DenseAnalytical() instead of BK.DenseAnalyticalInplace()? I am getting:
UndefVarError: `DenseAnalyticalInplace` not defined in `BifurcationKit` and I have the latest version (0.4.4)
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [November 26, 2024, 2:15pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/5 "2024-11-26T14:15:00Z")
Sorry, discard this. It is a new faster jacobian that will be added soon
---
### Author:  [@mmplac14](https://discourse.julialang.org/u/mmplac14)
#### Post date: [November 26, 2024, 2:26pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/6 "2024-11-26T14:26:45Z")
Ok! Can you send me the exact code you used? I am trying what you said and I am failing to converge
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [November 26, 2024, 2:30pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/7 "2024-11-26T14:30:27Z")
```julia
br_po2 = @time continuation(br, 3, opts_br_po,
PeriodicOrbitOCollProblem(50,5; jacobian = BK.DenseAnalytical(), meshadapt = true);
verbosity = 0,
plot = true,
alg = PALC(tangent = Bordered()),
linear_algo = COPBLS(),
normC = norminf,
callback_newton = BK.cbMaxNorm(1e2), #limit residual to avoid Inf or NaN
)
```
---
### Author:  [@mmplac14](https://discourse.julialang.org/u/mmplac14)
#### Post date: [November 26, 2024, 2:33pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/8 "2024-11-26T14:33:29Z")
Should the second param be 4 instead of 3?
`AssertionError: The provided index does not refer to a Hopf Point`
---
### Author:  [@mmplac14](https://discourse.julialang.org/u/mmplac14)
#### Post date: [November 26, 2024, 2:40pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/9 "2024-11-26T14:40:53Z")
If I run the below, I get:
┌ Error: Failure to converge with given tolerances.
│ We reached the smallest value [dsmin] valid for ds, namely 0.0001.
│ Stopping continuation at continuation step 374.
```julia
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)
print("PERIODIC ORBITS")
# continuation parameters
opts_br_po = ContinuationPar(p_min =-100.0, p_max = 400.0, max_steps = 3500, newton_options = NewtonPar(tol = 1e-12, verbose = true), plot_every_step = 150, detect_bifurcation = 0)
br_po2 = @time continuation(br, 4, opts_br_po,
PeriodicOrbitOCollProblem(50,5; jacobian = BK.DenseAnalytical(), meshadapt = true);
verbosity = 0,
plot = true,
alg = PALC(tangent = Bordered()),
linear_algo = COPBLS(),
normC = norminf,
callback_newton = BK.cbMaxNorm(1e2), #limit residual to avoid Inf or NaN
)
plot(br, br_po, br_po2)
plot!(br_po2, vars=(:param, :min))
```
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [November 26, 2024, 3:02pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/10 "2024-11-26T15:02:05Z")
> [@mmplac14](#):
>
> Should the second param be 4 instead of 3?
yes
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [November 26, 2024, 3:04pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/11 "2024-11-26T15:04:51Z")
> [@mmplac14](#):
>
> Failure to converge with given tolerances.
same as above decrease tol in newtonPar for periodicorbits
---
### Author:  [@mmplac14](https://discourse.julialang.org/u/mmplac14)
#### Post date: [November 26, 2024, 4:11pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/12 "2024-11-26T16:11:26Z")
Thanks, I have managed to reproduce the figure! I now want to colour the branches based on stability. However, br\_po.stable and br\_po\_2.stable contain `nothing`.
I have also noticed that the branches contain many fold points, with `guess` status. Is this related to the above?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [November 26, 2024, 5:33pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/13 "2024-11-26T17:33:46Z")
No, it is because it is close to the homoclinic point. Their distance is 1e-4 !
---
### Author:  [@photor](https://discourse.julialang.org/u/photor)
#### Post date: [November 27, 2024, 2:40am UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/14 "2024-11-27T02:40:35Z")
very interesting work
---
### Author:  [@mmplac14](https://discourse.julialang.org/u/mmplac14)
#### Post date: [November 27, 2024, 8:43am UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/15 "2024-11-27T08:43:44Z")
Is there another way of determining the stability of the periodic orbit branches? Or is it impossible because of the homoclinic orbit?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [November 27, 2024, 5:11pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/16 "2024-11-27T17:11:11Z")
It is because the tolerance on the stability is a bit too high. If you look at
You can remediate to these spurious bifurcations by looking at `br_po2.eig`, there should always be the zero Floquet exponent. But close to the homoclinic, it is difficult to get it precisely. This can cause spurious detection of bifurcatons.
You can change this by specifying a tolerance for stability:
```julia
opts_br_po = ContinuationPar(opts_br, p_min =-100.0, p_max = 400.0, max_steps = 3000, dsmax = 0.3, newton_options = NewtonPar(tol = 1e-10, verbose = false), plot_every_step = 150, detect_bifurcation = 3, tol_stability = 1e-4)
```
I will remediate to that in the near future
---
### Author:  [@mmplac14](https://discourse.julialang.org/u/mmplac14)
#### Post date: [November 28, 2024, 9:23am UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/17 "2024-11-28T09:23:50Z")
Thanks! I am getting the stability now but I still have way more special points than I should. Is there also a way of remediating this?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [November 28, 2024, 1:21pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/18 "2024-11-28T13:21:44Z")
if you zoom on the homoclinic point, it gives:
If you want to better resolve this zone you can increase the descretization. You can also try a different method to compute the periodic orbit
```julia
using OrdinaryDiffEq
prob_ode = ODEProblem(jansenrit!, z0, (0,1.), par)
br_po2 = continuation(br, 3, opts_br_po,
# PeriodicOrbitOCollProblem(50,5; jacobian = BK.DenseAnalyticalInplace(), meshadapt = true);
ShootingProblem(6, prob_ode, Rodas5(); abstol = 1e-10, reltol = 1e-8, parallel = true);
verbosity = 1,
plot = true,
alg = PALC(tangent = Bordered()),
normC = norminf,
callback_newton = BK.cbMaxNorm(1e2),
)
```
---
### Author:  [@mmplac14](https://discourse.julialang.org/u/mmplac14)
#### Post date: [November 28, 2024, 2:47pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/19 "2024-11-28T14:47:21Z")
And how do I increase the discretisation? Is that N in PeriodicOrbitOCollProblem?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [November 28, 2024, 3:22pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/20 "2024-11-28T15:22:41Z")
yes
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [November 28, 2024, 3:23pm UTC](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096/21 "2024-11-28T15:23:22Z")
but I would not worry too much about those bifurcations, you are very close to the homoclinic point
[Next page](https://discourse.julialang.org/t/bifurcationkit-detecting-fold-bifurcation/123096.md?page=2)