# Marking user points / events along a cxdimension-one curve of fold bifurcation
**URL:** https://discourse.julialang.org/t/marking-user-points-events-along-a-cxdimension-one-curve-of-fold-bifurcation/114325
**Category:** Modelling & Simulations
**Tags:** dynamical-systems
**Created:** [May 15, 2024, 9:54pm UTC](https://discourse.julialang.org/t/marking-user-points-events-along-a-cxdimension-one-curve-of-fold-bifurcation/114325 "2024-05-15T21:54:44Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@BigEpsilon.jl](https://discourse.julialang.org/u/BigEpsilon.jl)
#### Post date: [May 15, 2024, 9:54pm UTC](https://discourse.julialang.org/t/marking-user-points-events-along-a-cxdimension-one-curve-of-fold-bifurcation/114325/1 "2024-05-15T21:54:44Z")
This post is concerning BifurcationKit.jl
I doing a bifurcation analysis in two parameters, and want to be able to mark a point along a fold-curve, similar to how it is on in AUTO-07p.
I have gotten this working for a one-parameter continuation using a ContinuousEvent function. However, the same method does not detect user-defined points along a codimesion-one curve of fold bifurcation. My continuation parameters and continuation call is as follows
opts = ContinuationPar(
p\_min = -3.0,
p\_max = 0.0,
max\_steps=5000,
dsmin=1e-15,
ds=1e-4,
dsmax=2e-4,
detect\_bifurcation=0,
nev=4,
detect\_fold = false,
detect\_event = 1,
)
sn\_codim2 = continuation(branch, indFold, (@lens \_.eta), J=JAC,
opts;
bothside=true,
jacobian = :dense,
jacobian\_ma = :AutoDiffDenseAnalytical,
bdlinsolver = MatrixBLS(),
update\_minaug\_every\_step = 1,
detect\_codim2\_bifurcation = 0,
start\_with\_eigen = true,
detect\_loop=1,
event = BK.ContinuousEvent(1,
(iter, state) → (getp(state)+Target))
)
)
I want it to detect a point at “Target”. It simply does not detect any event during the continuation run.
Please help? \>.\<
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 16, 2024, 6:49am UTC](https://discourse.julialang.org/t/marking-user-points-events-along-a-cxdimension-one-curve-of-fold-bifurcation/114325/2 "2024-05-16T06:49:02Z")
Hi,
What version of BifurcationKit are you using?
For example, the option `jacobian_ma = :AutoDiffDenseAnalytical` is not valid. Also, `AutoDiffDenseAnalytical()` should be used instead but this is for PeriodicOrbits so I dont understand your code.
`jacobian = :dense` is strange to, it is not in the arguments of `continuation` for a codim 2.
Lastly, do you have a MWE of the problem so we can investigate this?
---
### Author:  [@BigEpsilon.jl](https://discourse.julialang.org/u/BigEpsilon.jl)
#### Post date: [May 16, 2024, 10:28am UTC](https://discourse.julialang.org/t/marking-user-points-events-along-a-cxdimension-one-curve-of-fold-bifurcation/114325/3 "2024-05-16T10:28:35Z")
Hello,
I am using version “BifurcationKit v0.3.3” according to Julia. I agree my code does not make much sense, I am new to BifurcationKit.jl and was/am fumbling a bit.
I have removed the incorrect “jacobian\_ma " and “jacobian” arguments that you pointed out. However, I now run into an additional problem – whenever I do not set " jacobian\_ma” to anything, even when incorrect, I get the following error:
```julia
ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{BifurcationKit.var"#352#353"{BifurcationKit.FoldMAProblem{…}, @NamedTuple{…}}, Float64}, Float64, 5})
Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:207
(::Type{T})(::T) where T<:Number
@ Core boot.jl:792
Float64(::IrrationalConstants.Fourπ)
@ IrrationalConstants ~/.julia/packages/IrrationalConstants/vp5v4/src/macro.jl:112
...
Stacktrace:
[1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{BifurcationKit.var"#352#353"{…}, Float64}, Float64, 5})
@ Base ./number.jl:7
[2] setindex!
@ ./array.jl:1024 [inlined]
[3] JAC(u::SubArray{ForwardDiff.Dual{…}, 1, Vector{…}, Tuple{…}, true}, p::@NamedTuple{mu::ForwardDiff.Dual{…}, eta::Float64, k2::Float64})
@ Main ~/Library/CloudStorage/OneDrive-TheUniversityofAuckland/Git/psweep_forced_julia/BailieKeane_2024_nonauto/Test:40
[4] jacobian(pb::BifFunction{…}, x::SubArray{…}, p::@NamedTuple{…})
@ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Problems.jl:72
[5] jacobian(pb::BifurcationProblem{…}, x::SubArray{…}, p::@NamedTuple{…})
@ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Problems.jl:234
[6] (::FoldProblemMinimallyAugmented{…})(x::SubArray{…}, p::ForwardDiff.Dual{…}, params::@NamedTuple{…})
@ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:35
[7] (::FoldProblemMinimallyAugmented{…})(x::Vector{…}, params::@NamedTuple{…})
@ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:48
[8] (::BifurcationKit.var"#352#353"{BifurcationKit.FoldMAProblem{…}, @NamedTuple{…}})(z::Vector{ForwardDiff.Dual{…}})
@ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:172
[9] vector_mode_dual_eval!(f::BifurcationKit.var"#352#353"{…}, cfg::ForwardDiff.JacobianConfig{…}, x::Vector{…})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/apiutils.jl:24
[10] vector_mode_jacobian(f::BifurcationKit.var"#352#353"{…}, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:125
[11] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{…}, Float64, 5, Vector{…}}, ::Val{true})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:0
[12] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{…}, Float64, 5, Vector{…}})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:19
[13] jacobian
@ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:172 [inlined]
[14] _newton(prob::BifurcationKit.FoldMAProblem{…}, x0::Vector{…}, p0::@NamedTuple{…}, options::NewtonPar{…}; normN::typeof(norminf), callback::typeof(BifurcationKit.cb_default), kwargs::@Kwargs{…})
@ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Newton.jl:88
[15] _newton
@ ~/.julia/packages/BifurcationKit/5kgAI/src/Newton.jl:62 [inlined]
[16] newton
@ ~/.julia/packages/BifurcationKit/5kgAI/src/Newton.jl:138 [inlined]
[17] iterate(it::ContIterable{…}; _verbosity::Int64)
@ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:271
[18] iterate
@ ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:250 [inlined]
[19] continuation(it::ContIterable{…})
@ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:474
[20] continuation(prob::BifurcationKit.FoldMAProblem{…}, alg::PALC{…}, contparams::ContinuationPar{…}; linear_algo::BorderingBLS{…}, bothside::Bool, kwargs::@Kwargs{…})
@ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:539
[21] continuation_fold(prob::BifurcationProblem{…}, alg::PALC{…}, foldpointguess::BorderedArray{…}, par::@NamedTuple{…}, lens1::Setfield.PropertyLens{…}, lens2::Setfield.PropertyLens{…}, eigenvec::Vector{…}, eigenvec_ad::Vector{…}, options_cont::ContinuationPar{…}; update_minaug_every_step::Int64, normC::typeof(norminf), bdlinsolver::MatrixBLS{…}, bdlinsolver_adjoint::MatrixBLS{…}, jacobian_ma::Symbol, compute_eigen_elements::Bool, usehessian::Bool, kind::BifurcationKit.FoldCont, record_from_solution::Nothing, kwargs::@Kwargs{…})
@ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:497
[22] continuation_fold
@ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:314 [inlined]
[23] continuation_fold(prob::BifurcationProblem{…}, br::ContResult{…}, ind_fold::Int64, lens2::Setfield.PropertyLens{…}, options_cont::ContinuationPar{…}; alg::PALC{…}, normC::typeof(norminf), nev::Int64, start_with_eigen::Bool, bdlinsolver::MatrixBLS{…}, bdlinsolver_adjoint::MatrixBLS{…}, kwargs::@Kwargs{…})
@ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:551
[24] continuation_fold
@ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:511 [inlined]
[25] #continuation#348
@ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/codim2.jl:237 [inlined]
[26] top-level scope
@ ~/Library/CloudStorage/OneDrive-TheUniversityofAuckland/Git/psweep_forced_julia/BK_2024_nonauto/Test:195
[27] include(fname::String)
@ Base.MainInclude ./client.jl:489
[28] top-level scope
@ REPL[49]:1
in expression starting at /Users/johnbailie/Library/CloudStorage/OneDrive-TheUniversityofAuckland/Git/psweep_forced_julia/BK_2024_nonauto/Test:195
Some type information was truncated. Use `show(err)` to see complete types.
```
On the other hand, when I set the “jacobian\_ma = :a” for example, the correct branch of fold bifurcation is computed. I am honestly really lost on this, I believed I followed the documentation correctly?
Any assistance would be really appreciated.
FYI, this is my code
```julia
function JAC(u, p)
m = p[1]
eta = p[2]
k2 = p[3]
xn, yn, yl, xb = u
### Set main constants
SV=10^6;
sig=2.1*10^4*SV;
alphaT=1.7*10^(-4);
alphaS=0.8*10^(-3);
ep=0.011;
TaN=7;
TaL=25
T0=2.65
S0 =35;
SB =34.538;
yb = alphaS/alphaT * (SB - S0)/(TaN - T0);
g=3.17*10^-8;
VN=7.2106*10^15 #7*10^16;
VL=6.3515 * 10^16
VB=25 * VN
phi=alphaT*(TaN-T0);
xl=(TaL-T0)/(TaN-T0);
#
TildeVL=VL/VN;
TildeVB=VB/VN;
K=5.456*SV/sig;
ps=phi*((yn-xn)-(yl-xl));
d = g*VN/sig
# Functions
k2n = k2;
k2l = k2/3
k1 = 10^-2*SV/sig;
J = zeros(4,4)
# FIRST ROW
J[1,1] = (1/ep)*((0.5*k1*xb - 0.5*k2n*xb - 0.5*k1*xn + 0.5*k2n*xn)*sech((eta + xn - yn)/ep)^2 + phi^2*(-xb*xl - xl^2 + xb*xn + 3*xl*xn - 2*xn^2 +
xb*yl + xl*yl - 2*xn*yl - xb*yn - xl*yn + 2*xn*yn)*sech((phi*(xl - xn - yl + yn))/ep)^2 + ep*(-d - K - 0.5*k1 -
0.5*k2n + (-0.5*k1 + 0.5*k2n)*tanh((eta + xn - yn)/ep) + phi*(-xb - 3*xl + 4*xn + 2*yl - 2*yn)*tanh((phi*(xl - xn - yl + yn))/ep)))
J[1,2] = (1/ep)*((-0.5*k1*xb + 0.5*k2n*xb + 0.5*k1*xn - 0.5*k2n*xn)*sech((eta + xn - yn)/ep)^2 + phi*sech((phi*(xl - xn - yl + yn))/ep)^2*(phi*(xl^2
+ xn*(2*xn + 2*yl - 2*yn) + xl*(-3*xn - yl + yn) + xb*(xl - xn - yl + yn)) + ep*(0.5*xb + 0.5*xl - xn)*sinh((2*phi*(xl - xn - yl + yn))/ep)))
J[1,3] = -((phi*(xb + xl - 2*xn)*sech((phi*(xl - xn - yl + yn))/ep)^2*(2*phi*(xl - xn - yl + yn) + ep*sinh((2*phi*(xl - xn - yl + yn))/ep)))/(2*ep))
J[1,4] = 0.5*k1 + 0.5*k2n + (0.5*k1 - 0.5*k2n)*tanh((eta + xn - yn)/ep) + phi*(xl - xn - yl + yn)*tanh((phi*(xl - xn - yl + yn))/ep)
# SECOND ROW
J[2,1] = (1/ep)*((0.5*k1*yb - 0.5*k2n*yb - 0.5*k1*yn + 0.5*k2n*yn)*sech((eta + xn - yn)/ep)^2 + phi*sech((phi*(xl - xn - yl
+ yn))/ep)^2*(phi*(-xl*yb + xn*yb - xl*yl + xn*yl + yb*yl + yl^2 + 2*xl*yn - 2*xn*yn - yb*yn - 3*yl*yn + 2*yn^2) + ep*(-0.5*yb
- 0.5*yl + yn)*sinh((2*phi*(xl - xn - yl + yn))/ep)))
J[2,2] = (1/ep)*((-0.5*k1*yb + 0.5*k2n*yb + 0.5*k1*yn - 0.5*k2n*yn)*sech((eta + xn - yn)/ep)^2 +
phi^2*(xl*yb - xn*yb + xl*yl - xn*yl - yb*yl - yl^2 - 2*xl*yn + 2*xn*yn + yb*yn + 3*yl*yn -
2*yn^2)*sech((phi*(xl - xn - yl + yn))/ep)^2 + ep*(-K - 0.5*k1 - 0.5*k2n + (-0.5*k1 + 0.5*k2n)*tanh((eta + xn - yn)/ep) +
phi*(-2*xl + 2*xn + yb + 3*yl - 4*yn)*tanh((phi*(xl - xn - yl + yn))/ep)))
J[2,3] = K + (1/(2*ep))*phi*sech((phi*(xl - xn - yl + yn))/ep)^2*(2*phi*(yb + yl - 2*yn)*(-xl + xn + yl - yn) +
ep*(-yb + yn)*sinh((2*phi*(xl - xn - yl + yn))/ep)) + phi*(xl - xn - 2*yl + 2*yn)*tanh((phi*(xl - xn - yl + yn))/ep)
J[2,4] = 0.0
# THIRD ROW
J[3,1] = -((phi*(yb - 2*yl + yn)*sech((phi*(xl - xn - yl + yn))/ep)^2*(2*phi*(xl - xn - yl + yn) + ep*sinh((2*phi*(xl - xn - yl + yn))/ep)))/(2*ep*TildeVL))
J[3,2] = (1/TildeVL)*(K + (1/(2*ep))*phi*sech((phi*(xl - xn - yl + yn))/ep)^2*(2*phi*(yb - 2*yl + yn)*(xl - xn - yl + yn)
+ ep*(yb - yl)*sinh((2*phi*(xl - xn - yl + yn))/ep)) + phi*(xl - xn - 2*yl + 2*yn)*tanh((phi*(xl - xn - yl + yn))/ep))
J[3,3] = (1/(ep*TildeVL))*((-0.5*k1*yb + 0.5*k2l*yb + 0.5*k1*yl - 0.5*k2l*yl)*sech((eta + xl - yl)/ep)^2 + phi^2*(-xl*yb + xn*yb + 2*xl*yl - 2*xn*yl + yb*yl -
2*yl^2 - xl*yn + xn*yn - yb*yn + 3*yl*yn - yn^2)*sech((phi*(xl - xn - yl + yn))/ep)^2 + ep*(-K - 0.5*k1 - 0.5*k2l + (-0.5*k1
+ 0.5*k2l)*tanh((eta + xl - yl)/ep) + phi*(-2*xl + 2*xn - yb + 4*yl - 3*yn)*tanh((phi*(xl - xn - yl + yn))/ep)))
J[3,4] = 0.0
# FOURTH ROW
J[4,1] = (1/(ep*TildeVB))*((-0.5*k1*xb + 0.5*k2n*xb + 0.5*k1*xn - 0.5*k2n*xn)*sech((eta + xn - yn)/ep)^2 + phi^2*(-xl^2 + xl*(yl - yn)
+ xn*(xn + yl - yn) + xb*(2*xl - 2*xn - 2*yl + 2*yn))*sech((phi*(xl - xn - yl + yn))/ep)^2 + ep*(0.5*k1 + 0.5*k2n + (0.5*k1
- 0.5*k2n)*tanh((eta + xn - yn)/ep) + phi*(2*xb - 2*xn - yl + yn)*tanh((phi*(xl - xn - yl + yn))/ep)))
J[4,2] = (1/(ep*TildeVB))*((0.5*k1*xb - 0.5*k2n*xb - 0.5*k1*xn + 0.5*k2n*xn)*sech((eta + xn - yn)/ep)^2
+ phi*sech((phi*(xl - xn - yl + yn))/ep)^2*(phi*(xl^2 - xn^2 - xl*yl - xn*yl + xb*(-2*xl + 2*xn
+ 2*yl - 2*yn) + xl*yn + xn*yn) + ep*(-xb + 0.5*xl + 0.5*xn)*sinh((2*phi*(xl - xn - yl + yn))/ep)))
J[4,3] = (1/(ep*TildeVB))*((0.5*k1*xb - 0.5*k2l*xb - 0.5*k1*xl + 0.5*k2l*xl)*sech((eta + xl - yl)/ep)^2 + phi*sech((phi*(xl - xn - yl
+ yn))/ep)^2*(phi*(-xl^2 + xn^2 + xl*yl + xn*yl - xl*yn - xn*yn + xb*(2*xl - 2*xn - 2*yl + 2*yn)) + ep*(xb
- 0.5*xl - 0.5*xn)*sinh((2*phi*(xl - xn - yl + yn))/ep)))
J[4,4] = (1/TildeVB)*(-k1 - 0.5*k2l - 0.5*k2n + (-0.5*k1 + 0.5*k2l)*tanh((eta + xl - yl)/ep)
+ (-0.5*k1 + 0.5*k2n)*tanh((eta + xn - yn)/ep) + phi*(-2*xl + 2*xn + 2*yl - 2*yn)*tanh((phi*(xl - xn - yl + yn))/ep))
# println(J)
return J
end
function ODEHandle2!(du, u, p)
m = p[1]
eta = p[2]
k2 = p[3]
xn, yn, yl, xb = u
### Set main constants
SV=10^6;
sig=2.1*10^4*SV;
alphaT=1.7*10^(-4);
alphaS=0.8*10^(-3);
ep=0.011;
TaN=7;
TaL=25
T0=2.65
S0 =35;
SB =34.538;
yb = alphaS/alphaT * (SB - S0)/(TaN - T0);
g=3.17*10^-8;
VN=7.2106*10^15 #7*10^16;
VL=6.3515 * 10^16
VB=25 * VN
phi=alphaT*(TaN-T0);
xl=(TaL-T0)/(TaN-T0);
#
TildeVL=VL/VN;
TildeVB=VB/VN;
K=5.456*SV/sig;
ps=phi*((yn-xn)-(yl-xl));
d = g*VN/sig
# Functions
k2n = k2;
k2l = k2/3
k1 = 10^-2*SV/sig;
kappaN=k1+0.5*(k2n-k1)*(1+tanh((yn-xn-eta)/ep));
kappaL=k1+0.5*(k2l-k1)*(1+tanh((yl-xl-eta)/ep));
PsiPlus=1/2*(1+tanh(ps/ep))*ps ;
PsiMinus=1/2*(1-tanh(ps/ep))*ps ;
Psi=ps*tanh(ps/ep) ;
# Set up function output
du[1] = -d*(xn-1)+(K+Psi)*(xl-xn)+(PsiMinus-PsiPlus-kappaN)*(xn-xb)
du[2] = m+(K+Psi)*(yl-yn)+(PsiMinus-PsiPlus-kappaN)*(yn-yb)
du[3] = 1/TildeVL*(-m-(K+Psi)*(yl-yn)+(PsiMinus-PsiPlus-kappaL)*(yl-yb))
du[4] = 1/TildeVB*((kappaN-PsiMinus+PsiPlus)*(xn-xb)+(kappaL-PsiMinus+PsiPlus)*(xl-xb))
du
end
# options
opts = ContinuationPar(
p_min = -3.0,
p_max = 0.01,
detect_bifurcation = 3,
nev = 4,
max_steps = 50000,
dsmin = 1e-15,
ds = 1e-3,
dsmax = 2e-3,
n_inversion = 8,
max_bisection_steps =150,
detect_event = 0
)
# Equilibria
z0 = [1.242167, -1.260954, 0.2613593, 3.190049]
# Define the bifurcation problem
prob = BifurcationProblem(
ODEHandle2!,
z0,
par_tm,
(@lens _.mu),
J=JAC)
# Codimension-one continuation
br = continuation(
prob,
PALC(),
J=JAC,
opts;
normC = norminf,
bothside=true,
detect_bifurcation=3)
# Codimension-two continuation from the branch point at index 2
indFold = 2
sn_codim2 = continuation(br, indFold, (@lens _.eta), J=JAC,
opts;
bothside=true,
normC = norminf,
bdlinsolver = MatrixBLS(),
update_minaug_every_step = 1,
detect_codim2_bifurcation = 2,
detect_loop=1
)
```
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 16, 2024, 2:48pm UTC](https://discourse.julialang.org/t/marking-user-points-events-along-a-cxdimension-one-curve-of-fold-bifurcation/114325/4 "2024-05-16T14:48:33Z")
Hi,
I took the freedom to move this to n issue to BifurcationKit, see here:
>
>
> Hello,
>
> I am using version “BifurcationKit v0.3.3” according to Julia. I agree… my code does not make much sense, I am new to BifurcationKit.jl and was/am fumbling a bit.
>
> I have removed the incorrect “jacobian\_ma " and “jacobian” arguments that you pointed out. However, I now run into an additional problem – whenever I do not set " jacobian\_ma” to anything, even when incorrect, I get the following error:
>
> \`\`\`julia
> ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{BifurcationKit.var"#352#353"{BifurcationKit.FoldMAProblem{…}, @NamedTuple{…}}, Float64}, Float64, 5})
>
> Closest candidates are:
> (::Type{T})(::Real, ::RoundingMode) where T\<:AbstractFloat
> @ Base rounding.jl:207
> (::Type{T})(::T) where T\<:Number
> @ Core boot.jl:792
> Float64(::IrrationalConstants.Fourπ)
> @ IrrationalConstants ~/.julia/packages/IrrationalConstants/vp5v4/src/macro.jl:112
> ...
>
> Stacktrace:
> \[1\] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{BifurcationKit.var"#352#353"{…}, Float64}, Float64, 5})
> @ Base ./number.jl:7
> \[2\] setindex!
> @ ./array.jl:1024 \[inlined\]
> \[3\] JAC(u::SubArray{ForwardDiff.Dual{…}, 1, Vector{…}, Tuple{…}, true}, p::@NamedTuple{mu::ForwardDiff.Dual{…}, eta::Float64, k2::Float64})
> @ Main ~/Library/CloudStorage/OneDrive-TheUniversityofAuckland/Git/psweep\_forced\_julia/BailieKeane\_2024\_nonauto/Test:40
> \[4\] jacobian(pb::BifFunction{…}, x::SubArray{…}, p::@NamedTuple{…})
> @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Problems.jl:72
> \[5\] jacobian(pb::BifurcationProblem{…}, x::SubArray{…}, p::@NamedTuple{…})
> @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Problems.jl:234
> \[6\] (::FoldProblemMinimallyAugmented{…})(x::SubArray{…}, p::ForwardDiff.Dual{…}, params::@NamedTuple{…})
> @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:35
> \[7\] (::FoldProblemMinimallyAugmented{…})(x::Vector{…}, params::@NamedTuple{…})
> @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:48
> \[8\] (::BifurcationKit.var"#352#353"{BifurcationKit.FoldMAProblem{…}, @NamedTuple{…}})(z::Vector{ForwardDiff.Dual{…}})
> @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:172
> \[9\] vector\_mode\_dual\_eval!(f::BifurcationKit.var"#352#353"{…}, cfg::ForwardDiff.JacobianConfig{…}, x::Vector{…})
> @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/apiutils.jl:24
> \[10\] vector\_mode\_jacobian(f::BifurcationKit.var"#352#353"{…}, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…})
> @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:125
> \[11\] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{…}, Float64, 5, Vector{…}}, ::Val{true})
> @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:0
> \[12\] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{…}, Float64, 5, Vector{…}})
> @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:19
> \[13\] jacobian
> @ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:172 \[inlined\]
> \[14\] \_newton(prob::BifurcationKit.FoldMAProblem{…}, x0::Vector{…}, p0::@NamedTuple{…}, options::NewtonPar{…}; normN::typeof(norminf), callback::typeof(BifurcationKit.cb\_default), kwargs::@Kwargs{…})
> @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Newton.jl:88
> \[15\] \_newton
> @ ~/.julia/packages/BifurcationKit/5kgAI/src/Newton.jl:62 \[inlined\]
> \[16\] newton
> @ ~/.julia/packages/BifurcationKit/5kgAI/src/Newton.jl:138 \[inlined\]
> \[17\] iterate(it::ContIterable{…}; \_verbosity::Int64)
> @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:271
> \[18\] iterate
> @ ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:250 \[inlined\]
> \[19\] continuation(it::ContIterable{…})
> @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:474
> \[20\] continuation(prob::BifurcationKit.FoldMAProblem{…}, alg::PALC{…}, contparams::ContinuationPar{…}; linear\_algo::BorderingBLS{…}, bothside::Bool, kwargs::@Kwargs{…})
> @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/Continuation.jl:539
> \[21\] continuation\_fold(prob::BifurcationProblem{…}, alg::PALC{…}, foldpointguess::BorderedArray{…}, par::@NamedTuple{…}, lens1::Setfield.PropertyLens{…}, lens2::Setfield.PropertyLens{…}, eigenvec::Vector{…}, eigenvec\_ad::Vector{…}, options\_cont::ContinuationPar{…}; update\_minaug\_every\_step::Int64, normC::typeof(norminf), bdlinsolver::MatrixBLS{…}, bdlinsolver\_adjoint::MatrixBLS{…}, jacobian\_ma::Symbol, compute\_eigen\_elements::Bool, usehessian::Bool, kind::BifurcationKit.FoldCont, record\_from\_solution::Nothing, kwargs::@Kwargs{…})
> @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:497
> \[22\] continuation\_fold
> @ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:314 \[inlined\]
> \[23\] continuation\_fold(prob::BifurcationProblem{…}, br::ContResult{…}, ind\_fold::Int64, lens2::Setfield.PropertyLens{…}, options\_cont::ContinuationPar{…}; alg::PALC{…}, normC::typeof(norminf), nev::Int64, start\_with\_eigen::Bool, bdlinsolver::MatrixBLS{…}, bdlinsolver\_adjoint::MatrixBLS{…}, kwargs::@Kwargs{…})
> @ BifurcationKit ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:551
> \[24\] continuation\_fold
> @ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/MinAugFold.jl:511 \[inlined\]
> \[25\] #continuation#348
> @ ~/.julia/packages/BifurcationKit/5kgAI/src/codim2/codim2.jl:237 \[inlined\]
> \[26\] top-level scope
> @ ~/Library/CloudStorage/OneDrive-TheUniversityofAuckland/Git/psweep\_forced\_julia/BK\_2024\_nonautoBailieKeane\_2024\_nonauto/Test:195
> \[27\] include(fname::String)
> @ Base.MainInclude ./client.jl:489
> \[28\] top-level scope
> @ REPL\[49\]:1
> in expression starting at /Users/johnbailie/Library/CloudStorage/OneDrive-TheUniversityofAuckland/Git/psweep\_forced\_julia/BK\_2024\_nonautoBailieKeane\_2024\_nonauto/Test:195
> Some type information was truncated. Use \`show(err)\` to see complete types.
> \`\`\`
>
> On the other hand, when I set the “jacobian\_ma = :a” for example, the correct branch of fold bifurcation is computed. I am honestly really lost on this, I believed I followed the documentation correctly?
>
> Any assistance would be really appreciated.
>
> FYI, this is my code
>
> \`\`\`julia
> function JAC(u, p)
> m = p\[1\]
> eta = p\[2\]
> k2 = p\[3\]
> xn, yn, yl, xb = u
>
> ### Set main constants
> SV=10^6;
> sig=2.1\*10^4\*SV;
> alphaT=1.7\*10^(-4);
> alphaS=0.8\*10^(-3);
> ep=0.011;
> TaN=7;
> TaL=25
> T0=2.65
> S0 =35;
> SB =34.538;
> yb = alphaS/alphaT \* (SB - S0)/(TaN - T0);
> g=3.17\*10^-8;
> VN=7.2106\*10^15 #7\*10^16;
> VL=6.3515 \* 10^16
> VB=25 \* VN
> phi=alphaT\*(TaN-T0);
> xl=(TaL-T0)/(TaN-T0);
> #
> TildeVL=VL/VN;
> TildeVB=VB/VN;
> K=5.456\*SV/sig;
> ps=phi\*((yn-xn)-(yl-xl));
> d = g\*VN/sig
>
> # Functions
> k2n = k2;
> k2l = k2/3
> k1 = 10^-2\*SV/sig;
>
> J = zeros(4,4)
>
> # FIRST ROW
> J\[1,1\] = (1/ep)\*((0.5\*k1\*xb - 0.5\*k2n\*xb - 0.5\*k1\*xn + 0.5\*k2n\*xn)\*sech((eta + xn - yn)/ep)^2 + phi^2\*(-xb\*xl - xl^2 + xb\*xn + 3\*xl\*xn - 2\*xn^2 +
> xb\*yl + xl\*yl - 2\*xn\*yl - xb\*yn - xl\*yn + 2\*xn\*yn)\*sech((phi\*(xl - xn - yl + yn))/ep)^2 + ep\*(-d - K - 0.5\*k1 -
> 0.5\*k2n + (-0.5\*k1 + 0.5\*k2n)\*tanh((eta + xn - yn)/ep) + phi\*(-xb - 3\*xl + 4\*xn + 2\*yl - 2\*yn)\*tanh((phi\*(xl - xn - yl + yn))/ep)))
>
> J\[1,2\] = (1/ep)\*((-0.5\*k1\*xb + 0.5\*k2n\*xb + 0.5\*k1\*xn - 0.5\*k2n\*xn)\*sech((eta + xn - yn)/ep)^2 + phi\*sech((phi\*(xl - xn - yl + yn))/ep)^2\*(phi\*(xl^2
> + xn\*(2\*xn + 2\*yl - 2\*yn) + xl\*(-3\*xn - yl + yn) + xb\*(xl - xn - yl + yn)) + ep\*(0.5\*xb + 0.5\*xl - xn)\*sinh((2\*phi\*(xl - xn - yl + yn))/ep)))
>
> J\[1,3\] = -((phi\*(xb + xl - 2\*xn)\*sech((phi\*(xl - xn - yl + yn))/ep)^2\*(2\*phi\*(xl - xn - yl + yn) + ep\*sinh((2\*phi\*(xl - xn - yl + yn))/ep)))/(2\*ep))
>
> J\[1,4\] = 0.5\*k1 + 0.5\*k2n + (0.5\*k1 - 0.5\*k2n)\*tanh((eta + xn - yn)/ep) + phi\*(xl - xn - yl + yn)\*tanh((phi\*(xl - xn - yl + yn))/ep)
>
>
> # SECOND ROW
> J\[2,1\] = (1/ep)\*((0.5\*k1\*yb - 0.5\*k2n\*yb - 0.5\*k1\*yn + 0.5\*k2n\*yn)\*sech((eta + xn - yn)/ep)^2 + phi\*sech((phi\*(xl - xn - yl
> + yn))/ep)^2\*(phi\*(-xl\*yb + xn\*yb - xl\*yl + xn\*yl + yb\*yl + yl^2 + 2\*xl\*yn - 2\*xn\*yn - yb\*yn - 3\*yl\*yn + 2\*yn^2) + ep\*(-0.5\*yb
> - 0.5\*yl + yn)\*sinh((2\*phi\*(xl - xn - yl + yn))/ep)))
>
> J\[2,2\] = (1/ep)\*((-0.5\*k1\*yb + 0.5\*k2n\*yb + 0.5\*k1\*yn - 0.5\*k2n\*yn)\*sech((eta + xn - yn)/ep)^2 +
> phi^2\*(xl\*yb - xn\*yb + xl\*yl - xn\*yl - yb\*yl - yl^2 - 2\*xl\*yn + 2\*xn\*yn + yb\*yn + 3\*yl\*yn -
> 2\*yn^2)\*sech((phi\*(xl - xn - yl + yn))/ep)^2 + ep\*(-K - 0.5\*k1 - 0.5\*k2n + (-0.5\*k1 + 0.5\*k2n)\*tanh((eta + xn - yn)/ep) +
> phi\*(-2\*xl + 2\*xn + yb + 3\*yl - 4\*yn)\*tanh((phi\*(xl - xn - yl + yn))/ep)))
>
>
> J\[2,3\] = K + (1/(2\*ep))\*phi\*sech((phi\*(xl - xn - yl + yn))/ep)^2\*(2\*phi\*(yb + yl - 2\*yn)\*(-xl + xn + yl - yn) +
> ep\*(-yb + yn)\*sinh((2\*phi\*(xl - xn - yl + yn))/ep)) + phi\*(xl - xn - 2\*yl + 2\*yn)\*tanh((phi\*(xl - xn - yl + yn))/ep)
>
> J\[2,4\] = 0.0
>
>
> # THIRD ROW
> J\[3,1\] = -((phi\*(yb - 2\*yl + yn)\*sech((phi\*(xl - xn - yl + yn))/ep)^2\*(2\*phi\*(xl - xn - yl + yn) + ep\*sinh((2\*phi\*(xl - xn - yl + yn))/ep)))/(2\*ep\*TildeVL))
>
> J\[3,2\] = (1/TildeVL)\*(K + (1/(2\*ep))\*phi\*sech((phi\*(xl - xn - yl + yn))/ep)^2\*(2\*phi\*(yb - 2\*yl + yn)\*(xl - xn - yl + yn)
> + ep\*(yb - yl)\*sinh((2\*phi\*(xl - xn - yl + yn))/ep)) + phi\*(xl - xn - 2\*yl + 2\*yn)\*tanh((phi\*(xl - xn - yl + yn))/ep))
>
> J\[3,3\] = (1/(ep\*TildeVL))\*((-0.5\*k1\*yb + 0.5\*k2l\*yb + 0.5\*k1\*yl - 0.5\*k2l\*yl)\*sech((eta + xl - yl)/ep)^2 + phi^2\*(-xl\*yb + xn\*yb + 2\*xl\*yl - 2\*xn\*yl + yb\*yl -
> 2\*yl^2 - xl\*yn + xn\*yn - yb\*yn + 3\*yl\*yn - yn^2)\*sech((phi\*(xl - xn - yl + yn))/ep)^2 + ep\*(-K - 0.5\*k1 - 0.5\*k2l + (-0.5\*k1
> + 0.5\*k2l)\*tanh((eta + xl - yl)/ep) + phi\*(-2\*xl + 2\*xn - yb + 4\*yl - 3\*yn)\*tanh((phi\*(xl - xn - yl + yn))/ep)))
>
> J\[3,4\] = 0.0
>
>
> # FOURTH ROW
> J\[4,1\] = (1/(ep\*TildeVB))\*((-0.5\*k1\*xb + 0.5\*k2n\*xb + 0.5\*k1\*xn - 0.5\*k2n\*xn)\*sech((eta + xn - yn)/ep)^2 + phi^2\*(-xl^2 + xl\*(yl - yn)
> + xn\*(xn + yl - yn) + xb\*(2\*xl - 2\*xn - 2\*yl + 2\*yn))\*sech((phi\*(xl - xn - yl + yn))/ep)^2 + ep\*(0.5\*k1 + 0.5\*k2n + (0.5\*k1
> - 0.5\*k2n)\*tanh((eta + xn - yn)/ep) + phi\*(2\*xb - 2\*xn - yl + yn)\*tanh((phi\*(xl - xn - yl + yn))/ep)))
>
> J\[4,2\] = (1/(ep\*TildeVB))\*((0.5\*k1\*xb - 0.5\*k2n\*xb - 0.5\*k1\*xn + 0.5\*k2n\*xn)\*sech((eta + xn - yn)/ep)^2
> + phi\*sech((phi\*(xl - xn - yl + yn))/ep)^2\*(phi\*(xl^2 - xn^2 - xl\*yl - xn\*yl + xb\*(-2\*xl + 2\*xn
> + 2\*yl - 2\*yn) + xl\*yn + xn\*yn) + ep\*(-xb + 0.5\*xl + 0.5\*xn)\*sinh((2\*phi\*(xl - xn - yl + yn))/ep)))
>
> J\[4,3\] = (1/(ep\*TildeVB))\*((0.5\*k1\*xb - 0.5\*k2l\*xb - 0.5\*k1\*xl + 0.5\*k2l\*xl)\*sech((eta + xl - yl)/ep)^2 + phi\*sech((phi\*(xl - xn - yl
> + yn))/ep)^2\*(phi\*(-xl^2 + xn^2 + xl\*yl + xn\*yl - xl\*yn - xn\*yn + xb\*(2\*xl - 2\*xn - 2\*yl + 2\*yn)) + ep\*(xb
> - 0.5\*xl - 0.5\*xn)\*sinh((2\*phi\*(xl - xn - yl + yn))/ep)))
>
> J\[4,4\] = (1/TildeVB)\*(-k1 - 0.5\*k2l - 0.5\*k2n + (-0.5\*k1 + 0.5\*k2l)\*tanh((eta + xl - yl)/ep)
> + (-0.5\*k1 + 0.5\*k2n)\*tanh((eta + xn - yn)/ep) + phi\*(-2\*xl + 2\*xn + 2\*yl - 2\*yn)\*tanh((phi\*(xl - xn - yl + yn))/ep))
>
> # println(J)
> return J
> end
>
>
> function ODEHandle2!(du, u, p)
> m = p\[1\]
> eta = p\[2\]
> k2 = p\[3\]
> xn, yn, yl, xb = u
> ### Set main constants
> SV=10^6;
> sig=2.1\*10^4\*SV;
> alphaT=1.7\*10^(-4);
> alphaS=0.8\*10^(-3);
> ep=0.011;
> TaN=7;
> TaL=25
> T0=2.65
> S0 =35;
> SB =34.538;
> yb = alphaS/alphaT \* (SB - S0)/(TaN - T0);
> g=3.17\*10^-8;
> VN=7.2106\*10^15 #7\*10^16;
> VL=6.3515 \* 10^16
> VB=25 \* VN
> phi=alphaT\*(TaN-T0);
> xl=(TaL-T0)/(TaN-T0);
> #
> TildeVL=VL/VN;
> TildeVB=VB/VN;
> K=5.456\*SV/sig;
> ps=phi\*((yn-xn)-(yl-xl));
> d = g\*VN/sig
>
> # Functions
> k2n = k2;
> k2l = k2/3
> k1 = 10^-2\*SV/sig;
>
> kappaN=k1+0.5\*(k2n-k1)\*(1+tanh((yn-xn-eta)/ep));
> kappaL=k1+0.5\*(k2l-k1)\*(1+tanh((yl-xl-eta)/ep));
> PsiPlus=1/2\*(1+tanh(ps/ep))\*ps ;
> PsiMinus=1/2\*(1-tanh(ps/ep))\*ps ;
> Psi=ps\*tanh(ps/ep) ;
> # Set up function output
> du\[1\] = -d\*(xn-1)+(K+Psi)\*(xl-xn)+(PsiMinus-PsiPlus-kappaN)\*(xn-xb)
> du\[2\] = m+(K+Psi)\*(yl-yn)+(PsiMinus-PsiPlus-kappaN)\*(yn-yb)
> du\[3\] = 1/TildeVL\*(-m-(K+Psi)\*(yl-yn)+(PsiMinus-PsiPlus-kappaL)\*(yl-yb))
> du\[4\] = 1/TildeVB\*((kappaN-PsiMinus+PsiPlus)\*(xn-xb)+(kappaL-PsiMinus+PsiPlus)\*(xl-xb))
>
> du
> end
>
> \# options
> opts = ContinuationPar(
> p\_min = -3.0,
> p\_max = 0.01,
> detect\_bifurcation = 3,
> nev = 4,
> max\_steps = 50000,
> dsmin = 1e-15,
> ds = 1e-3,
> dsmax = 2e-3,
> n\_inversion = 8,
> max\_bisection\_steps =150,
> detect\_event = 0
> )
>
> \# Equilibria
> z0 = \[1.242167, -1.260954, 0.2613593, 3.190049\]
>
> \# Define the bifurcation problem
> prob = BifurcationProblem(
> ODEHandle2!,
> z0,
> par\_tm,
> (@lens \_.mu),
> J=JAC)
>
> \# Codimension-one continuation
> br = continuation(
> prob,
> PALC(),
> J=JAC,
> opts;
> normC = norminf,
> bothside=true,
> detect\_bifurcation=3)
>
>
> \# Codimension-two continuation from the branch point at index 2
> indFold = 2
> sn\_codim2 = continuation(br, indFold, (@lens \_.eta), J=JAC,
> opts;
> bothside=true,
> normC = norminf,
> bdlinsolver = MatrixBLS(),
> update\_minaug\_every\_step = 1,
> detect\_codim2\_bifurcation = 2,
> detect\_loop=1
> )
> \`\`\`
Can you answer my questions there please?