Incorrect bifurcation diagram bifurcationkit.jl

[Sergey_Novak]

I apologize for the many versions vf and parametres

function sigma(x)
    return 1.0 / ( 1.0 + exp( -10.0 * ( x  - ( - 0.25 ) ) ) )
end

function HR(u, p)
    
    a, b, c, d, s, xr, r,  I, vs, k1, k2, el_link  = p
    x1, y1, z1, x2, y2, z2 = u
    
    du1 = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k1 * ( x1 - vs ) * sigma(x2) + el_link * ( x2 - x1 )
    du2 = c - d * x1 ^2 - y1
    du3 = r * ( s * ( x1 - xr ) - z1 )
    
    du4 = y2 + b * x2 ^ 2 - a * x2 ^3 - z2 + I - k2 * ( x2 - vs ) * sigma(x1) + el_link * ( x1 - x2 )
    du5 = c - d * x2 ^2 - y2
    du6 = r * ( s * ( x2 - xr ) - z2 )
    return [du1, du2, du3,
            du4, du5, du6]
end

param =  (a = 1.0, b = 3.0 , c = 1.0, d = 5.0, s = 5.0, xr = -1.6, r = 0.01,  I = 4.0, vs = 2.0, k1 = -0.17, k2 = -0.17 , k = 0.0)

[rveltz]

@Sergey_Novak

opt_newton = NewtonPar(tol = 1e-16, maxIter = 100)

You are asking too much here, make it

opt_newton = NewtonPar(tol = 1e-12, maxIter = 10)

You can easily debug these errors by going verbose like opt_newton = NewtonPar(tol = 1e-12, maxIter = 10, verbose = true)

Also, you cannot pass recordFromSolution to codim 2 (you could but you have to dispatch it on BorderedArrays), hence

hp_codim2 = continuation(br, 1, (@lens _.k1),
	ContinuationPar(opts_br, pMin = -0.5, pMax = 0.5,
		ds = -0.001, dsmax = 0.05) ;
	normC = norminf,
	# detection of codim 2 bifurcations with bisection
	detectCodim2Bifurcation = 2,
	# tell to start the Hopf problem using eigen elements: compute left eigenvector
	startWithEigen = true,
	# we update the Hopf problem at every continuation step
	updateMinAugEveryStep = 1,
	# compute both sides of the initial condition
	bothside = true,
	)

[Sergey_Novak]

Thank you! It work and i get next result
is it possible to additionally make sure that the bifurcation is found correctly?

[rveltz]

you can do getNormalForm(hp_codim2, 2)

[Sergey_Novak]

Thank you
I get the following result, and as I understand it, it is incorrect because real numbers of eigenvalues not zero, perhaps this is due to accuracy

[rveltz]

you are talking about G21 and G32? I did not take the real parts, see Bautin · Bifurcation Analysis in Julia.

For codim 1,2 I take the normal form for the vector field and see if my normal form computation gives the same. See hopf test, transcritical, pitchfork D6, cusp, bogdanov-takens, bautin

You can have access to more than what is given by show(nf) by looking at the fields of nf where nf is the normal form nf = getNormalForm(br, ind::Int).

In a nutshell, I try to test as much as I can to test the code against analytics although the above tests are not perfect because they do not involve the computatation of the center manifold since it is the whole space. But that’s already a test ;D

I am happy to know of a bug here because it will become one less in the package

[Sergey_Novak]

Okay, thank you. I wonder what this may indeed be a bautin bifurcation. I test critical parametres and see next attractor. the Lyapunov spectrum and the eigenvalues indicate that this is indeed a bautin bifurcation if of course believe the results of the calculation. Before critical value i have not zero real numbers but i have zero largest lyapunov exponent.
for calculate fixed point and spectrum uses DynamicalSystems.jl

[rveltz]

You are making a link between unrelated behaviors to my mind. It is more sophisticated than just a bauti and likely involve periodic orbits and torus.

[Sergey_Novak]

perhaps this is so, I hardly understand what is happening in the system, because there is only one state of equilibrium, which practically does not affect the dynamics in any way
I’ll attach another projection. perhaps you know what it is

[Sergey_Novak]

I take parameter before critical value and get this. could it be unrelated to the state of equilibrium?

[rveltz]

I’d compute the curve of periodic orbits from the Hopf point