increase compat bound for SciMLBase

improve doc strings improve show method of bifurcation points
bifurcationkit · Oct 1, 2023 · 3070390 · 3070390
1 parent 8dcc4df
commit 3070390
Show file tree

Hide file tree

Showing 12 changed files with 44 additions and 79 deletions.
diff --git a/Project.toml b/Project.toml
@@ -41,10 +41,10 @@ LinearMaps = "2.7, ^3.0"
 OrdinaryDiffEq = "^6.33"
 Parameters = "0.12.3"
 RecipesBase = "1.0, 1.1"
-RecursiveArrayTools = "2.4, 2.8, ^2.9"
-Requires = "1.1.2"
+RecursiveArrayTools = "2.3, 2.4, 2.8, ^2.9"
 Reexport = "1.2.2"
-SciMLBase = "^1.42.1"
+Requires = "1.1.2"
+SciMLBase = "2.0.7"
 Setfield = "0.5.0, 0.7.0, 0.8.0, ^1.1"
 StructArrays = "0.4, ^0.6"
 julia = "1.5"

diff --git a/src/ContParameters.jl b/src/ContParameters.jl
@@ -45,10 +45,7 @@ Returns a variable containing parameters to affect the `continuation` algorithm
     # parameters for arclength continuation
     dsmin::T    = 1e-3
     dsmax::T    = 1e-1
-    @assert dsmax >= dsmin "You must provide a valid interval (ordered) for ds"
-    ds::T        = 1e-2; @assert dsmax >= abs(ds) >= dsmin
-    @assert dsmin > 0 "The interval for ds must be positive"
-    @assert dsmax > 0
+    ds::T        = 1e-2
 
     # parameters for continuation
     a::T    = 0.5 # aggressiveness factor
@@ -86,14 +83,15 @@ Returns a variable containing parameters to affect the `continuation` algorithm
     # handling event detection
     detect_event::Int64 = 0               # event location
     tol_param_bisection_event::T = 1e-16  # tolerance on value of parameter
+    detect_loop::Bool = false                # detect if the branch loops
 
+    @assert dsmax >= abs(ds) >= dsmin > 0 "You must provide a valid interval (ordered) for ds"
     @assert p_max >= p_min "You must provide a valid interval [p_min, p_max]"
     @assert iseven(n_inversion) "The option `n_inversion` number must be even"
     @assert 0 <= detect_bifurcation <= 3 "The option `detect_bifurcation` must belong to {0,1,2,3}"
     @assert 0 <= detect_event <= 2 "The option `detect_event` must belong to {0,1,2}"
     @assert (detect_bifurcation > 1 && detect_event == 0) || (detect_bifurcation <= 1 && detect_event >= 0)  "One of these options must be disabled detect_bifurcation = $detect_bifurcation and detect_event = $detect_event"
     @assert tol_bisection_eigenvalue >= 0 "The option `tol_bisection_eigenvalue` must be positive"
-    detect_loop::Bool = false                # detect if the branch loops
     @assert plot_every_step > 0 "plot_every_step must be positive. You can turn off plotting by passing plot = false to `continuation`"
     @assert ~(detect_bifurcation > 1 && save_eig_every_step > 1) "We must at least save all eigenvalues for detection of bifurcation points. Please use save_eig_every_step = 1 or detect_bifurcation = 1."
 end

diff --git a/src/Utils.jl b/src/Utils.jl
@@ -129,13 +129,12 @@ $(SIGNATURES)
 This function extracts the indices of the blocks composing the matrix A which is a M x M Block matrix where each block N x N has the same sparsity.
 """
 function get_blocks(A::SparseMatrixCSC, N, M)
-    I,J,K = findnz(A)
+    I, J, K = findnz(A)
     out = [Vector{Int}() for i in 1:M+1, j in 1:M+1];
     for k in eachindex(I)
         m, l = div(I[k]-1, N), div(J[k]-1, N)
         push!(out[1+m, 1+l], k)
     end
-    res = [length(m) for m in out]
     out
 end
 ####################################################################################################

diff --git a/src/continuation/Palc.jl b/src/continuation/Palc.jl
@@ -434,7 +434,7 @@ function newton_palc(iter::AbstractContinuationIterable,
     while (step < max_iterations) && (res > tol) && line_step && compute
         # dFdp = (F(x, p + ϵ) - F(x, p)) / ϵ)
         copyto!(dFdp, residual(prob, x, set(par, paramlens, p + ϵ)))
-            minus!(dFdp, res_f); rmul!(dFdp, one(T) / ϵ)
+        minus!(dFdp, res_f); rmul!(dFdp, one(T) / ϵ)
 
         # compute jacobian
         J = jacobian(prob, x, set(par, paramlens, p))

diff --git a/src/periodicorbit/BifurcationPoints.jl b/src/periodicorbit/BifurcationPoints.jl
@@ -53,7 +53,8 @@ function Base.show(io::IO, pd::PeriodDoublingPO)
     else
         if ~pd.prm
             println("├─ ", get_lens_symbol(pd.nf.lens)," ≈ $(pd.nf.p).")
-            println(io, "├─ Normal form:\n├\t∂τ = 1 + a⋅ξ²\n├\t∂ξ = c⋅ξ³\n└─ ", pd.nf.nf)
+            println(io, "├─ Normal form:\n├\t∂τ = 1 + a⋅ξ²\n├\t∂ξ = c⋅ξ³")
+            println(io,"├─── a = ", pd.nf.nf.a, "\n└─── c = ", pd.nf.nf.b3)
         else
             show(io, pd.nf)
         end
@@ -115,7 +116,7 @@ type(bp::NeimarkSackerPO) = type(bp.nf)
 function Base.show(io::IO, ns::NeimarkSackerPO)
     printstyled(io, ns.nf.type, " - ",type(ns), color=:cyan, bold = true)
     println(io, " bifurcation point of periodic orbit\n┌─ ", get_lens_symbol(ns.nf.lens)," ≈ $(ns.p).")
-    println(io, "├─ Frequency θ ≈ ", abs(ns.ω))
+    println(io, "├─ Frequency θ ≈ ", ns.ω)
     println(io, "├─ Period at the periodic orbit T ≈ ", abs(ns.T))
     println(io, "├─ Second frequency of the bifurcated torus ≈ ", abs(2pi/ns.ω))
     if ns.prm

diff --git a/src/periodicorbit/NormalForms.jl b/src/periodicorbit/NormalForms.jl
@@ -295,11 +295,13 @@ function period_doubling_normal_form(pbwrap::WrapPOColl,
     B(u, p, du1, du2)      = d2F(coll.prob_vf, u, p, du1, du2)
     C(u, p, du1, du2, du3) = d3F(coll.prob_vf, u, p, du1, du2, du3)
 
-    _rand(n, r = 2) = r .* (rand(n) .- 1/2) # centered uniform random variables
+    _rand(n, r = 2) = r .* (rand(n) .- 1/2)        # centered uniform random variables
     local ∫(u,v) = BifurcationKit.∫(coll, u, v, 1) # define integral with coll parameters
 
     # we first compute the PD floquet eigenvector (for μ = -1)
     # we use an extended linear system for this
+    #########
+    # compute v1
     jac = jacobian(pbwrap, pd.x0, par)
     J = copy(jac.jacpb)
     nj = size(J, 1)
@@ -310,10 +312,6 @@ function period_doubling_normal_form(pbwrap::WrapPOColl,
     J[end-N:end-1, 1:N] .= I(N)
     J[end-N:end-1, end-N:end-1] .= I(N)
 
-    # the transpose problem is a bit different. 
-    # transposing J means the boundary condition is wrong
-    # however it seems Prop A.1 says the opposite
-
     rhs = zeros(nj); rhs[end] = 1;
     k = J  \ rhs; k = k[1:end-1]; k ./= norm(k) #≈ ker(J)
     l = J' \ rhs; l = l[1:end-1]; l ./= norm(l)
@@ -361,7 +359,7 @@ function period_doubling_normal_form(pbwrap::WrapPOColl,
     # for this, we generate the linear problem analytically
     # note that we could obtain the same by modifying inplace 
     # the previous linear problem Jψ
-    Jψ = analytical_jacobian(coll, pd.x0, par; _transpose = true, ρF = -1, ρI = 0)
+    Jψ = analytical_jacobian(coll, pd.x0, par; _transpose = true, ρF = -1)
     Jψ[end-N:end-1, 1:N] .= -I(N)
     Jψ[end-N:end-1, end-N:end-1] .= I(N)
     # build the extended linear problem
@@ -410,7 +408,8 @@ function period_doubling_normal_form(pbwrap::WrapPOColl,
         @warn "The value h₂[end] should be zero. We found $(h₂[end])"
     end
 
-    # computation of c. We need B(t, v₁(t), h₂(t))
+    # computation of c. 
+    # we need B(t, v₁(t), h₂(t))
     for i=1:size(Bₛ, 2)
         Bₛ[:,i]  .= B(u₀ₛ[:,i], par, v₁ₛ[:,i], h₂ₛ[:,i])
     end
@@ -420,9 +419,10 @@ function period_doubling_normal_form(pbwrap::WrapPOColl,
     c = 1/(3T) * ∫( v₁★ₛ, Cₛ ) + 
                  ∫( v₁★ₛ, Bₛ ) -
          2a₁/T * ∫( v₁★ₛ, Aₛ )
+                    @debug "" ∫( v₁★ₛ, Bₛ ) 2a₁/T * ∫( v₁★ₛ, Aₛ )
 
-    nf = (a = a₁, b3 = c) # keep b3 for PD-codim 2
-    @debug "" a₁ c c/a₁
+    nf = (a = a₁, b3 = c, h₂ₛ, ψ₁★ₛ, v₁ₛ) # keep b3 for PD-codim 2
+    @debug "" a₁ c
     return PeriodDoublingPO(pd.x0, T, v₁, v₁★, (@set pd.nf = nf), coll, false)
 end
 
@@ -525,7 +525,7 @@ function neimark_sacker_normal_form(pbwrap::WrapPOColl,
     # get the eigenvalue
     eigRes = br.eig
     λₙₛ = eigRes[bifpt.idx].eigenvals[bifpt.ind_ev]
-    ωₙₛ = imag(λₙₛ)
+    ωₙₛ = abs(imag(λₙₛ))
 
     if bifpt.x isa NamedTuple
         # the solution is mesh adapted, we need to restore the mesh.
@@ -572,41 +572,6 @@ function neimark_sacker_normal_form(pbwrap,
     return NeimarkSackerPO(bifpt.x, period, bifpt.param, ωₙₛ, nothing, nothing, ns0, pbwrap, true)
 end
 
-function neimark_sacker_normal_form(pbwrap::WrapPOColl,
-                                    br,
-                                    ind_bif::Int;
-                                    verbose = false,
-                                    nev = length(eigenvalsfrombif(br, ind_bif)),
-                                    newton_options = br.contparams.newton_options,
-                                    prm = true,
-                                    detailed = true,
-                                    kwargs_nf...)
-     # first, get the bifurcation point parameters
-     verbose && println("━"^53*"\n──▶ Period-Doubling normal form computation")
-     bifpt = br.specialpoint[ind_bif]
-     bptype = bifpt.type
-     par = setparam(br, bifpt.param)
-     period = getperiod(pbwrap.prob, bifpt.x, par)
-
-     if bifpt.x isa NamedTuple
-        # the solution is mesh adapted, we need to restore the mesh.
-        pbwrap = deepcopy(pbwrap)
-        updateMesh!(pbwrap.prob, bifpt.x._mesh )
-        bifpt = @set bifpt.x = bifpt.x.sol
-    end
-    ns0 =  NeimarkSacker(bifpt.x, bifpt.param, ωₙₛ, pars, getlens(br), nothing, nothing, nothing, :none)
-    if ~detailed || ~prm
-        # method based on Iooss method
-        return neimark_sacker_normal_form(pbwrap, ns0; detailed, verbose, nev, kwargs_nf...)
-    end
-    if prm # method based on Poincare Return Map (PRM)
-        # newton parameter
-        return neimark_sacker_normal_form_prm(pbwrap, ns0, newton_options; verbose, nev, kwargs_nf...)
-    end
-    return nothing
-
-end
-
 function neimark_sacker_normal_form_prm(pbwrap::WrapPOColl,
                                     ns0::NeimarkSacker,
                                     optn::NewtonPar;
@@ -615,7 +580,7 @@ function neimark_sacker_normal_form_prm(pbwrap::WrapPOColl,
                                     verbose = false,
                                     lens = getlens(pbwrap),
                                     kwargs_nf...)
-    @debug "methode PRM"
+    @debug "method PRM"
     coll = pbwrap.prob
     N, m, Ntst = size(coll)
     pars = ns0.params
@@ -664,9 +629,9 @@ function neimark_sacker_normal_form(pbwrap::WrapPOColl,
                                         nev = 3,
                                         verbose = false,
                                         lens = getlens(pbwrap),
+                                        _NRMDEBUG = false, # normalise to compare to ApproxFun
                                         kwargs_nf...)
-    _NRM = false # normalise to compare to ApproxFun
-    @warn "method IOOSS, NRM = $_NRM"
+    @warn "method IOOSS, NRM = $_NRMDEBUG"
 
     # based on the article
     # Kuznetsov, Yu. A., W. Govaerts, E. J. Doedel, and A. Dhooge. “Numerical Periodic Normalization for Codim 1 Bifurcations of Limit Cycles.” SIAM Journal on Numerical Analysis 43, no. 4 (January 2005): 1407–35. https://doi.org/10.1137/040611306.
@@ -687,15 +652,15 @@ function neimark_sacker_normal_form(pbwrap::WrapPOColl,
     C(u, p, du1, du2, du3) = TrilinearMap((dx1, dx2, dx3) -> d3F(coll.prob_vf, u, p, dx1, dx2, dx3))(du1, du2, du3)
 
     _plot(x; k...) = (_sol = get_periodic_orbit(coll, x, 1);display(plot(_sol.t, _sol.u'; k...)))
-    _rand(n, r = 2) = r .* (rand(n) .- 1/2) # centered uniform random variables
+    _rand(n, r = 2) = r .* (rand(n) .- 1/2)        # centered uniform random variables
     local ∫(u,v) = BifurcationKit.∫(coll, u, v, 1) # define integral with coll parameters
 
     #########
     # compute v1
     # we first compute the NS floquet eigenvector
     # we use an extended linear system for this
      # J = D  -  T*A(t) + iθ/T
-    θ = ns.ω
+    θ = abs(ns.ω)
     J = analytical_jacobian(coll, ns.x0, par; ρI = Complex(0,-θ/T), 𝒯 = ComplexF64)
 
     nj = size(J, 1)
@@ -716,7 +681,7 @@ function neimark_sacker_normal_form(pbwrap::WrapPOColl,
     v₁ ./= sqrt(∫(vr, vr))
     v₁ₛ = get_time_slices(coll, vcat(v₁,1))
 
-                if _NRM;v₁ₛ .*= (-0.46220415773497325 + 0.2722705470750184im)/v₁ₛ[1,1];end
+                if _NRMDEBUG;v₁ₛ .*= (-0.46231553686750715 - 0.27213798536704986im)/v₁ₛ[1,1];end
                 # re-scale the eigenvector
                 v₁ₛ ./= sqrt(∫(v₁ₛ, v₁ₛ))
                 v₁ = vec(v₁ₛ)
@@ -783,7 +748,7 @@ function neimark_sacker_normal_form(pbwrap::WrapPOColl,
     v₁★  = @view vr[1:end-1]
     v₁★ₛ = get_time_slices(coll, vcat(v₁★,1))
     v₁★ₛ ./= conj(∫(v₁★ₛ, v₁ₛ))
-                if _NRM; v₁★ₛ .*= (1.0371208296352463 + 4.170902638152008im)/v₁★ₛ[1,1];end
+                if _NRMDEBUG; v₁★ₛ .*= (1.0371208296352463 + 4.170902638152008im)/v₁★ₛ[1,1];end
                 # re-scale the eigenvector
     v₁★ₛ ./= conj(∫(v₁★ₛ, v₁ₛ))
     v₁★ = vec(v₁★ₛ)
@@ -833,7 +798,7 @@ function neimark_sacker_normal_form(pbwrap::WrapPOColl,
     h₁₁ ./= 2Ntst # this seems necessary to have something comparable to ApproxFun
     h₁₁ₛ = get_time_slices(coll, h₁₁)
                 # _plot(real(vcat(vec(h₁₁ₛ),1)),label="h11")
-                @info abs(∫( ϕ₁★ₛ, h₁₁ₛ))
+                @debug "" abs(∫( ϕ₁★ₛ, h₁₁ₛ))
     if abs(∫( ϕ₁★ₛ, h₁₁ₛ)) > 1e-10
         @warn "The integral ∫(coll,ϕ₁★ₛ, h₁₁ₛ) should be zero. We found $(∫( ϕ₁★ₛ, h₁₁ₛ ))"
     end
@@ -851,21 +816,21 @@ function neimark_sacker_normal_form(pbwrap::WrapPOColl,
 
     d = (1/T) * ∫( v₁★ₛ, Cₛ ) + 2 * ∫( v₁★ₛ, Bₛ )
 
-                @debug "h20*v1b" d  (1/T) * ∫( v₁★ₛ, Cₛ )     ∫( v₁★ₛ, Bₛ )
+                @debug "B(h11, v1)" d  (1/(2T)) * ∫( v₁★ₛ, Cₛ )     2*∫( v₁★ₛ, Bₛ )
 
     for i = 1:size(u₀ₛ, 2)
         Bₛ[:,i] .= B(u₀ₛ[:,i], par, h₂₀ₛ[:,i], conj(v₁ₛ[:,i]))
         Aₛ[:,i] .= A(u₀ₛ[:,i], par, v₁ₛ[:,i])
     end
-                @debug "h20*v1b" d   ∫( v₁★ₛ, Bₛ )
+                @debug "B(h20, v1b)" d   ∫( v₁★ₛ, Bₛ )
     d +=  ∫( v₁★ₛ, Bₛ )
     d = d/2
                 @debug ""  -a₁/T * ∫( v₁★ₛ, Aₛ ) + im * θ * a₁/T^2   im * θ * a₁/T^2
     d += -a₁/T * ∫( v₁★ₛ, Aₛ ) + im * θ * a₁/T^2
 
 
 
-    nf = (a = a₁, d, h₁₁ₛ, ϕ₁★ₛ, v₁★ₛ, h₂₀ₛ, _NRM) # keep b3 for ns-codim 2
+    nf = (a = a₁, d, h₁₁ₛ, ϕ₁★ₛ, v₁★ₛ, h₂₀ₛ, _NRMDEBUG) # keep b3 for ns-codim 2
     return NeimarkSackerPO(ns.x0, T, ns.p, θ, v₁, v₁★, (@set ns.nf = nf), coll, false)
 end
 

diff --git a/src/periodicorbit/PeriodicOrbitCollocation.jl b/src/periodicorbit/PeriodicOrbitCollocation.jl
@@ -156,7 +156,7 @@ Here are some useful methods you can apply to `pb`
 
 - `length(pb)` gives the total number of unknowns
 - `size(pb)` returns the triplet `(N, m, Ntst)`
-- `getmesh(pb)` returns the mesh `0 = τ0 < ... < τNtst+1 = 1`. This is useful because this mesh is born to vary by automatic mesh adaptation
+- `getmesh(pb)` returns the mesh `0 = τ0 < ... < τNtst+1 = 1`. This is useful because this mesh is born to vary during automatic mesh adaptation
 - `get_mesh_coll(pb)` returns the (static) mesh `0 = σ0 < ... < σm+1 = 1`
 - `get_times(pb)` returns the vector of times (length `1 + m * Ntst`) at the which the collocation is applied.
 - `generate_solution(pb, orbit, period)` generate a guess from a function `t -> orbit(t)` which approximates the periodic orbit.
@@ -166,7 +166,7 @@ Here are some useful methods you can apply to `pb`
 You will see below that you can evaluate the residual of the functional (and other things) by calling `pb(orbitguess, p)` on an orbit guess `orbitguess`. Note that `orbitguess` must be of size 1 + N * (1 + m * Ntst) where N is the number of unknowns in the state space and `orbitguess[end]` is an estimate of the period ``T`` of the limit cycle.
 
 # Constructors
-- `PeriodicOrbitOCollProblem(Ntst::Int, m::Int; kwargs)` creates an empty functional with `Ntst`and `m`.
+- `PeriodicOrbitOCollProblem(Ntst::Int, m::Int; kwargs)` creates an empty functional with `Ntst` and `m`.
 
 Note that you can generate this guess from a function using `generate_solution`.
 
@@ -528,7 +528,8 @@ Compute the jacobian of the problem defining the periodic orbits by orthogonal c
 @views function analytical_jacobian!(J,
                                     coll::PeriodicOrbitOCollProblem,
                                     u::AbstractVector,
-                                    pars; _transpose::Bool = false,
+                                    pars; 
+                                    _transpose::Bool = false,
                                     ρD = 1,
                                     ρF = 1,
                                     ρI = 0)

diff --git a/src/periodicorbit/PeriodicOrbitUtils.jl b/src/periodicorbit/PeriodicOrbitUtils.jl
@@ -75,14 +75,15 @@ function modify_po_finalise(prob, kwargs, updateSectionEveryStep)
     return _finsol2
 end
 
-# version specific to collocation. Handles mesh adaptation
+# version specific to collocation. Handle mesh adaptation
 function modify_po_finalise(prob::PeriodicOrbitOCollProblem, kwargs, updateSectionEveryStep)
     _finsol = get(kwargs, :finalise_solution, nothing)
     _finsol2 = (z, tau, step, contResult; kF...) ->
         begin
             # we first check that the continuation step was successful
             # if not, we do not update the problem with bad information
-            success = converged(get(kF, :state, nothing))
+            state = get(kF, :state, nothing)
+            success = isnothing(state) ? false : converged(state)
             # mesh adaptation
             if success && prob.meshadapt
                 oldsol = _copy(z)

diff --git a/src/periodicorbit/codim2/PeriodicOrbitCollocation.jl b/src/periodicorbit/codim2/PeriodicOrbitCollocation.jl
@@ -40,6 +40,7 @@ function continuation(br::AbstractResult{Tkind, Tprob},
     biftype = br.specialpoint[ind_bif].type
 
     # options to detect codim2 bifurcations
+    compute_eigen_elements = options_cont.detect_bifurcation > 0
     _options_cont = detect_codim2_parameters(detect_codim2_bifurcation, options_cont; kwargs...)
 
     if biftype == :bp
@@ -88,7 +89,6 @@ function continuation_coll_fold(br::AbstractResult{Tkind, Tprob},
         start_with_eigen = start_with_eigen,
         bdlinsolver = FloquetWrapperBLS(bdlinsolver),
         kind = FoldPeriodicOrbitCont(),
-        # detect_codim2_bifurcation = detect_codim2_bifurcation, # not necessary
         kwargs...
         )
 end
@@ -105,7 +105,7 @@ function continuation_coll_pd(br::AbstractResult{Tkind, Tprob},
     bifpt = br.specialpoint[ind_bif]
     biftype = bifpt.type
 
-    @assert biftype == :pd "We continue only PD points of Periodic orbits for now"
+    @assert biftype == :pd "Please open an issue on BifurcationKit website"
 
     pdpointguess = pd_point(br, ind_bif)
 

diff --git a/src/periodicorbit/codim2/StandardShooting.jl b/src/periodicorbit/codim2/StandardShooting.jl
@@ -190,7 +190,6 @@ function continuation_sh_fold(br::AbstractResult{Tkind, Tprob},
         br, ind_bif, lens2,
         options_cont;
         start_with_eigen = start_with_eigen,
-        # detect_codim2_bifurcation = detect_codim2_bifurcation, # not necessary
         bdlinsolver = FloquetWrapperBLS(bdlinsolver),
         kind = FoldPeriodicOrbitCont(),
         kwargs...)

diff --git a/src/periodicorbit/codim2/codim2.jl b/src/periodicorbit/codim2/codim2.jl
@@ -95,7 +95,7 @@ end
 
 function correct_bifurcation(contres::ContResult{<: Union{FoldPeriodicOrbitCont, PDPeriodicOrbitCont, NSPeriodicOrbitCont}})
     if contres.prob.prob isa FoldProblemMinimallyAugmented
-        conversion = Dict(:bp => :R1, :hopf => :foldNS, :fold => :cusp, :nd => :nd, :pd => :foldpd)
+        conversion = Dict(:bp => :R1, :hopf => :foldNS, :fold => :cusp, :nd => :nd, :pd => :foldpd, :bt => :R1)
     elseif contres.prob.prob isa PeriodDoublingProblemMinimallyAugmented
         conversion = Dict(:bp => :foldFlip, :hopf => :pdNS, :pd => :R2,)
     elseif contres.prob.prob isa NeimarkSackerProblemMinimallyAugmented

diff --git a/test/stuartLandauCollocation.jl b/test/stuartLandauCollocation.jl
@@ -227,11 +227,12 @@ prob_ana = BifurcationProblem(idvf, zeros(N), par_hopf, (@lens _.r) ; J = (x,p)
 prob_col = BK.PeriodicOrbitOCollProblem(Ntst, m; prob_vf = prob_ana, N = N, ϕ = rand(N*( 1 + m * Ntst)), xπ = rand(N*( 1 + m * Ntst)))
 _ci = BK.generate_solution(prob_col, t->cos(t) .* ones(N), 2pi);
 Jcofd = ForwardDiff.jacobian(z->prob_col(z, par_sl), _ci);
-Jco = BK.analytical_jacobian(prob_col, _ci, par_sl);
+Jco = BK.analytical_jacobian(prob_col, _ci, par_sl); # 0.006155 seconds (21.30 k allocations: 62.150 MiB)
 @test norminf(Jcofd - Jco) < 1e-15
 
 # same but with Stuart-Landau vector field
 N = 2
+@assert N == 2 "must be the dimension of the SL"
 Ntst = 20
 m = 4
 prob_col = BK.PeriodicOrbitOCollProblem(Ntst, m; prob_vf = probsl, N = N, ϕ = rand(N*( 1 + m * Ntst)), xπ = rand(N*( 1 + m * Ntst)))