Incorrect bifurcation diagram bifurcationkit.jl

Sergey_Novak · September 14, 2022, 8:04am

This is my first time using this package and I need to get a bifurcation diagram for the following system

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

a = 1.0; b = 3.0; c = 1.0; d = 5.0;
xr = -1.6; r = 0.01; s = 5.; I = 4.; xv = 2.;
k1= -0.17; k2 = -0.17;
k = 0.0

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

I following example from documentation and use next setting for NewtonPar and ContinuationPar

prob_HR = BifurcationProblem(HR, zeros(6), par, (@lens _.k),
                                J = (x, p) -> ForwardDiff.jacobian(z -> HR(z, p), x))
# newton options
opt_newton = NewtonPar(tol = 1e-11, maxIter = 10000)
# continuation options
opts_br = ContinuationPar(
                        pMax = 0.5, pMin = 0.0,
                        detectBifurcation = 3, nev = 6,
                        newtonOptions = opt_newton, maxSteps = 10000, nInversion = 4,
                        tolBisectionEigenvalue = 1e-11, dsminBisection = 1e-11)

bdiag = bifurcationdiagram(prob_HR, PALC(), 3,
	(args...) -> setproperties(opts_br; pMin = 0.0, pMax = 0.5, nInversion = 4, detectBifurcation = 3, maxBisectionSteps=100, newtonOptions = opt_newton);
	recordFromSolution = (x, p) -> norminf(x),
	xwnormC = norminf)

In result i get

In short, I know that there are bifurcations with a change in the parameter `k’. I have a diagram of the spectrum of Lyapunov exponents, projections of phase space with different values of k, and they show that there is a transition from chaos to regularity.

I understand that the diagram is most likely not correct and due to the wrong choice of settings. How to choose the right parameters and are there any recommendations for this?

I also attach a spectrum diagram and a projection of the phase space so that it is clearly visible that the diagram does not look quite correct. (The blue and red trajectories correspond to two different attractors)

RayleighLord · September 14, 2022, 8:54am

I have not used this package much, but I have a question of what you are trying to compute.

Have you identified a stationary solution from which you want to compute the new branches that may emerge? What is your approach to the problem?

Usually, to compute bifurcation diagrams, we have found a solution which is the baseline when constructing the bifurcation diagram. Although I think that BifurcationKit.jl has some implementations to help you construct the bifurcation diagram even for disconnected branches, you may have to think a little bit a about your problem before giving it to this package.

Sergey_Novak · September 14, 2022, 9:30am

I didn’t see any information about this in the documentation examples, I didn’t look at all the examples, of course, only what was necessary.
I have problems even when I try to use the code from this page for my system
https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/ode/lorenz84/#Extended-Lorenz-84-model-(codim-2-BT/ZH-aBS)

dawbarton · September 14, 2022, 10:03pm

BifurcationKit.jl uses numerical continuation to produce the bifurcation diagram and numerical continuation can only find connected branches of solutions. In this case, your run of BifurcationKit finds an initial branch of unstable steady-states which doesn’t undergo any bifurcations. It might find other branches if you vary the initial conditions. That said, the code you’ve used will only find steady-state solutions though (equilibria) - you need a slightly different code to find/track periodic orbits. (BifurcationKit cannot track chaotic attractors directly.)

The code I used to check that the branch of solutions computed is unstable is

using NonlinearSolve
using ForwardDiff
prob = NonlinearProblem(HR, zeros(6), par)  # root finding
sol = solve(prob, NewtonRaphson())  # find the steady-state
J = ForwardDiff.jacobian(u -> HR(u, par), sol.u)  # find it's linearised stability
isstable = all(real.(eigen(J).values) .< 0)  # check all eigenvalues have negative real part

Given the plots you’ve made of the stable orbits you see from simulation, you might be best to do a brute-force bifurcation diagram, at least to start with. That is, use Monte-Carlo simulation over the parameter range you are interested in and plot some summary statistics against the bifurcation parameter.

rveltz · September 15, 2022, 6:40am

As said above, I dont think there is anything wrong. BK just picked up the closest branch of solutions to your guess and continued it for k in [0, 0.5] as you asked. If you want the solutions to be saved, you can use ContinuationPar(saveSolEveryStep = 1) and the solutions will be in br.sol; this way it will be easier to check.

Also, the computed branch is unstable because

julia> bdiag.γ.stable
9-element Vector{Bool}:
 0
 0
 0
 0
 0
 0
 0
 0
 0

The eigenelements are in bdiag.γ.eig.

Usually, it is recommanded to use continuation with the option verbosity = 3 to get a feeling of the system and then use bifurcationdiagram.

BifurcationKit uses AD, so you can just do

prob_HR = BifurcationProblem(HR, zeros(6), par, (@lens _.k))

You should not use

because the number of Newton iterations is supposed to be small. Hence:

# newton options
opt_newton = NewtonPar(tol = 1e-11, maxIter = 10)

but then I understand why you did this, because newton does not converge to your guess in 10 iterations.
I changed the initial guess to ones .Nevertheless, if you ask to go to k=2, it finds a Hopf bifurcation.

prob_HR = BifurcationProblem(HR, ones(6), par, (@lens _.k);
            recordFromSolution = (x, p) -> (n = norminf(x), x = x[1]))
# newton options
opt_newton = NewtonPar(tol = 1e-11, maxIter = 10)

# continuation options
opts_br = ContinuationPar(
                        pMax = 0.5, pMin = 0.0,
                        detectBifurcation = 3, nev = 6,
                        newtonOptions = opt_newton, maxSteps = 10000, nInversion = 4,
                        tolBisectionEigenvalue = 1e-11, dsminBisection = 1e-11)

br = continuation(prob_HR, PALC(), setproperties(opts_br;  pMax = 1.5, newtonOptions = opt_newton);
    verbosity =3,
    normC = norminf)

and gives

Predictor:  Secant
 ┌─ Number of points: 16
 ├─ Curve of EquilibriumCont
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter k starts at 0.0, ends at 1.5
 ├─ Algo: PALC
 └─ Special points:

If `br` is the name of the branch,
ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`

- #  1,     hopf at k ≈ +1.23964658 ∈ (+1.23080774, +1.23964658), |δp|=9e-03, [converged], δ = (-2, -2), step =  13, eigenelements in eig[ 14], ind_ev =   4
- #  2, endpoint at k ≈ +1.50000000,                                                                     step =  15

Still the whole branch is unstable

Now, you can use deflated continuation to see if there are other equilibria:

alg = DefCont(deflationOperator = DeflationOperator(2, dot, 1., [ones(6)]), perturbSolution = (x,p,id) -> (x  .+ 0.1 .* rand(length(x))), maxIterDefOp = 10)

brdc = continuation(prob_HR, alg,
	ContinuationPar(opts_br, dsmin = 1e-6,ds = 0.001, pMax = 3.5, maxSteps = 20000, newtonOptions = NewtonPar(verbose = false, maxIter = 10), plotEveryStep = 400),
	; plot=true, verbosity = 0,
	callbackN = BifurcationKit.cbMaxNorm(1e7) # reject newton step if residual too large
	)

and the result seems to be no: you cannot be sure even with Deflated Continuation although it is a good indication.

Hence, I would say that chaos unfolds from the Hopf point, the bifurcated branch of periodic orbits PO must undergo some cascade of bifurcations. You can follow the branch of PO in BifurcationKit to investigate what is going on

Sergey_Novak · September 15, 2022, 7:16am

Thank you!
I got the same result in matcont that you have. I have configure the parameters a little differently and got a picture that doesn’t match anything.

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
a = 1.0; b = 3.0; c = 1.0; d = 5.0;
xr = -1.6; r = 0.01; s = 5.0; I = 4.; xv = 2.;
k1= -0.17; k2 = -0.17;
k = 0.0

par = (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 , el_link = 0.0)

prob_HR = BifurcationProblem(HR, zeros(6), par, (@lens _.el_link),
                                J = (x, p) -> ForwardDiff.jacobian(z -> HR(z, p), x))

# newton options
opt_newton = NewtonPar(tol = 1e-11, maxIter = 10000)
# continuation options
opts_br = ContinuationPar(
                        pMax = 1.0, pMin = -1.0,
                        detectBifurcation = 3, nev = 6,
                        newtonOptions = opt_newton, maxSteps = 10000, nInversion = 4,
                        tolBisectionEigenvalue = 1e-11, dsminBisection = 1e-11)

bdiag = bifurcationdiagram(prob_HR, PALC(), 3,
	(args...) -> setproperties(opts_br; pMin = -1.0, pMax = 1.0, nInversion = 4, detectBifurcation = 3, maxBisectionSteps=100, newtonOptions = opt_newton);
	recordFromSolution = (x, p) -> norminf(x),
	xwnormC = norminf)

Then I used as z0 equilibrium state coordinates obtained in the Dynamical Systems.jl and slightly changed the code

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

function HR!(dz, z, p, t)
    
    @unpack a, b, c, d, s, xr, r,  I, vs, k1, k2, el_link  = p
    x1, y1, z1, x2, y2, z2 = z
    
    dz[1] = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k1 * ( x1 - vs ) * sigma(x2) + el_link * ( x2 - x1 )
    dz[2] = c - d * x1 ^2 - y1
    dz[3] = r * ( s * ( x1 - xr ) - z1 )
    
    dz[4] = y2 + b * x2 ^ 2 - a * x2 ^3 - z2 + I - k2 * ( x2 - vs ) * sigma(x1) + el_link * ( x1 - x2 )
    dz[5] = c - d * x2 ^2 - y2
    dz[6] = r * ( s * ( x2 - xr ) - z2 )
    dz
end

HR(z, p) = HR!(similar(z), z, p, 0)
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 , el_link = 0.0)
z0 = [-0.738229,  -1.72491,  4.30886, -0.738229,  -1.72491,  4.30886]
prob_ = BifurcationProblem(HR, z0, param, (@lens _.el_link);
                            recordFromSolution = (x, p) -> (x1 = x[1], y1 = x[2], z1 = x[3],
                                                            x2 = x[4], y2 = x[5], z2 = x[6]))
opts_br = ContinuationPar(pMin = 0.0, pMax = 1.0,
	# parameters to have a smooth result
	ds = 0.04, dsmax = 0.05,
	# this is to detect bifurcation points precisely with bisection
	detectBifurcation = 3,
	# Optional: bisection options for locating bifurcations
	nInversion = 8, maxBisectionSteps = 25, nev = 6,
    detectFold = true);
br = continuation(prob_, PALC(tangent=Bordered()), opts_br; plot = true, normC = norminf, bothside = true)
br1 = continuation(prob_, PALC(), opts_br, detectCodim2Bifurcation = 2;
	plot = false, verbosity = 3, normC = norminf, bothside = true)

and get result above. What could be the reason for such a difference. Am I using the methods incorrectly or is it because of their configuration parameters?

rveltz · September 15, 2022, 8:01am

Euh … what do you mean? the continuation parameter is now el_link

Can you be more explicit? what do you want to do with br1?

Sergey_Novak · September 15, 2022, 8:07am

It’s the same parameter I just renamed it. I was checking what result br1 would give. For br i used PALC(tangent=Bordered()) and for br1 PALC() and wanted to compare

rveltz · September 15, 2022, 8:20am

Sergey_Novak:

br = continuation(prob_, PALC(tangent=Bordered()), opts_br; plot = true, normC = norminf, bothside = true)
br1 = continuation(prob_, PALC(), opts_br, detectCodim2Bifurcation = 2;
	plot = false, verbosity = 3, normC = norminf, bothside = true)

I get the same results, maybe restart Julia (Note that the option detectCodim2Bifurcation = 2 is not used, it is for a different context)

Sergey_Novak · September 15, 2022, 8:33am

I restart Julia this is not help. I can drop notebook jupyter with code. Also i want ask question about why can not find bifurcations at the parameter when a chaotic attractor is destroyed, and the destruction occurs due to fusion with another attractor. This is at least visible visually. I don’t think the problem is with the integrator, because with it I got the same results as in the article

rveltz · September 15, 2022, 8:47am

Can you try master ] add BifurcationKit#master please? This is the code

using Revise, BifurcationKit, Setfield, LinearAlgebra, Plots, Parameters

norminf(x) = norm(x, Inf)

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

function HR!(dz, z, p, t)

    @unpack a, b, c, d, s, xr, r,  I, vs, k1, k2, el_link  = p
    x1, y1, z1, x2, y2, z2 = z

    dz[1] = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k1 * ( x1 - vs ) * sigma(x2) + el_link * ( x2 - x1 )
    dz[2] = c - d * x1 ^2 - y1
    dz[3] = r * ( s * ( x1 - xr ) - z1 )

    dz[4] = y2 + b * x2 ^ 2 - a * x2 ^3 - z2 + I - k2 * ( x2 - vs ) * sigma(x1) + el_link * ( x1 - x2 )
    dz[5] = c - d * x2 ^2 - y2
    dz[6] = r * ( s * ( x2 - xr ) - z2 )
    dz
end

HR(z, p) = HR!(similar(z), z, p, 0)
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 , el_link = 0.0)
z0 = [-0.738229,  -1.72491,  4.30886, -0.738229,  -1.72491,  4.30886]
prob_ = BifurcationProblem(HR, z0, param, (@lens _.el_link);
                            recordFromSolution = (x, p) -> (x1 = x[1], y1 = x[2], z1 = x[3],
                                                            x2 = x[4], y2 = x[5], z2 = x[6]))
opts_br = ContinuationPar(pMin = 0.0, pMax = 1.0,
	# parameters to have a smooth result
	ds = 0.04, dsmax = 0.05,
	# this is to detect bifurcation points precisely with bisection
	detectBifurcation = 3,
	# Optional: bisection options for locating bifurcations
	nInversion = 8, maxBisectionSteps = 25, nev = 6,
    detectFold = true);

br = continuation(prob_, PALC(tangent=Bordered()), opts_br; plot = true, normC = norminf, bothside = true)
br1 = continuation(prob_, PALC(), opts_br, detectCodim2Bifurcation = 2;
		plot = false, verbosity = 3, normC = norminf, bothside = true)

plot(br)
plot(br1)

The Hopf bifurcation point is at k = 0.61, this is the right part on your lyapunov panel. Then there is destruction of chaotic attractor around k=0.15 but this is not associated with a bifurcation of an equilibrium. Presumably, the branch of periodic orbits (that you need to follow with BifurcationKit) will grow for k smaller than k=0.61 and undergoes a cascade of bifurcations when it reaches k=0.15

Sergey_Novak · September 15, 2022, 9:00am

Do you mean to run this code by installing the package from the master? If so, I got the following

rveltz · September 15, 2022, 9:06am

yes, it seems the same to me. You can try show(br) and show(br1) to be sure

Sergey_Novak · September 15, 2022, 9:09am

Yes, you are right. However I don’t understand why my result is different from yours

rveltz · September 15, 2022, 9:10am

I dont understand this. What is different now?

Sergey_Novak · September 15, 2022, 9:12am

I about this

rveltz · September 15, 2022, 9:13am

Ah I dont know, the vector fields must be different from your two implementations.

Sergey_Novak · September 15, 2022, 11:24am

I understood the difference. The order of transmission of parameters. Thank you

Sergey_Novak · September 15, 2022, 12:09pm

Can i ask how correctly calculate codim 2 biffurcation map. I use example from documentation

prob_HR = BifurcationProblem(HR, ones(6), param, (@lens _.k);
            recordFromSolution = (x, p) -> (n = norminf(x), x = x[1]))
# newton options
opt_newton = NewtonPar(tol = 1e-16, maxIter = 100)

# continuation options
opts_br = ContinuationPar(
                        pMax = 1.0, pMin = 0.0,
                        detectBifurcation = 3, nev = 6,
                        newtonOptions = opt_newton, maxSteps = 100, nInversion = 4,
                        tolBisectionEigenvalue = 1e-11, dsminBisection = 1e-11)

br = continuation(prob_HR, PALC(), setproperties(opts_br; newtonOptions = opt_newton);
    verbosity =3,
    normC = norminf)
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 save the first component for plotting
	recordFromSolution = (x, p) -> (n = norminf(x), x = x[1]),
	# we update the Hopf problem at every continuation step
	updateMinAugEveryStep = 1,
	# compute both sides of the initial condition
	bothside = true,
	)

with this code i have error

AssertionError: Newton failed to converge. Required for the computation of the initial tangent.

rveltz · September 15, 2022, 12:30pm

sorry I am lost with all your versions. can you paste the vector field and paraemeters?