# BifurcationKit.jl incorrect automatic branch switching at simple branch point
**URL:** https://discourse.julialang.org/t/bifurcationkit-jl-incorrect-automatic-branch-switching-at-simple-branch-point/78721
**Category:** General Usage
**Created:** [March 30, 2022, 6:25am UTC](https://discourse.julialang.org/t/bifurcationkit-jl-incorrect-automatic-branch-switching-at-simple-branch-point/78721 "2022-03-30T06:25:10Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@sceptri](https://discourse.julialang.org/u/sceptri)
#### Post date: [March 30, 2022, 6:25am UTC](https://discourse.julialang.org/t/bifurcationkit-jl-incorrect-automatic-branch-switching-at-simple-branch-point/78721/1 "2022-03-30T06:25:10Z")
Hello, I have a simple 2D system and I tried continuation from simple branch points
```julia
using Plots, BifurcationKit, Setfield, ForwardDiff, Parameters
const BK = BifurcationKit
function f(x,p)
@unpack a,g,γ,r,d = p
M,C = x
T = 1 - M - C
dM = @. a * M * C - (g * M) / (M + T) + γ * M * T
dC = @. r * T * C - d * C - a * M * C
return [dM, dC]
end
jet = BK.getJet(f; matrixfree = false)
# Starting variables
start = [0.2, 0.6]
# Default values for parameters
parameters = (a = 0.1, g = 0.01, γ = 0.8, r = 1., d = 0.44)
# Options for Newton-Krylov solver and continuation
optnewton = NewtonPar(tol = 1e-12, maxIter = 100)
options = ContinuationPar(pMin = 0.01, pMax = 1., dsmin = 0.0001, dsmax = 0.05, ds = 0.01,
maxSteps = 100, nev = 30, nInversion = 8, detectBifurcation = 3, newtonOptions = optnewton,
dsminBisection = 1e-7, maxBisectionSteps = 100)
# Solving for equilibrium from starting point using Newton-Krylon solver
solution, = newton(f, start, parameters, optnewton)
solution[2]
br, = continuation(f, solution, parameters, (@lens _.g), options,
plot = true,
recordFromSolution = (x,p) -> (C = x[2], M = x[1])
)
br1, = continuation(jet..., br, 1, setproperties(options;ds = 0.0001, maxSteps = 30), recordFromSolution = (x,p) -> (C = x[2], M = x[1]), bothside = true, plot = true)
br2, = continuation(jet..., br1, 2, setproperties(options; maxSteps = 20, ds = 0.0001), recordFromSolution = (x,p) -> (C = x[2], M = x[1]), bothside = true,plot = true)
plot(br, br1; branchlabel = ["br", "br1"], legend = :topright, ylims = (0,0.8), vars = (:param, :C))
plot(br, br1; branchlabel = ["br", "br1"], legend = :topright, ylims = (0,1), vars = (:param, :M))
```
The first continuation gives expected results - it picks the second branch and continues on it.
The second aBS should result in a horizontal line C = 0, but it stays on the same branch as it was.
Is there something I’m doing wrong or is it a bug in BifurcationKit.jl?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [March 30, 2022, 8:15am UTC](https://discourse.julialang.org/t/bifurcationkit-jl-incorrect-automatic-branch-switching-at-simple-branch-point/78721/2 "2022-03-30T08:15:53Z")
Hi,
You are right, something is wrong with the predictor from the transcitical bifurcation. In case you need it, this works:
```julia
parameters2 = @set parameters.g=br1.specialpoint[2].param + 0.01
solution2 = br1.specialpoint[2].x
br3, = continuation(f, solution2, parameters2, (@lens _.g), options,
plot = true,
recordFromSolution = (x,p) -> (C = x[2], M = x[1]),
plotSolution = (x, p; k...) -> begin
plot!(br1, subplot = 1)
end,
)
```
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [March 30, 2022, 6:43pm UTC](https://discourse.julialang.org/t/bifurcationkit-jl-incorrect-automatic-branch-switching-at-simple-branch-point/78721/3 "2022-03-30T18:43:21Z")
You can use deflation to help it converge to another branch:
```julia
br1, = continuation(jet..., br, 1, setproperties(options;ds = 0.0001, maxSteps = 30), recordFromSolution = (x,p) -> (C = x[2], M = x[1]), bothside = true, plot = true)
br2, = continuation(jet..., br1, 2, setproperties(options; maxSteps = 20, ds = 0.0001), recordFromSolution = (x,p) -> (C = x[2], M = x[1]), bothside = true,plot = true, usedeflation = true)
plot(br,br1,br2)
```
However, I am puzzled by the result from bifurcation theory:
```julia
# newton refinement of transcritical point
solt, = newton(jet[1:2]..., br1, 2) #[0.50909 2.4830e-16]
# jacobian at bifurcation point
J = jet[2](solt.u, @set parameters.g = solt.p)
vals, vecs = eigen(J)
J * vecs[:, 2]
# the singular eigenvector
e0 = vecs[:, 2] #[0.8068742471699807 -0.5907232425204519]
```
Essentially, the predictor is `solt + α e0` with α small. There is no way that this is close to [\* 0]. Something is off 🤔 from a theoritical point of view
---
### Author:  [@sceptri](https://discourse.julialang.org/u/sceptri)
#### Post date: [March 31, 2022, 10:15am UTC](https://discourse.julialang.org/t/bifurcationkit-jl-incorrect-automatic-branch-switching-at-simple-branch-point/78721/4 "2022-03-31T10:15:05Z")
Thanks a lot for your replies!
This whole model is taken from [The effect of fishing on hysteresis in Caribbean coral reefs | SpringerLink](https://link.springer.com/article/10.1007/s12080-010-0102-0) (more or less what they call the original model). As stated, it’s the 3D model and equality T + M + C = 1 holds, where T,M,C are non-negative.
But for educational purposes, we simplified the model such that T = 1 - M - C. I’m not sure if it could result in something fishy like this. I will try and talk about it with my lecturer to see what’s going on!
Once again, thanks for BifurcationKit and your replies 😀
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [March 31, 2022, 10:32am UTC](https://discourse.julialang.org/t/bifurcationkit-jl-incorrect-automatic-branch-switching-at-simple-branch-point/78721/5 "2022-03-31T10:32:55Z")
Also, I forgot that you should use the jet for more precision:
```julia
br, = continuation(jet[1], jet[2], solution, parameters, (@lens _.g), options,
plot = true,
recordFromSolution = (x,p) -> (C = x[2], M = x[1])
)
```
---
### Author:  [@sceptri](https://discourse.julialang.org/u/sceptri)
#### Post date: [April 3, 2022, 2:01pm UTC](https://discourse.julialang.org/t/bifurcationkit-jl-incorrect-automatic-branch-switching-at-simple-branch-point/78721/6 "2022-04-03T14:01:13Z")
Hello, I was playing around with it and noticed, that while newton refinement returns eigenvalues/vectors in this order
```julia
# newton refinement of transcritical point
solt, = newton(jet[1:2]..., br1, 2) #[0.50909 2.4830e-16]
# jacobian at bifurcation point
J = jet[2](solt.u, @set parameters.g = solt.p)
vals, vecs = eigen(J)
# vals
# 2-element Vector{Float64}:
# -0.407272727272727
# -2.0146434165322638e-17
# vecs
# 2×2 Matrix{Float64}:
# -1.0 0.806874
# 8.84403e-16 -0.590723
```
In the continuation itself, they are in a different order
```julia
tcb2 = br1.specialpoint[2] # bifurcation point in question
real(br1.eig[tcb2.idx].eigenvals)
# eigen values
# 2-element Vector{Float64}:
# -1.2239298102576833e-5
# -0.4072538456186035
real(br1.eig[tcb2.idx].eigenvec)
# eigen vectors
# 2×2 Matrix{Float64}:
# 0.80688 -1.0
# -0.590715 6.57802e-5
```
I have no clue if this can be to blame, but nonetheless, I find it interesting.
More so, if you look and the second eigenvector of the continuation, it would result in the correct branch switch.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [April 3, 2022, 8:35pm UTC](https://discourse.julialang.org/t/bifurcationkit-jl-incorrect-automatic-branch-switching-at-simple-branch-point/78721/7 "2022-04-03T20:35:06Z")
No this is fine. BK sorts the eigenvalues by decreasing real part. See [here](https://github.com/rveltz/BifurcationKit.jl/blob/master/src/EigSolver.jl#L23)
There is something I dont get at the theory level.
---
### Author:  [@sceptri](https://discourse.julialang.org/u/sceptri)
#### Post date: [April 4, 2022, 6:12am UTC](https://discourse.julialang.org/t/bifurcationkit-jl-incorrect-automatic-branch-switching-at-simple-branch-point/78721/8 "2022-04-04T06:12:30Z")
Oh sorry, I did not know that. I’m still new to all of this.
Thanks for your replies!
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [April 6, 2022, 8:46pm UTC](https://discourse.julialang.org/t/bifurcationkit-jl-incorrect-automatic-branch-switching-at-simple-branch-point/78721/9 "2022-04-06T20:46:02Z")
I opened an issue [here](https://github.com/rveltz/BifurcationKit.jl/issues/55).
I think I know what is going on. Nothing is wrong in the computations. It is just that I assumed that the non trivial solution must correspond to the nontrivial solution of the normal form whereas it could well be, as in this case, that it is the trivial one. I will think about it.