🟡 Period doubling in Lur'e problem (PD aBS)

• 🟡 Period doubling in Lur'e problem (PD aBS)
• • Periodic orbits with orthogonal collocation
  • Periodic orbits with Parallel Standard Shooting
  • Branch of periodic orbits with Trapeze

The following model is an adaptive control system of Lur’e type. It is an example from the MatCont library.

\[\left\{\begin{array}{l} \dot{x}=y \\ \dot{y}=z \\ \dot{z}=-\alpha z-\beta y-x+x^{2} \end{array}\right.\]

The model is interesting because there is a period doubling bifurcation and we want to show the branch switching capabilities of BifurcationKit.jl in this case. We provide 3 different ways to compute the periodic orbits and highlight their pro / cons.

It is easy to encode the ODE as follows

using Revise, Plotsimport BifurcationKit as BKimport BifurcationKit: @opticrecordFromSolution(x, p; k...) = (u1 = x[1], u2 = x[2])function lur!(dz, u, p, t = 0)	(;α, β) = p	x, y, z = u	dz[1] = y	dz[2] =	z	dz[3] = -α * z - β * y - x + x^2	dzend# bifurcation problemprob_bif = BK.ODEBifProblem(lur!, zeros(3), (α = -1.0, β = 1.0), (@optic _.α);    record_from_solution = recordFromSolution)

We first compute the branch of equilibria

# continuation optionsopts_br = BK.ContinuationPar(p_min = -1.4, p_max = 1.8, nev = 3)# computation of the branchbr = BK.continuation(prob_bif, BK.PALC(), opts_br)scene = plot(br)

With detailed information:

br

 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 26
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter α starts at -1.0, ends at 1.8
 ├─ Algo: PALC [Secant]
 └─ Special points:

- #  1,     hopf at α ≈ +1.00420578 ∈ (+0.99978637, +1.00420578), |δp|=4e-03, [converged], δ = (-2, -2), step =  19
- #  2, endpoint at α ≈ +1.80000000,                                                                     step =  25

We note the Hopf bifurcation point which we shall investigate now.

Periodic orbits with orthogonal collocation

We compute the branch of periodic orbits from the Hopf bifurcation point. We rely on a the state of the art method for computing periodic orbits of ODE: orthogonal collocation.

We first define a plotting function and a record function which are used for all cases below:

# plotting functionfunction plotPO(x, p; iter, state, k...)	xtt = BK.get_periodic_orbit(p.prob, x, BK.getparams(iter, state))	plot!(xtt.t, xtt[1,:]; markersize = 2, k...)	plot!(xtt.t, xtt[2,:]; k...)	plot!(xtt.t, xtt[3,:]; legend = false, k...)end# record functionfunction recordPO(x, p; iter, state, k...)	xtt = BK.get_periodic_orbit(p.prob, x, BK.getparams(iter, state))	period = BK.getperiod(p.prob, x, BK.getparams(iter, state))	mn, mx = extrema(xtt[1, :])	return (;max = mx, min = mn, period)end

recordPO (generic function with 1 method)

Continuation of periodic orbits from the Hopf point:

# continuation parametersopts_po_cont = BK.ContinuationPar(opts_br, dsmax = 0.03, dsmin = 1e-4, max_steps = 80, tol_stability = 1e-4, plot_every_step = 20)br_po = BK.continuation(	br, 1, opts_po_cont,	BK.Collocation(40, 4);	plot = true,	# linear_algo = BK.COPBLS(),	plot_solution = plotPO,	normC = BK.norminf)scene = plot(br, br_po)

Note that you can compute the PD normal form

BK.get_normal_form(br_po, 1)

Period-Doubling bifurcation point of periodic orbit
├─ Period = 6.364577072676575 -> 12.72915414535315
├─ Problem: Collocation
├─ α ≈ 0.6277212757673908
├─ type: ├─ Normal form (Iooss):
├	∂τ = 1 + a₀₁⋅δp + a₂⋅ξ²
├	∂ξ =  ξ⋅(c₁₁⋅δp + c₃⋅ξ²)
├─── a₀₁ = 0.030754360129198768
├─── a₂  = 0.0032149109236260733
├─── c₁₁ = -1.2973132584441993
└─── c₃  = -0.3019436508509669

We provide Automatic Branch Switching from the PD point and computing the bifurcated branch is as simple as:

# aBS from PDbr_po_pd = BK.continuation(deepcopy(br_po), 1, BK.ContinuationPar(br_po.contparams, max_steps = 100, dsmax = 0.02, plot_every_step = 10, ds = 0.005);	plot_solution = plotPO,	normC = BK.norminf,	)scene = plot(br_po, br_po_pd, title = "Collocation based")

Periodic orbits with Parallel Standard Shooting

We use a different method to compute periodic orbits: we rely on a fixed point of the flow. To compute the flow, we use OrdinaryDiffEq.jl. This way of computing periodic orbits should be less precise than the previous one. We rely on parallel multiple shooting. Finally, please note the close similarity to the code of the previous part. As before, we first rely on Hopf aBS.

import OrdinaryDiffEq as ODE# ODE problem for using DifferentialEquationsprob_ode = ODE.ODEProblem(lur!, prob_bif.u0, (0, 1.), prob_bif.params; abstol = 1e-12, reltol = 1e-10)# continuation parametersopts_po_cont = BK.ContinuationPar(opts_br, dsmax = 0.03, newton_options = BK.NewtonPar(tol = 1e-10), tol_stability = 1e-5, n_inversion = 8, max_steps = 100)br_po = BK.continuation(	br, 1, opts_po_cont,	BK.Shooting(5, prob_ode, ODE.Vern9(); parallel = true, abstol = 1e-10, reltol = 1e-8);	# parallel shooting functional with several sections	# plot = true,	record_from_solution = recordPO,	plot_solution = plotPO,	# limit the residual, useful to help DifferentialEquations	callback_newton = BK.cbMaxNorm(10),	normC = BK.norminf)plot(br_po)

Note that you can compute the PD normal form

BK.get_normal_form(br_po, 1; verbose = true)

Period-Doubling bifurcation point of periodic orbit
├─ Period = 6.364072923581334 -> 12.728145847162669
├─ Problem: Shooting
SuperCritical - Period-Doubling ┌─ Normal form:
├	 x ─▶ x⋅(a⋅δα - 1 + c⋅x²)
├─ a = 8.834936357763215
└─ c = 73.89114583600356

We provide Automatic Branch Switching from the PD point and computing the bifurcated branch is as simple as:

# aBS from PDbr_po_pd = BK.continuation(deepcopy(br_po), 1,	BK.ContinuationPar(br_po.contparams, max_steps = 10, ds = -0.015);	record_from_solution = recordPO,	normC = BK.norminf,	# verbosity = 3, plot = true,	callback_newton = BK.cbMaxNorm(10),	)scene = plot(br, br_po, br_po_pd, title = "Shooting based")

Branch of periodic orbits with Trapeze

We use finite differences to discretize the problem for finding periodic orbits. We appeal to automatic branch switching from the Hopf point as follows

# continuation parametersopts_po_cont = BK.ContinuationPar(dsmax = 0.02, dsmin = 1e-4, p_max = 1.1, max_steps = 80, tol_stability = 1e-4)br_po = BK.continuation(	br, 1, opts_po_cont,	BK.Trapeze(M = 120);	record_from_solution = recordPO,	plot_solution = (x, p; iter, state, k...) -> begin		plotPO(x, p; iter, state, k...)		## plot previous branch		plot!(br, subplot=1, putbifptlegend = false)		end,	normC = BK.norminf)scene = plot(br, br_po)

Two period doubling bifurcations were detected. We shall now compute the branch of periodic orbits from these PD points. We do not provide Automatic Branch Switching for Trapezoid method as we do not have yet the PD normal form computed for Trapeze. Hence, it takes some trial and error to find the ampfactor of the PD branch.

This is like in MatCont and Auto-07p here...

# aBS from PDbr_po_pd = BK.continuation(deepcopy(br_po), 1, BK.ContinuationPar(br_po.contparams, max_steps = 70);	# plot = true,	use_normal_form = false,	ampfactor = .2, δp = -0.005,	plot_solution = (x, p; iter, state, k...) -> begin		plotPO(x, p; iter, state, k...)		# add previous branch		plot!(br_po; legend=false, subplot=1)	end,	record_from_solution = recordPO,	normC = BK.norminf	)Scene = title!("")

plot(br, br_po, br_po_pd, title = "Trapezoid based")