The Public API

periodic_orbit(ds::DynamicalSystem, alg::PeriodicOrbitFinder [, ig::InitialGuess]) → PeriodicOrbit

Try to find a single periodic orbit of the dynamical system ds using the algorithm alg and optionally given some InitialGuess ig which defaults to InitialGuess(ds). Return the result as a PeriodicOrbit.

Depending on alg, it is not guaranteed that a periodic orbit will be found given ds, ig. If one is not found, nothing is returned instead.

For more details on the periodic orbit detection algorithm, see the documentation the alg.

periodic_orbits(ds::DynamicalSystem, alg::PeriodicOrbitFinder [, igs]::Vector{InitialGuess} = InitialGuess(ds)) → Vector{PeriodicOrbit}

Try to find multiple periodic orbits of the dynamical system ds using the algorithm alg and optionally given a Vector of initial guesses which defaults to [InitialGuess(ds)]. given some initial guesses igs. Return the result as a Vector of PeriodicOrbit.

This function exists because some algorithms optimize/specialize on detecting multiple periodic orbits.

PeriodicOrbitFinder

Supertype for all the periodic orbit detection algorithms. Each of the concrete subtypes of PeriodicOrbitFinder represents one given algorithm for detecting periodic orbits. This subtype includes all the necessary metaparameters for the algorithm to work.

InitialGuess

A structure that contains an initial guess for a periodic orbit detection algorithms.

• u0::AbstractArray{<:Real} - guess of a point in the periodic orbit
• T::Union{Real, Nothing} - guess of period of the orbit. Some algorithms do not require the period guess to be given, in which case T is set to nothing.

PeriodicOrbit

A structure that contains information about a periodic orbit.

• points::StateSpaceSet - points of the periodic orbit. This container always holds the complete orbit (in the sense of being continuously sampled with some sampling Δt for continuous time systems).
• T::Real - the period of the orbit
• stable::Union{Bool, Missing} - local stability of the periodic orbit. If the stability is unknown (not computed automatically by an algorithm), this is set to missing.

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.

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. The algorithm we use is NonlinearSolve.LevenbergMarquardt().

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.

postability(ds::CoreDynamicalSystem, po::PeriodicOrbit [, jac]) → new_po

Determine the local stability of the periodic orbit po using the jacobian rule jac. Return a new periodic orbit for which po.stable is set to true if the periodic orbit is stable or false if it is unstable.

For discrete-time systems, the stability is determined using eigenvalues of the jacobian of po.T-th iterate of the dynamical system ds at the point po.points[1]. If the maximum absolute value of the eigenvalues is less than 1, the periodic orbit is marked as stable.

For continuous-time systems, the stability is determined by the Floquet multipliers of the monodromy matrix. If the maximum absolute value of the Floquet multipliers is less than 1 (while neglecting the multiplier which is always 1), the periodic orbit is marked as stable.

The default value of jacobian rule jac is obtained via automatic differentiation.

minimal_period(ds::DynamicalSystem, po::PeriodicOrbit; kw...) → minT_po

Compute the minimal period of the periodic orbit po of the dynamical system ds. Return the periodic orbit minT_po with the minimal period. In the literature, minimal period is also called prime, principal or fundamental period.

Keyword arguments

• atol = 1e-4 : After stepping the point u0 for a time T, it must return to atol neighborhood of itself to be considered periodic.
• maxiter = 40 : Maximum number of Poincare map iterations. Continuous-time systems only. If the number of Poincare map iterations exceeds maxiter, but the point u0 has not returned to atol neighborhood of itself, the original period po.T is returned.
• Δt = missing : The time step between points in the trajectory minT_po.points. If Δt is missing, then Δt=minT_po.T/100 is used. Continuous-time systems only.

Description

For discrete systems, a valid period would be any natural multiple of the minimal period. Hence, all natural divisors of the period po.T are checked as a potential period. A point u0 of the periodic orbit po is iterated n times and if the distance between the initial point u0 and the final point is less than atol, the period of the orbit is n.

For continuous systems, a point u0 of the periodic orbit is integrated for a very short time. The resulting point u1 is used to create a normal vector a=(u1-u0) to a hyperplane perpendicular to the trajectory at u0. A Poincare map is created using this hyperplane. Using the Poincare map, the hyperplane crossings are checked. Time of the first crossing that is within atol distance of the initial point u0 is the minimal period. At most maxiter crossings are checked.

podistance(po1::PeriodicOrbit, po2::PeriodicOrbit, [, distance]) → Real

Compute the distance between two periodic orbits po1 and po2. Periodic orbits po1 and po2 and the dynamical system ds all have to be either discrete-time or continuous-time. Distance between the periodic orbits is computed using the given distance function distance. The default distance function is StrictlyMinimumDistance(true, Euclidean()) which finds the minimal Euclidean distance between any pair of points where one point belongs to po1 and the other to po2. For other options of the distance function, see StateSpaceSets.set_distance. Custom distance function can be provided as well.

poequal(po1::PeriodicOrbit, po2::PeriodicOrbit; kwargs...) → true/false

Return true if the periodic orbits po1 and po2 are approximately equal in period and in location.

Keyword arguments

• Tthres=1e-3 : difference in periods of the periodic orbits must be less than this threshold
• dthres=1e-3 : distance between periodic orbits must be less than this threshold
• distance : distance function used to compute the distance between the periodic orbits

Distance between the orbits is computed using podistance with distance.

uniquepos(pos::Vector{<:PeriodicOrbit}; atol=1e-6) → Vector{PeriodicOrbit}

Return a vector of unique periodic orbits from the vector pos of periodic orbits. By unique we mean that the distance between any two periodic orbits in the vector is greater than atol. To see details about the distance function, see podistance.

Keyword arguments

• atol : minimal distance between two periodic orbits for them to be considered unique.

isdiscretetime(po::PeriodicOrbit) → true/false

Return true if the periodic orbit belongs to a discrete-time dynamical system, false if it belongs to a continuous-time dynamical system. This simple function only checks whether the period is an integer or not.

isdiscretetime(ds::DynamicalSystem) → true/false

Return true if ds operates in discrete time, or false if it is in continuous time. This is information deduced from the type of ds.