test PD NF with mesh adaptation

add ampfactor for predictor of bifurcations of PO
correct show method of PD
correct warning for large amplitude in IG of PO from Hopf

src/BifurcationPoints.jl

@@ -211,7 +211,7 @@ end
 function Base.show(io::IO, bp::PeriodDoubling)
     printstyled(io, bp.type, " - Period-Doubling ", color=:cyan, bold = true)
     println("bifurcation point at ", get_lens_symbol(bp.lens), " ≈ $(bp.p)")
-    println(io, "┌─ Normal form:\n├\t x⋅(a⋅δp - x + c⋅x³)")
+    println(io, "┌─ Normal form:\n├\t x⋅(a⋅δp - 1 + c⋅x²)")
     if ~isnothing(bp.nf)
         println(io,"├─ a = ", bp.nf.a)
         println(io,"└─ c = ", bp.nf.b3)

src/periodicorbit/BifurcationPoints.jl

@@ -54,8 +54,8 @@ function Base.show(io::IO, pd::PeriodDoublingPO)
         if ~pd.prm
             println("├─ ", get_lens_symbol(pd.nf.lens)," ≈ $(pd.nf.p)")
             println("├─ type: ", "$(pd.nf.type)")
-            println(io, "├─ (Iooss) Normal form:\n├\t∂τ = 1 + a₀⋅δp + a⋅ξ²\n├\t∂ξ = c₀⋅δp⋅ξ + c⋅ξ³")
-            println(io,"├─── a = ", pd.nf.nf.a, "\n└─── c = ", pd.nf.nf.b3)
+            println(io, "├─ (Iooss) Normal form:\n├\t∂τ = 1 + a₀⋅δp + a⋅ξ²\n├\t∂ξ = ξ⋅(c₀⋅δp + c⋅ξ²)")
+            println(io, "├─── a = ", pd.nf.nf.a, "\n└─── c = ", pd.nf.nf.b3)
         else
             show(io, pd.nf)
         end

src/periodicorbit/NormalForms.jl

@@ -59,7 +59,7 @@ function branch_normal_form(pbwrap,
     ζ = geteigenvector(br.contparams.newton_options.eigsolver, br.eig[bifpt.idx].eigenvecs, bifpt.ind_ev)
     # we normalize it by the sup norm because it could be too small/big in L2 norm
     # TODO: user defined scaleζ
-    ζ ./= norm(ζ, Inf)
+    ζ ./= norminf(ζ)
     verbose && println("Done!")
 
     # compute the full eigenvector
@@ -942,18 +942,18 @@ function predictor(nf::PeriodDoublingPO{ <: PeriodicOrbitTrapProblem}, δp, ampf
     orbitguess = vec(orbitguess_c[:,1:2:end])
     # we append twice the period
     orbitguess = vcat(orbitguess, 2nf.T)
-    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = pb)
+    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = pb, ampfactor = ampfactor)
 end
 
 function predictor(nf::BranchPointPO{ <: PeriodicOrbitTrapProblem}, δp, ampfactor)
     orbitguess = copy(nf.po)
     orbitguess[1:end-1] .+= ampfactor .*  nf.ζ
-    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = nf.prob)
+    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = nf.prob, ampfactor = ampfactor)
 end
 
 function predictor(nf::NeimarkSackerPO, δp, ampfactor)
     orbitguess = copy(nf.po)
-    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = nf.prob)
+    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = nf.prob, ampfactor = ampfactor)
 end
 ####################################################################################################
 function predictor(nf::PeriodDoublingPO{ <: PeriodicOrbitOCollProblem }, δp, ampfactor)
@@ -975,26 +975,27 @@ function predictor(nf::PeriodDoublingPO{ <: PeriodicOrbitOCollProblem }, δp, am
     orbitguess = vcat(orbitguess, 2nf.T)
 
     # no need to change pbnew.cache
-    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = pbnew)
+    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = pbnew, ampfactor = ampfactor)
 end
 ####################################################################################################
 function predictor(nf::PeriodDoublingPO{ <: ShootingProblem }, δp, ampfactor)
     pbnew = deepcopy(nf.prob)
+    pnew = nf.nf.p + δp
     ζs = nf.ζ
     orbitguess = copy(nf.po)[1:end-1] .+ ampfactor .* ζs
     orbitguess = vcat(orbitguess, copy(nf.po)[1:end-1] .- ampfactor .* ζs, nf.po[end])
 
     @set! pbnew.M = 2nf.prob.M
     @set! pbnew.ds = _duplicate(pbnew.ds) ./ 2
     orbitguess[end] *= 2
-    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = pbnew)
+    return (orbitguess = orbitguess, pnew = pnew, prob = pbnew, ampfactor = ampfactor)
 end
 
 function predictor(nf::BranchPointPO{ <: ShootingProblem }, δp, ampfactor)
     ζs = nf.ζ
     orbitguess = copy(nf.po)
     orbitguess[1:length(ζs)] .+= ampfactor .* ζs
-    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = nf.prob)
+    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = nf.prob, ampfactor = ampfactor)
 end
 ####################################################################################################
 function predictor(nf::PeriodDoublingPO{ <: PoincareShootingProblem }, δp, ampfactor)
@@ -1006,12 +1007,12 @@ function predictor(nf::PeriodDoublingPO{ <: PoincareShootingProblem }, δp, ampf
     orbitguess = copy(nf.po) .+ ampfactor .* ζs
     orbitguess = vcat(orbitguess, orbitguess .- ampfactor .* ζs)
 
-    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = pbnew)
+    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = pbnew, ampfactor = ampfactor)
 end
 
 function predictor(nf::BranchPointPO{ <: PoincareShootingProblem}, δp, ampfactor)
     ζs = nf.ζ
     orbitguess = copy(nf.po)
     orbitguess .+= ampfactor .* ζs
-    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = nf.prob)
+    return (orbitguess = orbitguess, pnew = nf.nf.p + δp, prob = nf.prob, ampfactor = ampfactor)
 end

src/periodicorbit/PeriodicOrbitCollocation.jl

@@ -561,7 +561,7 @@ Compute the jacobian of the problem defining the periodic orbits by orthogonal c
         mul!(ϕj, ϕc[:, rg], ∂L)
         # put the jacobian of the vector field
         for l in 1:m
-            if ~_transpose
+            if _transpose == false
                 J0 .= jacobian(coll.prob_vf, pj[:,l], pars)
             else
                 J0 .= transpose(jacobian(coll.prob_vf, pj[:,l], pars))

src/periodicorbit/PeriodicOrbits.jl

@@ -363,7 +363,7 @@ function continuation(br::AbstractBranchResult, ind_bif::Int,
                     nev = length(eigenvalsfrombif(br, ind_bif)),
                     kwargs...)
     # compute the normal form of the branch point
-    verbose = get(kwargs, :verbosity, 0) > 1 ? true : false
+    verbose = get(kwargs, :verbosity, 0) > 1
     verbose && (println("──▶ Considering bifurcation point:"); _show(stdout, br.specialpoint[ind_bif], ind_bif))
 
     cb = get(kwargs, :callback_newton, cb_default)
@@ -392,8 +392,8 @@ function continuation(br::AbstractBranchResult, ind_bif::Int,
             "\n├─── phase        ϕ = ", ϕ / pi, "⋅π",
             "\n├─ Method = \n", probPO, "\n")
 
-    if pred.amp < 0.1
-        @error "The amplitude for the predictor for the first periodic orbit on the bifurcated branch is not small $(pred.amp). You can either decrease `ds`, or specify how far `δp` from the bifurcation point you want the branch of periodic orbits to start. Alternatively, you can specify a multiplicative factor `ampfactor` to be applied to the predictor"
+    if pred.amp > 0.1
+        @warn "The amplitude of the first periodic orbit on the bifurcated branch obtained by the predictor is not small $(pred.amp). You can either decrease `ds`, or specify how far `δp` from the bifurcation point you want the branch of periodic orbits to start. Alternatively, you can specify a multiplicative factor `ampfactor` to be applied to the predictor amplitude."
     end
 
     M = get_mesh_size(probPO)
@@ -487,7 +487,7 @@ function continuation(br::AbstractResult{PeriodicOrbitCont, Tprob},
             ind_bif::Int,
             _contParams::ContinuationPar;
             alg = br.alg,
-            δp = 0.1, ampfactor = 1,
+            δp = _contParams.ds, ampfactor = 1,
             usedeflation = false,
             linear_algo = nothing,
             detailed = false,
@@ -499,15 +499,6 @@ function continuation(br::AbstractResult{PeriodicOrbitCont, Tprob},
     @assert bptype in (:pd, :bp) "Branching from $(bptype) not possible yet."
     @assert abs(bifpt.δ[1]) == 1 "Only simple bifurcation points are handled"
 
-    verbose = get(kwargs, :verbosity, 0) > 0
-    verbose && printstyled(color = :green, "━"^55*
-            "\n┌─ Start branching from $(bptype) point to periodic orbits.\n├─ Bifurcation type = ", bifpt.type, " [PRM = $(prm)]",
-            "\n├─── bif. param  p0 = ", bifpt.param,
-            "\n├─── period at bif. = ", getperiod(br.prob.prob, bifpt.x, setparam(br, bifpt.param)),
-            "\n├─── new param    p = ", bifpt.param + δp, ", p - p0 = ", δp,
-            "\n├─── amplitude p.o. = ", ampfactor,
-            "\n")
-
     _linear_algo = isnothing(linear_algo) ? BorderingBLS(_contParams.newton_options.linsolver) : linear_algo
 
     # we copy the problem for not mutating the one passed by the user. This is a AbstractPeriodicOrbitProblem.
@@ -519,6 +510,15 @@ function continuation(br::AbstractResult{PeriodicOrbitCont, Tprob},
     newp = pred.pnew  # new parameter value
     pbnew = pred.prob # modified problem
 
+    verbose = get(kwargs, :verbosity, 0) > 0
+    verbose && printstyled(color = :green, "━"^55*
+            "\n┌─ Start branching from $(bptype) point to periodic orbits.\n├─ Bifurcation type = ", bifpt.type, " [PRM = $(prm & (pbnew isa PeriodicOrbitOCollProblem))]",
+            "\n├─── bif. param  p0 = ", bifpt.param,
+            "\n├─── period at bif. = ", getperiod(br.prob.prob, bifpt.x, setparam(br, bifpt.param)),
+            "\n├─── new param    p = ", newp, ", p - p0 = ", newp - bifpt.param,
+            "\n└─── amplitude p.o. = ", pred.ampfactor,
+            "\n")
+
     # a priori, the following do not overwrite the options in br
     # hence the results / parameters in br are kept intact
     if pb isa AbstractShootingProblem

test/testLure.jl

@@ -1,4 +1,4 @@
-# using Revise, Plots
+# using Revise, Plots, AbbreviatedStackTraces
 using Parameters, LinearAlgebra, Test
 using BifurcationKit, Test
 const BK = BifurcationKit
@@ -26,6 +26,20 @@ bothside = true, normC = norminf)
 
 # plot(br)
 ####################################################################################################
+function plotPO(x, p; k...)
+	xtt = get_periodic_orbit(p.prob, x, p.p)
+	plot!(xtt.t, xtt[1,:]; markersize = 2, marker = :d, k...)
+	plot!(xtt.t, xtt[2,:]; k...)
+	plot!(xtt.t, xtt[3,:]; legend = false, k...)
+end
+
+# record function
+function recordPO(x, p)
+	xtt = get_periodic_orbit(p.prob, x, p.p)
+	period = getperiod(p.prob, x, p.p)
+	return (max = maximum(xtt[1,:]), min = minimum(xtt[1,:]), period = period)
+end
+####################################################################################################
 # newton parameters
 optn_po = NewtonPar(tol = 1e-8,  max_iterations = 25)
 
@@ -39,7 +53,8 @@ br_po = continuation(
         jacobian = :Dense);
         ampfactor = 1., δp = 0.01,
         verbosity = 0, plot = false,
-        record_from_solution = (x, p) -> (xtt=reshape(x[1:end-1],3,Mt); return (max = maximum(xtt[1,:]), min = minimum(xtt[1,:]), period = x[end])),
+        record_from_solution = recordPO,
+        plot_solution = plotPO,
         finalise_solution = (z, tau, step, contResult; prob = nothing, kwargs...) -> begin
                 return z.u[end] < 40
                 true
@@ -75,35 +90,38 @@ br_po_pd = continuation(br_po, 1, setproperties(br_po.contparams, detect_bifurca
 # plot(br, br_po, br_po_pd, xlims=(0.0,1.5))
 ####################################################################################################
 # case of collocation
-opts_po_cont = ContinuationPar(dsmax = 0.03, ds= -0.0001, dsmin = 1e-4, p_max = 1.8, p_min=-5., max_steps = 120, newton_options = (@set optn_po.tol = 1e-11), nev = 3, tol_stability = 1e-4, detect_bifurcation = 3, plot_every_step = 20, save_sol_every_step = 1, n_inversion = 8)
+opts_po_cont = ContinuationPar(dsmax = 0.03, ds= -0.0001, dsmin = 1e-4, p_max = 1.8, p_min=-0.9, max_steps = 120, newton_options = (@set optn_po.tol = 1e-11), nev = 3, tol_stability = 1e-4, detect_bifurcation = 3, plot_every_step = 20, save_sol_every_step = 1, n_inversion = 8)
 
-br_po = continuation(
-    br, 2, opts_po_cont,
-    BK.PeriodicOrbitOCollProblem(40, 4);
-    alg = PALC(tangent = Bordered()),
-    ampfactor = 1., δp = 0.01,
-    # usedeflation = true,
-    # verbosity = 2,    plot = true,
-    normC = norminf)
+for meshadapt in (false, true)
+    br_po = continuation(
+        br, 2, opts_po_cont,
+        PeriodicOrbitOCollProblem(40, 4; meshadapt, K = 200);
+        alg = PALC(),
+        ampfactor = 1., δp = 0.01,
+        record_from_solution = recordPO,
+        plot_solution = plotPO,
+        # verbosity = 2,    plot = true,
+        normC = norminf)
 
-# test normal forms
-for _ind in (1,)
-    if length(br_po.specialpoint) >=1 && br_po.specialpoint[_ind].type ∈ (:bp, :pd, :ns)
-        println("")
-        for prm in (true, false)
-            pt = get_normal_form(br_po, _ind; verbose = true, prm)
-            show(pt)
+    # test normal forms
+    for _ind in (1,)
+        if length(br_po.specialpoint) >=1 && br_po.specialpoint[_ind].type ∈ (:bp, :pd, :ns)
+            println("")
+            for prm in (true, false)
+                pt = get_normal_form(br_po, _ind; verbose = true, prm)
+                show(pt)
+            end
         end
     end
-end
 
-pd = get_normal_form(br_po, 1; verbose = false, prm = false)
-@test pd.nf.nf.b3 ≈ -0.30509421737255177 rtol=1e-4 # reference value computed with ApproxFun
-@test pd.nf.nf.a  ≈ 0.020989802220981707 rtol=1e-3 # reference value computed with ApproxFun
+    pd = get_normal_form(br_po, 1; verbose = false, prm = false)
+    @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
+end
 ####################################################################################################
 using OrdinaryDiffEq
 
-probsh = ODEProblem(lur!, copy(z0), (0., 1000.), par_lur; abstol = 1e-10, reltol = 1e-8)
+probsh = ODEProblem(lur!, copy(z0), (0., 1000.), par_lur; abstol = 1e-12, reltol = 1e-10)
 
 optn_po = NewtonPar(tol = 1e-12, max_iterations = 25)
 
@@ -113,10 +131,10 @@ opts_po_cont = ContinuationPar(dsmax = 0.02, ds= -0.001, dsmin = 1e-4, max_steps
 br_po = continuation(
     br, 2, opts_po_cont,
     ShootingProblem(15, probsh, Rodas5P(); parallel = false, update_section_every_step = 1);
-    ampfactor = 1., δp = 0.0051,
+    # ampfactor = 1., δp = 0.0051,
     # verbosity = 3,    plot = true,
-    record_from_solution = (x, p) -> (return (max = getmaximum(p.prob, x, @set par_lur.β = p.p), period = getperiod(p.prob, x, @set par_lur.β = p.p))),
-    # plot_solution = plotSH,
+    record_from_solution = recordPO,
+    plot_solution = plotPO,
     # finalise_solution = (z, tau, step, contResult; prob = nothing, kwargs...) -> begin
     #         BK.haseigenvalues(contResult) && Base.display(contResult.eig[end].eigenvals)
     #         return z.u[end] < 30 && length(contResult.specialpoint) < 3
@@ -130,8 +148,8 @@ show(br_po)
 # plot(br, br_po)
 # plot(br_po, vars=(:param, :period))
 
-@test br_po.specialpoint[1].param ≈ 0.6273246 rtol = 1e-4
-@test br_po.specialpoint[2].param ≈ 0.5417461 rtol = 1e-4
+@test br_po.specialpoint[1].param ≈ 0.63030057 rtol = 1e-4
+@test br_po.specialpoint[2].param ≈ -0.63030476 rtol = 1e-4
 
 # test showing normal form
 for _ind in (1,3)
@@ -146,25 +164,26 @@ 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,
+    # verbosity = 0, plot = false,
     usedeflation = false,
-    ampfactor = .3, δp = -0.005,
-    record_from_solution = (x, p) -> (return (max = BK.getmaximum(p.prob, x, @set par_lur.β = p.p), period = getperiod(p.prob, x, @set par_lur.β = p.p))),
+    ampfactor = .1, δp = -0.005,
+    record_from_solution = recordPO,
     normC = norminf,
     callback_newton = BK.cbMaxNorm(10),
     )
-@test br_po_pd.sol[1].x[end] ≈ 16.956 rtol = 1e-4
 
 # plot(br_po, br_po_pd)
 #######################################
 @info "testLure Poincare"
-opts_po_cont_ps = @set opts_po_cont.newton_options.tol = 1e-7
+opts_po_cont_ps = @set opts_po_cont.newton_options.tol = 1e-9
 @set opts_po_cont_ps.dsmax = 0.0025
 br_po = continuation(br, 2, opts_po_cont_ps,
     PoincareShootingProblem(2, probsh, Rodas4P(); parallel = false, reltol = 1e-6, update_section_every_step = 1, jacobian = BK.AutoDiffDenseAnalytical());
-    ampfactor = 1., δp = 0.0051, #verbosity = 3,plot=true,
+    ampfactor = 1., δp = 0.0051, 
+    # verbosity = 3, plot=true,
     callback_newton = BK.cbMaxNorm(10),
-    record_from_solution = (x, p) -> (return (max = getmaximum(p.prob, x, @set par_lur.β = p.p), period = getperiod(p.prob, x, @set par_lur.β = p.p))),
+    record_from_solution = recordPO,
+    plot_solution = plotPO,
     normC = norminf)
 
 # plot(br_po, br)
@@ -181,13 +200,13 @@ for _ind in (1,)
 end
 
 # aBS from PD
-br_po_pd = BK.continuation(br_po, 1, setproperties(br_po.contparams, detect_bifurcation = 3, max_steps = 1, ds = 0.01, plot_every_step = 1);
-    verbosity = 3, 
-    # plot = true,
-    ampfactor = .3, δp = -0.005,
+br_po_pd = BK.continuation(br_po, 1, setproperties(br_po.contparams, detect_bifurcation = 0, max_steps = 3, ds = -0.01, plot_every_step = 1);
+    # verbosity = 3, plot = true,
+    ampfactor = .1, δp = -0.005,
     normC = norminf,
     callback_newton = BK.cbMaxNorm(10),
-    record_from_solution = (x, p) -> (return (max = getmaximum(p.prob, x, @set par_lur.β = p.p), period = getperiod(p.prob, x, @set par_lur.β = p.p))),
+    record_from_solution = recordPO,
+    plot_solution = plotPO,
     )
 
 # plot(br_po_pd, br_po)