Nonlinear laser model
The following model is taken from [Pusuluri]:
\[\left\{\begin{aligned} \dot{\beta} & =-\sigma \beta+g p_{23}, \\ \dot{p}_{21} & =-p_{21}-\beta p_{31}+a D_{21}, \\ \dot{p}_{23} & =-p_{23}+\beta D_{23}-a p_{31}, \\ \dot{p}_{31} & =-p_{31}+\beta p_{21}+a p_{23}, \\ \dot{D}_{21} & =-b\left(D_{21}-D_{21}^0\right)-4 a p_{21}-2 \beta p_{23}, \\ \dot{D}_{23} & =-b\left(D_{23}-D_{23}^0\right)-2 a p_{21}-4 \beta p_{23}, \end{aligned}\right.\]
It is easy to encode the ODE as follows
using Revise, Plots
using Parameters, Setfield, LinearAlgebra, Test, ForwardDiff
using BifurcationKit, Test
using HclinicBifurcationKit
const BK = BifurcationKit
recordFromSolution(x, p) = (D₂₃ = x[6], β = x[1],)
####################################################################################################
function OPL!(dz, u, p, t)
@unpack b, σ, g, a, D₂₁⁰, D₂₃⁰ = p
β, p₂₁, p₂₃, p₃₁, D₂₁, D₂₃ = u
dz[1] = -σ * β + g * p₂₃
dz[2] = -p₂₁ - β * p₃₁ + a * D₂₁
dz[3] = -p₂₃ + β * D₂₃ - a * p₃₁
dz[4] = -p₃₁ + β * p₂₁ + a * p₂₃
dz[5] = -b * (D₂₁ - D₂₁⁰) - 4a * p₂₁ - 2β * p₂₃
dz[6] = -b * (D₂₃ - D₂₃⁰) - 2a * p₂₁ - 4β * p₂₃
dz
end
OPL(z, p) = OPL!(similar(z), z, p, 0)
par_OPL = (b = 1.2, σ = 2.0, g=50., a = 1., D₂₁⁰ = -1., D₂₃⁰ = 0.)
z0 = zeros(6)
prob = BK.BifurcationProblem(OPL, z0, par_OPL, (@lens _.a); record_from_solution = recordFromSolution)WARNING: using Plots.wrap in module BifurcationKit conflicts with an existing identifier.We first compute the branch of equilibria
opts_br = ContinuationPar(p_min = -1., p_max = 8., ds = 0.001, dsmax = 0.06, n_inversion = 6, detect_bifurcation = 3, max_bisection_steps = 25, nev = 6, plot_every_step = 20, max_steps = 100, save_sol_every_step = 1, detect_fold = true)
opts_br = @set opts_br.newton_options.verbose = false
br = continuation(prob, PALC(tangent = Secant()), opts_br;
bothside = false, normC = norminf)
plot(br, plotfold=true)
br2 = continuation(re_make(prob; u0 = [1.6931472491037485, -0.17634826359471437, 0.06772588996414994, -0.23085768742546342, -0.5672243219935907, -0.09634826359471438]), PALC(tangent = Secant()), opts_br; verbosity = 0, bothside = false, normC = norminf)
scene = plot(br, br2)Codimension 2 bifurcations
sn_br = continuation(br, 1, (@lens _.b), ContinuationPar(opts_br, detect_bifurcation = 1, save_sol_every_step = 1, max_steps = 80) ;
alg = PALC(),
verbosity = 0,
detect_codim2_bifurcation = 2,
start_with_eigen = true,
update_minaug_every_step = 1,
bothside = true,
)
hopf_br = continuation(br, 2, (@lens _.b), ContinuationPar(opts_br, detect_bifurcation = 1, save_sol_every_step = 1, max_steps = 140),
detect_codim2_bifurcation = 2,
start_with_eigen = true,
update_minaug_every_step = 1,
bothside = true,
)
hopf_br2 = continuation(br2, 1, (@lens _.b), ContinuationPar(opts_br, detect_bifurcation = 1, save_sol_every_step = 1, max_steps = 140),
detect_codim2_bifurcation = 2,
start_with_eigen = true,
update_minaug_every_step = 1,
bothside = true,
)
plot(sn_br, vars = (:a, :b), branchlabel = "SN", )
plot!(hopf_br, branchlabel = "Hopf", vars = (:a, :b))
plot!(hopf_br2, branchlabel = "Hopf", vars = (:a, :b))
ylims!(0,1.5)Branch of homoclinic orbits with Orthogonal Collocation
# plotting function
function plotPO(x, p; k...)
xtt = BK.get_periodic_orbit(p.prob, x, set(getparams(p.prob), BK.getlens(p.prob), p.p))
plot!(xtt.t, xtt[1,:]; markersize = 2, k...)
plot!(xtt.t, xtt[6,:]; k...)
scatter!(xtt.t, xtt[1,:]; markersize = 1, legend = false, k...)
end
# record function
function recordPO(x, p)
xtt = BK.get_periodic_orbit(p.prob, x, set(getparams(p.prob), BK.getlens(p.prob), p.p))
period = BK.getperiod(p.prob, x, p.p)
return (max = maximum(xtt[6,:]), min = minimum(xtt[6,:]), period = period, )
end
# newton parameters
optn_po = NewtonPar(verbose = true, tol = 1e-8, max_iterations = 25)
# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.05, ds= 0.001, dsmin = 1e-4, p_max = 6.8, p_min=-5., max_steps = 100, newton_options = (@set optn_po.tol = 1e-8), detect_bifurcation = 0, plot_every_step = 3, save_sol_every_step=1,)
br_coll = continuation(
# br, 2,
br2, 1,
opts_po_cont,
PeriodicOrbitOCollProblem(20, 4; meshadapt = true, update_section_every_step = 2);
ampfactor = 1., δp = 0.0015,
verbosity = 2, plot = true,
alg = PALC(tangent = Bordered()),
record_from_solution = recordPO,
plot_solution = (x, p; k...) -> begin
plotPO(x, p; k...)
## add previous branch
plot!(br, subplot=1, putbifptlegend = false)
plot!(br2, subplot=1, putbifptlegend = false)
end,
finalise_solution = (z, tau, step, contResult; prob = nothing, kwargs...) -> begin
# limit the period
return z.u[end] < 150
true
end,
normC = norminf) ┌─ Curve type: PeriodicOrbitCont from Hopf bifurcation point.
├─ Number of points: 101
├─ Type of vectors: Vector{Float64}
├─ Parameter a starts at 4.269257131580585, ends at 4.259422939031913
├─ 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, fold at a ≈ +4.27074884 ∈ (+4.27074884, +4.27074884), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 9, eigenelements in eig[ 9], ind_ev = 0
- # 2, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 35, eigenelements in eig[ 35], ind_ev = 0
- # 3, fold at a ≈ +4.25942315 ∈ (+4.25942315, +4.25942315), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 36, eigenelements in eig[ 36], ind_ev = 0
- # 4, fold at a ≈ +4.25942313 ∈ (+4.25942313, +4.25942313), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 37, eigenelements in eig[ 37], ind_ev = 0
- # 5, fold at a ≈ +4.25942315 ∈ (+4.25942315, +4.25942315), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 38, eigenelements in eig[ 38], ind_ev = 0
- # 6, fold at a ≈ +4.25942313 ∈ (+4.25942313, +4.25942313), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 39, eigenelements in eig[ 39], ind_ev = 0
- # 7, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 40, eigenelements in eig[ 40], ind_ev = 0
- # 8, fold at a ≈ +4.25942313 ∈ (+4.25942313, +4.25942313), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 41, eigenelements in eig[ 41], ind_ev = 0
- # 9, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 42, eigenelements in eig[ 42], ind_ev = 0
- # 10, fold at a ≈ +4.25942313 ∈ (+4.25942313, +4.25942313), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 43, eigenelements in eig[ 43], ind_ev = 0
- # 11, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 44, eigenelements in eig[ 44], ind_ev = 0
- # 12, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 45, eigenelements in eig[ 45], ind_ev = 0
- # 13, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 46, eigenelements in eig[ 46], ind_ev = 0
- # 14, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 47, eigenelements in eig[ 47], ind_ev = 0
- # 15, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 48, eigenelements in eig[ 48], ind_ev = 0
- # 16, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 49, eigenelements in eig[ 49], ind_ev = 0
- # 17, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 50, eigenelements in eig[ 50], ind_ev = 0
- # 18, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 51, eigenelements in eig[ 51], ind_ev = 0
- # 19, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 52, eigenelements in eig[ 52], ind_ev = 0
- # 20, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 53, eigenelements in eig[ 53], ind_ev = 0
- # 21, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 54, eigenelements in eig[ 54], ind_ev = 0
- # 22, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 55, eigenelements in eig[ 55], ind_ev = 0
- # 23, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 56, eigenelements in eig[ 56], ind_ev = 0
- # 24, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 57, eigenelements in eig[ 57], ind_ev = 0
- # 25, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 58, eigenelements in eig[ 58], ind_ev = 0
- # 26, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 59, eigenelements in eig[ 59], ind_ev = 0
- # 27, fold at a ≈ +4.25942315 ∈ (+4.25942315, +4.25942315), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 60, eigenelements in eig[ 60], ind_ev = 0
- # 28, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 61, eigenelements in eig[ 61], ind_ev = 0
- # 29, fold at a ≈ +4.25942315 ∈ (+4.25942315, +4.25942315), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 62, eigenelements in eig[ 62], ind_ev = 0
- # 30, fold at a ≈ +4.25942315 ∈ (+4.25942315, +4.25942315), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 63, eigenelements in eig[ 63], ind_ev = 0
- # 31, fold at a ≈ +4.25942315 ∈ (+4.25942315, +4.25942315), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 64, eigenelements in eig[ 64], ind_ev = 0
- # 32, fold at a ≈ +4.25942315 ∈ (+4.25942315, +4.25942315), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 65, eigenelements in eig[ 65], ind_ev = 0
- # 33, fold at a ≈ +4.25942317 ∈ (+4.25942317, +4.25942317), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 75, eigenelements in eig[ 75], ind_ev = 0
- # 34, fold at a ≈ +4.25942314 ∈ (+4.25942314, +4.25942314), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 76, eigenelements in eig[ 76], ind_ev = 0
- # 35, fold at a ≈ +4.25942315 ∈ (+4.25942315, +4.25942315), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 77, eigenelements in eig[ 77], ind_ev = 0
- # 36, fold at a ≈ +4.25942302 ∈ (+4.25942302, +4.25942302), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 82, eigenelements in eig[ 82], ind_ev = 0
- # 37, fold at a ≈ +4.25942309 ∈ (+4.25942309, +4.25942309), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 91, eigenelements in eig[ 91], ind_ev = 0
- # 38, fold at a ≈ +4.25942307 ∈ (+4.25942307, +4.25942307), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 92, eigenelements in eig[ 92], ind_ev = 0
- # 39, fold at a ≈ +4.25942308 ∈ (+4.25942308, +4.25942308), |δp|=-1e+00, [ guess], δ = ( 0, 0), step = 93, eigenelements in eig[ 93], ind_ev = 0
- # 40, endpoint at a ≈ +4.25942292, step = 101
probhom, solh = generate_hom_problem(
setproperties(br_coll.prob.prob, meshadapt=true, K = 100),
br_coll.sol[end].x,
BK.setparam(br_coll, br_coll.sol[end].p),
BK.getlens(br_coll);
update_every_step = 4,
verbose = true,
# ϵ0 = 1e-7, ϵ1 = 1e-7, # WORK BEST
# ϵ0 = 1e-8, ϵ1 = 1e-5, # maxT = 70,
t0 = 0., t1 = 120.,
# freeparams = ((@lens _.T), (@lens _.ϵ1),)
# freeparams = ((@lens _.T), (@lens _.ϵ0)),
freeparams = ((@lens _.ϵ0), (@lens _.ϵ1)), # WORK BEST
# freeparams = ((@lens _.T),),
testOrbitFlip = false,
testInclinationFlip = false
)
#####
_sol = get_homoclinic_orbit(probhom, solh, BK.getparams(probhom);)
scene = plot(_sol.t, _sol[:,:]', marker = :d, markersize = 1)optn_hom = NewtonPar(verbose = true, tol = 1e-10, max_iterations = 5)
optc_hom = ContinuationPar(newton_options = optn_hom, ds = -0.0001, dsmin = 1e-5, plot_every_step = 10, max_steps = 100, detect_bifurcation = 0, detect_event = 2, save_sol_every_step = 1, p_min = -1.01)
br_hom_c = continuation(
deepcopy(probhom), solh, (@lens _.b),
PALC(tangent = Bordered()),
# PALC(),
# MoorePenrose(),
setproperties(optc_hom, max_steps = 100, save_sol_every_step = 1, dsmax = 4e-2, plot_every_step = 10, p_max = 1.5);
verbosity = 1, plot = true,
# callback_newton = BK.cbMaxNorm(1e1),
# bothside = true,
normC = norminf,
plot_solution = (x,p;k...) -> begin
𝐇𝐨𝐦 = p.prob
par0 = set(BK.getparams(𝐇𝐨𝐦), BK.getlens(𝐇𝐨𝐦), x.x[end][1])
par0 = set(par0, (@lens _.b), p.p)
sol = get_homoclinic_orbit(𝐇𝐨𝐦, x, par0)
m = (𝐇𝐨𝐦.bvp isa PeriodicOrbitOCollProblem && 𝐇𝐨𝐦.bvp.meshadapt) ? :d : :none
plot!(sol.t, sol[:,:]',subplot=3, markersize = 1, marker=m)
end,
)
plot(sn_br, vars = (:a, :b), branchlabel = "SN", )
plot!(hopf_br, branchlabel = "AH₀", vars = (:a, :b))
plot!(hopf_br2, branchlabel = "AH₁₂", vars = (:a, :b))
plot!(br_hom_c, branchlabel = "H₀", vars = (:a, :b))
ylims!(0,1.5)Branch of homoclinic orbits with Multiple Shooting
using DifferentialEquations
probsh = ODEProblem(OPL!, copy(z0), (0., 1000.), par_OPL; abstol = 1e-12, reltol = 1e-10)
# newton parameters
optn_po = NewtonPar(verbose = true, tol = 1e-8, max_iterations = 25)
# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.075, ds= -0.001, dsmin = 1e-4, p_max = 6.8, p_min=-5., max_steps = 130, newton_options = (@set optn_po.tol = 1e-8), tol_stability = 1e-4, detect_bifurcation = 0, plot_every_step = 10, save_sol_every_step=1)
br_sh = continuation(
# br, 2,
br2, 1,
opts_po_cont,
ShootingProblem(8, probsh, Rodas5P(); parallel = true, abstol = 1e-13, reltol = 1e-11);
ampfactor = 1., δp = 0.0015,
verbosity = 2, plot = true,
record_from_solution = recordPO,
callback_newton = BK.cbMaxNorm(1e0),
plot_solution = (x, p; k...) -> begin
plotPO(x, p; k...)
## add previous branch
# plot!(br, subplot=1, putbifptlegend = false)
end,
finalise_solution = (z, tau, step, contResult; prob = nothing, kwargs...) -> begin
# limit the period
return z.u[end] < 100
true
end,
normC = norminf)
_sol = get_periodic_orbit(br_sh.prob.prob, br_sh.sol[end].x, BK.setparam(br_sh, br_sh.sol[end].p))
plot(_sol.t, _sol[:,:]')# homoclinic
probhom, solh = generate_hom_problem(
br_sh.prob.prob, br_sh.sol[end].x,
BK.setparam(br_sh, br_sh.sol[end].p),
BK.getlens(br_sh);
verbose = true,
update_every_step = 4,
# ϵ0 = 1e-6, ϵ1 = 1e-5,
t0 = 75, t1 = 25,
# freeparams = ((@lens _.T), (@lens _.ϵ1),)
# freeparams = ((@lens _.T), (@lens _.ϵ0)),
# freeparams = ((@lens _.ϵ0), (@lens _.ϵ1)), # WORK BEST
freeparams = ((@lens _.T),),
)
_sol = get_homoclinic_orbit(probhom, solh, BK.getparams(probhom); saveat=.1)
plot(plot(_sol[1,:], _sol[2,:]), plot(_sol.t, _sol[:,:]'))
optn_hom = NewtonPar(verbose = true, tol = 1e-9, max_iterations = 7)
optc_hom = ContinuationPar(newton_options = optn_hom, ds = -1e-4, dsmin = 1e-6, dsmax = 1e-3, plot_every_step = 1, max_steps = 10, detect_bifurcation = 0, save_sol_every_step = 1)
br_hom_sh = continuation(
deepcopy(probhom), solh, (@lens _.b),
PALC(),
setproperties(optc_hom, max_steps = 600, save_sol_every_step = 1, dsmax = 12e-2, plot_every_step = 3, p_max = 7., detect_event = 2, a = 0.9);
verbosity = 3, plot = true,
callback_newton = BK.cbMaxNorm(1e0),
normC = norminf,
plot_solution = (x,p;k...) -> begin
𝐇𝐨𝐦 = p.prob
par0 = set(BK.getparams(𝐇𝐨𝐦), BK.getlens(𝐇𝐨𝐦), x.x[end][1])
par0 = set(par0, (@lens _.b), p.p)
sol = get_homoclinic_orbit(𝐇𝐨𝐦, x, par0)
m = (𝐇𝐨𝐦.bvp isa PeriodicOrbitOCollProblem && 𝐇𝐨𝐦.bvp.meshadapt) ? :d : :none
plot!(sol.t, sol[:,:]',subplot=3, markersize = 1, marker=m)
end,
)
_sol = get_homoclinic_orbit(probhom, br_hom_sh.sol[end].x, BK.setparam(br_hom_sh, br_hom_sh.sol[end].p); saveat=.1)
plot(_sol.t, _sol[:,:]')plot(sn_br, vars = (:a, :b), branchlabel = "SN", )
plot!(hopf_br, branchlabel = "AH₀", vars = (:a, :b))
plot!(hopf_br2, branchlabel = "AH₁₂", vars = (:a, :b))
plot!(br_hom_c, branchlabel = "Hc₀", vars = (:a, :b))
plot!(br_hom_sh, branchlabel = "Hsh₀", vars = (:a, :b), linewidth = 3)
ylims!(0,1.5)References
- Pusuluri
Pusuluri, K, H G E Meijer, and A L Shilnikov. “Homoclinic Puzzles and Chaos in a Nonlinear Laser Model,” n.d.