🟡 Colpitts–type oscillator

• 🟡 Colpitts–type oscillator
• • Curve of Hopf bifurcations
  • Periodic orbits with Trapezoid method
  • Periodic orbits with Multiple Standard Shooting
  • Branching from period doubling bifurcation
  • References

In this tutorial, we show how to study parametrized DAEs like:

\[M(\mu)\dot x = G(\mu,x).\]

In particular, we detect a Hopf bifurcation and compute the periodic orbit branching from it using a multiple standard shooting method.

The following DAE model is taken from ^[Rabier]:

\[\left(\begin{array}{cccc} -\left(C_{1}+C_{2}\right) & C_{2} & 0 & 0 \\ C_{2} & -C_{2} & 0 & 0 \\ C_{1} & 0 & 0 & 0 \\ 0 & 0 & L & 0 \end{array}\right)\left(\begin{array}{c} \dot{x}_{1} \\ \dot{x}_{2} \\ \dot{x}_{3} \\ \dot{x}_{4} \end{array}\right)=\left(\begin{array}{c} R^{-1}\left(x_{1}-V\right)+I E\left(x_{1}, x_{2}\right) \\ x_{3}+I C\left(x_{1}, x_{2}\right) \\ -x_{3}-x_{4} \\ -\mu+x_{2} \end{array}\right)\]

It is easy to encode the DAE as follows. The mass matrix is defined next.

using Revise, Plotsimport LinearAlgebra as LAimport BifurcationKit as BKimport BifurcationKit: @optic, @set, @reset# function to record information from the solutionrecordFromSolution(x, p; k...) = (u1 = BK.norminf(x), x1 = x[1], x2 = x[2], x3 = x[3], x4 = x[4])# vector fieldf(x, p) = p.Is * (exp(p.q * x) - 1)IE(x1, x2, p) = -f(x2, p) + f(x1, p) / p.αFIC(x1, x2, p) = f(x2, p)/ p.αR - f(x1, p)function Colpitts!(dz, z, p, t = 0)	(;C1, C2, L, R, Is, q, αF, αR, V, μ) = p	x1, x2, x3, x4 = z	dz[1] = (x1 - V) / R + IE(x1, x2, p)	dz[2] =	x3 + IC(x1, x2, p)	dz[3] = -x3-x4	dz[4] = -μ+x2	dzend# parameter valuespar_Colpitts = (C1 = 1.0, C2 = 1.0, L = 1.0, R = 1/4., Is = 1e-16, q = 40., αF = 0.99, αR = 0.5, μ = 0.5, V = 6.)# initial conditionz0 = [0.9957,0.7650,19.81,-19.81]# mass matrixBe = [-(par_Colpitts.C1+par_Colpitts.C2) par_Colpitts.C2 0 0;par_Colpitts.C2 -par_Colpitts.C2 0 0;par_Colpitts.C1 0 0 0; 0 0 par_Colpitts.L 0]# we group the differentials togetherprob = BK.ODEBifProblem(Colpitts!, z0, par_Colpitts, (@optic _.μ); record_from_solution = recordFromSolution)dae_problem = BK.DAEMassBifProblem(prob, Be)

We first compute the branch of equilibria. But we need a generalized eigenvalue solver for this.

opts_br = BK.ContinuationPar(p_min = 0.4, p_max = 0.8, ds = 0.01, dsmax = 0.02, n_inversion = 4)br = BK.continuation(dae_problem, BK.AutoSwitch(; alg = BK.PALC(tangent = BK.Bordered())), opts_br; normC = BK.norminf)scene = plot(br, vars = (:param, :x1))

Curve of Hopf bifurcations

BK.get_normal_form(br, 1; start_with_eigen = Val(false))

SuperCritical - Hopf bifurcation point at μ ≈ 0.765053790158384.
Frequency ω ≈ 8.108242584323488
Period of the periodic orbit ≈ 0.7749133356379253
Normal form z⋅(iω + a⋅δμ + b⋅|z|²):
┌─ a = -73.83147885834487 - 45.805447919198585im
└─ b = -1851.2524033591308 - 818.780448920847im

br_hopf = BK.continuation(br, 1, (@optic _.αR),		BK.ContinuationPar(BK.getcontparams(br), p_max = 41.5, p_min = 0.1);			start_with_eigen = false,			bothside = true,			callback_newton = BK.cbMaxNorm(1)			)scene = plot(br_hopf)

Periodic orbits with Trapezoid method

function recordPO(u, p; iter, state, k...)    outt = BK.get_periodic_orbit(p.prob, u, BK.getparams(iter, state))    m = maximum(outt[1,:])    return (s = m, period = BK.getperiod(p.prob, u, BK.getparams(iter, state)))endfunction PlotPO(x, p; iter, state, k...)    outt = BK.get_periodic_orbit(p.prob, x, BK.getparams(iter, state))    plot!(outt.t, outt[2, :], subplot = 3)    plot!(br, vars = (:param, :x1), subplot = 1)end# we lower the tolerance of newton for the periodic orbitsoptnpo = BK.NewtonPar(br.contparams.newton_options; tol = 1e-10)@reset optnpo.eigsolver = BK.DefaultEig()opts_po_cont = BK.ContinuationPar(dsmin = 0.0001, dsmax = 0.005, ds= -0.003, p_min = 0.2, max_steps = 50, newton_options = optnpo, tol_stability = 1e-3)# automatic branching from the Hopf pointbr_po = BK.continuation(br, 1, opts_po_cont,    BK.Trapeze(M = 200, jacobian = BK.Dense(), massmatrix = Be),    δp = -0.001,    start_with_eigen = Val(false),    record_from_solution = recordPO,    plot_solution = PlotPO,	callback_newton = BK.cbMaxNorm(1),    normC = BK.norminf    )scene = plot(br_po)

Periodic orbits with Multiple Standard Shooting

We use shooting to compute periodic orbits: we rely on a fixed point of the flow. To compute the flow, we use DifferentialEquations.jl.

Thanks to ^[Lamour], we can just compute the Floquet coefficients to get the nonlinear stability of the periodic orbit. Two period doubling bifurcations are detected.

Note that the Hopf normal form for a DAE with a constant mass matrix is now supported (using start_with_eigen = Val(false), see the page Differential-Algebraic Equations (DAE)). The Automatic Branch Switching from the Hopf bifurcation point is therefore based on the correct normal form; for this model the bifurcation is supercritical.

import OrdinaryDiffEq as ODE# this is the ODEProblem used with `DiffEqBase.solve`prob_dae = ODE.ODEFunction(Colpitts!; mass_matrix = Be)probFreez_ode = ODE.ODEProblem(prob_dae, z0, (0, 1), par_Colpitts)# we lower the tolerance of newton for the periodic orbitsoptnpo = BK.NewtonPar(br.contparams.newton_options; tol = 1e-9)@reset optnpo.eigsolver = BK.DefaultEig()opts_po_cont = BK.ContinuationPar(dsmin = 0.0001, dsmax = 0.005, ds= -0.001, p_min = 0.2, max_steps = 60, newton_options = optnpo, tol_stability = 1e-3)# automatic branching from the Hopf pointbr_po_sh = BK.continuation(br, 1, opts_po_cont,    BK.Shooting(10, probFreez_ode, ODE.Rodas5P(); reltol = 1e-9, abstol = 1e-11, parallel = true);    δp = -0.001,    start_with_eigen = Val(false),    record_from_solution = recordPO,    plot_solution = PlotPO,    # the newton callback is used to reject residual > 1    # this is to avoid numerical instabilities from DE.jl    callback_newton = BK.cbMaxNorm(1.0),    normC = BK.norminf)plot(br_po_sh)

Branching from period doubling bifurcation

br_po_sh2 = BK.continuation(deepcopy(br_po_sh), 1,    BK.ContinuationPar(opts_po_cont, max_steps = 15);    record_from_solution = recordPO,    use_normal_form = false, ampfactor = 0.01,    δp = -0.002,    callback_newton = BK.cbMaxNorm(1.0),    normC = BK.norminf)plot(br, vars = (:param, :x1));plot!(br_po_sh2, br_po_sh)

Let us show an example of periodic solution on the PD branch:

sol = BK.get_periodic_orbit(br_po_sh2, 3)plot(sol.t, sol[2, :])

References