Skip to content

Commit

Permalink
Make bautin_to_subcritical_lc work with SpecialPoint
Browse files Browse the repository at this point in the history
  • Loading branch information
tkf committed Jul 27, 2018
1 parent 5ca63d7 commit 00a3058
Showing 1 changed file with 30 additions and 1 deletion.
31 changes: 30 additions & 1 deletion src/codim2lc/factories.jl
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,9 @@ function FoldLimitCycleProblem(point::AbstractSpecialPoint,

flc_points = special_points(lc_solver,
Codim1LimitCycle.PointTypes.saddle_node)
if isempty(flc_points)
error("Cannot find fold bifurcation of the limit cycle.")
end
flc_point = resolve_point(flc_points[1], lc_solver)
if length(flc_points) > 1
param2 = get(solver.prob.param_axis2, lc_solver.prob.de_prob.p)
Expand Down Expand Up @@ -114,14 +117,40 @@ 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

coeff_1 = first_lyapunov_coefficient(prob_cache, point.u1, point.J1)
if coeff_1 > 0 # u1 is subcritical
u_sc = point.u1
Expand Down

0 comments on commit 00a3058

Please sign in to comment.