Skip to content

Commit

Permalink
complete update to [email protected]
Browse files Browse the repository at this point in the history
  • Loading branch information
rveltz committed Nov 27, 2024
1 parent c96b3ca commit 8440537
Show file tree
Hide file tree
Showing 5 changed files with 37 additions and 33 deletions.
2 changes: 0 additions & 2 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"

[compat]
BifurcationKit = "^0.3.6, ^0.4.4"
Expand All @@ -20,7 +19,6 @@ ForwardDiff = "^0.10"
Parameters = "^0.12"
RecursiveArrayTools = "^2.3, ^2.4, ^2.8, ^2.9, ^3"
SciMLBase = "^1, ^2"
Setfield = "0.5.0, 0.7.0, 0.8.0, ^1.1"
julia = "1.9"

[extras]
Expand Down
16 changes: 11 additions & 5 deletions src/HomProblemPBC.jl
Original file line number Diff line number Diff line change
Expand Up @@ -410,7 +410,7 @@ function BK.continuation(prob_vf,
T0 = 2/exp(ζ)
f(t) = sech* t) - ζ
pb = BifurcationProblem((t,p) -> [f(t[1])], [T0/ϵ], nothing)
solT = newton(pb, NewtonPar(tol = 1e-8, verbose = false))
solT = BK.solve(pb, Newton(), NewtonPar(tol = 1e-8, verbose = false))
@assert BK.converged(solT) "Newton iteration failed in determining the half return time T"
T = min(solT.u[1], maxT)
printstyled(color = :magenta,"├─ T = ", T, "\n")
Expand All @@ -431,18 +431,24 @@ function BK.continuation(prob_vf,

# define initial guess for newton
prob_vf = re_make(prob_vf, params = pars)
xflow, bvp = initBVPforPBC(bvp, prob_vf, Hom; N = N, T = T, ϵ = ϵ)
xflow, bvp = initBVPforPBC(bvp, prob_vf, Hom; N, T, ϵ)
bvp = BK.set_params_po(bvp, pars)

# define problem for Homoclinic functional
J = BK.jacobian(prob_vf, xsaddle, pars)
𝐇𝐨𝐦 = HomoclinicHyperbolicProblemPBC(bvp, lens1, length(xsaddle), copy(J); ϵ0 = ϵ0hom, ϵ1 = ϵ1hom, T = Thom, freeparams = freeparams, update_every_step = update_every_step, testOrbitFlip = test_orbit_flip,
testInclinationFlip = test_inclination_flip)
𝐇𝐨𝐦 = HomoclinicHyperbolicProblemPBC(bvp, lens1, length(xsaddle), copy(J);
ϵ0 = ϵ0hom,
ϵ1 = ϵ1hom,
T = Thom,
freeparams = freeparams,
update_every_step = update_every_step,
testOrbitFlip = test_orbit_flip,
testInclinationFlip = test_inclination_flip)

@assert BK.getparams(𝐇𝐨𝐦) == pars "Errors with setting the parameters. Please an issue on the website of BifurcationKit."

printstyled("──> convergence to saddle point:\n", color = :magenta)
solsaddle = BK.newton(BifurcationProblem((x,p) -> getF(𝐇𝐨𝐦,x,p), xsaddle, pars), NewtonPar(verbose = true), norm = x->norm(x,Inf))
solsaddle = BK.solve(BifurcationProblem((x,p) -> getF(𝐇𝐨𝐦,x,p), xsaddle, pars), Newton(), NewtonPar(verbose = true), norm = x->norm(x,Inf))
if BK.converged(solsaddle)
xsaddle .= solsaddle.u
end
Expand Down
2 changes: 1 addition & 1 deletion src/StandardShooting.jl
Original file line number Diff line number Diff line change
Expand Up @@ -206,7 +206,7 @@ function generate_hom_problem(sh::ShootingProblem,
ind_saddle = argmin(norm(BK.vf(sh.flow, solpo(t), pars)) for t in time)
xsaddle = solpo(time[ind_saddle])
tsaddle = time[ind_saddle]
solsaddle = BK.newton(BifurcationProblem((x,p) -> BK.vf(sh.flow,x,p),xsaddle,pars), NewtonPar(verbose = true), norm = x->norm(x,Inf))
solsaddle = BK.solve(BifurcationProblem((x,p) -> BK.vf(sh.flow,x,p),xsaddle,pars), Newton(), NewtonPar(verbose = true), norm = x->norm(x,Inf))
if BK.converged(solsaddle)
xsaddle .= solsaddle.u
end
Expand Down
24 changes: 12 additions & 12 deletions test/generateHom.jl
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
# using Revise, Plots, AbbreviatedStackTraces
using Setfield, LinearAlgebra, Test, ForwardDiff
using LinearAlgebra, Test, ForwardDiff
using BifurcationKit, Test
using HclinicBifurcationKit
const BK = BifurcationKit

recordFromSolution(x, p) = (x = x[1], y = x[2])
recordFromSolution(x, p; k...) = (x = x[1], y = x[2])
####################################################################################################
function freire!(dz, u, p, t)
(;ν, β, A₃, B₃, r, ϵ) = p
Expand All @@ -19,7 +19,7 @@ freire(z, p) = freire!(similar(z), z, p, 0)
par_freire == -0.75, β = -0.1, A₃ = 0.328578, B₃ = 0.933578, r = 0.6, ϵ = 0.01)
z0 = [0.7,0.3,0.1]
z0 = zeros(3)
prob = BK.BifurcationProblem(freire, z0, par_freire, (@lens _.β); record_from_solution = recordFromSolution)
prob = BK.BifurcationProblem(freire, z0, par_freire, (@optic _.β); record_from_solution = recordFromSolution)

opts_br = ContinuationPar(p_min = -1.4, p_max = 2.8, ds = 0.001, dsmax = 0.05, n_inversion = 6, detect_bifurcation = 3, max_bisection_steps = 25, nev = 3, max_steps = 2000)
br = continuation(prob, PALC(tangent = Bordered()), opts_br; verbosity = 0,
Expand All @@ -36,7 +36,7 @@ function plotPO(x, p; k...)
end

# record function
function recordPO(x, p)
function recordPO(x, p; k...)
xtt = BK.get_periodic_orbit(p.prob, x, @set par_freire.β = p.p)
period = BK.getperiod(p.prob, x, p.p)
return (period = period, max = maximum(xtt[1,:]), min = minimum(xtt[1,:]))
Expand Down Expand Up @@ -76,10 +76,10 @@ probhom, solh = generate_hom_problem(
update_every_step = 4,
ϵ0 = 1e-9,
ϵ1 = 1e-6,
# freeparams = ((@lens _.T), (@lens _.ϵ1),)
# freeparams = ((@lens _.T), (@lens _.ϵ0)),
freeparams = ((@lens _.ϵ0), (@lens _.ϵ1)),
# freeparams = ((@lens _.T),),
# freeparams = ((@optic _.T), (@optic _.ϵ1),)
# freeparams = ((@optic _.T), (@optic _.ϵ0)),
freeparams = ((@optic _.ϵ0), (@optic _.ϵ1)),
# freeparams = ((@optic _.T),),
)

show(probhom)
Expand Down Expand Up @@ -125,10 +125,10 @@ probhom, solh = generate_hom_problem(
update_every_step = 4,
ϵ0 = 7e-8,
ϵ1 = 8e-8,
# freeparams = ((@lens _.T), (@lens _.ϵ1),)
# freeparams = ((@lens _.T), (@lens _.ϵ0)),
freeparams = ((@lens _.ϵ0), (@lens _.ϵ1)), # WORK BEST
# freeparams = ((@lens _.T),),
# freeparams = ((@optic _.T), (@optic _.ϵ1),)
# freeparams = ((@optic _.T), (@optic _.ϵ0)),
freeparams = ((@optic _.ϵ0), (@optic _.ϵ1)), # WORK BEST
# freeparams = ((@optic _.T),),
)

solh.x[2] .=0
Expand Down
26 changes: 13 additions & 13 deletions test/testHomoclinic.jl
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
# using Revise
# using Plots
using Test
using BifurcationKit, LinearAlgebra, Setfield, ForwardDiff, HclinicBifurcationKit
using BifurcationKit, LinearAlgebra, ForwardDiff, HclinicBifurcationKit
const BK = BifurcationKit

####################################################################################################
Expand All @@ -13,21 +13,21 @@ function Fbt!(dx, x, p, t=0)
end
Fbt(x,p) = Fbt!(similar(x),x,p)
par = (β1 = -0.01, β2 = -0.1, a = 1., b = -1.)
prob = BK.BifurcationProblem(Fbt, [0.01, 0.01], par, (@lens _.β1))
prob = BK.BifurcationProblem(Fbt, [0.01, 0.01], par, (@optic _.β1))
opt_newton = NewtonPar(tol = 1e-9, max_iterations = 40, verbose = false)
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.01, ds = 0.01, p_max = 0.5, p_min = -0.5, detect_bifurcation = 3, nev = 2, newton_options = opt_newton, max_steps = 100, n_inversion = 8, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9, save_sol_every_step = 1)

br = continuation(prob, PALC(), opts_br; bothside = true, verbosity = 0)

sn_codim2 = continuation(br, 2, (@lens _.β2), ContinuationPar(opts_br, detect_bifurcation = 1, save_sol_every_step = 1, max_steps = 40) ;
sn_codim2 = continuation(br, 2, (@optic _.β2), ContinuationPar(opts_br, detect_bifurcation = 1, save_sol_every_step = 1, max_steps = 40) ;
detect_codim2_bifurcation = 2,
update_minaug_every_step = 1,
)
@test sn_codim2.specialpoint[1].type == :bt
@test sn_codim2.specialpoint[1].param 0 atol = 1e-6
@test length(unique(sn_codim2.BT)) == length(sn_codim2)

hopf_codim2 = continuation(br, 3, (@lens _.β2), ContinuationPar(opts_br, detect_bifurcation = 1, save_sol_every_step = 1, max_steps = 40, max_bisection_steps = 25) ; plot = false, verbosity = 0,
hopf_codim2 = continuation(br, 3, (@optic _.β2), ContinuationPar(opts_br, detect_bifurcation = 1, save_sol_every_step = 1, max_steps = 40, max_bisection_steps = 25) ; plot = false, verbosity = 0,
detect_codim2_bifurcation = 2,
update_minaug_every_step = 1,
bothside = true,
Expand Down Expand Up @@ -102,11 +102,11 @@ end
br_hom = continuation(prob, btpt,
PeriodicOrbitOCollProblem(20, 4; meshadapt = true),
PALC(tangent = Bordered()),
setproperties(optc_hom, max_steps = 200, save_sol_every_step = 1, dsmax = 3e-1, plot_every_step = 3, detect_event=2);
ContinuationPar(optc_hom, max_steps = 200, save_sol_every_step = 1, dsmax = 3e-1, plot_every_step = 3, detect_event=2);
amplitude = 5e-3, ϵ0 = 1e-3,
# freeparams = ((@lens _.T), (@lens _.ϵ0)),
freeparams = ((@lens _.ϵ0), (@lens _.ϵ1)), # WORK BEST
# freeparams = ((@lens _.T),),
# freeparams = ((@optic _.T), (@optic _.ϵ0)),
freeparams = ((@optic _.ϵ0), (@optic _.ϵ1)), # WORK BEST
# freeparams = ((@optic _.T),),
verbosity = 0, plot = false,
callback_newton = BK.cbMaxNorm(1e1),
plot_solution = plotHom,
Expand All @@ -131,12 +131,12 @@ br_hom = continuation(prob, btpt,
# PALC(tangent = Bordered()),
PALC(),
# MoorePenrose(),
setproperties(optc_hom, max_steps = 10, save_sol_every_step = 1, dsmax = 9e-1, plot_every_step = 1, detect_event=0);
ContinuationPar(optc_hom, max_steps = 10, save_sol_every_step = 1, dsmax = 9e-1, plot_every_step = 1, detect_event=0);
amplitude = 5e-3, ϵ0 = 2e-4,
# freeparams = ((@lens _.T), (@lens _.ϵ0)),
freeparams = ((@lens _.T), (@lens _.ϵ1)),
# freeparams = ((@lens _.ϵ0), (@lens _.ϵ1)), # WORK BEST
# freeparams = ((@lens _.T),),
# freeparams = ((@optic _.T), (@optic _.ϵ0)),
freeparams = ((@optic _.T), (@optic _.ϵ1)),
# freeparams = ((@optic _.ϵ0), (@optic _.ϵ1)), # WORK BEST
# freeparams = ((@optic _.T),),
verbosity = 0, plot = false,
update_every_step = 1, # ne marche pas sinon. pb n est pas avec updatesection
plot_solution = plotHom,
Expand Down

0 comments on commit 8440537

Please sign in to comment.