Periodic orbits based on orthogonal collocation

• Periodic orbits based on orthogonal collocation
• • • Number of unknowns
    • Phase condition
  • Precision of the scheme
  • Discretization of the BVP and jacobian
  • Interpolation
  • Mesh adaptation
  • Encoding of the functional
  • Jacobians
  • Linear solvers
  • Floquet multipliers computation
  • Computation with newton
  • Continuation
  • Examples of optimized use
  • • Inplace jacobian
    • Optimized linear solver
    • Optimized Floquet computation
  • References

We compute Ntst time slices of a periodic orbit using orthogonal collocation. This is implemented in the structure Collocation.

The general method is very well exposed in ^[Dankowicz],^[Doedel] and we adopt the notations of ^[Dankowicz]. However our implementation is based on ^[Doedel] because it is more economical (less equations) when it enforces the continuity of the solution.

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

\[\dot x = T\cdot F(x),\ x(0)=x(1)\in\mathbb R^n.\]

We focus on the differential equality and consider a partition of the time domain

\[0=\tau_{1}<\cdots<\tau_{j}<\cdots<\tau_{N_{tst}+1}=1\]

where the points are referred to as mesh points. On each mesh interval $[\tau_j,\tau_{j+1}]$ for $j=1,\cdots,N_{tst}$, we define the affine transformation

\[\tau=\tau^{(j)}(\sigma):=\tau_{j}+\frac{(1+\sigma)}{2}\left(\tau_{j+1}-\tau_{j}\right), \sigma \in[-1,1].\]

The functions $x^{(j)}$ defined on $[-1,1]$ by $x^{(j)}(\sigma) \equiv x(\tau_j(\sigma))$ satisfies the following equation on $[-1,1]$:

\[\dot x^{(j)} = T\frac{\tau_{j+1}-\tau_j}{2}\cdot F(x^{(j)})\tag{$E_j$}\]

with the continuity equation $x^{(j+1)}(-1) = x^{(j)}(1)$.

We now aim at solving $(E_j)$ by using an approximation with a polynomial of degree $m$. Following ^[Dankowicz], we define a (uniform) partition:

\[-1=\sigma_{1}<\cdots<\sigma_{i}<\cdots<\sigma_{m+1}=1.\]

The points $\tau_{i,j} = \tau^{(i)}(\sigma_j)$ are called the base points: they serve as collocation points.

The associated $m+1$ Lagrange polynomials of degree $m$ are:

\[\mathcal{L}_{i}(\sigma):=\prod_{k=1, k \neq i}^{m+1} \frac{\sigma-\sigma_{k}}{\sigma_{i}-\sigma_{k}}, i=1, \ldots, m+1.\]

We then introduce the approximation $p_j$ of $x^{(j)}$:

\[\mathcal p_j(\sigma)\equiv \sum\limits_{k=1}^{m+1}\mathcal L_k(\sigma)x_{j,k}\]

and the problem to be solved at the collocation nodes $z_l$, $l=1,\cdots,m$:

\[\forall 1\leq l\leq m,\quad 1\leq j\leq N_{tst},\quad \dot p_j(z_l) = T\frac{\tau_{j+1}-\tau_j}{2}\cdot F(p_j(z_l))\tag{$E_j^2$}.\]

The nodes $(z_l)$ are associated with a Gauss–Legendre quadrature.

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

We have $p_j(\sigma_k) = x_{j, k}$

Number of unknowns

Putting the period unknown aside, we have to find the $x_{j,k}$ which gives $n\times N_{tst}\times (m+1)$ unknowns.

The equations $E_j^2$ provides $n\times N_{tst}\times m$ plus the $(N_{tst}-1)\times n$ equations for the continuity equations. This makes a total of $(N_{tst}-1)\times m\times n+n\times N_{tst}\times m = n[N_{tst}(m+1)-1]$ equations to which we add the $n$ equations for the periodic boundary condition. In total, we have

\[n\times N_{tst}\times (m+1)\]

equations which matches the number of unknowns.

Phase condition

To ensure uniqueness of the solution to the functional, we use the following phase condition

\[\frac{1}{T} \int_{0}^{T}\left\langle x(s), \dot x_0(s)\right\rangle d s =0\]

During continuation at step $k$, we use $\frac{1}{T} \int_{0}^{T}\left\langle x(s), \dot x_{k-1}(s)\right\rangle d s$

Precision of the scheme

It is known that the error between the true periodic orbit $\gamma^*$ and the approximated one $\gamma$ scales as

\[\sup\limits_i\|\gamma^*(t_i) - \gamma(t_i) \| = O(h^{2m+1}).\]

Discretization of the BVP and jacobian

We only focus on the differential part. Summing up, we obtained the following equations for the $x_{j,l}\in\mathbb R^n$:

\[\sum\limits_{k=1}^{m+1}\mathcal L_k'(z_l)x_{j,k} = F\left(\sum\limits_{k=1}^{m+1}\mathcal L_k(z_l)x_{j,k}\right)\]

The jacobian in the case $m=2$ is given by:

\[\begin{array}{llllllll} x_{0,0} & x_{0,1} & x_{1,0} & x_{1,1} & x_{2,0} & x_{2,1} & x_{3,0} &\quad \mathbf{T} \end{array} \\ \left(\begin{array}{llllllll} H_{0,0}^0 & H_{0,1}^0 & H_{1,0}^0 & & & & & * \\ H_{0,0}^1 & H_{0,1}^1 & H_{1,0}^1 & & & & & * \\ & & H_{1,0}^0 & H_{1,1}^0 & H_{2,0}^0 & & & * \\ & & H_{1,0}^1 & H_{1,1}^1 & H_{2,0}^1 & & & * \\ & & & & H_{2,0}^0 & H_{2,1}^0 & H_{3,0}^0 & * \\ & & & & H_{2,0}^1 & H_{2,1}^1 & H_{3,0}^1 & * \\ & & & & & & & * \\ -I & & & & & & I & * \\ * & * & * & * & * & * & * & * \end{array}\right)\]

where

\[H_{k,l}^{l_2} = \mathcal L'_{l_2,l}\cdot I_n - T\frac{\tau_{j+1}-\tau_j}{2}\cdot\mathcal L_{l_2,l}\cdot dF\left(x_{k,l}\right)\in\mathbb R^n.\]

Interpolation

struct POSolution{Tpb, Tx, Tp}

sol = BifurcationKit.POSolution(prob_coll, x)
sol(t)

Mesh adaptation

The goal of this method^[Russell] is to adapt the mesh $\tau_i$ in order to minimize the error. It is particularly helpful near homoclinic solutions where the period diverges. It can also be useful in order to use a smaller $N_{tst}$.

Encoding of the functional

The functional is encoded in the composite type Collocation. See the link for more information, in particular on how to access the underlying functional, its jacobian...

Jacobians

We provide many different jacobians to take advantage of the problem or the dimensionality. These jacobians are available through the argument jacobian in the constructor of Collocation. For example, you can pass jacobian = FullSparse(). Note that all the internal linear solvers and jacobians are set up automatically so you don't need to do anything. However, for the sake of explanation, we detail how this works.

1. DenseAnalytical()

The jacobian is computed with an analytical formula, it works for dense matrices. This is the default algorithm.

1. DenseAnalyticalInplace()

Same as 1. but caches more information to limit allocations. (Meant to become the default jacobian in the future).

1. AutoDiffDense()

The jacobian is computed with automatic differentiation, it works for dense matrices. Can be used for debugging.

1. FullSparse()

The jacobian is computed with an analytical formula, it works for sparse matrices.

1. FullSparseInplace()

The sparse jacobian is computed in place, limiting memory allocations, with an analytical formula when the sparsity of the jacobian of the vector field is constant. This is much faster than FullSparse().

Linear solvers

When the jacobian is dense and dimension is small ~O(10), you should use condensation of parameters COPBLS for performance.

1. DefaultLS() can be used for most jacobians
2. COPLS() or COPBLS() which is the method of condensation of parameters (COP) implemented in Auto-07p ^[Govaerts]. For this to be most efficient, the vector field must be written in non-allocating form. In this case, this is almost 10x faster than using DefaultLS() for low dimensional ODEs (say dim = O(10)). The performances are closing on to that of Auto-07p.
3. You can use bordered linear solvers in large dimensions to take advantage of the specific shape of the jacobian. See also Trapezoid method for additional information.

Floquet multipliers computation

We provide three methods to compute the Floquet coefficients.

1. The algorithm (Default) FloquetColl is based on the method of condensation of parameters (COP) described in ^[Doedel]. It is the fastest method. It can also be used for sparse jacobians.
2. The algorithm FloquetCollGEV is a simplified version of the procedure described in ^[Fairgrieve]. It boils down to solving a large generalized eigenvalue problem. There is clearly room for improvements here but this can be used to check the results of the previous method.

These methods allow to detect bifurcations of periodic orbits. It seems to work reasonably well for the tutorials considered here. However they may be imprecise^[Lust].

1. 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

newton(
    coll::Collocation,
    orbitguess,
    options::NewtonPar;
    kwargs...
) -> NonLinearSolution{_A, Tprob, Tres} where {_A, Tprob<:(BifurcationKit.PeriodicOrbitFunctionalColl{Tdisc, _A, _B, Nothing, Nothing} where {Tdisc<:Collocation, _A, _B}), Tres<:(Vector)}

We provide a simplified call to newton to locate the periodic orbits. newton will look for prob.jacobian in order to select the requested way to compute the jacobian.

The docs for this specific newton are located at newton.

Continuation

We refer to continuation for more information regarding the arguments. continuation will look for prob.jacobian in order to select the requested way to compute the jacobian.

Examples of optimized use

We show in this part how to optimize the use of collocation for small systems. We consider a small dynamical system which has a Hopf bifurcation and we compute the branch of periodic orbits from the Hopf bifurcation.

The basic code is as follows (see Tutorials).

using BifurcationKit, BenchmarkTools

@inbounds function auto_cat!(dz, z, p, t = 0)
    (;κ, μ, σ, δ) = p
    a,b,c = z
    ab2 = a*b^2
    dz[1] = μ *(κ+c) - ab2 - a
    dz[2] = (ab2 + a - b)/σ
    dz[3] = (b-c)/δ
    dz
end

prob = ODEBifProblem(auto_cat!, [1.,0,1], (κ = 65., μ = 0.01, σ = 5e-3, δ = 2e-2 ) , (@optic _.μ);
        record_from_solution = (x, p; k...) -> (E = x[1], x = x[2], u = x[3]),
        )

opts_br = ContinuationPar(p_min = -0.0, p_max = .19, ds = 0.04, dsmax = 0.1, n_inversion = 8, detect_bifurcation = 3, max_bisection_steps = 25, nev = 3)
br = continuation(prob, PALC(tangent=Bordered()), opts_br; plot = false, normC = norminf, bothside = true)

 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 95
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter μ starts at -0.0, ends at 0.19
 ├─ Algo: PALC [Bordered]
 └─ Special points:

- #  1, endpoint at μ ≈ -0.00000000,                                                                     step =   0
- #  2,     hopf at μ ≈ +0.01531118 ∈ (+0.01531102, +0.01531118), |δp|=2e-07, [converged], δ = ( 2,  2), step =  11
- #  3,     hopf at μ ≈ +0.17443578 ∈ (+0.17443577, +0.17443578), |δp|=7e-09, [converged], δ = (-2, -2), step =  85
- #  4, endpoint at μ ≈ +0.19000000,                                                                     step =  94

# collocation method to discretize the periodic orbit
po_method = Collocation(80, 4,
                jacobian = BifurcationKit.DenseAnalyticalInplace(),
                );

opts_po_cont = ContinuationPar(dsmax = 0.03,
                            ds= -0.001,
                            dsmin = 1e-4,
                            p_max = 0.8,
                            p_min=-5.,
                            max_steps = 280,


                            newton_options = NewtonPar(tol = 1e-9,  max_iterations = 4),
                            nev = 3,
                            tol_stability = 1e-3,
                            detect_bifurcation = 3,
                            plot_every_step = 20,
                            save_sol_every_step = 1,
                            detect_fold = false,
                            n_inversion = 6)

br_po = @benchmark continuation(
    br, 3, opts_po_cont,
    po_method,
    alg = PALC(),
    δp = 0.001,
    normC = norminf,
    # under-optimized Floquet computation. ~Slow but general.
    eigsolver = BifurcationKit.FloquetColl(small_n = false),
    )

BenchmarkTools.Trial: 1 sample with 1 evaluation per sample.
 Single result which took 8.186 s (2.07% GC) to evaluate,
 with a memory estimate of 14.40 GiB, over 2727983 allocations.

The slow running time comes from \ used to invert the collocation jacobian.

Inplace jacobian

The first optimization is to pass an inplace jacobian of the vector field. Right now, ForwardDiff.jl provides this in BifurcationProblem but it is not very optimized.

@inbounds function auto_cat_jac!(J, z, p, t = 0)
    (;κ, μ, σ, δ) = p
    a,b,c = z
    ab = 2*a*b
    J[1,1] = -1-(b^2)
    J[1,2] = -ab
    J[1,3] = μ
    J[2,1] = (1+b^2) / σ
    J[2,2] = (-1+ab) / σ
    # J[2,3] = 0
    # J[3,1] = 0
    J[3,2] = 1 / δ
    J[3,3] = -1 / δ
    J
end

prob = BifurcationProblem(auto_cat!, [1.,0,1], (κ = 65., μ = 0.01, σ = 5e-3, δ = 2e-2 ) , (@optic _.μ);
        record_from_solution = (x, p; k...) -> (E = x[1], x = x[2], u = x[3]),
        J! = auto_cat_jac!,
        )

br = continuation(prob, PALC(tangent=Bordered()), opts_br; plot = false, normC = norminf, bothside = true)

br_po = @benchmark continuation(
    br, 3, opts_po_cont,
    po_method,
    alg = PALC(),
    δp = 0.001,
    normC = norminf,
    # under-optimized Floquet computation. ~Slow but general.
    eigsolver = BifurcationKit.FloquetColl(small_n = false),
    )

BenchmarkTools.Trial: 1 sample with 1 evaluation per sample.
 Single result which took 7.968 s (1.77% GC) to evaluate,
 with a memory estimate of 14.19 GiB, over 167663 allocations.

This is not a large gain.

Optimized linear solver

The next optimization is to used the specific linear solver based on condensation of parameters:

br_po = @benchmark continuation(
    br, 3, opts_po_cont,
    po_method,
    alg = PALC(),
    δp = 0.001,
    linear_algo = BifurcationKit.COPBLS(),
    eigsolver = BifurcationKit.FloquetColl(small_n = false),
    normC = norminf,
    )

BenchmarkTools.Trial: 6 samples with 1 evaluation per sample.
 Range (min … max):  681.605 ms …    1.685 s  ┊ GC (min … max):  1.98% … 56.90%
 Time  (median):     689.301 ms               ┊ GC (median):     2.49%
 Time  (mean ± σ):   856.173 ms ± 406.350 ms  ┊ GC (mean ± σ):  20.31% ± 22.24%

  █                                                              
  █▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆ ▁
  682 ms           Histogram: frequency by time          1.69 s <

 Memory estimate: 2.26 GiB, allocs estimate: 165363.

This is massive gain, almost 10x faster!

Optimized Floquet computation

br_po = @benchmark continuation(
    br, 3, opts_po_cont,
    po_method,
    alg = PALC(),
    δp = 0.001,
    # verbosity = 3,
    linear_algo = BifurcationKit.COPBLS(),
    eigsolver = BifurcationKit.FloquetColl(small_n = true),
    normC = norminf,
    )

BenchmarkTools.Trial: 12 samples with 1 evaluation per sample.
 Range (min … max):  413.072 ms … 453.508 ms  ┊ GC (min … max): 0.00% … 0.72%
 Time  (median):     416.437 ms               ┊ GC (median):    0.67%
 Time  (mean ± σ):   418.791 ms ±  11.085 ms  ┊ GC (mean ± σ):  0.48% ± 0.36%

  ▄  ▁  █                                                        
  █▆▁█▁▁█▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆ ▁
  413 ms           Histogram: frequency by time          454 ms <

 Memory estimate: 128.22 MiB, allocs estimate: 40314.

This is massive gain (typically from 4.6s to 180ms on Apple M2), almost 4x faster, overall 40x faster than basic version. The allocations are quite minimal.

References