Finding periodic orbits with BifurcationKit.jl

Specific Domains Modelling & Simulations
, ,
1

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)
2

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

3

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

4

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?

5

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
6

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

1 Like