Extended Lorenz-84 model (codim 2 + BT/ZH aBS)

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. 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 detected codim 2 points to curves of Fold / Hopf points (This part is still "work in progress")

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 that are useful in the following.

using Revise, ForwardDiff, Parameters, Setfield, Plots, LinearAlgebra
using BifurcationKit
const BK = BifurcationKit

# sup norm
norminf(x) = norm(x, Inf)

Problem setting

We can now encode the vector field (E) in a function and use automatic differentiation to compute its various derivatives.

# vector field
function Lor(u, p)
    @unpack α,β,γ,δ,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]

Continuation and codim 1 bifurcations

Once the problem is set up, we can continue the state w.r.t. $F$ to and detect codim 1 bifurcations. This is achieved as follows:

# bifurcation problem
recordFromSolutionLor(x, p) = (X = x[1], Y = x[2], Z = x[3], U = x[4])
prob = BifurcationProblem(Lor, z0, setproperties(parlor; T=0.04, F=3.), (@lens _.F);
    recordFromSolution = recordFromSolutionLor)

# continuation options
opts_br = ContinuationPar(pMin = -1.5, pMax = 3.0, ds = 0.002, dsmax = 0.15,
	# Optional: bisection options for locating bifurcations
	nInversion = 6, maxBisectionSteps = 25,
	# number of eigenvalues
	nev = 4, maxSteps = 200)

# compute the branch of solutions
br = @time 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:

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 F ≈ +3.00000000,                                                                     step =   0
- #  2,     hopf at F ≈ +2.85996783 ∈ (+2.85986480, +2.85996783), |δp|=1e-04, [converged], δ = ( 2,  2), step =   1, eigenelements in eig[  2], ind_ev =   4
- #  3,     hopf at F ≈ +2.46723305 ∈ (+2.46720734, +2.46723305), |δp|=3e-05, [converged], δ = (-2, -2), step =   3, eigenelements in eig[  4], ind_ev =   4
- #  4,     hopf at F ≈ +1.61975642 ∈ (+1.61959602, +1.61975642), |δp|=2e-04, [converged], δ = ( 2,  2), step =   9, eigenelements in eig[ 10], ind_ev =   4
- #  5,       bp at F ≈ +1.54664839 ∈ (+1.54664837, +1.54664839), |δp|=1e-08, [converged], δ = (-1,  0), step =  11, eigenelements in eig[ 12], ind_ev =   4
- #  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, (@lens _.T), ContinuationPar(opts_br, pMax = 3.2, pMin = -0.1, detectBifurcation = 1, dsmin=1e-5, ds = -0.001, dsmax = 0.005, nInversion = 10, maxSteps = 130, maxBisectionSteps = 55) ; normC = norminf,
	# detection of codim 2 bifurcations with bisection
	detectCodim2Bifurcation = 2,
	# we update the Fold problem at every continuation step
	updateMinAugEveryStep = 1,
	startWithEigen = false,
	# we save the different components for plotting
	recordFromSolution = 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}
 ├─ Parameter T starts at 0.04, ends at -0.1
 ├─ 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,       bt at T ≈ +0.02094017 ∈ (+0.02094017, +0.02094017), |δp|=3e-11, [converged], δ = ( 0,  0), step =  12, eigenelements in eig[ 13], ind_ev =   0
- #  2,       zh at T ≈ +0.00012654 ∈ (+0.00012654, +0.00012654), |δp|=9e-10, [converged], δ = ( 0,  0), step =  29, eigenelements in eig[ 30], ind_ev =   0
- #  3,       zh at T ≈ -0.00012654 ∈ (-0.00012654, -0.00012654), |δp|=1e-11, [converged], δ = ( 0,  0), step =  32, eigenelements in eig[ 33], ind_ev =   0
- #  4,       bt at T ≈ -0.02094018 ∈ (-0.02094018, -0.02094017), |δp|=7e-09, [converged], δ = ( 0,  0), step =  49, eigenelements in eig[ 50], ind_ev =   0
- #  5, endpoint at T ≈ -0.10000000,                                                                     step =  81

For example, we can compute the following normal form

getNormalForm(sn_codim2, 1; nev = 4)

Bogdanov-Takens bifurcation point at (:F, :T) ≈ (1.4467167009620117, 0.020940169656439543).
Normal form (B, β1 + β2⋅B + b⋅A⋅B + a⋅A²)
Normal form coefficients:
 a = 0.21442335085970526
 b = 0.6065145515450641

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((@set br.alg.tangent = Bordered()), 3, (@lens _.T), ContinuationPar(opts_br, ds = -0.001, dsmax = 0.02, dsmin = 1e-4, nInversion = 6, detectBifurcation = 1) ; normC = norminf,
	# detection of codim 2 bifurcations with bisection
	detectCodim2Bifurcation = 2,
	# we update the Fold problem at every continuation step
	updateMinAugEveryStep = 1,
	# we save the different components for plotting
	recordFromSolution = 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: 229
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter T starts at 0.020940169778140574, ends at -0.11830225185208443
 ├─ 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 T ≈ +0.02094017,                                                                     step =   0
- #  2,       bt at T ≈ +0.02094017 ∈ (+0.02094017, +0.02094017), |δp|=1e-10, [converged], δ = ( 0,  0), step =   0, eigenelements in eig[  1], ind_ev =   0
- #  3,       gh at T ≈ +0.05019747 ∈ (+0.05019363, +0.05019747), |δp|=4e-06, [converged], δ = ( 0,  0), step =  19, eigenelements in eig[ 20], ind_ev =   0
- #  4,       hh at T ≈ +0.02627369 ∈ (+0.02627369, +0.02627462), |δp|=9e-07, [converged], δ = (-2, -2), step =  35, eigenelements in eig[ 36], ind_ev =   2
- #  5, endpoint at T ≈ -0.11849955,                                                                     step = 229

For example, we can compute the following normal form

getNormalForm(hp_codim2_1, 3; nev = 4)

Bautin bifurcation point at (:F, :T) ≈ (2.376359555701081, 0.05019747303584815).
ω = 0.6903672728792098
Second lyapunov coefficient l₂ = 0.1557880753426296
Normal form: i⋅ω⋅u + l₂⋅u⋅|u|⁴

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((@set sn_codim2.alg.tangent = Bordered()), 4, ContinuationPar(opts_br, ds = -0.001, dsmax = 0.02, dsmin = 1e-4,
	nInversion = 6, detectBifurcation = 1) ; normC = norminf,
	# detection of codim 2 bifurcations with bisection
	detectCodim2Bifurcation = 2,
	# we update the Fold problem at every continuation step
	updateMinAugEveryStep = 1,
	# we save the different components for plotting
	recordFromSolution = 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: 201
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter T starts at -0.026822590800602113, ends at 0.11463598241572592
 ├─ 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,       gh at T ≈ -0.05019669 ∈ (-0.05020126, -0.05019669), |δp|=5e-06, [converged], δ = ( 0,  0), step =  23, eigenelements in eig[ 24], ind_ev =   0
- #  2,       hh at T ≈ -0.02627323 ∈ (-0.02627485, -0.02627323), |δp|=2e-06, [converged], δ = (-2, -2), step =  26, eigenelements in eig[ 27], ind_ev =   2
- #  3, endpoint at T ≈ +0.11483750,                                                                     step = 201

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((@set sn_codim2.alg.tangent = Bordered()), 2, ContinuationPar(opts_br, ds = 0.001, dsmax = 0.02, dsmin = 1e-4, nInversion = 6, detectBifurcation = 1, maxSteps = 150) ;
	normC = norminf,
	detectCodim2Bifurcation = 2,
	updateMinAugEveryStep = 1,
	startWithEigen = true,
	recordFromSolution = recordFromSolutionLor,
	bothside = false,
	bdlinsolver = MatrixBLS(),
	)

plot(sn_codim2,vars=(:X, :U),)
	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", plotspecialpoints = false, legend = :topleft)
	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: 151
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter T starts at 0.0001265405285364602, ends at 0.41337175987047703
 ├─ 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,       gh at T ≈ +0.00012653 ∈ (+0.00012653, +0.00012653), |δp|=9e-13, [   guessL], δ = ( 0,  0), step =   1, eigenelements in eig[  2], ind_ev =   0
- #  2,       zh at T ≈ +0.00012654 ∈ (+0.00012653, +0.00012654), |δp|=1e-08, [converged], δ = (-1,  0), step =   2, eigenelements in eig[  3], ind_ev =   1
- #  3,       hh at T ≈ +0.02627399 ∈ (+0.02627369, +0.02627399), |δp|=3e-07, [converged], δ = ( 2,  2), step =  27, eigenelements in eig[ 28], ind_ev =   2
- #  4, endpoint at T ≈ +0.41476034,                                                                     step = 151

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.