Continuation of Hopf curve for ComponentArray state type

[rseydam]

I tried to continue a Hopf curve using a ComponentArray.jl state type. As far as I can tell, there seems to be a problem when setting up the augmented problem. Here is a minimal working example that uses the example of the Hopf curve continuation from the tutorials (CO oxidation). I changed the state type to use the ComponentArray type.

using Revise, Plots
using BifurcationKit
using ComponentArrays

# function to record information from a solution
recordCO(x, p; k...) = (x = x.x[1], y = x.y[1], s = x.s[1])

# vector field of the problem
function COm(u, p)
	(;q1,q2,q3,q4,q5,q6,k) = p

    x = get_tmp(u.x, u)[1]
    y = get_tmp(u.y, u)[1]
    s = get_tmp(u.s, u)[1]

	z = 1-x-y-s
	out = similar(u)
	out[1] = 2q1 * z^2 - 2q5 * x^2 - q3 * x * y
	out[2] = q2 * z - q6 * y - q3 * x * y
	out[3] = q4 * z - k * q4 * s
	out
end

# parameters used in the model
par_com = (q1 = 2.5, q2 = 0.6, q3 = 10., q4 = 0.0675, q5 = 1., q6 = 0.1, k = 0.4)

# initial condition
z0 = ComponentArray(x=[0.07], y=[0.2], s=[05])

# Bifurcation Problem
prob = BifurcationProblem(COm, z0, par_com, (@optic _.q2); record_from_solution = recordCO)

# continuation parameters
opts_br = ContinuationPar(p_max = 1.9, dsmax = 0.01)

# compute the branch of solutions
br = continuation(prob, PALC(), opts_br; normC = norminf)

# plot the branch
scene = plot(br, xlims = (0.8,1.8))

hp_codim2 = continuation(br, 1, (@optic _.k),
	ContinuationPar(opts_br, p_max = 2.8, ds = -0.001, dsmax = 0.025) ;
	normC = norminf,
	# detection of codim 2 bifurcations
	detect_codim2_bifurcation = 2,
	# compute both sides of the initial condition
	bothside = true,
	)

# plotting
scene = plot(scene, hp_codim2, vars = (:q2, :x), branchlabel = "Hopf")
plot!(scene, br, xlims = (0.6, 1.5))

The error I get says:

ERROR: type Array has no field x

A similar problem might also appear in the fold curve continuation but I haven't checked. Any ideas on how to tackle this problem are very much appreciated.

[rveltz]

Hi,

You can try

hp_codim2 = continuation(br, 1, (@optic _.k),
	ContinuationPar(opts_br, p_max = 2.8, ds = -0.001, dsmax = 0.025) ;
    verbosity = 2,
	normC = norminf,
	# detection of codim 2 bifurcations
	detect_codim2_bifurcation = 2,
	# compute both sides of the initial condition
	bothside = true,
	jacobian_ma = :minaug,
	)

The issue is that using :autodiff for jacobian_ma generates a vcat that does not combine well with ComponentArrays. If the solution is good enough for you, please close the issue.

[rseydam]

This solution works well to do the continuation of fold and hopf curves with ComponentArrays state type. When trying to use a PeriodicOrbitOCollProblem starting from the hopf point one finds the same issue. How can it be solved for this problem?

opts_po_cont = ContinuationPar(opts_br, ds= 0.001, dsmin = 1e-4, dsmax = 0.1,
	max_steps = 110, tol_stability = 1e-5)

br_pocoll = continuation(
	br, 1, opts_po_cont,
	PeriodicOrbitOCollProblem(50, 4; meshadapt = true);
	)

[rveltz]

Yeah, this is not really gona work with collocation. I need to save the periodic orbit like (x1,y1,s1,x2,y2,s2,..., period) and this is not a component array.

You can either

• use shooting
• try MTK instead

[rseydam]

I also tried the continuation of the fold curve which works as well. Unfortunately, I couldn't get the periodic orbit continuation to work (also with shooting). Would it be easier to use RecursiveArrayTools.jl instead of ComponentArrays.jl or is it even safer to use just an Array?

[rveltz]

It depends. For low dimensional problems, I would try ModelingToolkit. For large dimensional problems, probably RecursiveArrayTools.jl but beware of #112

I have never tried ComponentArrays.jl.

[rseydam]

Thanks a lot! I will close this issue as the original problem using ComponentArrays types has been solved.