From f890a3d43581f8cf367405cc2db1b4ffe7643ac1 Mon Sep 17 00:00:00 2001 From: Takafumi Arakaki Date: Wed, 27 Jun 2018 20:58:14 -0700 Subject: [PATCH] Switching between Hopf & saddle-node continuations --- src/base/base.jl | 3 +- src/base/contkind.jl | 2 + src/codim2/augmented_systems.jl | 3 + src/codim2/codim2.jl | 1 + src/codim2/diffeq.jl | 3 + src/codim2/switch.jl | 125 ++++++++++++++++++++++++++++++++ test/test_bazykin_85.jl | 36 ++++++++- 7 files changed, 171 insertions(+), 2 deletions(-) create mode 100644 src/codim2/switch.jl diff --git a/src/base/base.jl b/src/base/base.jl index 1a50097..787c5fb 100644 --- a/src/base/base.jl +++ b/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") diff --git a/src/base/contkind.jl b/src/base/contkind.jl index 662810a..627dd98 100644 --- a/src/base/contkind.jl +++ b/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) diff --git a/src/codim2/augmented_systems.jl b/src/codim2/augmented_systems.jl index 2e16622..3e9161c 100644 --- a/src/codim2/augmented_systems.jl +++ b/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() diff --git a/src/codim2/codim2.jl b/src/codim2/codim2.jl index f6d487d..3c2b8a1 100644 --- a/src/codim2/codim2.jl +++ b/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") diff --git a/src/codim2/diffeq.jl b/src/codim2/diffeq.jl index fe1c104..a38c13d 100644 --- a/src/codim2/diffeq.jl +++ b/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)] diff --git a/src/codim2/switch.jl b/src/codim2/switch.jl new file mode 100644 index 0000000..1561223 --- /dev/null +++ b/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 diff --git a/test/test_bazykin_85.jl b/test/test_bazykin_85.jl index 42e3f21..4e07bee 100644 --- a/test/test_bazykin_85.jl +++ b/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