# Finding periodic orbits with BifurcationKit.jl

Hi everyone! Currently I’m working on a problem where I have non-standard Hopf-bifurcation points (Hopf-Hopf or degenerate Hopf) which get an “nd” (not documented) mark on my diagrams. Since these are not identified as Hopf points, the Hopf continuation fails of course. The question is: can I still somehow find the periodic orbits (stable and unstable)?

``````using BifurcationKit
using Plots,LaTeXStrings
using Setfield
using LinearAlgebra
const BK=BifurcationKit

function skipsum(a,i)
x = a[1:2:end]
return sum(x) - a[i]
end

function att_din!(u,p)
g,w,n = p
n = Int(n)
du = similar(u)
for i in 1:2:(2*n)
du[i] = g*u[i]*(1 - (1 - w)*u[i+1] - w/(n-1)*skipsum(u,i))
du[i+1] = u[i] - u[i+1]
end
du
end

par = (g = 5,w = 0.1,n = 3)

# sup norm
norminf(x) = norm(x, Inf)

initconds = [[0.0,0.0,0.0,0.0,0.0,0.0],
ones(2*par.n),
[1/(1-par.w),1/(1-par.w),0.0,0.0,0.0,0.0],
[0.0,0.0,1/(1-par.w/2),1/(1-par.w/2),1/(1-par.w/2),1/(1-par.w/2)]]

line_colors = [:gray10,:red,:blue,:green]
labels = [L"\$x^{*}_0\$",L"\$x^{*}_1\$",L"\$x^{*}_2\$",L"\$x^{*}_3\$"]
pl = plot()

for (i,ic) in enumerate(initconds)
@show ic

prob = BifurcationProblem(att_din!, ic, par, (@lens _.w),
recordFromSolution = (x, p) -> norminf(x),
plotSolution = (x, p; k...) -> plot!(x; ylabel="solution", label="", k...))

# continuation options
opts_br = ContinuationPar(pMin = 0.05, pMax = 0.8, dsmax = 0.01,
# options to detect bifurcations
detectBifurcation = 3, nInversion = 8, maxBisectionSteps = 25,
# number of eigenvalues
nev = 6,
# maximum number of continuation steps
maxSteps = 1000,
# parameter theta, see `? continuation`. Setting this to a non
# default value helps passing the transcritical bifurcation
θ = 0.3)

br = continuation(prob, PALC(tangent=Bordered()), opts_br;
plot = false,normC = norminf)

if i == 2
@show br
end

plot!(br;
putspecialptlegend=true,
plotstability=true,
linewidthunstable=1,
linewidthstable=5,
plotcirclesbif = false,
tickfontsize=10,
lc = line_colors[i],
xguidefontsize=20,
yguidefontsize=20,
ms=6,
markerstrokewidth=0.01,
legendfontsize = 10,
legend=:topleft,
xticks=[0:0.1:0.8;],
xlabel=L"\$w\$",
dpi=300)

end
plot!(pl,ylabel=L"\$max(|x^{*}|)\$")
display(pl)

``````

There is no Hopf-Hopf bifurcation in codim 1 generically. Your system must have some symmetries imbedded which rise the eigenvalue multiplicity.

There is no periodic orbit automatic branch switching from Hopf-Hopf as of now. The normal form is not robust against higher order terms.

If you happen to know a periodic orbit close to this point, you can modify the function to perform branch switching

It seems like an equivariant Hopf bifurcation with symmetry group SO(N)?

Yes, my system has symmetries which is responsible for the eigenvalue multiplicity. More precisely, for a certain equilibrium point there will be N-1 pairs of purely imaginary eigenvalues at the bifurcation (for 2N dimensional system).

Could you elaborate on that a little bit? I’m interested in a method to find periodic orbits without having a Hopf point. Could that be done by passing an `orbitguess`?

You could work out the equivariant Hopf normal form. Your system (I think) has SO(N) symmetry if the particles are on a circle. Golubitsky, Martin, Ian Stewart, and David G. Schaeffer. *Singularities and Groups in Bifurcation Theory. Vol. 2 1988. See XVIII.4 page 387. Then you define predictors for the periodic orbit, see here which can be used for branch switching.

Could that be done by passing an `orbitguess`

Yes! see here for shooting, here for trapezoid and here for collocation. Basically, you have to call `continuation(prob, orbitguess , alg, contparams)` where `prob <: AbstractPeriodicOrbitProblem`. I realize that I have no tutorials for that and I am writing one. Here is an example below. I will provide helper functions to ease this

``````using Revise
using Test, ForwardDiff, Parameters, Setfield, Plots, LinearAlgebra
using BifurcationKit
const BK = BifurcationKit

norminf(x) = norm(x, Inf)
####################################################################################################
function Pop!(du, X, p, t = 0)
@unpack r,K,a,ϵ,b0,e,d, = p
x, y, u, v = X
p = a*x/(b0*(1+ϵ*u)+x)
du[1] = r*(1-x/K)*x-p*y
du[2] = e*p*y-d*y
s = u^2+v^2
du[3] = u-2pi*v - s*u
du[4] = 2pi*u+v - s*v
du
end
Pop(u,p) = Pop!(similar(u),u,p,0)

par_pop = ( K = 1., r = 6.28, a = 12.56, b0 = 0.25, e = 1., d = 6.28, ϵ = 0.2,)

z0 = [0.1,0.1,1,0]

prob = BK.BifurcationProblem(Pop, z0, par_pop, (@lens _.b0); recordFromSolution = (x, p) -> (x = x[1], y = x[2], u = x[3]))

opts_br = ContinuationPar(pMin = 0., pMax = 20.0, ds = 0.002, dsmax = 0.01, nInversion = 6, detectBifurcation = 3, maxBisectionSteps = 25, nev = 4, maxSteps = 2000)
@set! opts_br.newtonOptions.verbose = true

################################################################################
using DifferentialEquations
prob_de = ODEProblem(Pop!, z0, (0,200.), par_pop)
sol = solve(prob_de, Rodas4())
prob_de = ODEProblem(Pop!, sol.u[end], (0,3.), par_pop, reltol = 1e-7, abstol = 1e-10)
sol = solve(prob_de, Rodas4())

plot(sol)

################################################################################
argspo = (recordFromSolution = (x, p) -> begin
xtt = BK.getPeriodicOrbit(p.prob, x, @set par_pop.b0 = p.p)
return (max = maximum(xtt[1,:]),
min = minimum(xtt[1,:]),
period = getPeriod(p.prob, x, @set par_pop.b0 = p.p))
end,
plotSolution = (x, p; k...) -> begin
xtt = BK.getPeriodicOrbit(p.prob, x, @set par_pop.b0 = p.p)
plot!(xtt.t, xtt[1,:]; label = "x", k...)
plot!(xtt.t, xtt[2,:]; label = "y", k...)
# plot!(br; subplot = 1, putspecialptlegend = false)
end)
################################################################################
probtrap = PeriodicOrbitTrapProblem(M = 100, N = 4, prob_vf = prob, xπ = copy(sol[end]), ϕ = copy(sol[end]), updateSectionEveryStep = 0)
resize!(probtrap.ϕ, probtrap.N*probtrap.M)
resize!(probtrap.xπ, probtrap.N*probtrap.M)
period = 1.9929
ci = BK.generateSolution(probtrap, t -> sol(t*period/2/pi), period)
_sol = BK.getPeriodicOrbit(probtrap, ci, nothing)
probtrap.xπ .= 0
probtrap.ϕ .= reduce(vcat, [Pop(_sol.u[:,i],par_pop) for i=1:probtrap.M])

plot(sol)
probtrap(ci, prob.params) |> plot

solpo = newton(probtrap, ci, NewtonPar(verbose = true))
``````
2 Likes

I forgot but you can adapt this 1d Ginzburg-Landau equation (TW) · Bifurcation Analysis in Julia

1 Like