🟡 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 finite differences

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, Plots
import BifurcationKit as BK
import BifurcationKit: @optic

recordFromSolution(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
	dz
end

# bifurcation problem
prob_bif = BK.ODEBifProblem(lur!, zeros(3), (α = -1.0, β = 1.0), (@optic _.α);
    record_from_solution = recordFromSolution)

We first compute the branch of equilibria

# continuation options
opts_br = BK.ContinuationPar(p_min = -1.4, p_max = 1.8, nev = 3)

# computation of the branch
br = 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 function
function plotPO(x, p; k...)
	xtt = BK.get_periodic_orbit(p.prob, x, p.p)
	plot!(xtt.t, xtt[1,:]; markersize = 2, k...)
	plot!(xtt.t, xtt[2,:]; k...)
	plot!(xtt.t, xtt[3,:]; legend = false, k...)
end

# record function
function recordPO(x, p; k...)
	xtt = BK.get_periodic_orbit(p.prob, x, p.p)
	period = BK.getperiod(p.prob, x, p.p)
	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 parameters
opts_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.PeriodicOrbitOCollProblem(40, 4);
	plot = true,
	# linear_algo = BK.COPBLS(),
	record_from_solution = recordPO,
	plot_solution = (x, p; k...) -> begin
		plotPO(x, p; k...)
		# plot previous branch
		plot!(br, subplot = 1, putbifptlegend = false)
		end,
	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.364577072721098 -> 12.729154145442196
├─ Problem: PeriodicOrbitOCollProblem
├─ α ≈ 0.6277212755398428
├─ type: ├─ Normal form (Iooss):
├	∂τ = 1 + a₀₁⋅δp + a₂⋅ξ²
├	∂ξ =  ξ⋅(c₁₁⋅δp + c₃⋅ξ²)
├─── a₀₁ = 0.030754360020421895
├─── a₂  = 0.0032024913046831715
├─── c₁₁ = -1.2706988113477693
└─── c₃  = -0.2970122980010248

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

# aBS from PD
br_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 = (x, p; k...) -> begin
		plotPO(x, p; k...)
		## add previous branch
		plot!(br_po; subplot = 1)
	end,
	record_from_solution = recordPO,
	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 DifferentialEquations
prob_ode = ODE.ODEProblem(lur!, prob_bif.u0, (0, 1), prob_bif.params; abstol = 1e-12, reltol = 1e-10)

# continuation parameters
opts_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,
	# parallel shooting functional with 10 sections
	BK.ShootingProblem(15, prob_ode, ODE.Vern9(); parallel = true);
	# 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; autodiff = true)

Period-Doubling bifurcation point of periodic orbit
├─ Period = 6.364074622670103 -> 12.728149245340205
├─ Problem: ShootingProblem
SuperCritical - Period-Doubling ┌─ Normal form:
├	 x ─▶ x⋅(a⋅δα - 1 + c⋅x²)
├─ a = 74.47963240896509
└─ c = 73.6406345619471

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

# aBS from PD
br_po_pd = BK.continuation(deepcopy(br_po), 1,
	BK.ContinuationPar(br_po.contparams, max_steps = 10, ds = -0.005);
	plot = true, verbosity = 2,
	plot_solution = (x, p; k...) -> begin
		plotPO(x, p; k...)
		# add previous branch
		plot!(br_po; subplot = 1)
	end,
	record_from_solution = recordPO,
	autodiff_nf = true,
	normC = BK.norminf,
	callback_newton = BK.cbMaxNorm(10),
	)

scene = plot(br, br_po, br_po_pd, title = "Shooting based")

Branch of periodic orbits with finite differences

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 parameters
opts_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.PeriodicOrbitTrapProblem(M = 120);
	record_from_solution = recordPO,
	plot_solution = (x, p; k...) -> begin
		plotPO(x, p; 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 PeriodicOrbitTrapProblem. 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 PD
br_po_pd = BK.continuation(deepcopy(br_po), 1, BK.ContinuationPar(br_po.contparams, max_steps = 70);
	# plot = true,
	ampfactor = .2, δp = -0.005,
	plot_solution = (x, p; k...) -> begin
		plotPO(x, p; 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")