π‘ Period doubling in Lur'e problem (PD aBS)
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 this periodic orbits and highlight their pro / cons.
It is easy to encode the ODE as follows
using Revise, Plots
using BifurcationKit
const BK = BifurcationKit
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 = ODEBifProblem(lur!, zeros(3), (Ξ± = -1.0, Ξ² = 1.0), (@optic _.Ξ±);
record_from_solution = recordFromSolution)
We first compute the branch of equilibria
# continuation options
opts_br = ContinuationPar(p_min = -1.4, p_max = 1.8, dsmax = 0.01, max_steps = 1000)
# computation of the branch
br = continuation(prob_bif, PALC(), opts_br)
scene = plot(br)
With detailed information:
br
ββ Curve type: EquilibriumCont
ββ Number of points: 200
ββ Type of vectors: Vector{Float64}
ββ Parameter Ξ± starts at -1.0, ends at 1.8
ββ Algo: PALC
ββ Special points:
- # 1, hopf at Ξ± β +1.00022831 β (+0.99934442, +1.00022831), |Ξ΄p|=9e-04, [converged], Ξ΄ = (-2, -2), step = 142
- # 2, endpoint at Ξ± β +1.80000000, step = 199
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 = 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 = get_periodic_orbit(p.prob, x, p.p)
period = 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 = ContinuationPar(opts_br, dsmax = 0.03, ds = 0.01, dsmin = 1e-4, max_steps = 80, tol_stability = 1e-4, plot_every_step = 20)
br_po = continuation(
br, 1, opts_po_cont,
PeriodicOrbitOCollProblem(40, 4);
plot = true,
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 = norminf)
scene = plot(br, br_po)
Note that you can compute the PD normal form
get_normal_form(br_po, 1)
Period-Doubling bifurcation point of periodic orbit
ββ Period = 6.364079037215513 -> 12.728158074431025
ββ Problem: PeriodicOrbitOCollProblem
ββ Ξ± β 0.6302627694057382
ββ type: ββ Normal form (Iooss):
β βΟ = 1 + aβββ
Ξ΄p + aββ
ΞΎΒ²
β βΞΎ = ΞΎβ
(cβββ
Ξ΄p + cββ
ΞΎΒ²)
ββββ aββ = 0.03082683879264809
ββββ aβ = 0.003295249439029141
ββββ cββ = -1.2961053260562345
ββββ cβ = -0.30508490446797343
We provide Automatic Branch Switching from the PD point and computing the bifurcated branch is as simple as:
# aBS from PD
br_po_pd = continuation(deepcopy(br_po), 1, 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 = 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 DifferentialEquations.jl
. This way of computing periodic orbits should be more precise than the previous one. We use a particular instance called multiple shooting which is computed in parallel. This is an additional advantage compared to the previous method. Finally, please note the close similarity to the code of the previous part. As before, we first rely on Hopf aBS.
import DifferentialEquations as DE
# ODE problem for using DifferentialEquations
prob_ode = DE.ODEProblem(lur!, prob_bif.u0, (0., 1.), prob_bif.params; abstol = 1e-12, reltol = 1e-10)
# continuation parameters
# we decrease a bit the newton tolerance to help automatic branch switching from PD point
opts_po_cont = ContinuationPar(dsmax = 0.03, ds= -0.001, newton_options = NewtonPar(tol = 1e-10), tol_stability = 1e-5, n_inversion = 8, nev = 3, plot_every_step = 1)
br_po = continuation(
br, 1, opts_po_cont,
# parallel shooting functional with 10 sections
ShootingProblem(10, prob_ode, DE.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 = norminf)
plot(br_po)
Note that you can compute the PD normal form
get_normal_form(br_po, 1; autodiff = true)
Period-Doubling bifurcation point of periodic orbit
ββ Period = 6.364071768026619 -> 12.728143536053238
ββ Problem: ShootingProblem
SuperCritical - Period-Doubling ββ Normal form:
β x ββΆ xβ
(aβ
Ξ΄p - 1 + cβ
xΒ²)
ββ a = 120.2299284458693
ββ c = 75.64571159720495
We provide Automatic Branch Switching from the PD point and computing the bifurcated branch is as simple as:
# aBS from PD
br_po_pd = continuation(deepcopy(br_po), 1,
ContinuationPar(br_po.contparams, max_steps = 6, ds = -0.008);
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 = 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 = ContinuationPar(dsmax = 0.02, dsmin = 1e-4, p_max = 1.1, max_steps = 80, tol_stability = 1e-4)
br_po = continuation(
br, 1, opts_po_cont,
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 = 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 as we do not have the PD normal form computed for PeriodicOrbitTrapProblem
. Hence, it takes some trial and error to find the ampfactor
of the PD branch.
# aBS from PD
br_po_pd = continuation(deepcopy(br_po), 1, 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 = norminf
)
Scene = title!("")
plot(br, br_po, br_po_pd, title = "Trapezoid based")