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