Marking user points / events along a cxdimension-one curve of fold bifurcation

BigEpsilon.jl · May 15, 2024, 9:54pm

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? >.<

rveltz · May 16, 2024, 6:49am

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?

BigEpsilon.jl · May 16, 2024, 10:28am

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:

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

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
)

rveltz · May 16, 2024, 2:48pm

Hi,

I took the freedom to move this to n issue to BifurcationKit, see here:

github.com/bifurcationkit/BifurcationKit.jl

Marking user points / events along a cxdimension-one curve of fold bifurcation

opened 02:45PM - 16 May 24 UTC

rveltz

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?