double the mesh size at branching at PD points for trapeze

change all predictors to return po, this is needed for deflated newton

src/BifurcationKit.jl

@@ -159,7 +159,7 @@ module BifurcationKit
     export bifurcationdiagram, bifurcationdiagram!, Branch, BifDiagNode, get_branch, get_branches_from_BP
 
     # Periodic orbit computation
-    export generate_solution, getperiod, getamplitude, getmaximum, get_periodic_orbit, guess_from_hopf, generate_ci_problem
+    export generate_solution, getperiod, get_periodic_orbit, guess_from_hopf, generate_ci_problem
 
     # Periodic orbit computation based on Trapeze method
     export PeriodicOrbitTrapProblem, continuation_potrap

src/periodicorbit/NormalForms.jl

@@ -993,14 +993,31 @@ function predictor(nf::PeriodDoublingPO{ <: PeriodicOrbitTrapProblem},
                     override = false)
     pb = nf.prob
     M, N = size(pb)
+    # we update the problem by doubling M
+    if pb.mesh isa TimeMesh{Int64}
+        # we double the mesh size
+        pbnew = @set pb.mesh.ds = 2M
+    else
+        oldmesh = get_times(pb)
+        new_mesh = vcat(old_mesh[begin:end-1] /2, old_mesh ./2 .+ 1/2)
+        pbnew = @set pb.mesh.ds = new_mesh
+    end
+    @reset pbnew.M = 2M
+
     orbitguess0c = get_time_slices(pb, nf.po)
     ζc = reshape(nf.ζ, N, M)
     orbitguess_c = orbitguess0c .+ ampfactor .*  ζc
     orbitguess_c = hcat(orbitguess_c, orbitguess0c .- ampfactor .*  ζc)
-    orbitguess = vec(orbitguess_c[:,1:2:end])
+    orbitguess = vec(orbitguess_c)
     # we append twice the period
     orbitguess = vcat(orbitguess, 2nf.T)
-    return (;orbitguess, pnew = nf.nf.p + δp, prob = pb, ampfactor)
+    # we updat the phase condition
+    @reset pbnew.xπ = orbitguess[begin:end-1]
+    @reset pbnew.ϕ = circshift(orbitguess[begin:end-1], length(orbitguess))
+    # we need to duplicate the po as well in case deflation is used
+    po0 = get_time_slices(pb, nf.po)
+    po = vcat(vec(hcat(po0, po0)), nf.T) 
+    return (;orbitguess, pnew = nf.nf.p + δp, prob = pbnew, ampfactor, po)
 end
 
 function predictor(nf::BranchPointPO{ <: PeriodicOrbitTrapProblem},
@@ -1009,15 +1026,15 @@ function predictor(nf::BranchPointPO{ <: PeriodicOrbitTrapProblem},
                     override = false)
     orbitguess = copy(nf.po)
     orbitguess[begin:end-1] .+= ampfactor .* nf.ζ
-    return (;orbitguess, pnew = nf.nf.p + δp, prob = nf.prob, ampfactor)
+    return (;orbitguess, pnew = nf.nf.p + δp, prob = nf.prob, ampfactor, po = nf.po)
 end
 
 function predictor(nf::NeimarkSackerPO,
                     δp,
                     ampfactor;
                     override = false)
     orbitguess = copy(nf.po)
-    return (;orbitguess, pnew = nf.nf.p + δp, prob = nf.prob, ampfactor)
+    return (;orbitguess, pnew = nf.nf.p + δp, prob = nf.prob, ampfactor, po = nf.po)
 end
 ####################################################################################################
 function predictor(nf::PeriodDoublingPO{ <: PeriodicOrbitOCollProblem }, 
@@ -1068,8 +1085,11 @@ function predictor(nf::PeriodDoublingPO{ <: PeriodicOrbitOCollProblem },
     # we append the doubled period
     orbitguess = vcat(orbitguess, 2nf.T * time_factor)
 
+    # we double po in case of use of deflated newton
+    po = vcat(orbitguess0[begin:end-N], orbitguess0, nf.T)
+
     # no need to change pbnew.cache
-    return (;orbitguess, pnew = nf.nf.p + δp, prob = pbnew, ampfactor, δp, time_factor)
+    return (;orbitguess, pnew = nf.nf.p + δp, prob = pbnew, ampfactor, δp, time_factor, po)
 end
 ####################################################################################################
 """
@@ -1097,7 +1117,10 @@ function predictor(nf::PeriodDoublingPO{ <: ShootingProblem },
     @reset pbnew.ds = _duplicate(pbnew.ds) ./ 2
     orbitguess[end] *= 2
     updatesection!(pbnew, orbitguess, setparam(pbnew, pnew))
-    return (;orbitguess, pnew, prob = pbnew, ampfactor, δp)
+
+    po = copy(nf.po)[begin:end-1]
+    po = vcat(po, copy(nf.po)[begin:end-1], nf.po[end])
+    return (;orbitguess, pnew, prob = pbnew, ampfactor, δp, po)
 end
 
 """
@@ -1112,7 +1135,7 @@ function predictor(nf::BranchPointPO{ <: ShootingProblem },
     ζs = nf.ζ
     orbitguess = copy(nf.po)
     orbitguess[eachindex(ζs)] .+= ampfactor .* ζs
-    return (;orbitguess, pnew = nf.nf.p + δp, prob = nf.prob, ampfactor)
+    return (;orbitguess, pnew = nf.nf.p + δp, prob = nf.prob, ampfactor, po = nf.po)
 end
 ####################################################################################################
 function predictor(nf::PeriodDoublingPO{ <: PoincareShootingProblem }, 
@@ -1126,8 +1149,8 @@ function predictor(nf::PeriodDoublingPO{ <: PoincareShootingProblem },
     @reset pbnew.M = pbnew.section.M
     orbitguess = copy(nf.po) .+ ampfactor .* ζs
     orbitguess = vcat(orbitguess, orbitguess .- ampfactor .* ζs)
-
-    return (;orbitguess, pnew = nf.nf.p + δp, prob = pbnew, ampfactor)
+    po = vcat(nf.po, nf.po)
+    return (;orbitguess, pnew = nf.nf.p + δp, prob = pbnew, ampfactor, po)
 end
 
 function predictor(nf::BranchPointPO{ <: PoincareShootingProblem},
@@ -1137,5 +1160,5 @@ function predictor(nf::BranchPointPO{ <: PoincareShootingProblem},
     ζs = nf.ζ
     orbitguess = copy(nf.po)
     orbitguess .+= ampfactor .* ζs
-    return (;orbitguess, pnew = nf.nf.p + δp, prob = nf.prob, ampfactor)
+    return (;orbitguess, pnew = nf.nf.p + δp, prob = nf.prob, ampfactor, po = nf.po)
 end

src/periodicorbit/PeriodicOrbits.jl

@@ -580,12 +580,11 @@ function continuation(br::AbstractResult{PeriodicOrbitCont, Tprob},
         optn = _contParams.newton_options
         # find point on the first branch
         pbnew = set_params_po(pbnew, setparam(br, newp))
-        sol0 = newton(pbnew, bifpt.x, optn; kwargs...)
+        sol0 = newton(pbnew, pred.po, optn; kwargs...)
         @assert converged(sol0) "The first guess did not converge"
 
         # find the bifurcated branch using deflation
-        @assert pbnew isa AbstractPOFDProblem || pbnew isa ShootingProblem "Deflated newton is not available for your problem. Try Trapezoid / collocation method or ShootingProblem"
-        deflationOp = DeflationOperator(2, (x, y) -> dot(x[1:end-1], y[1:end-1]), one(eltype(orbitguess)), [sol0.u]; autodiff = true)
+        deflationOp = DeflationOperator(2, (x, y) -> dot(x[begin:end-1], y[begin:end-1]), one(eltype(orbitguess)), [sol0.u]; autodiff = true)
         verbose && println("\n──> Compute point on the bifurcated branch...")
         solbif = newton(pbnew, orbitguess, deflationOp,
             (@set optn.max_iterations = 10 * optn.max_iterations) ; kwargs...,)

test/testLure.jl

@@ -116,8 +116,15 @@ for meshadapt in (false, true)
     predictor(pd, 0.1, 1)
     pd = get_normal_form(br_po, 1; verbose = false, prm = false)
     predictor(pd, 0.1, 1)
-    # @test pd.nf.nf.b3 ≈ -0.30509421737255177 rtol=1e-3 # reference value computed with ApproxFun
+    @test pd.nf.nf.b3 ≈ -0.30509421737255177 rtol=1e-3 # reference value computed with ApproxFun
     # @test pd.nf.nf.a  ≈ 0.020989802220981707 rtol=1e-3 # reference value computed with ApproxFun
+
+    # aBS from PD
+    continuation(br_po, 1, setproperties(br_po.contparams, detect_bifurcation = 3, max_steps = 5, ds = 0.01, dsmax = 0.01, plot_every_step = 10);
+    ampfactor = .2, δp = -0.005,
+    usedeflation = true,
+    )
+
 end
 ####################################################################################################
 using OrdinaryDiffEq
@@ -166,7 +173,7 @@ end
 # aBS from PD
 br_po_pd = continuation(br_po, 1, setproperties(br_po.contparams, detect_bifurcation = 3, max_steps = 5, ds = 0.01, plot_every_step = 1, save_sol_every_step = 1);
     # verbosity = 0, plot = false,
-    usedeflation = false,
+    usedeflation = true,
     ampfactor = .1, δp = -0.005,
     record_from_solution = recordPO,
     normC = norminf,