ο‚›GitHub

🟑 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)
Example block output

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,
	# 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))
Example block output

with detailed information

sn_codim2
 β”Œβ”€ Curve type: FoldCont
 β”œβ”€ Number of points: 81
 β”œβ”€ Type of vectors: Vector{Float64}
 β”œβ”€ Parameters (:F, :T)
 β”œβ”€ Parameter T starts at 0.04, ends at -0.1
 β”œβ”€ Algo: PALC
 └─ Special points:

- #  1,       zh at T β‰ˆ +0.02094014 ∈ (+0.02094014, +0.02094020), |Ξ΄p|=6e-08, [converged], Ξ΄ = (-2, -2), step =  12
- #  2,       zh at T β‰ˆ +0.00008001 ∈ (+0.00008001, +0.00019454), |Ξ΄p|=1e-04, [   guessL], Ξ΄ = (-2, -2), step =  29
- #  3,       zh at T β‰ˆ -0.00037404 ∈ (-0.00037404, -0.00006728), |Ξ΄p|=3e-04, [    guess], Ξ΄ = ( 2,  2), step =  32
- #  4,       bt at T β‰ˆ -0.02094024 ∈ (-0.02094024, -0.02094012), |Ξ΄p|=1e-07, [converged], Ξ΄ = ( 1,  0), step =  48
- #  5, endpoint at T β‰ˆ -0.10000000,                                                                     step =  80

For example, we can compute the following normal form

get_normal_form(sn_codim2, 1; nev = 4)

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,
	# 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)
Example block output
hp_codim2_1
 β”Œβ”€ Curve type: HopfCont
 β”œβ”€ Number of points: 802
 β”œβ”€ Type of vectors: Vector{Float64}
 β”œβ”€ Parameters (:F, :T)
 β”œβ”€ Parameter T starts at -0.14706943151472426, ends at -0.15357975432787585
 β”œβ”€ Algo: PALC
 └─ Special points:

- #  1, endpoint at T β‰ˆ -0.14723397,                                                                     step =  -1
- #  2,       hh at T β‰ˆ +0.02509172 ∈ (+0.02509172, +0.02676438), |Ξ΄p|=2e-03, [    guess], Ξ΄ = ( 2,  2), step = 352
- #  3,       hh at T β‰ˆ +0.02541156 ∈ (+0.02541156, +0.02637641), |Ξ΄p|=1e-03, [   guessL], Ξ΄ = (-2, -2), step = 408
- #  4, endpoint at T β‰ˆ -0.15373784,                                                                     step = 802

For example, we can compute the following normal form

get_normal_form(hp_codim2_1, 3; nev = 4)
Hopf-Hopf bifurcation point at (:F, :T) β‰ˆ (2.5365165590507286, 0.025411558557093924).
Eigenvalues:
Ξ»1 = -0.0024602367795248936 + 1.1506787857735017im
Ξ»2 = 1.3877787807814457e-16 + 0.744042999195727im
(Ξ»1 = -0.0024602367795248936 + 1.1506787857735017im, Ξ»2 = 1.3877787807814457e-16 + 0.744042999195727im, G2100 = 0.4951781219456526 + 0.13313738617258042im, G0021 = -0.3090832788495159 + 0.21218934917386378im, G1110 = -0.2734499471893299 - 1.423043698979615im, G1011 = 0.5470796102245956 - 1.2834479968585375im, γ₁₁₀ = -0.09108399332666783 + 0.018467722649667125im, γ₁₀₁ = 2.5015469319268537 + 1.0639933831760022im, γ₂₁₀ = 0.06282321799322829 + 0.28482215793049337im, γ₂₀₁ = 0.2381109605359213 + 0.14474134227255864im, Ξ“ = ComplexF64[-0.09108399332666783 + 0.018467722649667125im 2.5015469319268537 + 1.0639933831760022im; 0.06282321799322829 + 0.28482215793049337im 0.2381109605359213 + 0.14474134227255864im], h₁₁₀₀ = ComplexF64[0.15639669587604282 + 0.0im, -0.4078006629310759 - 0.0im, 0.5266012749665344 + 0.0im, 1.7278257922595728 - 0.0im], h₀₀₁₁ = ComplexF64[0.1646051265999119 + 0.0im, 0.5989620461000718 - 0.0im, 0.2388557308231652 + 0.0im, 0.9837726476260619 - 0.0im], h₀₀₀₀₁₀ = [-0.01892251445158322, 0.003117575650809403, 0.004514013862905299, -0.22929844177494238], h₀₀₀₀₀₁ = [1.834387607942936, -0.30222421841992114, -0.4375978210180883, 0.24817999904667404], hβ‚‚β‚€β‚€β‚€ = ComplexF64[-0.133952644355269 + 0.21632751142487133im, -0.5806580497621784 + 0.3393970912126489im, 0.36995264621629254 + 0.5510187793704826im, -0.12350863974757664 - 0.06653434693519789im], hβ‚€β‚€β‚‚β‚€ = ComplexF64[-0.5515330897904689 - 1.2608664684096849im, -1.3215029793411448 - 0.7025328440759735im, -0.7258412756696685 + 1.5826985322670368im, 0.9644588375127134 - 0.26471965520978896im], ns1 = (dΟ‰1 = 0.19893382741252413, dΟ‰2 = -0.23218568465589606, Ξ± = [-4.154482947777161, -0.05229477586394897]), ns2 = (dΟ‰1 = -1.3508055370492358, dΟ‰2 = 0.9128097263194911, Ξ± = [-2.8900089750776874, 0.1134682097712925]))

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)
Example block output

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.026823959138524478, ends at 0.1506371396616536
 β”œβ”€ Algo: PALC
 └─ Special points:

- #  1,       gh at T β‰ˆ -0.05018460 ∈ (-0.05022113, -0.05018460), |Ξ΄p|=4e-05, [converged], Ξ΄ = ( 0,  0), step =  23
- #  2,       hh at T β‰ˆ -0.02626316 ∈ (-0.02631517, -0.02626316), |Ξ΄p|=5e-05, [converged], Ξ΄ = (-2, -2), step =  26
- #  3, endpoint at T β‰ˆ +0.15079809,                                                                     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)
Example block output

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 8.000500426234362e-5, ends at 0.6669650656671116
 β”œβ”€ Algo: PALC
 └─ Special points:

- #  1,       zh at T β‰ˆ +0.00012654 ∈ (+0.00012652, +0.00012654), |Ξ΄p|=2e-08, [converged], Ξ΄ = (-1,  0), step =   2
- #  2,       hh at T β‰ˆ +0.02627930 ∈ (+0.02624096, +0.02627930), |Ξ΄p|=4e-05, [converged], Ξ΄ = ( 2,  2), step =  27
- #  3, endpoint at T β‰ˆ +0.66778213,                                                                     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.