π‘ Extended Lorenz-84 model (codim 2 + BT/ZH aBS)
In this tutorial, we study the extended Lorenz-84 model which is also treated in MatCont [Kuznetsov]. This model is interesting because it features all codim 2 bifurcations of equilibria. It is thus convenient to test our algorithms.
After this tutorial, you will be able to
- detect codim 1 bifurcation Fold / Hopf / Branch point
- follow Fold / Hopf points and detect codim 2 bifurcation points
- branch from the codim 2 points to curves of Fold / Hopf points
The model is as follows
\[\left\{\begin{array}{l} \dot{X}=-Y^{2}-Z^{2}-\alpha X+\alpha F-\gamma U^{2} \\ \dot{Y}=X Y-\beta X Z-Y+G \\ \dot{Z}=\beta X Y+X Z-Z \\ \dot{U}=-\delta U+\gamma U X+T \end{array}\right.\tag{E}\]
We start with some imports:
using Revise, Plots
using BifurcationKit
const BK = BifurcationKit
Problem setting
We can now encode the vector field (E) in a function.
# vector field
function Lor(u, p)
(;Ξ±,Ξ²,Ξ³,Ξ΄,G,F,T) = p
X,Y,Z,U = u
[
-Y^2 - Z^2 - Ξ±*X + Ξ±*F - Ξ³*U^2,
X*Y - Ξ²*X*Z - Y + G,
Ξ²*X*Y + X*Z - Z,
-Ξ΄*U + Ξ³*U*X + T
]
end
# parameter values
parlor = (Ξ± = 1//4, Ξ² = 1, G = .25, Ξ΄ = 1.04, Ξ³ = 0.987, F = 1.7620532879639, T = .0001265)
# initial condition
z0 = [2.9787004394953343, -0.03868302503393752, 0.058232737694740085, -0.02105288273117459]
# bifurcation problem
recordFromSolutionLor(x, p; k...) = (X = x[1], Y = x[2], Z = x[3], U = x[4])
prob = BifurcationProblem(Lor, z0, (parlor..., T=0.04, F=3.), (@optic _.F);
record_from_solution = recordFromSolutionLor)
Continuation and codim 1 bifurcations
Once the problem is set up, we can continue the state w.r.t. $F$ and detect codim 1 bifurcations. This is achieved as follows:
# continuation options
opts_br = ContinuationPar(p_min = -1.5, p_max = 3.0, ds = 0.002, dsmax = 0.15,
# Optional: bisection options for locating bifurcations
n_inversion = 6,
# number of eigenvalues
nev = 4)
# compute the branch of solutions
br = continuation(prob, PALC(), opts_br;
normC = norminf,
bothside = true)
scene = plot(br, plotfold = false, markersize = 4, legend = :topleft)
With detailed information:
br
ββ Curve type: EquilibriumCont
ββ Number of points: 33
ββ Type of vectors: Vector{Float64}
ββ Parameter F starts at 3.0, ends at 3.0
ββ Algo: PALC
ββ Special points:
- # 1, endpoint at F β +3.00000000, step = 0
- # 2, hopf at F β +2.85996783 β (+2.85986480, +2.85996783), |Ξ΄p|=1e-04, [converged], Ξ΄ = ( 2, 2), step = 1
- # 3, hopf at F β +2.46723305 β (+2.46720734, +2.46723305), |Ξ΄p|=3e-05, [converged], Ξ΄ = (-2, -2), step = 3
- # 4, hopf at F β +1.61975642 β (+1.61959602, +1.61975642), |Ξ΄p|=2e-04, [converged], Ξ΄ = ( 2, 2), step = 9
- # 5, bp at F β +1.54664839 β (+1.54664837, +1.54664839), |Ξ΄p|=1e-08, [converged], Ξ΄ = (-1, 0), step = 11
- # 6, endpoint at F β +3.00000000, step = 32
Continuation of Fold points
We follow the Fold points in the parameter plane $(T,F)$. We tell the solver to consider br.specialpoint[5]
and continue it.
# function to record the current state
sn_codim2 = continuation(br, 5, (@optic _.T),
ContinuationPar(opts_br, p_max = 3.2, p_min = -0.1,
dsmin=1e-5, ds = -0.001, dsmax = 0.005) ;
normC = norminf,
# detection of codim 2 bifurcations with bisection
detect_codim2_bifurcation = 2,
# we save the different components for plotting
record_from_solution = recordFromSolutionLor,
)
scene = plot(sn_codim2, vars=(:X, :U), branchlabel = "Folds", ylims=(-0.5, 0.5))
with detailed information
sn_codim2
ββ Curve type: FoldCont
ββ Number of points: 82
ββ Type of vectors: Vector{Float64}
ββ Parameters (:F, :T)
ββ Parameter T starts at 0.04, ends at -0.1
ββ Algo: PALC
ββ Special points:
- # 1, bt at T β +0.02094014 β (+0.02094014, +0.02094020), |Ξ΄p|=6e-08, [converged], Ξ΄ = ( 0, 0), step = 12
- # 2, zh at T β +0.00012644 β (+0.00012644, +0.00012666), |Ξ΄p|=2e-07, [converged], Ξ΄ = ( 0, 0), step = 29
- # 3, zh at T β -0.00012655 β (-0.00012655, -0.00012644), |Ξ΄p|=1e-07, [converged], Ξ΄ = ( 0, 0), step = 32
- # 4, bt at T β -0.02094041 β (-0.02094041, -0.02093949), |Ξ΄p|=9e-07, [converged], Ξ΄ = ( 0, 0), step = 49
- # 5, endpoint at T β -0.10000000, step = 81
For example, we can compute the following normal form
get_normal_form(sn_codim2, 1; nev = 4)
Bogdanov-Takens bifurcation point at (:F, :T) β (1.4467165285461494, 0.020940139624099376).
Normal form (B, Ξ²1 + Ξ²2β
B + bβ
Aβ
B + aβ
AΒ²)
Normal form coefficients:
a = 0.2144232743197446
b = 0.6065142321261716
You can call various predictors:
- predictor(::BogdanovTakens, ::Val{:HopfCurve}, ds)
- predictor(::BogdanovTakens, ::Val{:FoldCurve}, ds)
- predictor(::BogdanovTakens, ::Val{:HomoclinicCurve}, ds)
Continuation of Hopf points
We follow the Hopf points in the parameter plane $(T,F)$. We tell the solver to consider br.specialpoint[3]
and continue it.
hp_codim2_1 = continuation(br, 3, (@optic _.T),
ContinuationPar(opts_br, ds = -0.001, dsmax = 0.02, dsmin = 1e-4) ;
normC = norminf,
# detection of codim 2 bifurcations with bisection
detect_codim2_bifurcation = 2,
# we save the different components for plotting
record_from_solution = recordFromSolutionLor,
# compute both sides of the initial condition
bothside = true,
)
plot(sn_codim2, vars=(:X, :U), branchlabel = "Folds")
plot!(hp_codim2_1, vars=(:X, :U), branchlabel = "Hopfs")
ylims!(-0.7,0.7);xlims!(1,1.3)
hp_codim2_1
ββ Curve type: HopfCont
ββ Number of points: 429
ββ Type of vectors: Vector{Float64}
ββ Parameters (:F, :T)
ββ Parameter T starts at 0.020940169755227573, ends at -0.1535694927465657
ββ Algo: PALC
ββ Special points:
- # 1, endpoint at T β +0.02094017, step = 0
- # 2, bt at T β +0.02094017 β (+0.02094017, +0.02094017), |Ξ΄p|=6e-11, [converged], Ξ΄ = ( 0, 0), step = 0
- # 3, gh at T β +0.05019751 β (+0.05019655, +0.05019751), |Ξ΄p|=1e-06, [converged], Ξ΄ = ( 0, 0), step = 19
- # 4, hh at T β +0.02627340 β (+0.02627340, +0.02627528), |Ξ΄p|=2e-06, [converged], Ξ΄ = (-2, -2), step = 35
- # 5, endpoint at T β -0.15372759, step = 429
For example, we can compute the following normal form
get_normal_form(hp_codim2_1, 3; nev = 4)
Bautin bifurcation point at (:F, :T) β (2.3763590366726284, 0.05019751302158527).
Ο = 0.6903670769045964
Second lyapunov coefficient lβ = 0.1557753180139064
Normal form: iβ
Οβ
z + lββ
zβ
|z|β΄
Continuation of Hopf points from the Bogdanov-Takens point
When we computed the curve of Fold points, we detected a Bogdanov-Takens bifurcation. We can branch from it to get the curve of Hopf points. This is done as follows:
hp_from_bt = continuation(sn_codim2, 4,
ContinuationPar(opts_br, ds = -0.001, dsmax = 0.02, dsmin = 1e-4) ;
normC = norminf,
# detection of codim 2 bifurcations with bisection
detect_codim2_bifurcation = 2,
# we save the different components for plotting
record_from_solution = recordFromSolutionLor,
)
plot(sn_codim2, vars=(:X, :U), branchlabel = "SN")
plot!(hp_codim2_1, vars=(:X, :U), branchlabel = "Hopf1")
plot!(hp_from_bt, vars=(:X, :U), branchlabel = "Hopf2")
ylims!(-0.7,0.75); xlims!(0.95,1.3)
with detailed information
hp_from_bt
ββ Curve type: HopfCont from BogdanovTakens bifurcation point.
ββ Number of points: 401
ββ Type of vectors: Vector{Float64}
ββ Parameters (:F, :T)
ββ Parameter T starts at -0.026824826717918287, ends at 0.15063717399994617
ββ Algo: PALC
ββ Special points:
- # 1, gh at T β -0.05018361 β (-0.05022014, -0.05018361), |Ξ΄p|=4e-05, [converged], Ξ΄ = ( 0, 0), step = 23
- # 2, hh at T β -0.02626030 β (-0.02631231, -0.02626030), |Ξ΄p|=5e-05, [converged], Ξ΄ = (-2, -2), step = 26
- # 3, endpoint at T β +0.15079813, step = 401
Continuation of Hopf points from the Zero-Hopf point
When we computed the curve of Fold points, we detected a Zero-Hopf bifurcation. We can branch from it to get the curve of Hopf points. This is done as follows:
hp_from_zh = continuation(sn_codim2, 2,
ContinuationPar(opts_br, ds = 0.001, dsmax = 0.02) ;
normC = norminf,
detect_codim2_bifurcation = 2,
record_from_solution = recordFromSolutionLor,
)
plot(hp_codim2_1, vars=(:X, :U), branchlabel = "Hopf")
plot!(hp_from_bt, vars=(:X, :U), branchlabel = "Hopf2")
plot!( hp_from_zh, vars=(:X, :U), branchlabel = "Hopf", legend = :topleft)
plot!(sn_codim2,vars=(:X, :U),)
ylims!(-0.7,0.75); xlims!(0.95,1.3)
with detailed information
hp_from_zh
ββ Curve type: HopfCont from BifurcationKit.ZeroHopf bifurcation point.
ββ Number of points: 401
ββ Type of vectors: Vector{Float64}
ββ Parameters (:F, :T)
ββ Parameter T starts at 0.0001264399489437422, ends at 0.666963445643025
ββ Algo: PALC
ββ Special points:
- # 1, gh at T β +0.00012660 β (+0.00012654, +0.00012660), |Ξ΄p|=6e-08, [converged], Ξ΄ = ( 0, 0), step = 1
- # 2, hh at T β +0.02627444 β (+0.02627324, +0.02627444), |Ξ΄p|=1e-06, [converged], Ξ΄ = ( 2, 2), step = 27
- # 3, endpoint at T β +0.66778051, step = 401
References
- Kuznetsov
Kuznetsov, Yu A., H. G. E. Meijer, W. Govaerts, and B. Sautois. βSwitching to Nonhyperbolic Cycles from Codim 2 Bifurcations of Equilibria in ODEs.β Physica D: Nonlinear Phenomena 237, no. 23 (December 2008): 3061β68.