rveltz
commented
3 years ago Hi,
As I told you on discourse, this is a bit "border line" for my package which targets PDEs. However, I am currently improving a lot the interface so ODE can be studied too in regimes where this is really difficult for PDE (like homoclinic orbits,...). Indeed, with these improvements, you will be able to use much larger M than before (for Shooting methods and FD ones). However, it brings another problem: my computation of the Floquet coefficients become unstable for large Ms (but I am working on it). So using detectBifurcation = 2 is likely to give spurious bifurcations for now.
I will soon announce these improvements and document them with an ODE application.
For the FD method, you can do:
# automatic branch switching from Hopf point
opt_po = NewtonPar(tol = 1e-5, verbose = true, maxIter = 15)
opts_po_cont = ContinuationPar(dsmin = 0.0001, dsmax = 0.1, ds = 0.001, pMax = 5.5, pMin=-2.1, maxSteps = 200, newtonOptions = opt_po, saveSolEveryStep = 2, plotEveryStep = 1, nev = 11, precisionStability = 1e-4, detectBifurcation = 0, dsminBisection = 1e-4, maxBisectionSteps = 10);
# number of time slices for the periodic orbit
M = 251
br_po, = continuation(
# arguments for branch switching from the first
# Hopf bifurcation point
jet..., br, 1,
# arguments for continuation
setproperties(opts_po_cont; ds = -0.001, dsmax = 0.1),
PeriodicOrbitTrapProblem(M = M);
# to find non trivial periodic orbits
usedeflation = true,
# to update the section during continuation
updateSectionEveryStep = 1,
# OPTIONAL parameters
# we want to jump on the new branch at phopf + δp
# ampfactor is a factor to increase the amplitude of the guess
δp = 0.003, ampfactor = 1.0,
# specific method for computing the Jacobian
linearPO = :Dense,
# regular options for continuation
verbosity = 3, plot = true,
printSolution = (x, p) -> begin
xtt=reshape(x[1:end-1],2,M);
return (max = maximum(xtt[1,:]), min = minimum(xtt[1,:]), period = x[end])
end,
plotSolution = (x,p; k...) -> begin
xtt = BK.getTrajectory(p.prob, x, p.p)
plot!(xtt.t, xtt.u[1,:]; label = "x", k...)
plot!(xtt.t, xtt.u[2,:]; label = "y", k...)
plot!(br, legend = false; subplot = 1)
end,
normC = norminf)It gives:
For the Shooting method, it is very similar:
probsh = ODEProblem(Fitz!, copy(sol0), (0., 1000.), par_Fitz; atol = 1e-10, rtol = 1e-9)
br_po, = continuation(
# arguments for branch switching from the first
# Hopf bifurcation point
jet..., br, 1,
# arguments for continuation
setproperties(opts_po_cont; ds = -0.001, dsmax = 0.31, detectBifurcation = 0),
ShootingProblem(25, par_Fitz, probsh, Rodas4(), parallel = true);
usedeflation = false,
# OPTIONAL parameters
# we want to jump on the new branch at phopf + δp
# ampfactor is a factor to increase the amplitude of the guess
δp = 0.001, ampfactor = 1.0,
# specific method for computing the Jacobian
linearPO = :autodiffDense,
updateSectionEveryStep = 2,
# regular options for continuation
verbosity = 3, plot = true,
printSolution = (x, p) -> begin
(max = getMaximum(p.prob, x, @set par_Fitz.I = p.p), period = getPeriod(p.prob, x, @set par_Fitz.I = p.p))
end,
plotSolution = (x,p; k...) -> begin
xtt = BK.getTrajectory(p.prob, x, @set par_Fitz.I = p.p)
plot!(xtt; k...)
plot!(br, legend = false; subplot = 1)
end,
normC = norminf)But it does not work that well, I think this is because of the section.
Best.
daybrown
It's really sensitive to the parameters (like δp).
I expect to get the following result in AUTO. Though I already have the result, I am looking for a pure Julia package to handle this kind of problem.
My question is how do I modify the code to get a similar result?