# Finding periodic orbits with BifurcationKit.jl
**URL:** https://discourse.julialang.org/t/finding-periodic-orbits-with-bifurcationkit-jl/93085
**Category:** Modelling & Simulations
**Tags:** question, package, diffeq
**Created:** [January 17, 2023, 11:21am UTC](https://discourse.julialang.org/t/finding-periodic-orbits-with-bifurcationkit-jl/93085 "2023-01-17T11:21:54Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@rusandris](https://discourse.julialang.org/u/rusandris)
#### Post date: [January 17, 2023, 11:21am UTC](https://discourse.julialang.org/t/finding-periodic-orbits-with-bifurcationkit-jl/93085/1 "2023-01-17T11:21:54Z")
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)?
```julia
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)
```
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 17, 2023, 12:09pm UTC](https://discourse.julialang.org/t/finding-periodic-orbits-with-bifurcationkit-jl/93085/2 "2023-01-17T12:09:43Z")
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](https://github.com/bifurcationkit/BifurcationKit.jl/blob/master/src/periodicorbit/PeriodicOrbits.jl#L327) to perform branch switching
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 18, 2023, 8:26am UTC](https://discourse.julialang.org/t/finding-periodic-orbits-with-bifurcationkit-jl/93085/3 "2023-01-18T08:26:28Z")
It seems like an equivariant Hopf bifurcation with symmetry group SO(N)?
---
### Author:  [@rusandris](https://discourse.julialang.org/u/rusandris)
#### Post date: [January 18, 2023, 8:45am UTC](https://discourse.julialang.org/t/finding-periodic-orbits-with-bifurcationkit-jl/93085/4 "2023-01-18T08:45:35Z")
> [@rveltz](#):
>
> There is no Hopf-Hopf bifurcation in codim 1 generically. Your system must have some symmetries imbedded which rise the eigenvalue multiplicity.
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).
> [@rveltz](#):
>
> If you happen to know a periodic orbit close to this point, you can modify the [function](https://github.com/bifurcationkit/BifurcationKit.jl/blob/master/src/periodicorbit/PeriodicOrbits.jl#L327) to perform branch switching
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`?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 18, 2023, 10:22am UTC](https://discourse.julialang.org/t/finding-periodic-orbits-with-bifurcationkit-jl/93085/5 "2023-01-18T10:22:26Z")
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](https://github.com/bifurcationkit/BifurcationKit.jl/blob/master/src/NormalForms.jl#L799) which can be used for branch switching.
> Could that be done by passing an `orbitguess`
Yes! see [here](https://github.com/bifurcationkit/BifurcationKit.jl/blob/master/src/periodicorbit/PeriodicOrbits.jl#L220) for shooting, [here](https://github.com/bifurcationkit/BifurcationKit.jl/blob/master/src/periodicorbit/PeriodicOrbitTrapeze.jl#L932) for trapezoid and [here](https://github.com/bifurcationkit/BifurcationKit.jl/blob/master/src/periodicorbit/PeriodicOrbitCollocation.jl#L519) 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
```julia
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))
```
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 19, 2023, 5:35pm UTC](https://discourse.julialang.org/t/finding-periodic-orbits-with-bifurcationkit-jl/93085/6 "2023-01-19T17:35:02Z")
I forgot but you can adapt this [1d Ginzburg-Landau equation (TW) · Bifurcation Analysis in Julia](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/tutorials/cgl1dwave/#Computation-of-the-travelling-wave)