Humphries model (codim 2, periodic orbit)

• Humphries model (codim 2, periodic orbit)
• • Continuation and codim 1 bifurcations
  • Continuation of Hopf point
  • Branch of periodic orbits
  • References

Consider the model ^[Hum] as an example of state-dependent delays

\[x^{\prime}(t)=-\gamma x(t)-\kappa_1 x\left(t-a_1-c x(t)\right)-\kappa_2 x\left(t-a_2-c x(t)\right)\]

Continuation and codim 1 bifurcations

We first instantiate the model

using Revise, DDEBifurcationKit, Parameters, Plots
using BifurcationKit
const BK = BifurcationKit

function humpriesVF(x, xd, p)
   @unpack κ1,κ2,γ,a1,a2,c = p
   [
      -γ * x[1] - κ1 * xd[1][1] - κ2 * xd[2][1]
   ]
end

function delaysF(x, par)
   [
      par.a1 + par.c * x[1],
      par.a2 + par.c * x[1],
   ]
end


pars = (κ1=0.,κ2=2.3,a1=1.3,a2=6,γ=4.75,c=1.)
x0 = zeros(1)

prob = SDDDEBifProblem(humpriesVF, delaysF, x0, pars, (@lens _.κ1))

┌─ State-dependent delays Bifurcation Problem with uType Vector{Float64}
├─ Inplace:  false
└─ Parameter: κ1

We then compute the branch

optn = NewtonPar(verbose = true, eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 13., p_min = 0., newton_options = optn, ds = -0.01, detect_bifurcation = 3, nev = 3, )
br = continuation(prob, PALC(), opts; verbosity = 0, bothside = true)

 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 100
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter κ1 starts at 13.0, ends at 0.0
 ├─ Algo: PALC
 └─ Special points:

If `br` is the name of the branch,
ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`

- #  1, endpoint at κ1 ≈ +13.00000000,                                                                      step =   0
- #  2,     hopf at κ1 ≈ +4.28020573 ∈ (+4.27136690, +4.28020573), |δp|=9e-03, [converged], δ = (-2, -2), step =  62, eigenelements in eig[ 63], ind_ev =   2
- #  3, endpoint at κ1 ≈ +0.00000000,                                                                      step =  99

and plot it

scene = plot(br)

Continuation of Hopf point

We follow the Hopf points in the parameter plane $(\kappa_1,\kappa_2)$. We tell the solver to consider br.specialpoint[2] and continue it.

brhopf = continuation(br, 2, (@lens _.κ2),
         setproperties(br.contparams, detect_bifurcation = 2, dsmax = 0.04, max_steps = 230, p_max = 5., p_min = -1.,ds = -0.02);
         verbosity = 0, plot = false,
         # we disable detection of Bautin bifurcation as the
         # Hopf normal form is not implemented for SD-DDE
         detect_codim2_bifurcation = 0,
         bothside = true,
         start_with_eigen = true)

scene = plot(brhopf, vars = (:κ1, :κ2))

Branch of periodic orbits

We compute the branch of periodic orbits from the Hopf bifurcation points using orthogonal collocation. We use a lot of time sections $N_{tst}=200$ to have enough precision to resolve the sophisticated branch of periodic solutions especially near the first Fold point around $\kappa_1\approx 10$.

# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.05, ds= 0.001, dsmin = 1e-4, p_max = 12., p_min=-5., max_steps = 3000,
	nev = 3, tol_stability = 1e-8, detect_bifurcation = 0, plot_every_step = 20, save_sol_every_step=1)
@set! opts_po_cont.newton_options.tol = 1e-9
@set! opts_po_cont.newton_options.verbose = true

# arguments for periodic orbits
args_po = (	record_from_solution = (x, p) -> begin
		xtt = BK.get_periodic_orbit(p.prob, x, nothing)
		_max = maximum(xtt[1,:])
		_min = minimum(xtt[1,:])
		return (amp = _max - _min,
				period = getperiod(p.prob, x, nothing))
	end,
	plot_solution = (x, p; k...) -> begin
		xtt = BK.get_periodic_orbit(p.prob, x, nothing)
		plot!(xtt.t, xtt[1,:]; label = "x", k...)
		plot!(br; subplot = 1, putspecialptlegend = false)
		end,
	normC = norminf)

probpo = PeriodicOrbitOCollProblem(200, 2; N = 1, jacobian = BK.AutoDiffDense())
br_pocoll = continuation(
	br, 2, opts_po_cont,
	probpo;
	alg = PALC(tangent = Bordered()),
	# regular continuation options
	verbosity = 2,	plot = true,
	args_po...,
	ampfactor = 1/0.467829783456199 * 0.1,
	δp = 0.01,
	callback_newton = BK.cbMaxNorm(10.0),
	)	

which gives

References