# π 1d Brusselator

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
```

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, Ξ²`

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`

.

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.

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.

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];
# enable threading
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)
```

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;
# enable threading
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.