Attractor Basins, Tipping Points
In this page we list several functions related with basins of attraction and tipping points. In the example 2D basins of higher dimensional system we try to apply every single function listed below, so check this for an example application of everything listed here!
Computing basins of attraction
ChaosTools.basins_of_attraction
— Functionbasins_of_attraction(grid::Tuple, ds::DynamicalSystem; kwargs...) -> basins, attractors
Compute an estimate of the basins of attraction of a dynamical system ds
on a partitioning of the state space given by grid
. The method has been inspired by the 2D grid approach developed by Nusse & Yorke [Yorke1997]. It works without knowledge of where attractors are; it identifies them automatically.
The dynamical system can be:
- An actual
DiscreteDynamicalSystem
orContinuousDynamicalSystem
. - A Poincaré map of
ContinuousDynamicalSystem
:poincaremap
. - A stroboscopic map, i.e. a periodically forced
ContinuousDynamicalSystem
(see examples for this particular application).
grid
is a tuple of ranges defining the grid of initial conditions, for example grid = (xg, yg)
where xg = yg = range(-5, 5; length = 100)
. The grid is not necessarilly of the same dimension as the state space, attractors can be found in lower dimensional projections.
The output basins
is an integer-valued array on the grid
, with its entries labelling which basin of attraction the given grid point belongs to. The output attractors
is a dictionary whose keys correspond to the attractor number and the values contains the points of the attractors found. Notice that for some attractors this list may be incomplete.
See also match_attractors!
, basin_fractions
, tipping_probabilities
.
Keyword Arguments
Δt
: Approximate time step of the integrator, which is1
for discrete systems. For continuous systems, an automatic value is calculated usingautomatic_Δt_basins
. See that function for more info.T
: Period of the stroboscopic map, in case of a continuous dynamical system with periodic time forcing. This argument is incompatible withΔt
.idxs = 1:length(grid)
: This vector selects the variables of the system that will define the subspace the dynamics will be projected into.complete_state = zeros(D-Dg)
: This argument allows setting the remaining variables of the dynamical system state on each initial conditionu
, withDg
the dimension of the grid. It can be either a vector of lengthD-Dg
, or a functionf(y)
that returns a vector of lengthD-Dg
given the projected initial condition on the gridy
.diffeq...
: Keyword arguments propagated tointegrator
.mx_chk_att = 2
: A parameter that sets the maximum checks of consecutives hits of an attractor before deciding the basin of the initial condition.mx_chk_hit_bas = 10
: Maximum check of consecutive visits of the same basin of attraction. This number can be increased for higher accuracy.mx_chk_fnd_att = 100
: Maximum check of unnumbered cell before considering we have an attractor. This number can be increased for higher accuracy.mx_chk_loc_att = 100
: Maximum check of consecutive cells marked as an attractor before considering that we have all the available pieces of the attractor.mx_chk_lost
: Maximum check of iterations outside the defined grid before we consider the orbit lost outside. This number can be increased for higher accuracy. It defaults to20
if no attractors are given (see discussion on refining basins), and to1000
if attractors are given.horizon_limit = 1e6
: If the norm of the integrator state reaches this limit we consider that the orbit diverges.show_progress = true
: By default a progress bar is shown using ProgressMeter.jl.attractors, ε
: See discussion on refining basins below.
Description
basins
has the following organization:
- The basins are coded in sequential order from 1 up to the number of attractors.
- If the trajectory diverges or converges to an attractor outside the defined grid it is numbered
-1
attractors
has the following organization:
- The keys of the dictionary correspond to the number of the attractor.
- The value associated to a key is a
Dataset
with the guessed location of the attractor on the state space.
The method starts by picking the first available initial condition on the grid not yet numbered. The dynamical system is then iterated until one of the following conditions happens:
- The trajectory hits a known attractor already numbered
mx_chk_att
consecutive times: the initial condition is numbered with the corresponding number. - The trajectory spends
mx_chk_lost
steps outside the defiend grid: the initial condition is set to -1. - The trajectory hits a known basin
mx_chk_hit_bas
times in a row: the initial condition belongs to that basin and is numbered accordingly. - The trajectory hits
mx_chk_fnd_att
times in a row an unnumbered cell: it is considered an attractor and is labelled with a new number.
Regarding performace, this method is at worst as fast as tracking the attractors. In most cases there is a signicative improvement in speed.
Notice that in the case we have to project the dynamics on a lower dimensional space, there are edge cases where the system may have two attractors that are close on the projected space but are far apart in another dimension. They could be collapsed or confused into the same attractor. This is a drawback of this method.
Refining basins of attraction
Sometimes one would like to be able to refine the found basins of attraction by recomputing basins_of_attraction
on a smaller, and more fine-grained, grid
. If however this new grid
does not contain the attractors, basins_of_attraction
would (by default) attribute the value -1
to all grid points. For these cases, an extra search clause can be provided by setting the keywords attractors, ε
. The attractors
is a dictionary mapping attractor IDs to Dataset
s (i.e., the same as the return value of basins_of_attraction
). The algorithm checks at each step whether the system state is ε
-close (Euclidean norm) to any of the given attractors, and if so it attributes the stating grid point to the basin of the close attractor. By default ε
is equal to the mean grid spacing.
A word of advice while using this method: in order to work properly, ε
should be about the size of a grid cell that has been used to compute the given attractors
. It is recomended to keep the same step size (i.e., use the same integrator) since it may have an influence in some cases. This algorithm is usually slower than the method with the attractors on the grid.
ChaosTools.automatic_Δt_basins
— Functionautomatic_Δt_basins(ds, grid; kwargs...) → Δt
Calculate an optimal Δt
value for basins_of_attraction
. This is done by evaluating the dynamic rule f
(vector field) at N
randomly chosen points of the grid. The average f
is then compared with the diagonal length of a grid cell and their ratio provides Δt
.
Keywords idxs, complete_state, diffeq
are exactly as in basins_of_attraction
, and the keyword N
is 5000
by default.
Notice that Δt
should not be too small which happens typically if the grid resolution is high. It is okay for basins_of_attraction
if the trajectory skips a few cells. But if Δt
is too small the default values for all other keywords such as mx_chk_hit_bas
need to be increased drastically.
Also, Δt
that is smaller than the internal step size of the integrator will cause a performance drop.
ChaosTools.match_attractors!
— Functionmatch_attractors!(b₋, a₋, b₊, a₊, [, method = :distance])
Attempt to match the attractors in basins/attractors b₊, a₊
with those at b₋, a₋
. b
is an array whose values encode the attractor ID, while a
is a dictionary mapping IDs to Dataset
s containing the attractors (e.g. output of basins_of_attraction
). Typically the +,- mean after and before some change of parameter for a system.
In basins_of_attraction
different attractors get assigned different IDs, however which attractor gets which ID is somewhat arbitrary, and computing the basins of the same system for slightly different parameters could label the "same" attractors (at the different parameters) with different IDs. match_attractors!
tries to "match" them by modifying the attractor IDs.
The modification of IDs is always done on the b, a
that have less attractors.
method
decides the matching process:
method = :overlap
matches attractors whose basins before and after have the most overlap (in pixels).method = :distance
matches attractors whose state space distance the smallest.
Final state sensitivity / fractal boundaries
ChaosTools.basins_fractal_dimension
— Functionbasins_fractal_dimension(basins; kwargs...) -> V_ε, N_ε ,d
Estimate the Fractal Dimension d
of the boundary between basins of attraction using the box-counting algorithm.
The output N_ε
is a vector with the number of the balls of radius ε
(in pixels) that contain at least two initial conditions that lead to different attractors. V_ε
is a vector with the corresponding size of the balls. The ouput d
is the estimation of the box-counting dimension of the boundary by fitting a line in the log.(N_ε)
vs log.(1/V_ε)
curve. However it is recommended to analyze the curve directly for more accuracy.
Keyword arguments
range_ε = 2:maximum(size(basins))÷20
is the range of sizes of the ball to test (in pixels).
Description
It is the implementation of the popular algorithm of the estimation of the box-counting dimension. The algorithm search for a covering the boundary with N_ε
boxes of size ε
in pixels.
ChaosTools.basin_entropy
— Functionbasin_entropy(basins, ε = 20) -> Sb, Sbb
This algorithm computes the basin entropy Sb
of the basins of attraction. First, the input basins
is divided regularly into n-dimensional boxes of side ε
(along all dimensions). Then Sb
is simply the average of the Gibbs entropy computed over these boxes. The function returns the basin entropy Sb
as well as the boundary basin entropy Sbb
. The later is the average of the entropy only for boxes that contains at least two different basins, that is, for the boxes on the boundary.
The basin entropy is a measure of the uncertainty on the initial conditions of the basins. It is maximum at the value log(n_att)
being n_att
the number of attractors. In this case the boundary is intermingled: for a given initial condition we can find another initial condition that lead to another basin arbitriraly close. It provides also a simple criterion for fractality: if the boundary basin entropy Sbb
is above log(2)
then we have a fractal boundary. It doesn't mean that basins with values below cannot have a fractal boundary, for a more precise test see basins_fractal_test
. An important feature of the basin entropy is that it allows comparisons between different basins using the same box size ε
.
ChaosTools.basins_fractal_test
— Functionbasins_fractal_test(basins; ε = 20, Ntotal = 1000) -> test_res, Sbb
This is an automated test to decide if the boundary of the basins has fractal structures. The bottom line is to look at the basins with a magnifier of size ε
at random in basins
. If what we see in the magnifier looks like a smooth boundary (in average) we decide that the boundary is smooth. If it is not smooth we can say that at the scale ε
we have structures, i.e., it is fractal.
In practice the algorithm computes the boundary basin entropy Sbb
basin_entropy
for Ntotal
random balls of radius ε
. If the computed value is equal to theoretical value of a smooth boundary (taking into account statistical errors and biases) then we decide that we have a smooth boundary. Notice that the response test_res
may depend on the chosen ball radius ε
. For larger size, we may observe structures for smooth boundary and we obtain a different answer.
The output test_res
is a symbol describing the nature of the basin and the output Sbb
is the estimated value of the boundary basin entropy with the sampling method.
[Puy2021] Andreu Puy, Alvar Daza, Alexandre Wagemakers, Miguel A. F. Sanjuán. A test for fractal boundaries based on the basin entropy. Commun Nonlinear Sci Numer Simulat, 95, 105588, 2021.
Keyword arguments
ε = 20
: size of the ball for the test of basin. The result of the test may change with the size.Ntotal = 1000
: number of balls to test in the boundary for the computation ofSbb
ChaosTools.uncertainty_exponent
— Functionuncertainty_exponent(basins; kwargs...) -> ε, N_ε ,α
Estimate the uncertainty exponent[Grebogi1983] of the basins of attraction. This exponent is related to the final state sensitivity of the trajectories in the phase space. An exponent close to 1
means basins with smooth boundaries whereas an exponent close to 0
represent complety fractalized basins, also called riddled basins.
The output N_ε
is a vector with the number of the balls of radius ε
(in pixels) that contain at least two initial conditions that lead to different attractors. The ouput α
is the estimation of the uncertainty exponent using the box-counting dimension of the boundary by fitting a line in the log.(N_ε)
vs log.(1/ε)
curve. However it is recommended to analyze the curve directly for more accuracy.
Keyword arguments
range_ε = 2:maximum(size(basins))÷20
is the range of sizes of the ball to test (in pixels).
Description
A phase space with a fractal boundary may cause a uncertainty on the final state of the dynamical system for a given initial condition. A measure of this final state sensitivity is the uncertainty exponent. The algorithm probes the basin of attraction with balls of size ε
at random. If there are a least two initial conditions that lead to different attractors, a ball is tagged "uncertain". f_ε
is the fraction of "uncertain balls" to the total number of tries in the basin. In analogy to the fractal dimension, there is a scaling law between, f_ε ~ ε^α
. The number that characterizes this scaling is called the uncertainty exponent α
.
Notice that the uncertainty exponent and the box counting dimension of the boundary are related. We have Δ₀ = D - α
where Δ₀
is the box counting dimension computed with basins_fractal_dimension
and D
is the dimension of the phase space. The algorithm first estimates the box counting dimension of the boundary and returns the uncertainty exponent.
Tipping points
ChaosTools.basin_fractions
— Functionbasin_fractions(basins::Array) → fs::Dict
Calculate the fraction of the basins of attraction encoded in basins
. The elements of basins
are integers, enumerating the attractor that the entry of basins
converges to. Return a dictionary that maps attractor IDs to their relative fractions.
In[Menck2013] the authors use these fractions to quantify the stability of a basin of attraction, and specifically how it changes when a parameter is changed.
ChaosTools.tipping_probabilities
— Functiontipping_probabilities(basins_before, basins_after) → P
Return the tipping probabilities of the computed basins before and after a change in the system parameters (or time forcing), according to the definition of[Kaszás2019].
The input basins
are integer-valued arrays, where the integers enumerate the attractor, e.g. the output of basins_of_attraction
.
Description
Let $\mathcal{B}_i(p)$ denote the basin of attraction of attractor $A_i$ at parameter(s) $p$. Kaszás et al[Kaszás2019] define the tipping probability from $A_i$ to $A_j$, given a parameter change in the system of $p_- \to p_+$, as
\[P(A_i \to A_j | p_- \to p_+) = \frac{|\mathcal{B}_j(p_+) \cap \mathcal{B}_i(p_-)|}{|\mathcal{B}_i(p_-)|}\]
where $|\cdot|$ is simply the volume of the enclosed set. The value of $P(A_i \to A_j | p_- \to p_+)$ is P[i, j]
. The equation describes something quite simple: what is the overlap of the basin of attraction of $A_i$ at $p_-$ with that of the attractor $A_j$ at $p_+$. If basins_before, basins_after
contain values of -1
, corresponding to trajectories that diverge, this is considered as the last attractor of the system in P
.
Discrete system example
using DynamicalSystems, PyPlot
function newton_map(dz, z, p, n)
z1 = z[1] + im*z[2]
dz1 = newton_f(z1, p[1])/newton_df(z1, p[1])
z1 = z1 - dz1
dz[1]=real(z1)
dz[2]=imag(z1)
return
end
newton_f(x, p) = x^p - 1
newton_df(x, p)= p*x^(p-1)
# dummy Jacobian function due to https://github.com/JuliaDiff/ForwardDiff.jl/issues/520
function newton_map_J(J,z0, p, n) end
ds = DiscreteDynamicalSystem(newton_map, [0.1, 0.2], [3.0], newton_map_J)
xg = yg = range(-1.5,1.5,length=400)
basins, attractors = basins_of_attraction((xg, yg), ds; show_progress = false)
basins
400×400 Matrix{Int16}:
1 1 1 1 1 1 1 1 1 1 1 1 1 … 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 … 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
⋮ ⋮ ⋮ ⋱ ⋮ ⋮
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 … 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
attractors
Dict{Int16, Dataset{2, Float64}} with 3 entries:
2 => 2-dimensional Dataset{Float64} with 1 points
3 => 2-dimensional Dataset{Float64} with 1 points
1 => 2-dimensional Dataset{Float64} with 1 points
Now let's plot this as a heatmap
# Set up some code for plotting attractors
function scatter_attractors(attractors)
for k ∈ keys(attractors)
x, y = columns(attractors[k])
scatter(x, y; color = "C$(k-1)", edgecolor = "white")
end
end
LC = matplotlib.colors.ListedColormap
cmap = LC([matplotlib.colors.to_rgb("C$k") for k in 0:3])
vmin = 1; vmax = 4
fig = figure()
pcolormesh(xg, yg, basins'; cmap, vmin, vmax)
scatter_attractors(attractors)
fig.tight_layout(pad=0.3); fig
Stroboscopic map example
This example targets periodically driven 2D continuous dynamical systems, like the Duffing oscillator:
using DynamicalSystems, PyPlot
ω=1.0; f = 0.2
ds = Systems.duffing([0.1, 0.25]; ω, f, d = 0.15, β = -1)
2-dimensional continuous dynamical system
state: [0.1, 0.25]
rule f: duffing_rule
in-place? false
jacobian: duffing_jacob
parameters: [1.0, 0.2, 0.15, -1.0]
For stroboscopic maps, we strongly recommend using a higher precision integrator from OrdinaryDiffEq.jl.
using OrdinaryDiffEq
diffeq = (alg = Tsit5(), reltol = 1e-9)
xg = yg = range(-2.2,2.2,length=200)
basins, attractors = basins_of_attraction((xg, yg), ds; T=2π/ω, diffeq, show_progress = false)
basins
200×200 Matrix{Int16}:
1 2 1 2 1 1 1 1 1 2 2 1 2 … 1 1 2 2 2 2 2 1 2 2 2 1
1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 2 2 1 1 2 2 2 2 1 2
1 1 2 1 1 1 1 2 1 2 2 1 1 1 1 1 1 2 1 1 1 2 1 2 2
1 1 1 2 2 2 2 1 1 2 1 2 2 2 2 2 1 1 1 1 2 1 1 2 2
2 2 1 2 2 2 1 2 2 2 2 2 2 2 1 1 1 2 1 1 1 1 1 2 2
1 2 2 2 2 2 2 2 1 1 1 1 2 … 1 1 1 1 1 1 1 1 2 1 1 1
2 2 2 1 2 2 1 1 1 2 2 2 1 1 1 1 1 2 1 1 1 1 1 2 2
2 2 2 1 1 2 2 2 1 1 1 1 2 2 1 1 1 1 1 2 1 1 1 2 1
1 2 2 2 1 2 1 2 1 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 2
2 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1 2 2 2 1 1 1 1 2 1
⋮ ⋮ ⋮ ⋱ ⋮ ⋮
2 1 2 1 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2
2 2 2 1 2 2 2 2 2 2 1 1 1 1 1 1 1 1 2 1 1 1 2 2 1
1 2 1 2 2 2 2 1 1 2 2 1 2 1 1 1 2 1 2 1 2 2 1 2 2
1 2 1 2 2 2 1 2 1 1 2 1 2 1 1 1 2 1 2 1 2 2 1 2 1
2 2 1 2 2 2 2 1 1 1 2 1 1 … 1 2 1 2 2 2 2 1 2 1 2 2
1 1 2 2 2 2 1 1 2 2 2 2 1 1 2 2 1 2 2 1 2 1 2 2 2
2 1 2 2 1 2 2 2 2 2 2 2 2 2 2 1 2 1 2 2 2 2 1 2 1
2 2 2 2 2 2 1 1 2 1 1 2 1 2 2 2 1 1 2 2 2 1 1 1 2
2 1 2 1 2 2 2 2 1 2 2 1 2 1 1 2 2 1 1 1 2 1 2 1 2
And visualize the result as a heatmap, scattering the found attractors via scatter.
fig = figure()
pcolormesh(xg, yg, basins'; cmap, vmin, vmax)
scatter_attractors(attractors)
fig.tight_layout(pad=0.3); fig
2D basins of higher dimensional system
In this section we will calculate the basins of attraction of the four-dimensional magnetic pendulum. We know that the attractors of this system are all individual fixed points on the (x, y) plane so we will only compute the basins there. See the 3D basins section for something even more complex.
Computing the basins
ds = Systems.magnetic_pendulum(d=0.2, α=0.2, ω=0.8, N=3)
xg = yg = range(-4, 4, length=150)
@time basins, attractors = basins_of_attraction((xg, yg), ds; diffeq = (reltol = 1e-9,), show_progress = false)
attractors
Dict{Int16, Dataset{4, Float64}} with 3 entries:
2 => 4-dimensional Dataset{Float64} with 1 points
3 => 4-dimensional Dataset{Float64} with 1 points
1 => 4-dimensional Dataset{Float64} with 2 points
Alright, so far so good, we found 3 attractors (the 3 magnets). Let's visualize this beauty now
fig = figure()
pcolormesh(xg, yg, basins'; cmap, vmin, vmax)
scatter_attractors(attractors)
fig.tight_layout(pad=0.3); fig
Computing the uncertainty exponent
Let's now calculate the uncertainty_exponent
for this system as well. The calculation is straightforward:
ε, f_ε, α = uncertainty_exponent(basins)
fig = figure()
plot(log.(ε), log.(f_ε))
plot(log.(ε), log.(ε) .* α; label = "α = $(round(α; digits=3))")
legend()
fig.tight_layout(pad=0.3); fig
The actual uncertainty exponent is the slope of the curve.
Computing the tipping probabilities
We will compute the tipping probabilities using the magnetic pendulum's example as the "before" state. For the "after" state we will change the γ
parameter of the third magnet to be so small, it's basin of attraction will virtually disappear.
ds = Systems.magnetic_pendulum(d=0.2, α=0.2, ω=0.8, N=3, γs = [1.0, 1.0, 0.1])
basins_after, attractors_after = basins_of_attraction(
(xg, yg), ds; diffeq = (reltol = 1e-9,), show_progress = false
)
match_attractors!(basins, attractors, basins_after, attractors_after)
fig = figure()
pcolormesh(xg, yg, basins_after'; vmin, vmax, cmap)
scatter_attractors(attractors_after)
fig.tight_layout(pad=0.3); fig
P = tipping_probabilities(basins, basins_after)
3×2 Matrix{Float64}:
0.504877 0.495123
0.456224 0.543776
0.561929 0.438071
As you can see P
has size 3×2, as after the change only 2 attractors have been identified in the system (3 still exist but our state space discretization isn't accurate enough to find the 3rd because it has such a small basin). Also, the first row of P
is 50% probability to each other magnet, as it should be due to the system's symmetry.
3D basins
To showcase the true power of basins_of_attraction
we need to use a system whose attractors span higher-dimensional space. An example is
ds = Systems.thomas_cyclical(b = 0.1665)
3-dimensional continuous dynamical system
state: [1.0, 0.0, 0.0]
rule f: thomas_rule
in-place? false
jacobian: thomas_jacob
parameters: [0.1665]
which, for this parameter, contains 5 coexisting attractors. 3 of them are entangled periodic orbits that span across all three dimensions, and the remaining 2 are fixed points.
To compute the basins we define a three-dimensional grid and call on it basins_of_attraction
.
# This computation takes about an hour
xg = yg = zg = range(-6.0, 6.0; length = 251)
basins, attractors = basins_of_attraction((xg, yg, zg), ds)
attractors
Dict{Int16, Dataset{3, Float64}} with 5 entries:
5 => 3-dimensional Dataset{Float64} with 1 points
4 => 3-dimensional Dataset{Float64} with 379 points
6 => 3-dimensional Dataset{Float64} with 1 points
2 => 3-dimensional Dataset{Float64} with 538 points
3 => 3-dimensional Dataset{Float64} with 537 points
1 => 3-dimensional Dataset{Float64} with 1 points
The basins of attraction are very complicated. We can try to visualize them by animating the 2D slices at each z value, to obtain:
Then, we visualize the attractors to obtain:
In the animation above, the scattered points are the attractor values the function basins_of_attraction
found by itself. Of course, for the periodic orbits these points are incomplete. Once the function's logic understood we are on an attractor, it stops computing. However, we also simulated lines, by evolving initial conditions colored appropriately with the basins output.
The animation was produced with the code:
using GLMakie
fig = Figure()
display(fig)
ax = fig[1,1] = Axis3(fig; title = "found attractors")
cmap = cgrad(:dense, 6; categorical = true)
for i in keys(attractors)
tr = attractors[i]
markersize = length(attractors[i]) > 10 ? 2000 : 6000
marker = length(attractors[i]) > 10 ? :circle : :rect
scatter!(ax, columns(tr)...; markersize, marker, transparency = true, color = cmap[i])
j = findfirst(isequal(i), bsn)
x = xg[j[1]]
y = yg[j[2]]
z = zg[j[3]]
tr = trajectory(ds, 100, SVector(x,y,z); Ttr = 100)
lines!(ax, columns(tr)...; linewidth = 1.0, color = cmap[i])
end
a = range(0, 2π; length = 200) .+ π/4
record(fig, "cyclical_attractors.mp4", 1:length(a)) do i
ax.azimuth = a[i]
end
Poincaré map example
In the previous example we saw that this system has periodic attractors. In the Poincaré map these periodic attractors become points. We can use the functionality of basins_of_attraction
and poincaremap
to find basins of attraction on the Poincaré surface of section.
ds = Systems.thomas_cyclical(b = 0.1665)
xg = yg = range(-6.0, 6.0; length = 100)
pmap = poincaremap(ds, (3, 0.0), 1e6;
rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9
)
basins, attractors = basins_of_attraction((xg,yg), pmap; show_progress = false)
attractors
Dict{Int16, Dataset{3, Float64}} with 3 entries:
2 => 3-dimensional Dataset{Float64} with 1 points
3 => 3-dimensional Dataset{Float64} with 3 points
1 => 3-dimensional Dataset{Float64} with 1 points
attractors[1]
3-dimensional Dataset{Float64} with 1 points
1.83899 -4.15575 2.07651e-11
Looks good so far, but let's plot it as well:
fig = figure()
pcolormesh(xg, yg, basins'; cmap, vmin, vmax)
scatter_attractors(attractors)
fig.tight_layout(pad=0.3); fig
This aligns perfectly with the video we produced above.
- Yorke1997H. E. Nusse and J. A. Yorke, Dynamics: numerical explorations Ch. 7, Springer, New York, 1997
- Daza2016A. Daza, A. Wagemakers, B. Georgeot, D. Guéry-Odelin and M. A. F. Sanjuán, Basin entropy: a new tool to analyze uncertainty in dynamical systems, Sci. Rep., 6, 31416, 2016.
- Grebogi1983C. Grebogi, S. W. McDonald, E. Ott and J. A. Yorke, Final state sensitivity: An obstruction to predictability, Physics Letters A, 99, 9, 1983
- Menck2013Menck, Heitzig, Marwan & Kurths. How basin stability complements the linear stability paradigm. Nature Physics, 9(2), 89–92
- Kaszás2019Kaszás, Feudel & Tél. Tipping phenomena in typical dynamical systems subjected to parameter drift. Scientific Reports, 9(1)