Make via-LC-fold strategy work with test_reparametrized_bautin.jl

src/codim2lc/factories.jl

@@ -2,21 +2,24 @@ using DiffEqBase: AbstractODEProblem, remake
 using StaticArrays: SVector
 using Setfield: get, set
 
-using ..Continuations: init, solve!, residual_jacobian
+using ..Continuations: init, solve!, residual_jacobian, ContinuationSolver,
+    ZeroNotFoundError, FindZeroInputError
 using ..BifurcationsBase: AbstractSpecialPoint, SpecialPoint, special_points
 using ..Codim1
-using ..Codim1: resolve_point
+using ..Codim1: resolving_points
 using ..Codim2
 using ..Codim2: Codim2Solver, DiffEqCodim2Problem, first_lyapunov_coefficient
-using ..Codim1LimitCycle: Codim1LCSolver
+using ..Codim1LimitCycle: Codim1LCSolver, diameter
 
 # See also:
 # * [[../codim2/problem.jl::function BifurcationProblem]]
 # * [[../codim2/switch.jl::function BifurcationProblem]]
 # * [[../codim1lc/factories.jl::function LimitCycleProblem]]
-function FoldLimitCycleProblem(point::AbstractSpecialPoint,
-                               solver::Codim2Solver;
-                               kwargs...)
+function FoldLimitCycleProblem(
+        point::AbstractSpecialPoint,
+        solver::Codim2Solver;
+        lc_solver_opts = [],
+        kwargs...)
     @assert point.point_type == Codim2.PointTypes.bautin
     @assert timekind(point) isa Continuous
 
@@ -29,42 +32,85 @@ function FoldLimitCycleProblem(point::AbstractSpecialPoint,
         lc_prob;
         start_from_nearest_root = true,
         max_branches = 0,
-        bidirectional_first_sweep = false,
         nominal_angle_rad = 2π * (1 / 360),  # TODO: should it be the default?
+        lc_solver_opts...
     )
     solve!(lc_solver)
 
-    flc_points = special_points(lc_solver,
-                                Codim1LimitCycle.PointTypes.saddle_node)
-    flc_point = resolve_point(flc_points[1], lc_solver)
-    if length(flc_points) > 1
+    PT = Codim1LimitCycle.PointTypes.saddle_node
+    point_candidates = special_points(lc_solver, PT)
+    if isempty(point_candidates)
+        error("Cannot find fold bifurcation of the limit cycle.")
+    end
+    sw1_points = special_points(lc_solver.sol.sweeps[1], PT)
+    sw2_points = special_points(lc_solver.sol.sweeps[2], PT)
+    sweep_index = -1
+    if length(sw2_points) == 0
+        @assert length(sw1_points) > 0
+        flc_point, = sw1_points
+        sweep_index = 1
+    elseif length(sw1_points) == 0
+        @assert length(sw2_points) > 0
+        flc_point, = sw2_points
+        sweep_index = 2
+    else
+        # Choose increasing sweep
+        cont_sweeps = as(lc_solver, ContinuationSolver).sol.sweeps
+        at_most = (n, u) -> u[1:max(length(u), n)]
+        ds1 = [diameter(u, lc_prob) for u in at_most(10, cont_sweeps[1].u)]
+        ds2 = [diameter(u, lc_prob) for u in at_most(10, cont_sweeps[2].u)]
+        mean_dd1 = mean(diff(ds1))
+        mean_dd2 = mean(diff(ds2))
+        if mean_dd1 < 0 && mean_dd2 < 0
+            @warn "Both sweeps have decreasing diameters."
+        end
+        if mean_dd1 > mean_dd2
+            flc_point, = sw1_points
+            sweep_index = 1
+        else
+            flc_point, = sw2_points
+            sweep_index = 2
+        end
+    end
+
+    # TODO?: resolve flc_point
+    u_flc = (flc_point.u0 .+ flc_point.u1) ./ 2
+    prob_cache = as(lc_solver.cache, ContinuationCache).prob_cache
+    _H, J_flc = residual_jacobian(u_flc, prob_cache)
+
+    if length(point_candidates) > 1
         param2 = get(solver.prob.param_axis2, lc_solver.prob.de_prob.p)
         @warn """
-        $(length(flc_points)) fold bifurcations are found.
-        Using the first one found:
-            $(solver.prob.param_axis1) : $(flc_point.u[end])
+        $(length(point_candidates)) fold bifurcations are found.
+        Using the first point found in sweep $sweep_index:
+            $(solver.prob.param_axis1) : $(u_flc[end])
             $(solver.prob.param_axis2) : $(param2) (fixed)
+            diameter : $(diameter(u_flc, lc_prob))
         """
     end
 
-    n = length(lc_prob.xs0)
-    vals, vecs = eig(@view flc_point.J[1:n, 1:n])
-    _, idx0 = findmin(abs.(real.(vals)))  # real closest to zero
-    @assert all(x -> imag(x) ≈ 0, vecs[:, idx0])
+    n = length(lc_prob.xs0) + 1  # +1 to include period
+    vals, vecs = eig(@view J_flc[1:n, 1:n])
+    _, idx0 = findmin(abs.(vals))  # closest to zero
+    opts = as(solver, ContinuationSolver).opts
+    if any(x -> abs(imag(x)) > opts.atol, @view vecs[:, idx0])
+        @warn "Nearest-to-null vector is complex (eigval: $(vals[idx0]))"
+    end
     v0 = normalize(real(@view vecs[:, idx0]))
 
-    xs0 = reshape(flc_point.u[1:length(lc_prob.xs0)], size(lc_prob.xs0))
-    l0 = flc_point.u[length(lc_prob.xs0) + 1]
-    vs0 = reshape(v0, size(lc_prob.xs0))
-    t0 = SVector(flc_point.u[end],
+    xs0 = reshape(u_flc[1:length(lc_prob.xs0)], size(lc_prob.xs0))
+    l0 = u_flc[length(lc_prob.xs0) + 1]
+    vs0 = reshape(v0[1:end-1], size(lc_prob.xs0))
+    dl0 = v0[end]
+    t0 = SVector(u_flc[end],
                  get(solver.prob.param_axis2, lc_solver.prob.de_prob.p))
 
     return FoldLimitCycleProblem(
         lc_prob;
         xs0 = xs0,
         l0 = l0,
         vs0 = vs0,
-        dl0 = 0.0,  # ???
+        dl0 = dl0,
         t0 = t0,
         param_axis1 = solver.prob.param_axis1,
         param_axis2 = solver.prob.param_axis2,
@@ -114,14 +160,43 @@ function bautin_to_subcritical_lc(
         solver::Codim2Solver;
         kwargs...)
 
-    point :: SpecialPointInterval  # TODO: support SpecialPoint
     @assert point.point_type == Codim2.PointTypes.bautin
     @assert timekind(point) isa Continuous
 
     solver.prob :: DiffEqCodim2Problem
     de_prob = solver.prob.de_prob :: AbstractODEProblem
     prob_cache = as(solver.cache, ContinuationCache).prob_cache
 
+    if point isa SpecialPoint
+        point = let
+            if point.sweep.value == nothing
+                error("""
+                `point` is not associated with a sweep.
+                point = $(point)
+                 """)
+            end
+            # See: [[../base/solver.jl::SpecialPointInterval]]
+            sweep = as(point.sweep.value, ContinuationSweep)
+            u0 = sweep.u[point.point_index - 1]
+            u1 = sweep.u[point.point_index]
+            _H, J1 = residual_jacobian(u1, prob_cache)
+            SpecialPointInterval(
+                timekind(point),
+                point.point_type,
+                point.point_index,
+                u0,
+                u1,
+                J1,
+                WeakRef(sweep),
+            )
+        end
+    else
+        point :: SpecialPointInterval
+    end
+
+    # TODO: Specify target l1 (> 0) value and step into it using zero
+    # finder (so that hopefully it's not too close to/far away from
+    # the Bautin point).
     coeff_1 = first_lyapunov_coefficient(prob_cache, point.u1, point.J1)
     if coeff_1 > 0  # u1 is subcritical
         u_sc = point.u1

src/continuations/zero_point.jl

@@ -27,6 +27,26 @@ u1 = $(e.u1)
 f = $(e.f)
 """)
 
+@with_kw struct ZeroNotFoundError <: Exception
+    msg = "Zero not found"
+    u::AbstractVector
+    f::Real
+    h::Real
+end
+
+function Base.showerror(io::IO, e::ZeroNotFoundError)
+    println(io, e.msg)
+    println(io, "h = ", e.h)
+    println(io, "f = ", e.f)
+    println(io, "u =")
+    show(IOContext(io,
+                   :displaysize => displaysize(io),
+                   :limit => true),
+         MIME("text/plain"),
+         e.u')
+    println(io)
+end
+
 function find_zero!(cache, opts, f, u0, u1, direction)
     prob_cache = cache.prob_cache
     H = cache.H
@@ -66,7 +86,7 @@ function find_zero!(cache, opts, f, u0, u1, direction)
             end
             v = w
         end
-        error("corrector failed")
+        throw(ZeroNotFoundError(msg="Corrector failed", u=v, f=fu, h=h))
 
         @label corrector_success
         tJv = tangent(L, Q)
@@ -87,5 +107,5 @@ function find_zero!(cache, opts, f, u0, u1, direction)
         end
         tJ = tJv
     end
-    error("zero not found")
+    throw(ZeroNotFoundError(u=u, f=fu, h=h))
 end

test/test_reparametrized_bautin.jl

@@ -70,7 +70,8 @@ using Bifurcations.Examples.Reparametrization: orig_p
     β_bautin = codim2_points[1].u[end-1:end]
     @test all(@. abs(β_bautin) < 1e-6)
 
-    flc_prob = FoldLimitCycleProblem(
+    @info "Creating FoldLimitCycleProblem..."
+    @time flc_prob = FoldLimitCycleProblem(
         codim2_points[1],
         hopf_solver;
         num_mesh = 20,
@@ -122,7 +123,7 @@ using Bifurcations.Examples.Reparametrization: orig_p
     @show maximum(@. abs(4 * flc_β₁ + flc_β₂^2))
     @test all(@. abs(4 * flc_β₁ + flc_β₂^2) < 5e-2)
     @test maximum(flc_β₂) > 2
-    @test minimum(flc_β₂) > -1e-3
+    @test minimum(flc_β₂) > -5e-2
 end
 
 end  # module