incorrect automatic branch switching at simple branch point

[rveltz]

This is a duplicate of the mistake on discourse posted by @sceptri

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?

[rveltz]

This is a duplicate of the mistake on discourse posted by @sceptri

[rveltz]

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

[rveltz]

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:

# 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

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

[rveltz]

I understand the reason for the mistake now. I will try to fix it asap

[rveltz]

solved with f59fb7c