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
— Methodnewton(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
— Methodcontinuation(
prob,
orbitguess,
alg,
_contParams;
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 Gorbitguess
a guess for the periodic orbit whereorbitguess[end]
is an estimate of the period of the orbit. It could be a vector of sizeN * M + 1
whereM
is the number of time slices,N
is the dimension of the phase space. This must be compatible with the numbersN, M
inprob
.alg
continuation algorithmcontParams
same as for the regularcontinuation
method
Keyword arguments
linear_algo
same as incontinuation
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 ofdG
. 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 matrixdG
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 matrixdG
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 invertdG
using a bordered linear solver. - For
jacobian = :BorderedSparseInplace
, this is the same as for:BorderedLU
but the cyclic matrixdG
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 fordG
: note that a preconditioner is very likely required here because of the cyclic shape ofdG
which affects negatively the convergence properties of GMRES. - For
jacobian = :BorderedMatrixFree
, a matrix free linear solver is used but forJc
only (see docs): it means thatoptions.linsolver
is used to invertJc
. These two Matrix-Free options thus expose different part of the jacobiandG
in order to use specific preconditioners. For example, an ILU preconditioner onJc
could remove the constraints indG
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.