correct erroneous Floquet coefficients for odd Ntst in case of colloc…

…ation

src/NormalForms.jl

@@ -238,7 +238,7 @@ function predictor(bp::Union{Pitchfork, PitchforkMap}, ds::T; verbose = false, a
 
     # we need to find the type, supercritical or subcritical
     dsfactor = b1 * b3 < 0 ? T(1) : T(-1)
-    if 1==1
+    if true
         # we solve b1 * ds + b3 * amp^2 / 6 = 0
         amp = ampfactor * sqrt(-6abs(ds) * dsfactor * b1 / b3)
         pnew = bp.p + abs(ds) * dsfactor

src/codim2/NormalForms.jl

@@ -1016,7 +1016,7 @@ function zero_hopf_normal_form(_prob,
         tol_ev = max(1e-10, 10abs(imag(_λ[_ind0])))
         # imaginary eigenvalue iω1
         _ind2 = [ii for ii in eachindex(_λ) if ((abs(imag(_λ[ii])) > tol_ev) & (ii != _ind0))]
-        verbose && (@info "EV" _λ _ind2)
+        verbose && (@info "Eigenvalue :" _λ _ind2)
         _indIm = argmin(abs(real(_λ[ii])) for ii in _ind2)
         λI = _λ[_ind2[_indIm]]
         q1 = geteigenvector(optionsN.eigsolver, _ev, _ind2[_indIm])

src/periodicorbit/Floquet.jl

@@ -459,9 +459,13 @@ end
     # this removes the internal unknowns of each mesh interval
     # this matrix is diagonal by blocks and each block is the L Matrix
     # which makes the corresponding J block upper triangular
-    # the following matrix collects the LU factorizations by blocks
+    # the following matrix 𝐅𝐬 collects the LU factorizations by blocks
+    # recall that if F = lu(A) then
+    # F.L * F.U = A * F.P
+    # (F.P⁻¹ * F.L) * F.U = A
+    # hence 𝐅𝐬⁻¹ = (P⁻¹ * L)⁻¹ = L⁻¹ * P
     𝐅𝐬 = Matrix{𝒯}(LinearAlgebra.I(size(J, 1)))
-    rg = 1:nbcoll # range
+    rg = 1:nbcoll
     for k = 1:Ntst
         F = lu(J[rg, rg .+ n])
         𝐅𝐬[rg, rg] .= (F.P \ F.L)
@@ -487,11 +491,14 @@ end
         r2  = r2 .+ m * n
     end
 
+    # in theory, it should be multiplied by (-1)ᴺᵗˢᵗ
+    factor = iseven(Ntst) ? 1 : -1
+
     # floquet multipliers
     vals, vecs = eigen(M)
 
     nev = min(n, nev)
-    logvals = log.(Complex.(vals))
+    logvals = @. log(Complex(factor * vals))
     I = sortperm(logvals, by = real, rev = true)[1:nev]
 
     # floquet exponents

src/periodicorbit/PeriodicOrbitCollocation.jl

@@ -626,7 +626,16 @@ Compute the jacobian of the problem defining the periodic orbits by orthogonal c
     J[end, end] = -phase / period
     return J
 end
-analytical_jacobian(coll::PeriodicOrbitOCollProblem, u::AbstractVector, pars; 𝒯 = eltype(u), k...) = analytical_jacobian!(zeros(𝒯, length(coll)+1, length(coll)+1), coll, u, pars; k...)
+
+analytical_jacobian(coll::PeriodicOrbitOCollProblem, 
+                            u, 
+                            pars; 
+                            𝒯 = eltype(u), 
+                            k...) = analytical_jacobian!(zeros(𝒯, length(coll)+1, length(coll)+1), 
+                                                        coll, 
+                                                        u, 
+                                                        pars; 
+                                                        k...)
 
 function analytical_jacobian_sparse(coll::PeriodicOrbitOCollProblem,
                                     u::AbstractVector,
@@ -1181,8 +1190,9 @@ end
     nbcoll = N * m
     nⱼ = size(J, 1)
     𝒯 = eltype(coll)
+    In = I(N)
 
-    P = Matrix{𝒯}(LinearAlgebra.I(nⱼ))
+    𝐅𝐬 = Matrix{𝒯}(I(nⱼ))
     blockⱼ = zeros(𝒯, nbcoll, nbcoll)
     rg = 1:nbcoll
     for _ in 1:Ntst
@@ -1203,11 +1213,10 @@ end
     rN = 1:N
 
     # solving for the external variables
-    In = I(N)
     Jext = zeros(Ntst*N+N+1, Ntst*N+N+1)
     Jext[end-N:end-1,end-N:end-1] .= In
     Jext[end-N:end-1,1:N] .= -In
-    Jext[end, end] = J[end,end]
+    Jext[end, end] = J[end, end]
     rhs_ext = 𝒯[]
 
     # we solve for the external unknowns