# Periodic orbits based on Trapezoidal rule

- Periodic orbits based on Trapezoidal rule

The Trapezoid method allows to compute periodic orbits by discretizing time using Finite Differences based on a trapezoidal rule . The method is implemented in the structure `PeriodicOrbitTrapProblem`

. The general method is very well exposed in ^{[Uecker]},^{[Lust]} and we adopt the notations of the first reference.

We look for periodic orbits as solutions $(x(0),T)$ of

\[M_a\dot x = T\cdot F(x),\ x(0)=x(1)\in\mathbb R^n\tag{1}\]

where $M_a$ is a mass matrix (default is the identity one).

In order to have a unique solution, we need to remove the phase freedom. This is done by imposing a *phase* condition

\[\frac1T\int_0^T\langle x(s)-x_\pi(s), \phi(s)\rangle ds\approx \frac{1}{m}\sum\limits_{i=1}^m\langle x_{i} - x_{\pi,i}, \phi_{i}\rangle = 0\]

for some $x_\pi,\phi$ which are chosen (wisely).

We note `m`

the number of time slices of the periodic orbit. By discretizing the above problem, we obtain

\[\begin{array}{l} 0= M_a\left(x_{j}-x_{j-1}\right)-\frac{h}{2} \left(F\left(x_{j}\right)+F\left(x_{j-1}\right)\right)\equiv G_j(x),\quad j=1,\cdots,m-1 \\ 0= x_m-x_1 \equiv G_m(x) \\ 0=\sum\limits_{i=1}^m\langle x_{i} - x_{\pi,i}, \phi_{i}\rangle=0 \end{array}\]

where $x_0=x_m$ and $h=T/m$. In view of the Newton method, we study the jacobian of the above system. The Jacobian *w.r.t.* $(x_0,T)$ is given by

\[\mathcal{J}=\left(\begin{array}{cc}{A_1} & {\partial_TG} \\ {\star} & {d}\end{array}\right)\tag{2}\]

where

\[A_{\gamma}:=\left(\begin{array}{ccccccc} {M_{1}} & {0} & {0} & {0} & {\cdots} & {-H_{1}} & {0} \\ {-H_{2}} & {M_{2}} & {0} & {0} & {\cdots} & {0} & {0} \\ {0} & {-H_{3}} & {M_{3}} & {0} & {\cdots} & {0} & {0} \\ {\vdots} & {\cdots} & {\ddots} & {\ddots} & {\ddots} & {\vdots} & {\vdots} \\ {0} & {\cdots} & {\cdots} & {\ddots} & {\ddots} & {0} & {0} \\ {0} & {\cdots} & {\cdots} & {0} & {-H_{m-1}} & {M_{m-1}} & {0} \\ {-\gamma I} & {0} & {\cdots} & {\cdots} & {\cdots} & {0} & {I} \end{array}\right)\]

with $M_i := M_a- \frac h2dF(x_i)$ and $H_i := M_a+\frac h2dF(x_{i-1})$.

We solve the linear equation $\mathcal J\cdot sol = rhs$ with a bordering strategy (*i.e.* the linear solver is a subtype of `<: AbstractBorderedLinearSolver`

) which in turn requires to solve $A_\gamma z=b$ where $z=(x,x_m)$. We also solve this equation with a bordering strategy but this time, it can be simplified as follows. If we write $b=(f,g)$, one gets $J_c x=f$ and $x_m=g+\gamma x_1$ where $x_1$ is the first time slice of $x$ and $J_c$ is the following **cyclic matrix**:

\[J_c:=\left(\begin{array}{ccccccc} {M_{1}} & {0} & {0} & {0} & {\cdots} & {-H_{1}} \\ {-H_{2}} & {M_{2}} & {0} & {0} & {\cdots} & {0} \\ {0} & {-H_{3}} & {M_{3}} & {0} & {\cdots} & {0} \\ {\vdots} & {\cdots} & {\ddots} & {\ddots} & {\ddots} & {\vdots} \\ {0} & {\cdots} & {\cdots} & {\ddots} & {\ddots} & {0} \\ {0} & {\cdots} & {\cdots} & {0} & {-H_{m-1}} & {M_{m-1}} \\ \end{array}\right)\]

Our code thus provides methods to invert $J_c$ and $A_\gamma$ using a sparse solver or a Matrix-Free one. A preconditioner can be used.

## Encoding of the functional

The functional is encoded in the composite type `PeriodicOrbitTrapProblem`

. See the link for more information, in particular on how to access the underlying functional, its jacobian and other matrices related to it like $A_\gamma, J_c$...

## Preconditioning

We strongly advise you to use a preconditioner to deal with the above linear problem. See 2d Ginzburg-Landau equation (finite differences, codim 2, Hopf aBS) for an example.

## Linear solvers

We provide many different linear solvers to take advantage of the formulations. These solvers are available through the argument `jacobian`

in the constructor of `PeriodicOrbitTrapProblem`

. For example, you can pass `jacobian = :FullLU`

. Note that all the internal solvers and Jacobians are set up automatically, you don't need to do anything. However, for the sake of explanation, we detail how this works.

### 1. FullLU

When using `jacobianPO = :FullLU`

, this triggers the computation of $\mathcal J$ as in (2) at each step of newton/continuation. The Jacobian matrix $\mathcal J$ is stored a SparseArray. This can be quite costly flow large $n$ (see (1)). This Jacobian is often used with the the linear solver `DefaultLS()`

.

### 2. FullSparseInplace

Same as `:FullLU`

but the Jacobian is allocated only once and updated inplace. This is much faster than `:FullLU`

but the sparsity pattern of `dF`

must be constant.

### 3. Dense

Same as `: FullSparseInplace`

above but the matrix `dG`

is dense. It is also updated inplace. This is useful to study ODE of small dimension.

### 4. FullMatrixFree

A matrix free linear solver is used for $\mathcal J$: note that a preconditioner is very likely required here because of the cyclic shape of $\mathcal J$ which affects negatively the convergence properties of iterative solvers. Note that $\mathcal J$ is never formed in this case.

### 5. BorderedLU

For `:BorderedLU`

, we take advantage of the bordered shape of the linear solver and use a LU decomposition to invert `dG`

using a bordered linear solver. More precisely, the bordered structure of $\mathcal J$ is stored using the internal structure `POTrapJacobianBordered`

. Then, $\mathcal J$ is inverted using the custom bordered linear solver `PeriodicOrbitTrapBLS`

which is based on the bordering strategy (see Bordered linear solvers (BLS)). This particular solver is based on an explicit formula which only requires to invert $A_\gamma$: this is done by the linear solver `AγLinearSolver`

. In a nutshell, we have:

`PeriodicOrbitTrapBLS = BorderingBLS(solver = AγLinearSolver(), check_precision = false)`

### 6. BorderedSparseInplace

Same as `:BorderedLU`

but the Jacobian is allocated only once and updated inplace. This is much faster than `:BorderedLU`

but the sparsity pattern of `dF`

must be constant.

### 7. BorderedMatrixFree

A matrix free linear solver is used but for $\mathcal J_c$ only: it means that `options.linsolver`

is used to invert $\mathcal J_c$.

These two Matrix-Free options, `:FullMatrixFree`

and `:BorderedMatrixFree`

, thus expose different part of the Jacobian $\mathcal J$ in order to use specific preconditioners. For example, an ILU preconditioner on $\mathcal J_c$ could remove the constraints in $\mathcal J$ and lead to poor convergence. Of course, for these last two methods, a preconditioner is likely be required.

## Floquet multipliers computation

### Default method

A **not very precise** algorithm for computing the Floquet multipliers is provided in the package. The method, dubbed Quick and Dirty (QaD), is not numerically very precise for large / small Floquet exponents because it relies on constructing the monodromy matrix.

Note that the computation of the eigenvalues can be iterative or direct based on the eigensolver passed in arguments.

It amounts to computing the eigenvalues of the monodromy matrix

\[\mathcal{M}=M_{1}^{-1} H_{1} M_{m-1}^{-1} H_{m-1} \cdots M_{2}^{-1} H_{2}.\]

The method allows, nevertheless, to detect bifurcations of periodic orbits. It seems to work reasonably well for the tutorials considered here. For more information, have a look at `FloquetQaD`

.

### Most precise method

The state of the art method is based on a Periodic Schur decomposition. It is available through the package PeriodicSchurBifurcationKit.jl. For more information, have a look at `FloquetPQZ`

.

## Computation with `newton`

We provide a simplified call to `newton`

to locate the periodic orbits. Compared to the regular `newton`

function, there is an additional option `linearalgo`

to select one of the many ways to deal with the above linear problem. The default solver `linearalgo`

is `:BorderedLU`

.

Have a look at the Periodic orbits based on Trapezoidal rule example for the Brusselator for a basic example and at 2d Ginzburg-Landau equation for a more advanced one.

The docs for this specific `newton`

are located at `newton`

.

## Computation with `newton`

and deflation

We also provide a simplified call to `newton`

to locate the periodic orbit with a deflation operator.

`BifurcationKit.newton`

— Method```
newton(probPO, orbitguess, defOp, options; kwargs...)
```

This function is similar to `newton(probPO, orbitguess, options, jacobianPO; kwargs...)`

except that it uses deflation in order to find periodic orbits different from the ones stored in `defOp`

. We refer to the mentioned method for a full description of the arguments. The current method can be used in the vicinity of a Hopf bifurcation to prevent the Newton-Krylov algorithm from converging to the equilibrium point.

## Continuation

Have a look at the Periodic orbits based on Trapezoidal rule example for the Brusselator. We refer to `continuation`

for more information regarding the arguments.

`BifurcationKit.continuation`

— Method```
continuation(
prob,
orbitguess,
alg,
_contParams;
894,
record_from_solution,
linear_algo,
kwargs...
)
```

This is the continuation routine for computing a periodic orbit using a functional G based on Finite Differences and a Trapezoidal rule.

**Arguments**

`prob::PeriodicOrbitTrapProblem`

encodes the functional G`orbitguess`

a guess for the periodic orbit where`orbitguess[end]`

is an estimate of the period of the orbit. It could be a vector of size`N * M + 1`

where`M`

is the number of time slices,`N`

is the dimension of the phase space. This must be compatible with the numbers`N, M`

in`prob`

.`alg`

continuation algorithm`contParams`

same as for the regular`continuation`

method

**Keyword arguments**

`linear_algo`

same as in`continuation`

Specify the choice of the jacobian (and linear algorithm), `jacobian`

must belong to `[:FullLU, :FullSparseInplace, :Dense, :DenseAD, :BorderedLU, :BorderedSparseInplace, :FullMatrixFree, :BorderedMatrixFree, :FullMatrixFreeAD]`

. This is used to select a way of inverting the jacobian `dG`

of the functional G.

- For
`jacobian = :FullLU`

, we use the default linear solver based on a sparse matrix representation of`dG`

. This matrix is assembled at each newton iteration. This is the default algorithm. - For
`jacobian = :FullSparseInplace`

, this is the same as for`:FullLU`

but the sparse matrix`dG`

is updated inplace. This method allocates much less. In some cases, this is significantly faster than using`:FullLU`

. Note that this method can only be used if the sparsity pattern of the jacobian is always the same. - For
`jacobian = :Dense`

, same as above but the matrix`dG`

is dense. It is also updated inplace. This option is useful to study ODE of small dimension. - For
`jacobian = :DenseAD`

, evaluate the jacobian using ForwardDiff - For
`jacobian = :BorderedLU`

, we take advantage of the bordered shape of the linear solver and use a LU decomposition to invert`dG`

using a bordered linear solver. - For
`jacobian = :BorderedSparseInplace`

, this is the same as for`:BorderedLU`

but the cyclic matrix`dG`

is updated inplace. This method allocates much less. In some cases, this is significantly faster than using`:BorderedLU`

. Note that this method can only be used if the sparsity pattern of the jacobian is always the same. - For
`jacobian = :FullMatrixFree`

, a matrix free linear solver is used for`dG`

: note that a preconditioner is very likely required here because of the cyclic shape of`dG`

which affects negatively the convergence properties of GMRES. - For
`jacobian = :BorderedMatrixFree`

, a matrix free linear solver is used but for`Jc`

only (see docs): it means that`options.linsolver`

is used to invert`Jc`

. These two Matrix-Free options thus expose different part of the jacobian`dG`

in order to use specific preconditioners. For example, an ILU preconditioner on`Jc`

could remove the constraints in`dG`

and lead to poor convergence. Of course, for these last two methods, a preconditioner is likely to be required. - For
`jacobian = :FullMatrixFreeAD`

, the evaluation map of the differential is derived using automatic differentiation. Thus, unlike the previous two cases, the user does not need to pass a Matrix-Free differential.

Note that by default, the method prints the period of the periodic orbit as function of the parameter. This can be changed by providing your `record_from_solution`

argument.