# π  1d Brusselator

Goal

The goal of this tutorial is to show similar computations as in the previous tutorial but without using the automatic branch switching tools. This is for the experienced used who wants to dive more in the internals of the package.

We look at the Brusselator in 1d (see [Lust]). The equations are as follows

\begin{aligned} \frac { \partial X } { \partial t } & = \frac { D _ { 1 } } { l ^ { 2 } } \frac { \partial ^ { 2 } X } { \partial z ^ { 2 } } + X ^ { 2 } Y - ( Ξ² + 1 ) X + Ξ± \\ \frac { \partial Y } { \partial t } & = \frac { D _ { 2 } } { l ^ { 2 } } \frac { \partial ^ { 2 } Y } { \partial z ^ { 2 } } + Ξ² X - X ^ { 2 } Y \end{aligned}

with Dirichlet boundary conditions

$$$\begin{array} { l } { X ( t , z = 0 ) = X ( t , z = 1 ) = Ξ± } \\ { Y ( t , z = 0 ) = Y ( t , z = 1 ) = Ξ² / Ξ± } \end{array}$$$

These equations have been introduced to reproduce an oscillating chemical reaction. There is an obvious equilibrium $(Ξ±, Ξ² / Ξ±)$. Here, we consider bifurcations with respect to the parameter $l$.

We start by writing the PDE

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

f1(u, v) = u * u * v

function Fbru!(f, x, p, t = 0)
@unpack Ξ±, Ξ², D1, D2, l = p
n = div(length(x), 2)
h2 = 1.0 / n^2
c1 = D1 / l^2 / h2
c2 = D2 / l^2 / h2

u = @view x[1:n]
v = @view x[n+1:2n]

# Dirichlet boundary conditions
f[1]   = c1 * (Ξ±	  - 2u[1] + u[2] ) + Ξ± - (Ξ² + 1) * u[1] + f1(u[1], v[1])
f[end] = c2 * (v[n-1] - 2v[n] + Ξ² / Ξ±)			 + Ξ² * u[n] - f1(u[n], v[n])

f[n]   = c1 * (u[n-1] - 2u[n] +  Ξ±   ) + Ξ± - (Ξ² + 1) * u[n] + f1(u[n], v[n])
f[n+1] = c2 * (Ξ² / Ξ±  - 2v[1] + v[2])			 + Ξ² * u[1] - f1(u[1], v[1])

for i=2:n-1
f[i] = c1 * (u[i-1] - 2u[i] + u[i+1]) + Ξ± - (Ξ² + 1) * u[i] + f1(u[i], v[i])
f[n+i] = c2 * (v[i-1] - 2v[i] + v[i+1])			  + Ξ² * u[i] - f1(u[i], v[i])
end
return f
end

Fbru(x, p, t = 0) = Fbru!(similar(x), x, p, t)

For computing periodic orbits, we will need a Sparse representation of the Jacobian:

function Jbru_sp(x, p)
@unpack Ξ±, Ξ², D1, D2, l = p
# compute the Jacobian using a sparse representation
n = div(length(x), 2)
π― = eltype(x)
h = 1.0 / n; h2 = h*h

c1 = D1 / p.l^2 / h2
c2 = D2 / p.l^2 / h2

u = @view x[1:n]
v = @view x[n+1:2n]

diag   = zeros(π―, 2n)
diagp1 = zeros(π―, 2n-1)
diagm1 = zeros(π―, 2n-1)

diagpn = zeros(π―, n)
diagmn = zeros(π―, n)

@. diagmn = Ξ² - 2 * u * v
@. diagm1[1:n-1] = c1
@. diagm1[n+1:end] = c2

@. diag[1:n]    = -2c1 - (Ξ² + 1) + 2 * u * v
@. diag[n+1:2n] = -2c2 - u * u

@. diagp1[1:n-1] = c1
@. diagp1[n+1:end] = c2

@. diagpn = u * u
return spdiagm(0 => diag, 1 => diagp1, -1 => diagm1, n => diagpn, -n => diagmn)
end
Tip

We could have used DiffEqOperators.jl like for the Swift-Hohenberg tutorial.

We shall now compute the equilibria and their stability.

n = 500

# parameters of the Brusselator model and guess for the stationary solution
par_bru = (Ξ± = 2., Ξ² = 5.45, D1 = 0.008, D2 = 0.004, l = 0.3)
sol0 = vcat(par_bru.Ξ± * ones(n), par_bru.Ξ² / par_bru.Ξ± * ones(n))

# bifurcation problem
probBif = BK.BifurcationProblem(Fbru, sol0, par_bru, (@lens _.l);
J = Jbru_sp,
plot_solution = (x, p; kwargs...) -> (plotsol(x; label="", kwargs... )),
record_from_solution = (x, p) -> x[div(n,2)])

For the eigensolver, we use a Shift-Invert algorithm (see Eigen solvers (Eig))

eigls = EigArpack(1.1, :LM)

We continue the trivial equilibrium to find the Hopf points

opt_newton = NewtonPar(eigsolver = eigls, tol = 1e-9)
opts_br_eq = ContinuationPar(dsmin = 0.001, dsmax = 0.1,
p_max = 1.9, nev = 21,
newton_options = opt_newton, max_steps = 1000,
# specific options for precise localization of Hopf points
n_inversion = 6)

br = continuation(probBif, PALC(),opts_br_eq, normC = norminf)
 ββ Curve type: EquilibriumCont
ββ Number of points: 18
ββ Type of vectors: Vector{Float64}
ββ Parameter l starts at 0.3, ends at 1.9
ββ Algo: PALC
ββ Special points:

- #  1,     hopf at l β +0.51206298 β (+0.51199393, +0.51206298), |Ξ΄p|=7e-05, [converged], Ξ΄ = ( 2,  2), step =   6
- #  2,     hopf at l β +1.02402486 β (+1.02388676, +1.02402486), |Ξ΄p|=1e-04, [converged], Ξ΄ = ( 2,  2), step =  10
- #  3,     hopf at l β +1.53598675 β (+1.53584864, +1.53598675), |Ξ΄p|=1e-04, [converged], Ξ΄ = ( 2,  2), step =  14
- #  4, endpoint at l β +1.90000000,                                                                     step =  17


We obtain the following bifurcation diagram with 3 Hopf bifurcation points

scene = plot(br)

## Normal form computation

We can compute the normal form of the Hopf points as follows

hopfpt = get_normal_form(br, 1)
SuperCritical - Hopf bifurcation point at l β 0.512062980959364.
Frequency Ο β 2.139470722402076
Period of the periodic orbit β 2.9367942460671688
Normal form zβ(iΟ + aβΞ΄p + bβ|z|Β²):
ββ a = 0.8785637258030531 + 0.5680545896552575im
ββ b = -0.0009377741298574548 + 0.0009393095359773213im


## Continuation of Hopf points

We use the bifurcation points guesses located in br.specialpoint to turn them into precise bifurcation points. For the second one, we have

# index of the Hopf point in br.specialpoint
ind_hopf = 2

# newton iterations to compute the Hopf point
hopfpoint = newton(br, ind_hopf; normN = norminf)
BK.converged(hopfpoint) && printstyled(color=:red, "--> We found a Hopf Point at l = ", hopfpoint.u.p[1], ", Ο = ", hopfpoint.u.p[2], ", from l = ", br.specialpoint[ind_hopf].param, "\n")
--> We found a Hopf Point at l = 1.0239851623501666, Ο = -2.1395092637754476, from l = 1.0240248633536086

We now perform a Hopf continuation with respect to the parameters l, Ξ²

Tip

You don't need to call newton first in order to use continuation.

optcdim2 = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, p_max = 6.5, p_min = 0.0, newton_options = opt_newton, detect_bifurcation = 0)
br_hopf = continuation(br, ind_hopf, (@lens _.Ξ²), optcdim2, normC = norminf, jacobian_ma = :minaug)
scene = plot(br_hopf)

## Continuation of periodic orbits (Finite differences)

Here, we perform continuation of periodic orbits branching from the Hopf bifurcation points.We need an educated guess for the periodic orbit which is given by guess_from_hopf:

# number of time slices
M = 51

l_hopf, Th, orbitguess2, hopfpt, vec_hopf = BK.guess_from_hopf(br, ind_hopf,
opts_br_eq.newton_options.eigsolver,
M, 2.7; phase = 0.25)

We wish to make two remarks at this point. The first is that an initial guess is composed of a space time solution and of the guess for the period Th of the solution. Note that the argument 2.7 is a guess for the amplitude of the orbit.

# orbit initial guess from guess_from_hopf, is not a vector, so we reshape it
orbitguess_f2 = reduce(vcat, orbitguess2)
orbitguess_f = vcat(vec(orbitguess_f2), Th) |> vec

The second remark concerns the phase 0.25 written above. To account for the additional unknown (i.e. the period), periodic orbit localisation using Finite Differences requires an additional constraint (see Periodic orbits based on Trapezoidal rule for more details). In the present case, this constraint is

$$$< u(0) - u_{hopf}, \phi> = 0$$$

where u_{hopf} is the equilibrium at the Hopf bifurcation and $\phi$ is real.(vec_hopf) where vec_hopf is the eigenvector. This is akin to a PoincarΓ© section. We do not put any constraint on $u(t)$ albeit this is possible (see Periodic orbits based on Trapezoidal rule.

The phase of the periodic orbit is set so that the above constraint is satisfied. We shall now use Newton iterations to find a periodic orbit.

Given our initial guess, we create a (family of) problem which encodes the functional associated to finding Periodic orbits based on Trapezoidal rule (see Periodic orbits based on Trapezoidal rule for more information):

poTrap = PeriodicOrbitTrapProblem(
probBif,    		   # pass the bifurcation problem
real.(vec_hopf),	 # used to set Ο, see the phase constraint
hopfpt.u,          # used to set uhopf, see the phase constraint
M, 2n;             # number of time slices
jacobian = :FullSparseInplace) # jacobian of PO functional

nothing #hide

To evaluate the functional at x, you call it like a function: poTrap(x, par) for the parameter par.

Using the functional for deflation, Fold of limit cycles...

The functional poTrap gives you access to the underlying methods to call a regular newton. For example the functional is x -> poTrap(x, par) at parameter par. The (sparse) Jacobian at (x,p) is computed like this poTrap(Val(:JacFullSparse), x, p) while the Matrix Free version is dx -> poTrap((x, p, dx). This also allows you to call the newton deflated method (see Deflated problems) or to locate Fold point of limit cycles see PeriodicOrbitTrapProblem. You can also use preconditioners. In the case of more computationally intense problems (like the 2d Brusselator), this might be mandatory as using LU decomposition for the linear solve will use too much memory. See also the example 2d Ginzburg-Landau equation (finite differences, codim 2, Hopf aBS)

For convenience, we provide a simplified newton / continuation methods for periodic orbits. One has just to pass a PeriodicOrbitTrapProblem.

# we use the linear solver LSFromBLS to speed up the computations
opt_po = NewtonPar(tol = 1e-10, verbose = true, max_iterations = 14, linsolver = BK.LSFromBLS())

# we set the parameter values
poTrap = @set poTrap.prob_vf.params = (@set par_bru.l = l_hopf + 0.01)

outpo_f = @time newton(poTrap, orbitguess_f, opt_po, normN = norminf)
BK.converged(outpo_f) && printstyled(color=:red, "--> T = ", outpo_f.u[end], ", amplitude = ", BK.amplitude(outpo_f.u, n, M; ratio = 2),"\n")
# plot of the periodic orbit
BK.plot_periodic_potrap(outpo_f.u, n, M; ratio = 2)

Finally, we can perform continuation of this periodic orbit using the specialized call continuationPOTrap

opt_po = @set opt_po.eigsolver = EigArpack(; tol = 1e-5, v0 = rand(2n))
opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.03, ds= 0.01,
p_max = 3.0, max_steps = 20,
newton_options = opt_po, nev = 5, tol_stability = 1e-8, detect_bifurcation = 0)

br_po = continuation(poTrap,
outpo_f.u, PALC(),
opts_po_cont;
verbosity = 2,	plot = true,
plot_solution = (x, p;kwargs...) -> heatmap!(reshape(x[1:end-1], 2*n, M)'; ylabel="time", color=:viridis, kwargs...),
normC = norminf)

Scene = title!("")

## Deflation for periodic orbit problems

Looking for periodic orbits branching of bifurcation points, it is very useful to use newton algorithm with deflation. We thus define a deflation operator (see previous example)

deflationOp = DeflationOperator(2, (x,y) -> dot(x[1:end-1], y[1:end-1]), 1.0, [zero(orbitguess_f)])

which allows to find periodic orbits different from orbitguess_f. Note that the dot product removes the last component, i.e. the period of the cycle is not considered during this particular deflation. We can now use

outpo_f = @time newton(poTrap, orbitguess_f, deflationOp, opt_po; normN = norminf)

## Floquet coefficients

A basic method for computing Floquet coefficients based on the eigenvalues of the monodromy operator is available (see FloquetQaD). It is precise enough to locate bifurcations. Their computation is triggered like in the case of a regular call to continuation:

opt_po = @set opt_po.eigsolver = DefaultEig()
opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.04, ds= -0.01, p_max = 3.0, max_steps = 200, newton_options = opt_po, nev = 5, tol_stability = 1e-6)
br_po = @time continuation(poTrap, outpo_f.u, PALC(),
opts_po_cont; verbosity = 3, plot = true,
plot_solution = (x, p;kwargs...) -> heatmap!(reshape(x[1:end-1], 2*n, M)'; ylabel="time", color=:viridis, kwargs...), normC = norminf)

A more complete diagram can be obtained combining the methods (essentially deflation and Floquet) described above. It shows the period of the periodic orbits as function of l. See example/brusselator.jl for more information.

Floquet multipliers computation

The computation of Floquet multipliers is necessary for the detection of bifurcations of periodic orbits (which is done by analyzing the Floquet exponents obtained from the Floquet multipliers). Hence, the eigensolver needs to compute the eigenvalues with largest modulus (and not with largest real part which is their default behavior). This can be done by changing the option which = :LM of the eigensolver. Nevertheless, note that for most implemented eigensolvers in the current Package, the proper option is set when the computation of Floquet multipliers is requested.

Performances

This example is clearly not optimized because we wanted to keep it simple. We can use a Matrix-Free version of the functional and preconditioners to speed this up. Floquet multipliers could also be computed in a Matrix-Free manner. See examples/brusselator.jl for more efficient methods. See also 2d Ginzburg-Landau equation (finite differences, codim 2, Hopf aBS) for a more advanced example where we introduce those methods.

## Continuation of periodic orbits (Standard Shooting)

Note that what follows is not really optimized on the DifferentialEquations.jl side. Indeed, we do not use automatic differentiation, we do not pass the sparsity pattern,...

We now turn to a different method based on the flow of the Brusselator. To compute this flow (time stepper), we need to be able to solve the differential equation (actually a PDE) associated to the vector field Fbru. We will show how to do this with an implicit method Rodas4P from DifferentialEquations.jl. Note that the user can pass its own time stepper but for convenience, we use the ones in DifferentialEquations.jl. More information regarding the shooting method is contained in Periodic orbits based on the shooting method.

n = 100

# different parameters to define the Brusselator model and guess for the stationary solution
par_bru = (Ξ± = 2., Ξ² = 5.45, D1 = 0.008, D2 = 0.004, l = 0.3)
sol0 = vcat(par_bru.Ξ± * ones(n), par_bru.Ξ²/par_bru.Ξ± * ones(n))
probBif = re_make(probBif, u0 = sol0)

eigls = EigArpack(1.1, :LM)
opts_br_eq = ContinuationPar(dsmin = 0.001, dsmax = 0.00615, ds = 0.0061, p_max = 1.9,
detect_bifurcation = 3, nev = 21, plot_every_step = 50,
newton_options = NewtonPar(eigsolver = eigls, tol = 1e-9), max_steps = 1060)

br = @time continuation(probBif, PALC(), opts_br_eq, normC = norminf)

We need to create a guess for the periodic orbit. We proceed as previously:

# number of time slices
M = 5

# index of the Hopf point in the branch br
ind_hopf = 1

l_hopf, Th, orbitguess2, hopfpt, vec_hopf = BK.guess_from_hopf(br, ind_hopf,
opts_br_eq.newton_options.eigsolver, M, 22*0.075)
#
orbitguess_f2 = reduce(hcat, orbitguess2)
orbitguess_f = vcat(vec(orbitguess_f2), Th) |> vec

Let us now initiate the Standard Shooting method. To this aim, we need to provide a guess of the periodic orbit at times $T/M_{sh}$ where $T$ is the period of the cycle and $M_{sh}$ is the number of slices along the periodic orbits. If $M_{sh} = 1$, this the Standard Simple Shooting and the Standard Multiple one otherwise. See ShootingProblem for more information.

dM = 2
orbitsection = Array(orbitguess_f2[:, 1:dM:M])

# the last component is an estimate of the period of the cycle.
initpo = vcat(vec(orbitsection), 3.0)

Finally, we need to build a problem which encodes the Shooting functional. This done as follows where we first create the time stepper. For performance reasons, we rely on SparseDiffTools

using DifferentialEquations, DiffEqOperators, SparseDiffTools, SparseArrays

FOde(f, x, p, t) = Fbru!(f, x, p)

u0 = sol0 .+ 0.01 .* rand(2n)

# parameter close to the Hopf bifurcation point
par_hopf = (@set par_bru.l = l_hopf + 0.01)

jac_prototype = Jbru_sp(ones(2n), @set par_bru.Ξ² = 0)
jac_prototype.nzval .= ones(length(jac_prototype.nzval))
_colors = matrix_colors(jac_prototype)
vf = ODEFunction(FOde; jac_prototype = jac_prototype, colorvec = _colors)
prob = ODEProblem(vf,  sol0, (0.0, 520.), par_bru)

We create the parallel standard shooting problem:

# this encodes the functional for the Shooting problem
probSh = ShootingProblem(
# we pass the ODEProblem encoding the flow and the time stepper
prob, Rodas4P(),

# this is for the phase condition, you can pass your own section as well
[orbitguess_f2[:,ii] for ii=1:dM:M];

parallel = true,

# these are options passed to the ODE time stepper
abstol = 1e-10, reltol = 1e-8,

# parameter axis
lens = (@lens _.l),

# parameters
par = par_hopf,

# jacobian of the periodic orbit functional
jacobian = BK.FiniteDifferencesMF())

We are now ready to call newton

ls = GMRESIterativeSolvers(reltol = 1e-7, N = length(initpo), maxiter = 100)
optn_po = NewtonPar(verbose = true, tol = 1e-9,  max_iterations = 20, linsolver = ls)
outpo = @time newton(probSh, initpo, optn_po; normN = norminf)
plot(initpo[1:end-1], label = "Init guess")
plot!(outpo.u[1:end-1], label = "sol")

which gives (note that we did not have a really nice guess...)

βββββββββββββββββββββββββββββββββββββββββββββββββββββββ
β Newton step         residual     linear iterations  β
βββββββββββββββ¬βββββββββββββββββββββββ¬βββββββββββββββββ€
β       0     β       1.9613e-01     β        0       β
β       1     β       5.6101e-02     β       45       β
β       2     β       1.0307e-01     β       49       β
β       3     β       4.1119e-03     β       48       β
β       4     β       8.0511e-03     β       49       β
β       5     β       3.8250e-02     β       48       β
β       6     β       9.8080e-03     β       49       β
β       7     β       2.1179e+01     β       53       β
β       8     β       2.0105e+00     β       36       β
β       9     β       2.0545e+00     β       49       β
β      10     β       4.8793e-01     β       49       β
β      11     β       4.8457e-02     β       46       β
β      12     β       2.3299e-02     β       49       β
β      13     β       1.6365e-02     β       48       β
β      14     β       1.3534e-04     β       49       β
β      15     β       1.4582e-05     β       48       β
β      16     β       1.5886e-08     β       49       β
β      17     β       1.7228e-11     β       49       β
βββββββββββββββ΄βββββββββββββββββββββββ΄βββββββββββββββββ
9.706977 seconds (7.61 M allocations: 13.964 GiB, 3.62% gc time)

and

Note that using Simple Shooting, the convergence is much faster. Indeed, running the code above with dM = 10 gives:

βββββββββββββββββββββββββββββββββββββββββββββββββββββββ
β Newton step         residual     linear iterations  β
βββββββββββββββ¬βββββββββββββββββββββββ¬βββββββββββββββββ€
β       0     β       6.1712e-03     β        0       β
β       1     β       3.4465e-03     β        6       β
β       2     β       1.0516e-01     β        8       β
β       3     β       7.4614e-03     β        6       β
β       4     β       1.6620e-03     β        7       β
β       5     β       3.9589e-04     β        7       β
β       6     β       4.3043e-05     β        8       β
β       7     β       1.7232e-06     β        8       β
β       8     β       8.0455e-09     β        8       β
β       9     β       3.9453e-11     β        8       β
βββββββββββββββ΄βββββββββββββββββββββββ΄βββββββββββββββββ
0.612070 seconds (217.21 k allocations: 523.069 MiB, 4.83% gc time)
Convergence and speedup

The convergence is much worse for the multiple shooting than for the simple one. This is reflected above in the number of linear iterations made during the newton solve. The reason for this is because of the cyclic structure of the jacobian which impedes GMRES from converging fast. This can only be resolved with an improved GMRES which we'll provide in the future.

Finally, we can perform continuation of this periodic orbit using a specialized version of continuation:

# note the eigensolver computes the eigenvalues of the monodromy matrix. Hence
# the dimension of the state space for the eigensolver is 2n
opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, p_max = 1.5,
max_steps = 500, newton_options = (@set optn_po.tol = 1e-7), nev = 25,
tol_stability = 1e-8, detect_bifurcation = 0)

br_po = @time continuation(probSh, outpo.u, PALC(),
opts_po_cont; verbosity = 2,
# specific bordered linear solver
linear_algo = MatrixFreeBLS(@set ls.N = ls.N+1),
plot = true,
plot_solution = (x, p; kwargs...) -> BK.plot_periodic_shooting!(x[1:end-1], length(1:dM:M); kwargs...),
record_from_solution = (u, p) -> u[end], normC = norminf)

We can observe that simple shooting is faster but the Floquet multipliers are less accurate than for multiple shooting. Also, when the solution is very unstable, simple shooting can have spurious branch switching. Finally, note the $0=\log 1$ eigenvalue of the monodromy matrix in the graph below.

## Continuation of periodic orbits (PoincarΓ© Shooting)

We now turn to another Shooting method, namely the PoincarΓ© one. We can provide this method thanks to the unique functionalities of DifferentialEquations.jl. More information is provided at PoincareShootingProblem and Periodic orbits based on the shooting method but basically, it is a shooting method between PoincarΓ© sections $\Sigma_i$ (along the orbit) defined by hyperplanes. As a consequence, the dimension of the unknowns is $M_{sh}\cdot(N-1)$ where $N$ is the dimension of the phase space. Indeed, each time slice lives in an hyperplane $\Sigma_i$. Additionally, the period $T$ is not an unknown of the method but rather a by-product. However, the method requires the time stepper to find when the flow hits an hyperplane $\Sigma_i$, something called event detection.

We show how to use this method, the code is very similar to the case of the Standard Shooting. We first define the functional for PoincarΓ© Shooting Problem

# sub-sampling factor of a initial guess for the periodic orbit
dM = 5

# vectors to define the hyperplanes Sigma_i
normals = [Fbru(orbitguess_f2[:,ii], par_hopf)/(norm(Fbru(orbitguess_f2[:,ii], par_hopf))) for ii = 1:dM:M]
centers = [orbitguess_f2[:,ii] for ii = 1:dM:M]

# functional to hold the Poincare Shooting Problem
probHPsh = PoincareShootingProblem(
# ODEProblem, ODE solver used to compute the flow
prob, Rodas4P(),

# parameters for the PoincarΓ© sections
normals, centers;

parallel = true,

# Parameters passed to the ODE solver
abstol = 1e-10, reltol = 1e-8,

# parameter axis
lens = (@lens _.l),

# parameters
par = par_hopf,

# jacobian of the periodic orbit functional
jacobian = BK.FiniteDifferencesMF())

Let us now compute an initial guess for the periodic orbit, it must live in the hyperplanes $\Sigma_i$. Fortunately, we provide projections on these hyperplanes.

# projection of the initial guess on the hyperplanes. We assume that the centers[ii]
# form the periodic orbit initial guess.
initpo_bar = reduce(vcat, BK.projection(probHPsh, centers))

We can now call continuation to get the first branch.

# eigen / linear solver
eig = EigKrylovKit(tol= 1e-12, xβ = rand(2n-1), dim = 40)
ls = GMRESIterativeSolvers(reltol = 1e-11, N = length(vec(initpo_bar)), maxiter = 500)

# newton options
optn = NewtonPar(verbose = true, tol = 1e-9,  max_iterations = 140, linsolver = ls)

# continuation options
opts_po_cont_floquet = ContinuationPar(dsmin = 0.0001, dsmax = 0.05, ds= 0.001,
p_max = 2.5, max_steps = 500, nev = 10,
tol_stability = 1e-5, detect_bifurcation = 3, plot_every_step = 3)
opts_po_cont_floquet = @set opts_po_cont_floquet.newton_options =
NewtonPar(linsolver = ls, eigsolver = eig, tol = 1e-9, verbose = true)

# continuation run
br_po = @time continuation(probHPsh, vec(initpo_bar), PALC(),
opts_po_cont_floquet; verbosity = 3,
linear_algo = MatrixFreeBLS(@set ls.N = ls.N+1),
plot = true,
plot_solution = (x, p; kwargs...) -> BK.plot!(x; label="", kwargs...),
normC = norminf)		

We also obtain the following information:

julia> br_po
Branch number of points: 41
Bifurcation points:
- #  1,    bp at p β 1.20987963 β (1.20128196, 1.20987963), |Ξ΄p|=9e-03, [converged], Ξ΄ = ( 1,  0), step =  21, eigenelements in eig[ 22], ind_ev =   1
- #  2,    ns at p β 1.78687615 β (1.77831727, 1.78687615), |Ξ΄p|=9e-03, [converged], Ξ΄ = ( 2,  2), step =  30, eigenelements in eig[ 31], ind_ev =   3
- #  3,    pd at p β 1.85103701 β (1.84676466, 1.85103701), |Ξ΄p|=4e-03, [converged], Ξ΄ = ( 1,  1), step =  31, eigenelements in eig[ 32], ind_ev =   4
- #  4,    ns at p β 1.87667870 β (1.86813520, 1.87667870), |Ξ΄p|=9e-03, [converged], Ξ΄ = ( 2,  2), step =  32, eigenelements in eig[ 33], ind_ev =   6

## References

• Lust

Numerical Bifurcation Analysis of Periodic Solutions of Partial Differential Equations, Lust, 1997.