# Fixed points & Periodicity

## Fixed points

`ChaosTools.fixedpoints`

— Function`fixedpoints(ds::DynamicalSystem, box, p = ds.p; kwargs...) → fp, eigs, stable`

Return all fixed points `fp`

of the given `ds`

that exist within the state space subset `box`

for parameter configuration `p`

. Fixed points are returned as a `Dataset`

. For convenience, a vector of the Jacobian eigenvalues of each fixed point, and whether the fixed points are stable or not, are also returned.

`box`

is an appropriate `IntervalBox`

from IntervalRootFinding.jl. E.g. for a 3D system it would be something like

```
v, z = -5..5, -2..2 # 1D intervals, can use `interval(-5, 5)` instead
box = v × v × z # `\times = ×`, or use `IntervalBox(v, v, z)` instead
```

Internally IntervalRootFinding.jl is used and as a result we are guaranteed to find all fixed points that exist in `box`

, regardless of stability. Since IntervalRootFinding.jl returns an interval containing a unique fixed point, we return the midpoint of the interval as the actual fixed point. Naturally, limitations inherent to IntervalRootFinding.jl apply here.

The output of `fixedpoints`

can be used in the BifurcationKit.jl as a start of a continuation process. See also `periodicorbits`

.

**Keywords**

`method = IntervalRootFinding.Krawczyk`

configures the root finding method, see the docs of IntervalRootFinding.jl for all posibilities.`tol = 1e-15`

is the root-finding tolerance.`warn = true`

throw a warning if no fixed points are found.

A rather simple example of the fixed points can be demonstrated using E.g., the Lorenz-63 system, whose fixed points can be easily calculated analytically to be the following three

\[(0,0,0) \\ \left( \sqrt{\beta(\rho-1)}, \sqrt{\beta(\rho-1)}, \rho-1 \right) \\ \left( -\sqrt{\beta(\rho-1)}, -\sqrt{\beta(\rho-1)}, \rho-1 \right) \\ \]

So, let's calculate

```
using DynamicalSystems
ρ, β = 30.0, 10/3
ds = Systems.lorenz(; ρ, β)
# Define the box within which to find fixed points:
x = -20..20
y = -20..20
z = 0.0..40
box = x × y × z
fp, eigs, stable = fixedpoints(ds, box)
fp
```

```
3-dimensional Dataset{Float64} with 3 points
9.83192 9.83192 29.0
-9.83192 -9.83192 29.0
3.23647e-17 3.23647e-17 3.43473e-16
```

and compare this with the analytic ones:

```
lorenzfp(ρ, β) = [
SVector(0, 0, 0.0),
SVector(sqrt(β*(ρ-1)), sqrt(β*(ρ-1)), ρ-1),
SVector(-sqrt(β*(ρ-1)), -sqrt(β*(ρ-1)), ρ-1),
]
lorenzfp(ρ, β)
```

```
3-element Vector{SVector{3, Float64}}:
[0.0, 0.0, 0.0]
[9.83192080250175, 9.83192080250175, 29.0]
[-9.83192080250175, -9.83192080250175, 29.0]
```

## Stable and Unstable Periodic Orbits of Maps

Chaotic behavior of low dimensional dynamical systems is affected by the position and the stability properties of the periodic orbits of a dynamical system.

Finding unstable (or stable) periodic orbits of a discrete mapping analytically rapidly becomes impossible for higher orders of fixed points. Fortunately there is a numeric algorithm due to Schmelcher & Diakonos which allows such a computation. Notice that even though the algorithm can find stable fixed points, it is mainly aimed at *unstable* ones.

The functions `periodicorbits`

and `lambdamatrix`

implement the algorithm:

`ChaosTools.periodicorbits`

— Function```
periodicorbits(ds::DiscreteDynamicalSystem,
o, ics [, λs, indss, singss]; kwargs...) -> FP
```

Find fixed points `FP`

of order `o`

for the map `ds`

using the algorithm due to Schmelcher & Diakonos^{[Schmelcher1997]}. `ics`

is a collection of initial conditions (container of vectors) to be evolved.

**Optional Arguments**

The optional arguments `λs, indss, singss`

*must be containers* of appropriate values, besides `λs`

which can also be a number. The elements of those containers are passed to: `lambdamatrix(λ, inds, sings)`

, which creates the appropriate $\mathbf{\Lambda}_k$ matrix. If these arguments are not given, a random permutation will be chosen for them, with `λ=0.001`

.

**Keyword Arguments**

`maxiters::Int = 100000`

: Maximum amount of iterations an i.c. will be iterated before claiming it has not converged.`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.`inftol = 10.0`

: If a state reaches`norm(state) ≥ inftol`

it is assumed that it has escaped to infinity (and is thus abandoned).`roundtol::Int = 4`

: The found fixed points are rounded to`roundtol`

digits before pushed into the list of returned fixed points`FP`

,*if*they are not already contained in`FP`

. This is done so that`FP`

doesn't contain duplicate fixed points (notice that this has nothing to do with`disttol`

). Turn this to`typemax(Int)`

to get the full precision of the algorithm.

**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)\]

with $f$ = `eom`

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

**Performance Notes**

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

`ChaosTools.lambdamatrix`

— Function`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 function `periodicorbits`

.

**Arguments**

`λ<:Real`

: the multiplier of the $C_k$ matrix, with`0<λ<1`

.`inds::Vector{Int}`

: The`i`

th entry of this vector gives the*row*of the nonzero element of the`i`

th column of $C_k$.`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.^{[Pingel2000]} 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^{[Diakonos1998]} for details.

`ChaosTools.lambdaperms`

— Function`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`

).

### Standard Map example

For example, let's find the fixed points of the `Systems.standardmap`

of order 2, 3, 4, 5, 6 and 8. We will use all permutations for the `signs`

but only one for the `inds`

. We will also only use one `λ`

value, and a 21×21 density of initial conditions.

First, initialize everything

```
using DynamicalSystems
ds = Systems.standardmap()
xs = range(0, stop = 2π, length = 11); ys = copy(xs)
ics = [SVector{2}(x,y) for x in xs for y in ys]
# All permutations of [±1, ±1]:
singss = lambdaperms(2)[2] # second entry are the signs
# I know from personal research I only need this `inds`:
indss = [[1,2]] # <- must be container of vectors!
λs = 0.005 # <- only this allowed to not be vector (could also be vector)
orders = [2, 3, 4, 5, 6, 8]
ALLFP = Dataset{2, Float64}[];
```

`Dataset{2, Float64}[]`

Then, do the necessary computations for all orders

```
for o in orders
FP = periodicorbits(ds, o, ics, λs, indss, singss)
push!(ALLFP, FP)
end
```

Plot the phase space of the standard map

```
using CairoMakie
iters = 1000
dataset = trajectory(ds, iters)
for x in xs
for y in ys
append!(dataset, trajectory(ds, iters, SVector{2}(x, y)))
end
end
fig = Figure()
ax = Axis(fig[1,1]; xlabel = L"\theta", ylabel = L"p",
limits = ((xs[1],xs[end]), (xs[1],xs[end]))
)
scatter!(ax, dataset[:, 1], dataset[:, 2]; markersize = 1, color = "black")
fig
```

and finally, plot the fixed points

```
markers = [:diamond, :utriangle, :rect, :pentagon, :hexagon, :circle]
for i in 1:6
FP = ALLFP[i]
o = orders[i]
scatter!(ax, columns(FP)...; marker=markers[i], color = Cycled(i),
markersize = 30 - 2i, strokecolor = "grey", strokewidth = 1, label = "order $o"
)
end
axislegend(ax)
fig
```

Okay, this output is great, and we can easily tell that it is correct for many reasons:

- It is the same fig. 12 of this publication.
- Fixed points of order $n$ are also fixed points of order $2n, 3n, 4n, ...$
- Besides fixed points of previous orders,
*original*fixed points of order $n$ come in (possible multiples of) $2n$-sized pairs (see e.g. order 5). This is a direct consequence of the Poincaré–Birkhoff theorem.

## Estimating the Period

The function `estimate_period`

offers ways for estimating the period (either exact for periodic timeseries, or approximate for near-periodic ones) of a given timeseries. We offer five methods to estimate periods, some of which work on evenly sampled data only, and others which accept any data. The figure below summarizes this:

`ChaosTools.estimate_period`

— Function`estimate_period(v::Vector, method, t=0:length(v)-1; kwargs...)`

Estimate the period of the signal `v`

, with accompanying time vector `t`

, using the given `method`

.

If `t`

is an AbstractArray, then it is iterated through to ensure that it's evenly sampled (if necessary for the algorithm). To avoid this, you can pass any `AbstractRange`

, like a `UnitRange`

or a `LinRange`

, which are defined to be evenly sampled.

**Methods requiring evenly sampled data**

These methods are faster, but some are error-prone.

`:periodogram`

or`:pg`

: Use the fast Fourier transform to compute a periodogram (power-spectrum) of the given data. Data must be evenly sampled.`:multitaper`

or`mt`

: The multitaper method reduces estimation bias by using multiple independent estimates from the same sample. Data tapers are then windowed and the power spectra are obtained. Available keywords follow:`nw`

is the time-bandwidth product, and`ntapers`

is the number of tapers. If`window`

is not specified, the signal is tapered with`ntapers`

discrete prolate spheroidal sequences with time-bandwidth product`nw`

. Each sequence is equally weighted; adaptive multitaper is not (yet) supported. If`window`

is specified, each column is applied as a taper. The sum of periodograms is normalized by the total sum of squares of`window`

.`:autocorrelation`

or`:ac`

: Use the autocorrelation function (AC). The value where the AC first comes back close to 1 is the period of the signal. The keyword`L = length(v)÷10`

denotes the length of the AC (thus, given the default setting, this method will fail if there less than 10 periods in the signal). The keyword`ϵ = 0.2`

(`\epsilon`

) means that`1-ϵ`

counts as "1" for the AC.`:yin`

: The YIN algorithm. An autocorrelation-based method to estimate the fundamental

period of the signal. See the original paper ^{[CheveigneYIN2002]} or the implementation `yin`

.

speech and music. The Journal of the Acoustical Society of America, 111(4), 1917-1930.

**Methods not requiring evenly sampled data**

These methods tend to be slow, but versatile and low-error.

`:lombscargle`

or`:ls`

: Use the Lomb-Scargle algorithm to compute a periodogram. The advantage of the Lomb-Scargle method is that it does not require an equally sampled dataset and performs well on undersampled datasets. Constraints have been set on the period, since Lomb-Scargle tends to have false peaks at very low frequencies. That being said, it's a very flexible method. It is extremely customizable, and the keyword arguments that can be passed to it are given in the documentation.`:zerocrossing`

or`:zc`

: Find the zero crossings of the data, and use the average difference between zero crossings as the period. This is a naïve implementation, with only linear interpolation; however, it's useful as a sanity check. The keyword`line`

controls where the "crossing point" is. It defaults to`mean(v)`

.

For more information on the periodogram methods, see the documentation of DSP.jl and LombScargle.jl.

`ChaosTools.yin`

— Function```
yin(sig::Vector, sr::Int; w_len::Int = 512, f_step::Int = 256, f0_min = 100, f0_max = 500,
harmonic_threshold = 0.1,
difference_function::Function = difference_function_original, kwargs...) -> F0s, frame_times
```

Estimates the fundamental frequency (F0) of the signal `sig`

using the YIN algorithm ^{[1]}. The signal `sig`

is a vector of points uniformly sampled at a rate `sr`

.

**Keyword arguments**

`w_len`

: size of the analysis window [samples == number of points]`f_step`

: size of the lag between two consecutive frames [samples == number of points]`f0_min`

: Minimum fundamental frequency that can be detected [linear frequency]`f0_max`

: Maximum fundamental frequency that can be detected [linear frequency]`harmonic_threshold`

: Threshold of detection. The algorithm returns the first minimum of the CMNDF function below this threshold.`diffference_function`

: The difference function to be used (by default`ChaosTools.difference_function_original`

).

**Description**

The YIN algorithm ^{[CheveigneYIN2002]} estimates the signal's fundamental frequency `F0`

by basically looking for the period `τ0`

which minimizes the signal's autocorrelation. This autocorrelation is calculated for signal segments (frames), composed of two windows of length `w_len`

. Each window is separated by a distance `τ`

, and the idea is that the distance which minimizes the pairwise difference between each window is considered to be the fundamental period `τ0`

of that frame.

More precisely, the algorithm first computes the cumulative mean normalized difference function (MNDF) between two windows of a frame for several candidate periods `τ`

ranging from `τ_min=sr/f0_max`

to `τ_max=sr/f0_min`

. The MNDF is defined as

\[d_t^\prime(\tau) = \begin{cases} 1 & \text{if} ~ \tau=0 \\ d_t(\tau)/\left[{(\frac 1 \tau) \sum_{j=1}^{\tau} d_{t}(j)}\right] & \text{otherwise} \end{cases}\]

where `d_t`

is the difference function:

\[d_t(\tau) = \sum_{j=1}^W (x_j - x_{j+\tau})^2\]

It then refines the local minima of the MNDF using parabolic (quadratic) interpolation. This is done by taking each minima, along with their first neighbor points, and finding the minimum of the corresponding interpolated parabola. The MNDF minima are substituted by the interpolation minima. Finally, the algorithm chooses the minimum with the smallest period and with a corresponding MNDF below the `harmonic threshold`

. If this doesn't exist, it chooses the period corresponding to the global minimum. It repeats this for frames starting at the first signal point, and separated by a distance `f_step`

(frames can overlap), and returns the vector of frequencies `F0=sr/τ0`

for each frame, along with the start times of each frame.

As a note, the physical unit of the frequency is 1/[time], where [time] is decided by the sampling rate `sr`

. If, for instance, the sampling rate is over seconds, then the frequency is in Hertz.

speech and music. The Journal of the Acoustical Society of America, 111(4), 1917-1930.

`ChaosTools.difference_function_original`

— Function`difference_function_original(x, W, τmax) -> df`

Computes the difference function of `x`

. `W`

is the window size, and `τmax`

is the maximum period. This corresponds to equation (6) in ^{[1]}:

\[d_t(\tau) = \sum_{j=1}^W (x_j - x_{j+\tau})^2\]

### Example

Here we will use a modified FitzHugh-Nagumo system that results in periodic behavior, and then try to estimate its period. First, let's see the trajectory:

```
using DynamicalSystems, CairoMakie
function FHN(u, p, t)
e, b, g = p
v, w = u
dv = min(max(-2 - v, v), 2 - v) - w
dw = e*(v - g*w + b)
return SVector(dv, dw)
end
g, e, b = 0.8, 0.04, 0.0
p0 = [e, b, g]
fhn = ContinuousDynamicalSystem(FHN, SVector(-2, -0.6667), p0)
T, Δt = 1000.0, 0.1
v = trajectory(fhn, T; Δt)[:, 1]
t = 0:Δt:T
fig, ax = lines(0:Δt:T, v)
fig
```

Examining the figure, one can see that the period of the system is around `91`

time units. To estimate it numerically let's use some of the methods:

`estimate_period(v, :autocorrelation, t)`

`91.0`

`estimate_period(v, :periodogram, t)`

`91.62720091627202`

`estimate_period(v, :zerocrossing, t)`

`91.07000000000001`

`estimate_period(v, :yin, t; sr=round(Int, (1/Δt)), f0_min=0.01)`

`91.07002491049154`

- Schmelcher1997P. Schmelcher & F. K. Diakonos, Phys. Rev. Lett.
**78**, pp 4733 (1997) - Pingel2000D. Pingel
*et al.*, Phys. Rev. E**62**, pp 2119 (2000) - Diakonos1998F. K. Diakonos
*et al.*, Phys. Rev. Lett.**81**, pp 4349 (1998) - CheveigneYIN2002De Cheveigné, A., & Kawahara, H. (2002). YIN, a fundamental frequency estimator for
- CheveigneYIN2002De Cheveigné, A., & Kawahara, H. (2002). YIN, a fundamental frequency estimator for