🟠 Lorenz-84 model, take 2

• 🟠 Lorenz-84 model, take 2
• • Two parameters curves of Fold / Hopf bifurcation
  • Fold bifurcations of periodic orbits from Bautin bifurcation
  • NS bifurcations of periodic orbits from Hopf-Hopf bifurcation
  • References

In this tutorial, we study the extended Lorenz-84 model which is also treated in MatCont ^[Kuznetsov]. We use this model to showcase the automatic branch switching procedure to

• Fold of periodic orbits from Bautin bifurcation point
• NS of periodic orbits from ZH bifurcation point
• NS of periodic orbits from HH bifurcation point.

As this model has been studied in this tutorial, we do not give much details and refer to the corresponding tutorial to get the 2 parameters curves of Hopf / Fold bifurcations.

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 recall the problem setting:

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

# vector field
function Lor(u, p, t = 0)
	@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

parlor = (α = 1//4, β = 1., G = .25, δ = 1.04, γ = 0.987, F = 1.7620532879639, T = .0001265)

z0 = [2.9787004394953343, -0.03868302503393752,  0.058232737694740085, -0.02105288273117459]

recordFromSolutionLor(x, p) = (u = BK.getVec(x);(X = u[1], Y = u[2], Z = u[3], U = u[4]))
prob = BK.BifurcationProblem(Lor, z0, parlor, (@lens _.F);
	record_from_solution = (x, p) -> (X = x[1], Y = x[2], Z = x[3], U = x[4]),)

opts_br = ContinuationPar(p_min = -1.5, p_max = 3.0, ds = 0.002, dsmax = 0.05, n_inversion = 6, nev = 4, max_steps = 200)
@set! opts_br.newton_options.tol = 1e-12
br = @time continuation(re_make(prob, params = setproperties(parlor;T=0.04,F=3.)),
	 	PALC(), opts_br;
		normC = norminf, bothside = true)

scene = plot(br, plotfold=false, markersize=4, legend=:topleft)

Two parameters curves of Fold / Hopf bifurcation

We follow the Fold points in the parameter plane $(T,F)$. We tell the solver to consider br.specialpoint[5] and continue it.

sn_codim2 = continuation(br, 5, (@lens _.T), ContinuationPar(opts_br, p_max = 3.2, p_min = -0.1, detect_bifurcation = 1, dsmin=1e-5, ds = -0.001, dsmax = 0.005, n_inversion = 10, max_steps = 130) ;
	# plot = true,
	normC = norminf,
	detect_codim2_bifurcation = 2,
	start_with_eigen = false,
	bothside = false,
	)

hp_codim2_1 = continuation(br, 3, (@lens _.T), ContinuationPar(opts_br, ds = -0.001, dsmax = 0.02, dsmin = 1e-4, n_inversion = 8, save_sol_every_step = 1, detect_bifurcation = 1) ;
	# plot = false, verbosity = 0,
	normC = norminf,
	detect_codim2_bifurcation = 2,
	start_with_eigen = true,
	bothside = true,
	)

plot(sn_codim2, vars=(:F, :T), branchlabel = "SN")
plot!(hp_codim2_1, vars=(:F, :T), branchlabel = "Hopf1", xlims = (1,2.7), ylims = (-0.06,0.06))

Fold bifurcations of periodic orbits from Bautin bifurcation

We compute the branch of Fold of periodic orbits from the Bautin bifurcation (labelled :gh) in the previous figure. In this tutorial, we focus on orthogonal collocation but standard shooting would do too.

opts_fold_po = ContinuationPar(hp_codim2_1.contparams, dsmax = 0.01, detect_bifurcation = 0, max_steps = 30, detect_event = 0, ds = 0.001, plot_every_step = 10)
# @set! opts_fold_po.newton_options.verbose = false
@set! opts_fold_po.newton_options.tol = 1e-8
fold_po = continuation(hp_codim2_1, 3, opts_fold_po,
		PeriodicOrbitOCollProblem(20, 3, meshadapt = false);
		normC = norminf,
		δp = 0.02,
		update_minaug_every_step = 0,
		jacobian_ma = :minaug,
		verbosity = 0, plot = false,
	)
plot!(fold_po, vars=(:F, :T), branchlabel = "Fold-PO")

NS bifurcations of periodic orbits from Hopf-Hopf bifurcation

When we computed the curve of Hopf points, we detected a Hopf-Hopf bifurcation. We can branch from it to get the curve of NS points. This is done as follows:

opts_ns_po = ContinuationPar(hp_codim2_1.contparams, dsmax = 0.02, detect_bifurcation = 1, max_steps = 20, ds = -0.001, detect_event = 0)
@set! opts_ns_po.newton_options.verbose = false
@set! opts_ns_po.newton_options.tol = 1e-9
@set! opts_ns_po.newton_options.max_iterations = 10
ns_po1 = continuation(hp_codim2_1, 4, opts_ns_po,
		PeriodicOrbitOCollProblem(20, 3, update_section_every_step = 1);
		detect_codim2_bifurcation = 0,
		normC = norminf,
		δp = 0.02,
		update_minaug_every_step = 1,
		# which of the 2 NS curves should we compute?
		whichns = 1,
		jacobian_ma = :minaug,
		)
plot!(ns_po1, vars=(:F, :T), branchlabel = "NS1")

ns_po2 = continuation(hp_codim2_1, 4, opts_ns_po,
		PeriodicOrbitOCollProblem(30, 3, update_section_every_step = 1);
		detect_codim2_bifurcation = 0,
		normC = norminf,
		δp = 0.02,
		update_minaug_every_step = 1,
		# which of the 2 NS curves should we compute?
		whichns = 2,
		jacobian_ma = :minaug,
		)
plot!(ns_po2, vars=(:F, :T), branchlabel = "NS2")

References