# Unable to reproduce a simple bifurcation diagram
**URL:** https://discourse.julialang.org/t/unable-to-reproduce-a-simple-bifurcation-diagram/124738
**Category:** Modelling & Simulations
**Tags:** question
**Created:** [January 13, 2025, 7:02pm UTC](https://discourse.julialang.org/t/unable-to-reproduce-a-simple-bifurcation-diagram/124738 "2025-01-13T19:02:06Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@vembha](https://discourse.julialang.org/u/vembha)
#### Post date: [January 13, 2025, 7:02pm UTC](https://discourse.julialang.org/t/unable-to-reproduce-a-simple-bifurcation-diagram/124738/1 "2025-01-13T19:02:06Z")
Hello,
New to `BifurcationKit.jl` here! I have a simple ODE-based model from an article (Fig. 2A in [https://www.pnas.org/doi/10.1073/pnas.1419162112](https://www.pnas.org/doi/10.1073/pnas.1419162112)), which is known to have three bifurcation points when varied over a parameter `m`. The figure is attached:
_The parameter is varied between 0 and 1, the dashed branches are unstable, and the solid branches are stable._
However, I am unable to reproduce the bifurcation diagram. The `PALC()` algorithm fails to identify the bifurcation points. It does not give out any errors. The code is below:
```julia
using BifurcationKit, Plots
# ODE system
function ConwayPerelson_2015(u, p)
(λ, d, β, αL, ρ, a, dL, δ, m, pr, c, λE, bE, KB, dE, KD, μ) = p
T, L, I, V, E = u # The variable V is plotted in the bifurcation diagram
[
λ - d*T - β*V*T,
αL*β*V*T + (ρ - a - dL)*L,
(1 - αL)*β*V*T - δ*I + a*L - m*E*I,
pr*I - c*V,
λE + bE*E*I/(KB + I) - dE*E*I/(KD + I) - μ*E
]
end
# Parameters
params = (
λ = 1E4, d = 0.01, β = 1.5E-8, δ = 1, pr = 2_000, c = 23,
a = 0.001, dL = 0.004, ρ = 0.0045, αL = 1E-6,
λE = 1, bE = 1, KB = 0.1, dE = 2, KD = 5, μ = 2, m = 0.42
)
# Initial guess
u0 = zeros(5)
# Bifurcation problem
prob = BifurcationProblem(ConwayPerelson_2015, u0, params, (@optic _.m);
record_from_solution = (x, p; k...) -> (V = x[4]))
opts = ContinuationPar(p_min = 0.0, p_max = 1.0)
# Solve
br = continuation(prob, PALC(), opts, bothside = true)
```
Can someone please help me out? Is there a way to ensure that the algorithm does not miss out on identifying the bifurcation points?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 13, 2025, 7:28pm UTC](https://discourse.julialang.org/t/unable-to-reproduce-a-simple-bifurcation-diagram/124738/2 "2025-01-13T19:28:25Z")
You made a tiny mistake (semicolon missing) in the vector field
```julia
(;λ, d, β, αL, ρ, a, dL, δ, m, pr, c, λE, bE, KB, dE, KD, μ) = p
```
You should also make the computation of bifurcation points more precise `opts = ContinuationPar(p_min = 0.0, p_max = 1.0, n_inversion=6)`
---
### Author:  [@vembha](https://discourse.julialang.org/u/vembha)
#### Post date: [January 13, 2025, 8:02pm UTC](https://discourse.julialang.org/t/unable-to-reproduce-a-simple-bifurcation-diagram/124738/3 "2025-01-13T20:02:38Z")
Thank you very much for the prompt reply, and I am sorry for the silly miss of the semicolon.
The algorithm works, but it is giving me only one of the bifurcation points.
FYI, the estimated bifurcation point is highlighted as red II in the figure in the question. But there are two more. This is despite modifying the `opts` to detect fold bifurcations.
```julia
opts = ContinuationPar(p_min = 0.0, p_max = 1.0, n_inversion = 6,
detect_bifurcation = 3, detect_fold = true, dsmin = 1E-4)
br = continuation(prob, PALC(), opts, bothside = true, normC = norminf)
```
Any idea how to ensure the algorithm finds all of them?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 13, 2025, 8:33pm UTC](https://discourse.julialang.org/t/unable-to-reproduce-a-simple-bifurcation-diagram/124738/4 "2025-01-13T20:33:52Z")
diagram = bifurcationdiagram(prob,PALC(), 2, ContinuationPar(opts, max\_steps = 10000))
gives the warning
```julia
┌ Error: Failure to converge with given tolerance = 1.0e-12.
│ Step = 9741
│ You can decrease the tolerance or pass a different norm using the argument `normC`.
│ We reached the smallest value [dsmin] valid for ds, namely 0.0001.
│ Stopping continuation at continuation step 9741.
```
but `plot(diagram)` yields:
The warning basically says that the newton tolerance are too high, hence:
```julia
@reset opts.newton_options.tol=1e-10
diagram = bifurcationdiagram(prob,PALC(), 2, ContinuationPar(opts, max_steps = 150000))
plot(diagram)
```
or:
```julia
plot(diagram, applytoY = x->log(1+abs(x)))
```
---
### Author:  [@vembha](https://discourse.julialang.org/u/vembha)
#### Post date: [January 14, 2025, 6:46am UTC](https://discourse.julialang.org/t/unable-to-reproduce-a-simple-bifurcation-diagram/124738/5 "2025-01-14T06:46:54Z")
Thank you very much! I now see where I fell short, and the package is super cool 🙌. I was able to recapitulate the figure:
Two quick questions, however:
1. A non-existent `hopf` point is being predicted, very close to the point where a simple fold bifurcation is known to exist (see the StackTrace below). So, I am guessing it is an error of the algorithm. Are there any settings that I could change to avoid such incorrect estimations?
```julia
[Bifurcation diagram]
┌─ From 1-th bifurcation point.
├─ Children number: 0
└─ Root (recursion level 2)
┌─ Curve type: EquilibriumCont from Transcritical bifurcation point.
├─ Number of points: 150001
├─ Type of vectors: Vector{Float64}
├─ Parameter m starts at 0.6087198085395702, ends at 0.22451213423608637
├─ Algo: PALC
└─ Special points:
- # 1, bp at m ≈ +0.40983773 ∈ (+0.40983773, +0.40983773), |δp|=2e-11, [converged], δ = ( 1, 0), step = 53
- # 2, bp at m ≈ +0.83909847 ∈ (+0.83909847, +0.83909847), |δp|=2e-15, [guess], δ = ( 1, 0), step = 7161
- # 3, hopf at m ≈ +0.83853493 ∈ (+0.83853493, +0.83853493), |δp|=2e-09, [converged], δ = (-2, -2), step = 8138
- # 4, endpoint at m ≈ +0.22450773, step = 150001
```
1. Is there a way to force the top stable branch to go all the way to `m=0.0`? I figured increasing `max_steps` does that, but it seems I would need a large value for it, probably because the top branch is being explored very slowly.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 14, 2025, 12:48pm UTC](https://discourse.julialang.org/t/unable-to-reproduce-a-simple-bifurcation-diagram/124738/6 "2025-01-14T12:48:45Z")
> [@vembha](#):
>
> Is there a way to force the top stable branch to go all the way to `m=0.0`? I figured increasing `max_steps` does that, but it seems I would need a large value for it, probably because the top branch is being explored very slowly.
You can use a larger `dsmax`:
```julia
diagram = bifurcationdiagram(prob,PALC(), 2, ContinuationPar(opts, max_steps = 150000, dsmax = 0.7))
```
The issue here is that the range of values of `ConwayPerelson_2015 ` is pretty large. I would be better numerically to rescale the model (log or other).
> A non-existent `hopf` point is being predicted, very close to the point where a simple fold bifurcation is known to exist (see the StackTrace below).
You can see
```julia
julia> diagram[2]
[Bifurcation diagram]
┌─ From 1-th bifurcation point.
├─ Children number: 0
└─ Root (recursion level 2)
┌─ Curve type: EquilibriumCont from Transcritical bifurcation point.
├─ Number of points: 150001
├─ Type of vectors: Vector{Float64}
├─ Parameter m starts at 0.6087056695593298, ends at 0.7526129563574361
├─ Algo: MoorePenrose
└─ Special points:
- # 1, bp at m ≈ +0.40983773 ∈ (+0.40983773, +0.40983773), |δp|=5e-12, [converged], δ = ( 1, 0), step = 225
- # 2, bp at m ≈ +0.83909847 ∈ (+0.83909847, +0.83909847), |δp|=9e-16, [converged], δ = ( 1, 0), step = 35762
- # 3, hopf at m ≈ +0.83853493 ∈ (+0.83853493, +0.83853493), |δp|=2e-09, [converged], δ = (-2, -2), step = 40644
- # 4, endpoint at m ≈ +0.75261202,
```
It seems the Hopf bifurcation is not a fluke. See
```julia
julia> eigenvals(diagram[2].γ, 40644)
5-element Vector{ComplexF64}:
-2.0506719532705528e-10 - 0.005389356755654721im
-2.0506719532705528e-10 + 0.005389356755654721im
-0.000499605123267488 + 0.0im
-2.8970614512243027 + 0.0im
-24.290528399220722 + 0.0im
julia> eigenvals(diagram[2].γ, 40645)
5-element Vector{ComplexF64}:
-1.2226521003841756e-7 - 0.005389880502352536im
-1.2226521003841756e-7 + 0.005389880502352536im
-0.0004996051977692838 + 0.0im
-2.8970638500412025 + 0.0im
-24.29052805667206 + 0.0im
julia> eigenvals(diagram[2].γ, 40643)
5-element Vector{ComplexF64}:
8.657536849168388e-8 - 0.005388984365665424im
8.657536849168388e-8 + 0.005388984365665424im
-0.000499605070280656 + 0.0im
-2.8970597457475535 + 0.0im
-24.29052864275146 + 0.0im
```
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 14, 2025, 12:49pm UTC](https://discourse.julialang.org/t/unable-to-reproduce-a-simple-bifurcation-diagram/124738/7 "2025-01-14T12:49:35Z")
You could use MTK to have analytical jacobian if you are unsure
---
### Author:  [@vembha](https://discourse.julialang.org/u/vembha)
#### Post date: [January 15, 2025, 6:35am UTC](https://discourse.julialang.org/t/unable-to-reproduce-a-simple-bifurcation-diagram/124738/8 "2025-01-15T06:35:04Z")
Yes. Thank you!
However, despite using MTK for using the symbolic jacobian, the hopf is being predicted. Need to understand why … Thanks again for the help! 🙌
```julia
@variables T(t) L(t) I(t) V(t) E(t)
@parameters λ, d, β, αL, ρ, a, dL, δ, m, pr, c, λE, bE, KB, dE, KD, μ
eqs = [
D(T) ~ λ - d*T - β*V*T,
D(L) ~ αL*β*V*T + (ρ - a - dL)*L,
D(I) ~ (1 - αL)*β*V*T - δ*I + a*L - m*E*I,
D(V) ~ pr*I - c*V,
D(E) ~ λE + bE*E*I/(KB + I) - dE*E*I/(KD + I) - μ*E
]
@mtkbuild odesys = ODESystem(eqs, t)
bif_par = m
plot_var = V
params = [
λ => 1E4, d => 0.01, β => 1.5E-8, δ => 1, pr => 2_000, c => 23,
a => 0.001, dL => 0.004, ρ => 0.0045, αL => 1E-6,
λE => 1, bE => 1, KB => 0.1, dE => 2, KD => 5, μ => 2, m => 0.0
]
u0 = zeros(5)
prob = BifurcationProblem(odesys, u0, params, bif_par;
plot_var = plot_var, jac = true)
```
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 15, 2025, 6:37am UTC](https://discourse.julialang.org/t/unable-to-reproduce-a-simple-bifurcation-diagram/124738/9 "2025-01-15T06:37:56Z")
Cool, maybe there is a small paper there if they did not see that.