Algorithms

OptimizedShooting(; kwargs...)

A shooting method (Dednam and Botha, 2014) that uses Levenberg-Marquardt optimization to find periodic orbits of continuous-time dynamical systems.

Keyword arguments

• Δt::Float64 = 1e-6: step between the points in the residual R. See below for details.
• n::Int64 = 2: n*dimension(ds) is the number of points in the residual R. See below for details.
• nonlinear_solve_kwargs = (reltol=1e-6, abstol=1e-6, maxiters=1000): keyword arguments to pass to the solve function from NonlinearSolve.jl. For details on the keywords see the respective package documentation.

Description

Let us consider the following continuous-time dynamical system

\[\frac{dx}{dt} = f(x, p, t)\]

Dednam and Botha (Dednam and Botha, 2014) suggest to minimize the residual $R$ defined as

\[R = (x(T)-x(0), x(T+\Delta t)-x(\Delta t), \dots, x(T+(n-1)\Delta t)-x((n-1)\Delta t))\]

where $T$ is unknown period of a periodic orbit and $x(t)$ is a solution at time $t$ given some unknown initial point. Initial guess of the period $T$ and the initial point is optimized by the Levenberg-Marquardt algorithm.

In our implementation, the keyword argument n corresponds to $n$ in the residual $R$. The keyword argument Δt corresponds to $\Delta t$ in the residual $R$.

Note that for the algorithm to converge to a periodic orbit, the initial guess has to be close to an existing periodic orbit.

SchmelcherDiakonos(; kwargs...)

Detect periodic orbits of ds <: DiscreteTimeDynamicalSystem using algorithm proposed by Schmelcher & Diakonos (Schmelcher and Diakonos, 1997).

Keyword arguments

• o : period of the periodic orbit
• λs : vector of λ parameters, see (Schmelcher and Diakonos, 1997) for details
• indss : vector of vectors of indices for the permutation matrix
• signss : vector of vectors of signs for the permutation matrix
• maxiters=10000 : maximum amount of iterations an initial guess will be iterated before claiming it has not converged
• inftol=10.0 : if a state reaches norm(state) ≥ inftol it is assumed that it has escaped to infinity (and is thus abandoned)
• disttol=1e-10 : distance tolerance. If the 2-norm of a previous state with the next one is ≤ disttol then it has converged to a fixed point

Description

The algorithm used can detect periodic orbits by turning fixed points of the original map ds to stable ones, through the transformation

\[\mathbf{x}_{n+1} = \mathbf{x}_n + \mathbf{\Lambda}_k\left(f^{(o)}(\mathbf{x}_n) - \mathbf{x}_n\right)\]

The index $k$ counts the various possible $\mathbf{\Lambda}_k$.

Performance notes

All initial guesses are evolved for all $\mathbf{\Lambda}_k$ which can very quickly lead to long computation times.

lambdamatrix(λ, inds::Vector{Int}, sings) -> Λk

Return the matrix $\mathbf{\Lambda}_k$ used to create a new dynamical system with some unstable fixed points turned to stable in the algorithm SchmelcherDiakonos.

Arguments

1. λ<:Real : the multiplier of the $C_k$ matrix, with 0 < λ < 1.
2. inds::Vector{Int} : The i-th entry of this vector gives the row of the nonzero element of the ith column of $C_k$.
3. sings::Vector{<:Real} : The element of the i-th column of $C_k$ is +1 if signs[i] > 0 and -1 otherwise (sings can also be Bool vector).

Calling lambdamatrix(λ, D::Int) creates a random $\mathbf{\Lambda}_k$ by randomly generating an inds and a signs from all possible combinations. The collections of all these combinations can be obtained from the function lambdaperms.

Description

Each element of inds must be unique such that the resulting matrix is orthogonal and represents the group of special reflections and permutations.

Deciding the appropriate values for λ, inds, sings is not trivial. However, in ref.(Pingel et al., 2000) there is a lot of information that can help with that decision. Also, by appropriately choosing various values for λ, one can sort periodic orbits from e.g. least unstable to most unstable, see (Diakonos et al., 1998) for details.

lambdaperms(D) -> indperms, singperms

Return two collections that each contain all possible combinations of indices (total of $D!$) and signs (total of $2^D$) for dimension D (see lambdamatrix).

DavidchackLai(; kwargs...)

Find periodic orbits fps of periods 1 to n+1 for the dynamical system ds using the algorithm propesed by Davidchack & Lai (Davidchack and Lai, 1999).

Keyword arguments

• n::Int64 : Periodic orbits of period up to n will be detected. Some (but not all) POs of period n+1 will be detected. Keyword argument n must be a positive integer.
• m::Int64 : Initial guesses will be used to find POs of period 1 to m. These periodic orbits will then be used to detect periodic orbits of periods from m+1 to n+1. Keyword argument m must be a positive integer between 1 and n.
• β = nothing: If it is nothing, then β(n) = 10*1.2^n. Otherwise can be a function that takes period n and return a number. It is a parameter mentioned in the paper(Davidchack and Lai, 1999).
• maxiters = nothing: If it is nothing, then initial condition will be iterated max(100, 4*β(p)) times (where p is the period of the periodic orbit) before claiming it has not converged. If it is an integer, then it is the maximum amount of iterations an initial condition will be iterated before claiming it has not converged.
• disttol = 1e-10: Distance tolerance. If norm(f^{n}(x)-x) < disttol where f^{n} is the n-th iterate of the dynamic rule f, then x is an n-periodic point.
• abstol = 1e-8: A detected periodic point isn't stored if it is in abstol neighborhood of some previously detected point. Distance is measured by euclidian norm. If you are getting duplicate periodic points, increase this value.

Description

The algorithm is an extension of Schmelcher & Diakonos(Schmelcher and Diakonos, 1997) implemented as SchmelcherDiakonos.

The algorithm can detect periodic orbits by turning fixed points of the original dynamical system ds to stable ones, through the transformation

\[\mathbf{x}_{n+1} = \mathbf{x}_{n} + [\beta |g(\mathbf{x}_{n})| C^{T} - J(\mathbf{x}_{n})]^{-1} g(\mathbf{x}_{n})\]

where

\[g(\mathbf{x}_{n}) = f^{n}(\mathbf{x}_{n}) - \mathbf{x}_{n}\]

and

\[J(\mathbf{x}_{n}) = \frac{\partial g(\mathbf{x}_{n})}{\partial \mathbf{x}_{n}}\]

The main difference between SchmelcherDiakonos and DavidchackLai is that the latter uses periodic points of previous period as seeds to detect periodic points of the next period. Additionally, SchmelcherDiakonos only detects periodic points of a given period, while davidchacklai detects periodic points of all periods up to n.

Important note

For low periods n circa less than 6, you should select m = n otherwise the algorithm won't detect periodic orbits correctly. For higher periods, you can select m as 6. We recommend experimenting with m as it may depend on the specific problem. Increase m in case the orbits are not being detected correctly.

Initial guesses for this algorithm can be selected as a uniform grid of points in the state space or subset of a chaotic trajectory.