Skip to content

Commit

Permalink
change coefficient a in PD normal form using Iooss method.
Browse files Browse the repository at this point in the history 
add prm option for NS continuation
add warning when using Bordered predictor with Fold continuation
  • Loading branch information
@rveltz
rveltz committed Jan 7, 2024
1 parent 9f79b7d commit 880b276
Show file tree
Hide file tree
Showing 11 changed files with 130 additions and 87 deletions.
105 changes: 78 additions & 27 deletions examples/COModel.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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")
49 changes: 15 additions & 34 deletions examples/TMModel.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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...)
Expand All @@ -52,79 +47,65 @@ 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.),
)

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...,
plot_solution = (x, p; k...) -> begin
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,
Expand Down
20 changes: 10 additions & 10 deletions examples/pd-1d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
6 changes: 3 additions & 3 deletions src/LinearBorderSolver.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
4 changes: 4 additions & 0 deletions src/codim2/MinAugFold.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
12 changes: 8 additions & 4 deletions src/continuation/MoorePenrose.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down
8 changes: 4 additions & 4 deletions src/events/EventDetection.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
5 changes: 3 additions & 2 deletions src/periodicorbit/NormalForms.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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₀ + δ)) .-
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit 880b276

Please sign in to comment.