diff --git a/examples/COModel.jl b/examples/COModel.jl index 63c00561..84c1dfd3 100644 --- a/examples/COModel.jl +++ b/examples/COModel.jl @@ -21,43 +21,101 @@ z0 = [0.001137, 0.891483, 0.062345] prob = BifurcationProblem(COm!, z0, par_com, (@lens _.q2); record_from_solution = (x, p) -> (x = x[1], y = x[2], s = x[3])) -opts_br = ContinuationPar(dsmax = 0.05, p_min = 0.5, p_max = 2.0, n_inversion = 6, detect_bifurcation = 3, max_bisection_steps = 25, nev = 3) +opts_br = ContinuationPar(dsmax = 0.015, dsmin=1e-4, ds=1e-4, p_min = 0.5, p_max = 2.0, n_inversion = 6, detect_bifurcation = 3, nev = 3) br = @time continuation(prob, PALC(), opts_br; - # plot = false, verbosity = 0, + plot = true, verbosity = 0, normC = norminf, bothside = true) -show(br) -plot(br, plotfold=false, markersize=4, legend=:topright, ylims=(0,0.16)) +BK.plot(br, plotfold=false, )#markersize=4, legend=:topright, ylims=(0,0.16)) #################################################################################################### -@set! opts_br.newton_options.verbose = false -@set! opts_br.newton_options.max_iterations = 10 -opts_br = @set opts_br.newton_options.tol = 1e-12 +# periodic orbits +function plotSolution(x, p; k...) + xtt = BK.get_periodic_orbit(p.prob, x, p.p) + plot!(xtt.t, xtt[1,:]; label = "", k...) + plot!(xtt.t, xtt[2,:]; label = "", k...) + plot!(xtt.t, xtt[3,:]; label = "", k...) + plot!(xtt.t, xtt[1,:]; label = "", marker =:d, markersize = 1.5, k...) + plot!(br; subplot = 1, putspecialptlegend = false, xlims = (1.02, 1.07)) +end -sn = newton(br, 3; options = opts_br.newton_options, bdlinsolver = MatrixBLS()) +function plotSolution(ax, x, p; ax1 = nothing, k...) + @info "plotsol Makie" + xtt = BK.get_periodic_orbit(p.prob, x, p.p) + lines!(ax1, br; putspecialptlegend = false) + lines!(ax, xtt.t, xtt[1,:]; k...) + lines!(ax, xtt.t, xtt[2,:]; k...) + lines!(ax, xtt.t, xtt[3,:]; k...) + scatter!(ax, xtt.t, xtt[1,:]; markersize = 1.5, k...) +end -hp = newton(br, 2; options = NewtonPar( opts_br.newton_options; max_iterations = 10),start_with_eigen=true) +args_po = ( record_from_solution = (x, p) -> begin + xtt = BK.get_periodic_orbit(p.prob, x, @set par_com.q2 = p.p) + return (max = maximum(xtt[1,:]), + min = minimum(xtt[1,:]), + period = getperiod(p.prob, x, @set par_com.q2 = p.p)) + end, + plot_solution = plotSolution, + normC = norminf) + +opts_po_cont = ContinuationPar(opts_br, dsmax = 1.95, ds= 2e-2, dsmin = 1e-6, p_max = 5., p_min=-5., +max_steps = 300, detect_bifurcation = 0, plot_every_step = 10) + +# @set! opts_po_cont.newton_options.verbose = false +# @set! opts_po_cont.newton_options.tol = 1e-11 +# @set! opts_po_cont.newton_options.max_iterations = 10 +# @set! opts_po_cont.newton_options.linesearch = true + +using DifferentialEquations +prob_ode = ODEProblem(COm!, copy(z0), (0., 1000.), par_com; abstol = 1e-11, reltol = 1e-9) + +function callbackCO(state; fromNewton = false, kwargs...) + # check that the solution is not too far + δ0 = 1e-3 + z0 = get(state, :z0, nothing) + p = get(state, :p, nothing) + if state.residual > 1 + @error "Reject Newton step, res too big!!" + return false + end + # abort of the δp is too large + if ~fromNewton && ~isnothing(z0) + # @show abs(p - z0.p) + return abs(p - z0.p) <= 2e-3 + end + return true +end -hpnf = get_normal_form(br, 2) +brpo = @time continuation(br, 5, opts_po_cont, + PeriodicOrbitOCollProblem(50, 4 ; update_section_every_step = 1, jacobian = BK.DenseAnalytical(), meshadapt = true, K = 5000, verbose_mesh_adapt = true); + # ShootingProblem(25, prob_ode, TaylorMethod(25); parallel = true; update_section_every_step = 1, jacobian = BK.AutoDiffDense()); + verbosity = 3, plot = true, + normC = norminf, + # alg = PALC(tangent = Bordered()), + # alg = PALC(), + # alg = MoorePenrose(tangent=PALC(tangent = Bordered()), method = BK.direct), + δp = 0.0005, + # linear_algo = DefaultLS(), + callback_newton = callbackCO, + args_po... + ) +BK.plot(br, brpo, putspecialptlegend = false, branchlabel = ["eq","max"]) +xlims!((1.037, 1.055)) +scatter!(br) +plot!(brpo.param, brpo.min, label = "min", xlims = (1.037, 1.055)) +#################################################################################################### sn_codim2 = continuation(br, 3, (@lens _.k), ContinuationPar(opts_br, p_max = 3.2, p_min = 0., detect_bifurcation = 0, dsmin=1e-5, ds = -0.001, dsmax = 0.05, n_inversion = 6, detect_event = 2, detect_fold = false) ; plot = true, verbosity = 3, normC = norminf, update_minaug_every_step = 1, start_with_eigen = true, # record_from_solution = (u,p; kw...) -> (x = u.u[1] ), - bothside=true, + bothside = true, bdlinsolver = MatrixBLS() ) -using Test -@test sn_codim2.specialpoint[2].printsol.k ≈ 0.971397 rtol = 1e-4 -@test sn_codim2.specialpoint[2].printsol.q2 ≈ 1.417628 rtol = 1e-4 -@test sn_codim2.specialpoint[4].printsol.k ≈ 0.722339 rtol = 1e-4 -@test sn_codim2.specialpoint[4].printsol.q2 ≈ 1.161199 rtol = 1e-4 - BK.plot(sn_codim2)#, real.(sn_codim2.BT), ylims = (-1,1), xlims=(0,2)) - BK.plot(sn_codim2, vars=(:q2, :x), branchlabel = "Fold", plotstability = false);plot!(br,xlims=(0.8,1.8)) hp_codim2 = continuation((@set br.alg.tangent = Bordered()), 2, (@lens _.k), ContinuationPar(opts_br, p_min = 0., p_max = 2.8, detect_bifurcation = 0, ds = -0.0001, dsmax = 0.08, dsmin = 1e-4, n_inversion = 6, detect_event = 2, detect_loop = true, max_steps = 50, detect_fold=false) ; plot = true, @@ -70,18 +128,11 @@ hp_codim2 = continuation((@set br.alg.tangent = Bordered()), 2, (@lens _.k), Con bothside = true, bdlinsolver = MatrixBLS()) -@test hp_codim2.branch[6].l1 |> real ≈ 33.15920 rtol = 1e-1 -@test hp_codim2.specialpoint[3].printsol.k ≈ 0.305879 rtol = 1e-3 -@test hp_codim2.specialpoint[3].printsol.q2 ≈ 0.924255 rtol = 1e-3 -@test hp_codim2.specialpoint[4].printsol.k ≈ 0.23248736 rtol = 1e-4 -@test hp_codim2.specialpoint[4].printsol.q2 ≈ 0.8913189828755895 rtol = 1e-4 - BK.plot(sn_codim2, vars=(:q2, :x), branchlabel = "Fold", plotcirclesbif = true) plot!(hp_codim2, vars=(:q2, :x), branchlabel = "Hopf",plotcirclesbif = true) plot!(br,xlims=(0.6,1.5)) -plot(sn_codim2, vars=(:k, :q2), branchlabel = "Fold") +BK.plot(sn_codim2, vars=(:k, :q2), branchlabel = "Fold") plot!(hp_codim2, vars=(:k, :q2), branchlabel = "Hopf",) -plot(hp_codim2, vars=(:q2, :x), branchlabel = "Hopf") -#################################################################################################### +BK.plot(hp_codim2, vars=(:q2, :x), branchlabel = "Hopf") diff --git a/examples/TMModel.jl b/examples/TMModel.jl index 9ed22609..0a156c71 100644 --- a/examples/TMModel.jl +++ b/examples/TMModel.jl @@ -19,18 +19,13 @@ par_tm = (α = 1.5, τ = 0.013, J = 3.07, E0 = -2.0, τD = 0.200, U0 = 0.3, τF z0 = [0.238616, 0.982747, 0.367876 ] prob = BifurcationProblem(TMvf!, z0, par_tm, (@lens _.E0); record_from_solution = (x, p) -> (E = x[1], x = x[2], u = x[3]),) -opts_br = ContinuationPar(p_min = -10.0, p_max = -0.9, ds = 0.04, dsmax = 0.125, n_inversion = 8, detect_bifurcation = 3, max_bisection_steps = 25, nev = 3) -br = continuation(prob, PALC(tangent = Bordered()), opts_br; plot = true, normC = norminf) +opts_br = ContinuationPar(p_min = -10.0, p_max = -0.9, ds = 0.04, dsmax = 0.1, n_inversion = 8, detect_bifurcation = 3, max_bisection_steps = 25, nev = 3) +br = continuation(prob, PALC(), opts_br; plot = true, normC = norminf) BK.plot(br, plotfold=false) #################################################################################################### -hopfpt = get_normal_form(br, 4) - -# newton parameters -optn_po = NewtonPar(verbose = false, tol = 1e-8, max_iterations = 9) - # continuation parameters -opts_po_cont = ContinuationPar(dsmax = 0.1, ds= 0.0005, dsmin = 1e-4, p_max = 0., p_min=-2., max_steps = 120, newton_options = (@set optn_po.tol = 1e-8), nev = 3, tol_stability = 1e-6, detect_bifurcation = 2, plot_every_step = 20, save_sol_every_step=1) +opts_po_cont = ContinuationPar(dsmax = 0.1, ds= 0.0005, dsmin = 1e-4, p_max = 0., p_min=-2., max_steps = 20, newton_options = NewtonPar(tol = 1e-10, max_iterations = 9), nev = 3, tol_stability = 1e-6, detect_bifurcation = 2, plot_every_step = 20) # arguments for periodic orbits function plotSolution(x, p; k...) @@ -52,8 +47,8 @@ args_po = ( record_from_solution = (x, p) -> begin normC = norminf) br_potrap = continuation(br, 4, opts_po_cont, - PeriodicOrbitTrapProblem(M = 250, jacobian = :Dense, update_section_every_step = 0); - verbosity = 2, plot = false, + PeriodicOrbitTrapProblem(M = 250, jacobian = :Dense, update_section_every_step = 1); + verbosity = 2, plot = true, args_po..., callback_newton = BK.cbMaxNorm(10.), ) @@ -61,18 +56,12 @@ br_potrap = continuation(br, 4, opts_po_cont, plot(br, br_potrap, markersize = 3) plot!(br_potrap.param, br_potrap.min, label = "") #################################################################################################### -# based on collocation -hopfpt = get_normal_form(br, 4) - -# newton parameters -optn_po = NewtonPar(verbose = false, tol = 1e-8, max_iterations = 10) - # continuation parameters -opts_po_cont = ContinuationPar(dsmax = 0.1, ds= 0.0001, dsmin = 1e-4, p_max = 0., p_min=-5., max_steps = 110, newton_options = (@set optn_po.tol = 1e-7), nev = 3, tol_stability = 1e-5, detect_bifurcation = 2, plot_every_step = 40, save_sol_every_step=1) +opts_po_cont = ContinuationPar(dsmax = 0.1, ds= 0.0001, dsmin = 1e-4, p_max = 0., p_min=-5., max_steps = 110, newton_options = NewtonPar(verbose = false, tol = 1e-10, max_iterations = 10), nev = 3, tol_stability = 1e-5, detect_bifurcation = 2, plot_every_step = 40, save_sol_every_step=1) br_pocoll = @time continuation( br, 4, opts_po_cont, - PeriodicOrbitOCollProblem(30, 5, update_section_every_step = 0, meshadapt = true); + PeriodicOrbitOCollProblem(30, 5, update_section_every_step = 1, meshadapt = true); verbosity = 2, plot = true, args_po..., @@ -80,51 +69,43 @@ br_pocoll = @time continuation( xtt = BK.get_periodic_orbit(p.prob, x, p.p) plot!(xtt.t, xtt[1,:]; label = "", marker =:d, markersize = 1.5, k...) plot!(br; subplot = 1, putspecialptlegend = false) - end, - callback_newton = BK.cbMaxNorm(1000.), ) plot(br, br_pocoll, markersize = 3) - # plot!(br_pocoll.param, br_pocoll.min, label = "") - # plot!(br, br_potrap, markersize = 3) - # plot!(br_potrap.param, br_potrap.min, label = "", marker = :d) #################################################################################################### # idem with Standard shooting -using DifferentialEquations#, TaylorIntegration +using DifferentialEquations # this is the ODEProblem used with `DiffEqBase.solve` -probsh = ODEProblem(TMvf!, copy(z0), (0., 1000.), par_tm; abstol = 1e-10, reltol = 1e-9) +prob_ode = ODEProblem(TMvf!, copy(z0), (0., 1000.), par_tm; abstol = 1e-11, reltol = 1e-9) -opts_po_cont = ContinuationPar(dsmax = 0.09, ds= -0.0001, dsmin = 1e-4, p_max = 0., p_min=-5., max_steps = 120, newton_options = NewtonPar(optn_po; tol = 1e-6, max_iterations = 7), nev = 3, tol_stability = 1e-8, detect_bifurcation = 2, plot_every_step = 10, save_sol_every_step=1) +opts_po_cont = ContinuationPar(dsmax = 0.09, ds= -0.0001, dsmin = 1e-4, p_max = 0., p_min=-5., max_steps = 120, newton_options = NewtonPar(tol = 1e-6, max_iterations = 7), nev = 3, tol_stability = 1e-8, detect_bifurcation = 2, plot_every_step = 10, save_sol_every_step=1) br_posh = @time continuation( br, 4, # arguments for continuation opts_po_cont, # this is where we tell that we want Standard Shooting - ShootingProblem(15, probsh, Rodas4P(), parallel = true, update_section_every_step = 1, jacobian = BK.AutoDiffDense(),); - # ShootingProblem(15, probsh, TaylorMethod(15), parallel = false); - ampfactor = 1.0, δp = 0.0005, - usedeflation = true, + ShootingProblem(15, prob_ode, Rodas4P(), parallel = true, update_section_every_step = 1, jacobian = BK.AutoDiffDense(),); linear_algo = MatrixBLS(), - verbosity = 2, plot = true, + verbosity = 2, + plot = true, args_po..., ) plot(br_posh, br, markersize=3) - # plot(br, br_potrap, br_posh, markersize=3) #################################################################################################### # idem with Poincaré shooting -opts_po_cont = ContinuationPar(dsmax = 0.02, ds= 0.001, dsmin = 1e-6, p_max = 0., p_min=-5., max_steps = 50, newton_options = NewtonPar(optn_po;tol = 1e-9, max_iterations=15), nev = 3, tol_stability = 1e-8, detect_bifurcation = 2, plot_every_step = 5, save_sol_every_step = 1) +opts_po_cont = ContinuationPar(dsmax = 0.02, ds= 0.001, p_max = 0., p_min=-5., max_steps = 50, newton_options = NewtonPar(tol = 1e-9, max_iterations=15), nev = 3, tol_stability = 1e-6, detect_bifurcation = 2, plot_every_step = 5) br_popsh = @time continuation( br, 4, # arguments for continuation opts_po_cont, # this is where we tell that we want Poincaré Shooting - PoincareShootingProblem(5, probsh, Rodas4P(); parallel = true, update_section_every_step = 1, jacobian = BK.AutoDiffDenseAnalytical(), abstol = 1e-10, reltol = 1e-9); + PoincareShootingProblem(5, prob_ode, Rodas4P(); parallel = true, update_section_every_step = 1, jacobian = BK.AutoDiffDenseAnalytical()); alg = PALC(tangent = Bordered()), ampfactor = 1.0, δp = 0.005, # usedeflation = true, diff --git a/examples/pd-1d.jl b/examples/pd-1d.jl index be5a748a..9ca18f79 100644 --- a/examples/pd-1d.jl +++ b/examples/pd-1d.jl @@ -62,19 +62,19 @@ using SparseDiffTools #################################################################################################### N = 100 - n = 2N - lx = 3pi /2 - X = LinRange(-lx,lx, N) - Δ = Laplacian(N, lx, :Neumann) - D = 0.08 - par_br = (η = 1.0, a = -1., b = -3/2., H = 3.0, D = D, C = -0.6, Δ = blockdiag(D*Δ, Δ)) +n = 2N +lx = 3pi /2 +X = LinRange(-lx,lx, N) +Δ = Laplacian(N, lx, :Neumann) +D = 0.08 +par_br = (η = 1.0, a = -1., b = -3/2., H = 3.0, D = D, C = -0.6, Δ = blockdiag(D*Δ, Δ)) - u0 = cos.(2X) - solc0 = vcat(u0, u0) +u0 = cos.(2X) +solc0 = vcat(u0, u0) probBif = BifurcationProblem(Fbr, solc0, par_br, (@lens _.C) ;J = Jbr, - record_from_solution = (x, p) -> norm(x, Inf), - plot_solution = (x, p; kwargs...) -> plot!(x[1:end÷2];label="",ylabel ="u", kwargs...)) + record_from_solution = (x, p) -> norm(x, Inf), + plot_solution = (x, p; kwargs...) -> plot!(x[1:end÷2];label="",ylabel ="u", kwargs...)) #################################################################################################### # eigls = DefaultEig() eigls = EigArpack(0.5, :LM) diff --git a/src/LinearBorderSolver.jl b/src/LinearBorderSolver.jl index 1e854cfb..022bd061 100644 --- a/src/LinearBorderSolver.jl +++ b/src/LinearBorderSolver.jl @@ -30,7 +30,7 @@ $(TYPEDFIELDS) """ @with_kw struct BorderingBLS{S <: Union{AbstractLinearSolver, Nothing}, Ttol} <: AbstractBorderedLinearSolver - "Linear solver used for the Bordering method." + "Linear solver for the Bordering method." solver::S = nothing "Tolerance for checking precision" @@ -167,7 +167,7 @@ function (lbs::BorderingBLS)(::Val{:Block}, J, b::NTuple{M, AbstractVector}, c:: u2 = δd \ (rhsb - cx1) # TODO USE mul! - u1 = x1 - x2_mat * u2 + u1 = x1 - x2_mat * u2 return u1, u2, cv, (its...) end @@ -344,7 +344,7 @@ This struct is used to provide the bordered linear solver based a matrix free o $(TYPEDFIELDS) """ struct MatrixFreeBLS{S <: Union{AbstractLinearSolver, Nothing}} <: AbstractBorderedLinearSolver - "Linear solver used to solve the extended linear system" + "Linear solver for solving the extended linear system" solver::S "What is the structure used to hold `(x, p)`. If `true`, this is achieved using `BorderedArray`. If `false`, a `Vector` is used." use_bordered_array::Bool diff --git a/src/codim2/MinAugFold.jl b/src/codim2/MinAugFold.jl index 79132e78..7a3b959b 100644 --- a/src/codim2/MinAugFold.jl +++ b/src/codim2/MinAugFold.jl @@ -331,6 +331,10 @@ function continuation_fold(prob, alg::AbstractContinuationAlgorithm, @assert lens1 != lens2 "Please choose 2 different parameters. You only passed $lens1" @assert lens1 == getlens(prob) + if alg isa PALC && alg.tangent isa Bordered + @warn "You selected the PALC continuation algorithm with Bordered predictor. The jacobian being singular on Fold points, this could lead to bad prediction and converge. If you have issues, try a different tangent predictor like Secant for example." + end + # options for the Newton Solver inherited from the ones the user provided options_newton = options_cont.newton_options diff --git a/src/continuation/MoorePenrose.jl b/src/continuation/MoorePenrose.jl index 71607a6e..7310cdeb 100644 --- a/src/continuation/MoorePenrose.jl +++ b/src/continuation/MoorePenrose.jl @@ -17,9 +17,9 @@ Additional information is available on the [website](https://bifurcationkit.gith $(TYPEDFIELDS) """ struct MoorePenrose{T, Tls <: AbstractLinearSolver} <: AbstractContinuationAlgorithm - "Tangent predictor, example `PALC()`" + "Tangent predictor, for example `PALC()`" tangent::T - "Use a direct linear solver. Can be BifurcationKit.direct, BifurcationKit.pInv or BifurcationKit.iterative" + "Moore Penrose linear solver. Can be BifurcationKit.direct, BifurcationKit.pInv or BifurcationKit.iterative" method::MoorePenroseLS "(Bordered) linear solver" ls::Tls @@ -32,7 +32,9 @@ internal_adaptation!(alg::MoorePenrose, swch::Bool) = internal_adaptation!(alg.t """ $(SIGNATURES) """ -function MoorePenrose(;tangent = PALC(), method = direct, ls = nothing) +function MoorePenrose(;tangent = PALC(), + method = direct, + ls = nothing) if ~(method == iterative) ls = isnothing(ls) ? DefaultLS() : ls else @@ -55,7 +57,9 @@ function Base.empty!(alg::MoorePenrose) alg end -function update(alg0::MoorePenrose, contParams::ContinuationPar, linearAlgo) +function update(alg0::MoorePenrose, + contParams::ContinuationPar, + linearAlgo) tgt = update(alg0.tangent, contParams, linearAlgo) alg = @set alg0.tangent = tgt if isnothing(linearAlgo) diff --git a/src/events/EventDetection.jl b/src/events/EventDetection.jl index 9085ac5b..d8484bd0 100644 --- a/src/events/EventDetection.jl +++ b/src/events/EventDetection.jl @@ -174,11 +174,11 @@ function locate_event!(event::AbstractEvent, iter, _state, verbose::Bool = true) if verbose printstyled(color=:red, "────> Found at p = ", getp(state), " ∈ $interval, \n\t\t\t δn = ", abs.(2 .* nsigns[end] .- n1 .- n2), ", from p = ", getp(_state), "\n") - printstyled(color=:blue, "─"^40*"\n────> Stopping reason:\n──────> isnothing(next) = ", isnothing(next), - "\n──────> |ds| < dsmin_bisection = ", abs(state.ds) < contParams.dsmin_bisection, + printstyled(color=:blue, "─"^40*"\n────> Stopping reason:\n──────> isnothing(next) = ", isnothing(next), + "\n──────> |ds| < dsmin_bisection = ", abs(state.ds) < contParams.dsmin_bisection, "\n──────> step >= max_bisection_steps = ", state.step >= contParams.max_bisection_steps, - "\n──────> n_inversion >= n_inversion = ", n_inversion >= contParams.n_inversion, - "\n──────> eventlocated = ", eventlocated == true, "\n") + "\n──────> n_inversion >= n_inversion = ", n_inversion >= contParams.n_inversion, + "\n──────> eventlocated = ", eventlocated == true, "\n") end diff --git a/src/periodicorbit/NormalForms.jl b/src/periodicorbit/NormalForms.jl index 812dadca..ff3af898 100644 --- a/src/periodicorbit/NormalForms.jl +++ b/src/periodicorbit/NormalForms.jl @@ -265,13 +265,13 @@ function period_doubling_normal_form(pbwrap::WrapPOColl, T = getperiod(coll, pd.x0, par) lens = getlens(coll) δ = getdelta(coll) + # identity matrix for collocation problem Icoll = analytical_jacobian(coll, pd.x0, par; ρD = 0, ρF = 0, ρI = -1/T) Icoll[:,end] .=0; Icoll[end,:] .=0 Icoll[end-N:end-1, 1:N] .= 0 Icoll[end-N:end-1, end-N:end-1] .= 0 - F(u, p) = residual(coll.prob_vf, u, p) # dₚF(u, p) = ForwardDiff.derivative(z -> residual(coll.prob_vf, u, set(p, lens, z)), get(par, lens)) dₚF(u, p) = (residual(coll.prob_vf, u, set(p, lens, p₀ + δ)) .- @@ -464,7 +464,8 @@ function period_doubling_normal_form(pbwrap::WrapPOColl, c₁₁ = ∫(v₁★ₛ, rhsₛ) - a₀₁ * ∫(v₁★ₛ, Aₛ) c₁₁ *= 2 - nf = (a = a₁, b3 = c, h₂ₛ, ψ₁★ₛ, v₁ₛ, a₀₁, c₁₁) # keep b3 for PD-codim 2 + # we want the parameter a, not the rescaled a₁ + nf = (a = a₁/T, b3 = c, h₂ₛ, ψ₁★ₛ, v₁ₛ, a₀₁, c₁₁) # keep b3 for PD-codim 2 newpd = @set pd.nf = nf @debug "[PD-NF-Iooss]" a₁ c if real(c) < 0 diff --git a/src/periodicorbit/PeriodicOrbits.jl b/src/periodicorbit/PeriodicOrbits.jl index 278751e1..b84d4437 100644 --- a/src/periodicorbit/PeriodicOrbits.jl +++ b/src/periodicorbit/PeriodicOrbits.jl @@ -262,7 +262,8 @@ Similar to [`continuation`](@ref) except that `probPO` is either a [`ShootingPro - `eigsolver` specify an eigen solver for the computation of the Floquet exponents, defaults to `FloquetQaD` $DocStrjacobianPOSh """ -function continuation(probPO::AbstractShootingProblem, orbitguess, +function continuation(probPO::AbstractShootingProblem, + orbitguess, alg::AbstractContinuationAlgorithm, contParams::ContinuationPar, linear_algo::AbstractBorderedLinearSolver; diff --git a/src/periodicorbit/codim2/MinAugNS.jl b/src/periodicorbit/codim2/MinAugNS.jl index 415054fc..c6155a23 100644 --- a/src/periodicorbit/codim2/MinAugNS.jl +++ b/src/periodicorbit/codim2/MinAugNS.jl @@ -227,6 +227,7 @@ function continuation_ns(prob, alg::AbstractContinuationAlgorithm, kind = NSCont(), usehessian = true, plot_solution = BifurcationKit.plot_solution(prob), + prm = false, kwargs...) where {𝒯b, vectype} @assert lens1 != lens2 "Please choose 2 different parameters. You only passed $lens1" @assert lens1 == getlens(prob) diff --git a/src/periodicorbit/codim2/codim2.jl b/src/periodicorbit/codim2/codim2.jl index a2484283..adcdc78f 100644 --- a/src/periodicorbit/codim2/codim2.jl +++ b/src/periodicorbit/codim2/codim2.jl @@ -96,11 +96,11 @@ end function correct_bifurcation(contres::ContResult{<: Union{FoldPeriodicOrbitCont, PDPeriodicOrbitCont, NSPeriodicOrbitCont}}) if contres.prob.prob isa FoldProblemMinimallyAugmented - conversion = Dict(:bp => :R1, :hopf => :foldNS, :fold => :cusp, :nd => :nd, :pd => :foldpd, :bt => :R1) + conversion = Dict(:bp => :R1, :hopf => :foldNS, :fold => :cusp, :nd => :nd, :pd => :foldpd, :bt => :R1, :zh => :R1) elseif contres.prob.prob isa PeriodDoublingProblemMinimallyAugmented conversion = Dict(:bp => :foldFlip, :hopf => :pdNS, :pd => :R2) elseif contres.prob.prob isa NeimarkSackerProblemMinimallyAugmented - conversion = Dict(:bp => :foldNS, :hopf => :nsns, :pd => :pdNS,) + conversion = Dict(:bp => :foldNS, :hopf => :nsns, :pd => :pdNS, :R1ch => :R1, :fold => :R1) else throw("Error! this should not occur. Please open an issue on the website of BifurcationKit.jl") end