BifurcationKit.jl incorrect automatic branch switching at simple branch point

General Usage
1

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[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?

2

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[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,
	)
1 Like
3

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

Screen Shot 2022-03-30 at 20.38.57
Screen Shot 2022-03-30 at 20.38.57962×648 26.7 KB

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[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 :thinking: from a theoritical point of view

1 Like
4

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 :grinning:

5

Also, I forgot that you should use the jet for more precision:

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

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

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.

7

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
8

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

Thanks for your replies!

9

I opened an issue here.

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.

1 Like