Parsing lens for codim2 bifurcation

Hi, I’m trying to perform a codimension2 continuation of a SN point, but I’m runnning into errors with parsing lens. I’ve continued the initial equilibrium solution for parameter A1, found some bifurcation points. Now I want to pick some of those points and continue them for different parameter pairs. I’m using Julia v1.11.4. The example code is 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] = A1*K1^n/(K1^n + y^n) - B1*x
    out[2] = A2*(m1*K2^n/(K2^n + z^n) + m2*K3^n/(K3^n + x^n)) - B2*y
    out[3] = A3*(m3*x^n/(K4^n + x^n) + m4*K5^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
)

The following error is generated

AssertionError: lens isa Int || lens isa AllOpticTypes

Stacktrace:
[1] BifurcationProblem(_F::typeof(ODE_3D), u0::Vector{Float64}, parms::@NamedTuple{A1::Float64, A2::Float64, A3::Float64, B1::Float64, B2::Float64, B3::Float64, m1::Float64, m2::Float64, m3::Float64, m4::Float64, K1::Float64, K2::Float64, K3::Float64, K4::Float64, K5::Float64, n::Float64}, lens::Tuple{PropertyLens{:A2}, PropertyLens{:A3}}; jvp::Nothing, vjp::Nothing, J::Symbol, J!::Nothing, Jᵗ::Nothing, d2F::Nothing, d3F::Nothing, issymmetric::Bool, record_from_solution::typeof(record3D), plot_solution::typeof(BifurcationKit.plot_default), delta::Float64, save_solution::typeof(BifurcationKit.save_solution_default), inplace::Bool, kwargs_jet::@Kwargs{})
@ BifurcationKit ~/.julia/packages/BifurcationKit/DXZuX/src/Problems.jl:282
[2] top-level scope
@ In[28]:2

Can anyone please help me setup this codim2 problem correctly. I’ve tried parsing lens in different ways but none have worked. I would also like to perform this continuation for other parameter pairs, so I appreciate any help.

The following steps would be
“br_codim2 = continuation(prob_codim2, PALC(tangent = BifurcationKit.Secant()), opts_codim2; verbosity = 2)”

Thank you

Hi,

That’s not how fold point are suppsed to be continued, it is more sophisticated. However, a simple interface for the user is provided. See 🟡 CO oxidation (codim 2) · Bifurcation Analysis in Julia

Hi, I followed that example in another step in the analysis to follow a fold in (A1, A2) - space. See

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

However, what I want to do is see what happens the fold point at A1 = 631.8055 in (A2, A3) - space. And I’m having thouble setting this problem up.

Regards,
Ielyaas

This is an interesting question. In case you have some time, you should open an issue on BifurcationKit.jl so we can find a proper solution. The proper solution would be to allow the user to provide two optics for continuation of bifurcation points, like

continuation(br, 1, (@optic _.A2), (@optic _.A3), ...)

for now, you can probably use the following hack to trick BK into thinking that the current parameter is A3. You have to provide the proper value for A1 as well and you can then call classical Fold continuation routine:

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

Thank you. I tried your suggestion, which didn’t produce any errors (phew!) but I didn’t find any bifurcation that parameter set in the given ranges. I’ll try a few other parameter pairs I’m interested in.

I’ll also open an issue in BifurcationKit.jl.

Thank you for suggestion.

And there should be ?

I’m not sure to be honest. I would have guessed that there’d be stability changes as A2 and A3 were varied but, there may not be for that bifurcation point. I’m trying others and other parameter pairs, but I’m currently running into convergence issues.