# BifurcationKit.jl incorrect automatic branch switching at simple branch point

Hello, I have a simple 2D system and I tried continuation from simple branch points

``````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

br, = continuation(f, solution, parameters, (@lens _.g), options,
plot = true,
recordFromSolution = (x,p) -> (C = x, M = x)
)

br1, = continuation(jet..., br, 1, setproperties(options;ds = 0.0001, maxSteps = 30), recordFromSolution = (x,p) -> (C = x, M = x), bothside = true, plot = true)
br2, = continuation(jet..., br1, 2, setproperties(options; maxSteps = 20, ds = 0.0001), recordFromSolution = (x,p) -> (C = x, M = x), 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?

Hi,

You are right, something is wrong with the predictor from the transcitical bifurcation. In case you need it, this works:

``````parameters2 = @set parameters.g=br1.specialpoint.param + 0.01
solution2 = br1.specialpoint.x
br3, = continuation(f, solution2, parameters2, (@lens _.g), options,
plot = true,
recordFromSolution = (x,p) -> (C = x, M = x),
plotSolution = (x, p; k...) -> begin
plot!(br1, subplot = 1)
end,
)
``````
1 Like

You can use deflation to help it converge to another branch:

``````br1, = continuation(jet..., br, 1, setproperties(options;ds = 0.0001, maxSteps = 30), recordFromSolution = (x,p) -> (C = x, M = x), bothside = true, plot = true)
br2, = continuation(jet..., br1, 2, setproperties(options; maxSteps = 20, ds = 0.0001), recordFromSolution = (x,p) -> (C = x, M = x), bothside = true,plot = true, usedeflation = true)
plot(br,br1,br2)
``````

However, I am puzzled by the result from bifurcation theory:

``````# newton refinement of transcritical point
solt, = newton(jet[1:2]..., br1, 2) #[ 0.50909 2.4830e-16]
# jacobian at bifurcation point
J = jet(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

1 Like

Thanks a lot for your replies!

This whole model is taken from The effect of fishing on hysteresis in Caribbean coral reefs | SpringerLink (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 Also, I forgot that you should use the jet for more precision:

``````br, = continuation(jet, jet, solution, parameters, (@lens _.g), options,
plot = true,
recordFromSolution = (x,p) -> (C = x, M = x)
)
``````
1 Like

Hello, I was playing around with it and noticed, that while newton refinement returns eigenvalues/vectors in this order

``````# newton refinement of transcritical point
solt, = newton(jet[1:2]..., br1, 2) #[ 0.50909 2.4830e-16]
# jacobian at bifurcation point
J = jet(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

``````tcb2 = br1.specialpoint # 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.

No this is fine. BK sorts the eigenvalues by decreasing real part. See here

There is something I dont get at the theory level.

1 Like

Oh sorry, I did not know that. I’m still new to all of this.