Performing codm2 continuation for multiple parameter pairs

[Ielyaas]

Hi All,

I've moved a discussion [(https://discourse.julialang.org/t/parsing-lens-for-codim2-bifurcation/127710)] from Julia discourse to here.

I want to continue a fold bifurcation from an equilibrium solution in a codim2 continuation for multiple different parameter pairs. So, in the example below I performed the initial continuation for for parameter A1, and found a few fold bifurcations, I would like to continue some of those fold bifurcations for other parameter pairs, for eg. (A2, A3). My initial try, which generated errors, is shown below

`# vector field
function ODE_3D(u, p)
(;A1,A2,A3,B1,B2,B3,
m1,m2,m3,m4,
K1,K2,K3,K4,K5,n) = p
x, y, z = u
out = similar(u)
out[1] = A1K1^n/(K1^n + y^n) - B1x
out[2] = A2*(m1K2^n/(K2^n + z^n) + m2K3^n/(K3^n + x^n)) - B2y
out[3] = A3(m3x^n/(K4^n + x^n) + m4K5^n/(K5^n + y^n)) - B3*z
out
end

Model parameters

params = (A1 = 0.0, A2 = 400.0, A3 = 450.0,
B1 = 1.0, B2 = 1.0, B3 = 1.0,
m1 = 1.0, m2 = 0.1, m3 = 1.0, m4 = 0.1,
K1 = 0.2, K2 = 0.2, K3 = 0.2, K4 = 0.2, K5 = 0.2, n = 2.0)

record3D(x, p; k...) = (x = x[1], y = x[2], z = x[3])

Initial conditions

u0 = [0.0, 400.0, 0.0];
prob = BifurcationProblem(ODE_3D, u0, params, (@optic _.A1); record_from_solution = record3D)

Continuation parameters

opts_br = ContinuationPar(
dsmin = 1e-6, # Minimum step size
dsmax = 1e-2, # Maximum step size
ds = 1e-3, # Initial step size
p_min = 0.0, # Minimum value of the parameter
p_max = 10000.0, # Maximum value of the parameter
max_steps = 1e9, # Maximum number of continuation steps
tol_stability = 1e-12, # Tolerance for stability computation
detect_fold = true, # Detect fold bifurcations
detect_bifurcation = 3 # Number of bifurcation detection methods
)

br = continuation(prob, PALC(), opts_br, bothside = true, normC = norminf)
┌─ Curve type: EquilibriumCont
├─ Number of points: 2493725
├─ Type of vectors: Vector{Float64}
├─ Parameter A1 starts at 0.0, ends at 10000.0
├─ Algo: PALC
└─ Special points:

• 1, endpoint at A1 ≈ +0.00000000, step = 0

• 2, bp at A1 ≈ +9872.77442260 ∈ (+9872.77442260, +9872.77442266), |δp|=6e-08, [converged], δ = ( 1, 0), step = 698179

• 3, bp at A1 ≈ +631.80555757 ∈ (+631.80555642, +631.80555757), |δp|=1e-06, [converged], δ = (-1, 0), step = 1351832

• 4, bp at A1 ≈ +2598.41066651 ∈ (+2598.41066602, +2598.41066651), |δp|=5e-07, [converged], δ = ( 1, 0), step = 1491233

• 5, bp at A1 ≈ +4.95752725 ∈ (+4.95752725, +4.95752726), |δp|=1e-08, [converged], δ = (-1, 0), step = 1676352

• 6, endpoint at A1 ≈ +10000.00000000,

Codim2 continuation

Extract the bifurcation points to continue

SN1 = br.specialpoint[3] # A1 = 631.80555757
SN2 = br.specialpoint[4] # A1 = 2598.41066651

p0 = setproperties(params, A1 = SN1.param)
u0 = SN1.x

lensA2 = (@optic _.A2)
lensA3 = (@optic _.A3)

Continuation in A2 vs A3 starting from the fold point

prob_codim2 = BifurcationProblem(
ODE_3D, u0, p0,
(lensA2, lensA3);
record_from_solution = record3D,
J = :finite
)
opts_codim2 = ContinuationPar(
dsmin = 1e-5,
dsmax = 1e-1,
ds = 1e-2,
max_steps = 5000,
p_min = 0.0,
p_max = 10000.0,
newton_options = NewtonPar(tol = 1e-8),
detect_bifurcation = 0,
tol_stability = 1e-10
)`

However, Romain Veltz suggested the following workaround which didn't generate any errors
`br2 = @set br.prob.lens = @optic _.A3
@reset br.prob.params.A1 = 631.8055

sn_codim2_A2 = continuation(br2, 2, (@optic _.A2),
ContinuationPar(opts_br, p_min = 0.0, p_max = 10000.0, dsmin = 1e-9, max_steps = 1e3, tol_stability = 1e-15);
normC = norminf,
bothside = true,
detect_codim2_bifurcation = 2,)`

I was wondering if there's a proper way to perform such a codim2 continuation, soomething like continuation(br, 1, (@optic _.A2), (@optic _.A3), ...).

Thank you

[rveltz]

see #229