Switching between Hopf & saddle-node continuations

src/base/base.jl

@@ -1,6 +1,7 @@
 module BifurcationsBase
 using Compat
-using ..Continuations: AbstractContinuationCache, AbstractProblemCache
+using ..Continuations: AbstractContinuationCache, AbstractProblemCache,
+    AbstractContinuationSolver
 include("timekind.jl")
 include("statekind.jl")
 include("contkind.jl")

src/base/contkind.jl

@@ -16,3 +16,5 @@ function contkind(cache::AbstractContinuationCache)  # TODO: trait
     prob = as(cache, ContinuationCache).prob_cache.prob  # TODO: interface
     return contkind(prob)
 end
+
+contkind(solver::AbstractContinuationSolver) = contkind(solver.prob)

src/codim2/augmented_systems.jl

@@ -45,3 +45,6 @@ function preferred_augsys(point)
     end
     return BackReferencingAS()
 end
+
+preferred_augsys(::HopfCont) = BackReferencingAS()
+preferred_augsys(::SaddleNodeCont) = NormalizingAS()

src/codim2/codim2.jl

@@ -27,6 +27,7 @@ include("problem.jl")
 include("diffeq.jl")
 include("solver.jl")
 include("analysis.jl")
+include("switch.jl")
 include("api.jl")
 include("display.jl")
 include("diagnosis.jl")

src/codim2/diffeq.jl

@@ -170,6 +170,9 @@ end
 ds_eigvec(prob::DiffEqCodim2Problem, u::AbstractArray) =
     _ds_eigvec(eltype(prob.v0), u::AbstractArray)
 
+ds_eigvec(::HopfCont, u::AbstractArray) = _ds_eigvec(Complex, u)
+ds_eigvec(::SaddleNodeCont, u::AbstractArray) = _ds_eigvec(Real, u)
+
 function _ds_eigvec(::Type{E}, u::AbstractArray) where {E <: Real}
     d = dims_from_augsys(length(u), E)
     return @view u[eigvec_range(d)]

src/codim2/switch.jl

@@ -0,0 +1,125 @@
+using ..Continuations: _normalize!
+using ..ArrayUtils: _eig
+
+other(::HopfCont) = SaddleNodeCont()
+other(::SaddleNodeCont) = HopfCont()
+
+orthogonalize(v, u) = v .- u .* ((u ⋅ v) / (u ⋅ u))
+
+
+function BifurcationProblem(point::AbstractSpecialPoint,
+                            solver::Codim2Solver;
+                            t_domain = solver.prob.t_domain,
+                            param_axis1 = solver.prob.param_axis1,
+                            param_axis2 = solver.prob.param_axis2,
+                            kwargs...)
+
+    cd2_prob = solver.prob :: DiffEqCodim2Problem
+    de_prob = cd2_prob.de_prob :: AbstractODEProblem
+
+    # Manually cast.  The fact that `resolved.u` is not of type xtype
+    # indicates type instability somewhere.  Until it is fixed, let's
+    # silently cast always.  After type stability is established, it
+    # has to be changed to emit warning before cast.
+    xtype = typeof(de_prob.u0)
+
+    resolved = resolve_point(point, solver)
+    N = length(de_prob.u0)
+    x0 = xtype(resolved.u[1:N])
+    t0 = SVector{2}(resolved.u[end-1:end])
+    @assert param_axis1 == solver.prob.param_axis1
+    @assert param_axis2 == solver.prob.param_axis2
+
+    @assert point.point_type in (PointTypes.bogdanov_takens,
+                                 PointTypes.fold_hopf)
+    if (contkind(solver) :: ContinuationKind) isa SaddleNodeCont
+        next_ckind = HopfCont()
+    else
+        next_ckind = SaddleNodeCont()
+    end
+
+    @assert timekind(point) isa Continuous
+
+    v0, w0 = next_eig(next_ckind, resolved)
+    @assert length(v0) == length(x0)
+    if next_ckind isa HopfCont
+        v0 = v0 .+ 0im
+    else
+        v0 = real.(v0)
+    end
+    v0 = cast_container(xtype, v0)
+
+    # Sanity checks:
+    if next_ckind isa SaddleNodeCont
+        @assert eltype(v0) <: Real
+    else
+        @assert eltype(v0) <: Complex
+    end
+    if xtype <: SVector
+        v0 :: SVector
+    end
+
+    return DiffEqCodim2Problem(
+        de_prob,
+        param_axis1,
+        param_axis2,
+        t_domain;
+        x0 = x0,
+        v0 = v0,
+        w0 = w0,
+        t0 = t0,
+        augmented_system = preferred_augsys(next_ckind),
+        kwargs...)
+end
+
+next_eig(next_ckind::ContinuationKind, point::SpecialPoint) =
+    next_eig(next_ckind, Val(point.point_type), point)
+
+# Switching to Hopf continuation via Bogdanov-Takens from saddle-node
+function next_eig(::HopfCont, ::Val{PointTypes.bogdanov_takens},
+                  point)
+    prev_ckind = SaddleNodeCont()
+    v_prev = ds_eigvec(prev_ckind, point.u)
+    vals, vecs = _eig(ds_jacobian(prev_ckind, point.J))
+
+    if eltype(vals) <: Complex
+        _, i0 = findmin(abs.(real.(vals)))
+        return _normalize!(vecs[:, i0]), imag(vals[i0])
+    end
+
+    i1, i2 = sortperm(vals, by=abs)
+    v1 = vecs[:, i1]
+    v2 = vecs[:, i2]
+    if abs(v1 ⋅ v_prev) > abs(v2 ⋅ v_prev)
+        v_img = v1
+    else
+        v_img = v2
+    end
+    v_next = _normalize!(@. v_prev + v_img * im)
+    return v_next, zero(real(eltype(v_prev)))
+end
+
+# Switching to saddle-node continuation via Bogdanov-Takens from Hopf
+function next_eig(::SaddleNodeCont, ::Val{PointTypes.bogdanov_takens},
+                  point)
+    prev_ckind = HopfCont()
+    v_prev = ds_eigvec(prev_ckind, point.u)
+    v_next = _normalize!(real.(v_prev))
+    return v_next, zero(eltype(v_next))
+end
+
+# Switching to Hopf (saddle-node) continuation via Fold-Hopf from
+# saddle-node (Hopf)
+function next_eig(next_ckind::ContinuationKind,
+                  ::Val{PointTypes.fold_hopf},
+                  point)
+    prev_ckind = other(next_ckind)
+    vals, vecs = _eig(ds_jacobian(prev_ckind, point.J))
+    i1, i2, i3 = sortperm(vals, by=abs∘real)
+    ir, _, ii = sort((i1, i2, i3), by=i -> abs(imag(vals[i])))
+    if next_ckind isa HopfCont
+        return _normalize!(vecs[:, ii]), imag(vals[ii])
+    else
+        return _normalize!(vecs[:, ir]), imag(vals[ir])
+    end
+end

test/test_bazykin_85.jl

@@ -2,12 +2,15 @@ module TestBazykin85
 include("preamble_plots.jl")
 
 using Bifurcations: Codim2, resolved_points
+using Bifurcations.BifurcationsBase: contkind, HopfCont
 using Bifurcations.Examples: Bazykin85
 
 @testset "smoke Bazykin85 codim-2" begin
     codim1_solver = init(Bazykin85.prob)
     solve!(codim1_solver)
 
+    hopf_solver1 = nothing
+    sn_solver1 = nothing
     point_list = sort!(special_points(codim1_solver), by=p->p.u0[end])
     for point in point_list[2:3]  # TODO: avoid manual indexing
         codim2_prob = BifurcationProblem(
@@ -37,8 +40,39 @@ using Bifurcations.Examples: Bazykin85
                                         Codim2.PointTypes.hopf_hopf]
             @test_nothrow resolved_points(codim2_solver, p)
         end
+
+        if contkind(codim2_solver) isa HopfCont
+            hopf_solver1 = codim2_solver
+        else
+            sn_solver1 = codim2_solver
+        end
     end
-end
 
+    @testset "switch to saddle-node" begin
+        @assert hopf_solver1 !== nothing
+        # Find the right-most Bogdanov-Takens bifurcation:
+        point = first(sort(
+            special_points(hopf_solver1,
+                           Codim2.PointTypes.bogdanov_takens);
+            by = p -> p.u0[end - 1],  # α
+            rev = true))
+        sn_prob = BifurcationProblem(point, hopf_solver1)
+        sn_solver2 = init(sn_prob)
+        @test_nothrow solve!(sn_solver2)
+    end
+
+    @testset "switch to hopf" begin
+        @assert sn_solver1 !== nothing
+        # Find the right-most Bogdanov-Takens bifurcation:
+        point = first(sort(
+            special_points(sn_solver1,
+                           Codim2.PointTypes.bogdanov_takens);
+            by = p -> p.u0[end - 1],  # α
+            rev = true))
+        hopf_prob = BifurcationProblem(point, sn_solver1)
+        hopf_solver2 = init(hopf_prob)
+        @test_nothrow solve!(hopf_solver2)
+    end
+end
 
 end  # module