Basins of Attraction

In DynamicalSystems.jl we provide performant methods for estimating basins of attraction of various attractors. The performance of these methods comes from a constraint on a 2D plane, as you will see below.

ChaosTools.basins_map2D — Function

basins_map2D(xg, yg, integ; kwargs...) → basins, attractors

Compute an estimate of the basins of attraction on a two-dimensional plane using a map of the plane onto itself according to the method of Nusse & Yorke^[Yorke1997]. The dynamical system should be a discrete two dimensional system such as:

Discrete 2D map.
2D poincaré map.
A 2D stroboscopic map.

For a higher-dimensional dynamical system, use basins_general.

integ is an istance of an integrator, not a DynamicalSystem. This includes the output of poincaremap. See documentation online for examples for all cases! xg, yg are 1-dimensional ranges that define the grid of the initial conditions to test. The output basins is a matrix on the grid (xg, yg), see below for details. The output attractors is a dictionary whose keys correspond to the attractor number and the values contains the points of the attractors found on the map. Notice that for some attractors this list may be incomplete.

Keyword Arguments

T : Period of the stroboscopic map, in case integ is an integrator of a 2D continuous dynamical system with periodic time forcing.
Ncheck : A parameter that sets the number of consecutives hits of an attractor before deciding the basin of the initial condition.

Description

This method identifies the attractors and their basins of attraction on the grid without prior knowledge about the system. At the end of a successfull computation the function returns a matrix coding the basins of attraction and a dictionary with all attractors found.

basins has the following organization:

The atractors points are even numbers in the matrix. For example, 2 and 4 refer to distinct attractors.
The basins are coded with odd numbers, 2n+1 corresponding the attractor 2n.
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 (even numbers).
The value associated to a key is a Dataset with the guessed location of the attractor on the plane.

The method starts by picking the first available initial condition on the plane not yet numbered. The dynamical system is then iterated until one of the following conditions happens:

The trajectory hits a known attractor already numbered: the initial condition is numbered with corresponding odd number.
The trajectory diverges or hits an attractor outside the defined grid: the initial condition is set to -1
The trajectory hits a known basin 10 times in a row: the initial condition belongs to that basin and is numbered accordingly.
The trajectory hits 60 times in a row an unnumbered cell: it is considered an attractor and is labelled with a even number.

Regarding performace, this method is at worst as fast as tracking the attractors. In most cases there is a signicative improvement in speed.

ChaosTools.basins_general — Function

basins_general(xg, yg, ds::DynamicalSystem; kwargs...) -> basin, attractors

Compute an estimate of the basins of attraction of a higher-dimensional dynamical system ds on a projection of the system dynamics on a two-dimensional plane.

Like basins_map2D, xg, yg are ranges defining the grid of initial conditions on the plane.

Refer to basins_map2D for detailed information on the computation and the structure of basin and attractors.

Keyword Arguments

dt = 1: Approxiamte time step of the integrator. It is recommended to use values ≥ 1.
idxs = 1:2: This vector selects the two variables of the system that will define the "plane" the dynamics will be projected into.
complete_state = zeros(D-2): This argument allows setting the remaining variables of the dynamical system state on each planar initial condition x, y. It can be either a vector of length D-2, or a function f(x, y) that returns a vector of length D-2.
Ncheck = 10: A parameter that sets the number of consecutives hits of an attractor before deciding the basin of the initial condition.
diffeq...: Keyword arguments propagated to integrator.

Stroboscopic map example

First define a dynamical system on the plane, for example with a stroboscopic map or Poincaré section. For example we can set up an dynamical system with a stroboscopic map using a periodically driven 2D continuous dynamical system, like the Duffing oscillator:

using DynamicalSystems
ω=1.; F = 0.2
ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1)
integ = integrator(ds; reltol=1e-8)

t: 0.0
u: 2-element SVector{2, Float64} with indices SOneTo(2):
 0.1
 0.25

Now we define the grid of ICs that we want to analyze and launch the procedure:

xg = range(-2.2,2.2,length=200)
yg = range(-2.2,2.2,length=200)
basin, attractors = basins_map2D(xg, yg, integ; T=2π/ω)
basin

200×200 Matrix{Int16}:
 3  5  3  5  3  3  5  3  3  5  5  3  5  …  3  5  5  5  5  5  3  3  5  5  5  5
 3  3  3  3  3  5  3  5  3  3  3  3  5     3  5  5  3  3  3  5  5  5  5  5  5
 3  3  5  3  3  3  3  5  3  5  5  3  3     3  3  5  5  3  3  5  3  5  5  3  5
 3  3  3  5  5  5  3  5  3  5  3  5  5     5  5  3  3  3  5  3  5  3  3  5  5
 5  5  3  3  5  5  3  5  5  5  5  5  5     3  3  3  3  5  3  3  3  3  3  5  3
 3  5  5  5  5  5  5  5  3  3  5  3  5  …  3  5  3  3  3  3  3  3  5  3  3  3
 5  5  5  3  5  5  3  5  5  5  5  5  3     3  3  3  5  5  5  3  3  3  3  3  5
 5  5  5  5  5  5  5  5  3  5  3  3  3     5  3  3  3  3  3  5  5  3  3  5  3
 5  5  5  5  3  5  5  5  5  3  5  3  3     3  5  5  3  3  3  3  5  5  5  3  5
 5  5  3  5  3  3  5  3  3  3  3  3  3     3  5  3  5  5  3  5  3  3  3  5  3
 ⋮              ⋮              ⋮        ⋱        ⋮              ⋮           
 3  3  5  3  3  5  5  5  3  3  3  3  3     3  3  3  3  3  3  3  3  3  5  5  5
 3  5  5  3  5  5  5  5  5  5  3  3  3     3  3  3  3  3  5  3  3  3  5  5  5
 5  3  3  3  5  5  5  3  3  3  5  3  5     3  3  3  5  3  5  3  5  5  3  5  3
 3  5  5  5  5  5  3  5  5  3  3  3  5     5  3  3  5  3  5  3  5  5  3  3  3
 5  5  3  3  5  5  5  3  3  3  5  3  3  …  3  5  3  5  5  3  5  5  3  3  5  5
 5  3  5  5  5  3  3  3  5  5  5  3  3     3  3  5  5  5  5  3  5  5  3  5  3
 5  3  5  5  3  5  5  5  5  5  3  5  5     5  3  3  5  5  5  5  5  3  3  5  5
 5  5  5  5  5  5  5  3  5  3  3  5  5     5  5  5  3  3  3  3  5  3  3  3  5
 5  3  5  3  5  5  5  5  3  3  5  3  5     5  3  5  3  3  3  3  5  3  5  3  5

using PyPlot
fig = figure()
pcolormesh(xg, yg, basin'; cmap = "Accent")
fig.tight_layout(pad=0.3); fig

Poincaré map example

A Poincaré map of a 3D continuous system is a 2D discrete system and can be directly passed into basins_map2D.

ds = Systems.rikitake(μ = 0.47, α = 1.0)
plane = (3, 0.0)
pmap = poincaremap(ds, (3, 0.), Tmax=1e6;
    idxs = 1:2, rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9
)

Iterator of the Poincaré map
 rule f:      rikitake_rule
 hyperplane:  (3, 0.0)
 selection:   [1, 2]

Once the Poincaré map has been created, we simply call basins_map2D

xg=range(-6.,6.,length=200)
yg=range(-6.,6.,length=200)
basin, attractors = basins_map2D(xg, yg, pmap)

fig = figure()
pcolormesh(xg, yg, basin'; cmap = "Accent")
fig.tight_layout(pad=0.3); fig

Discrete system example

The process to compute the attraction basins of a discrete 2D dynamical system is trivial, as one passes its integrator directly into basins_map2D

function newton_map(dz, z, p, n)
    z1 = z[1] + im*z[2]
    dz1 = f(z1, p[1])/df(z1, p[1])
    z1 = z1 - dz1
    dz[1]=real(z1)
    dz[2]=imag(z1)
    return
end
f(x, p) = x^p - 1
df(x, p)= p*x^(p-1)

# dummy Jacobian function to keep the initializator quiet
function newton_map_J(J,z0, p, n)
   return
end

ds = DiscreteDynamicalSystem(newton_map,[0.1, 0.2], [3.0], newton_map_J)
integ  = integrator(ds)

xg=range(-1.5,1.5,length=200)
yg=range(-1.5,1.5,length=200)

basin, attractors  = basins_map2D(xg, yg, integ)
fig = figure()
pcolormesh(xg, yg, basin'; cmap = "Accent")
fig.tight_layout(pad=0.3); fig

Basins in Higher Dimensions

When it is not so simple to define a 2D stroboscopic map or Poincaré map, which is the case in continuous dynamical systems of higher dimensionality, you can always try the general method basins_general. It is slower and may requires some tuning. The algorithm looks for attractors on a 2D grid. The initial conditions are set on this grid and all others variables are set to zero by default.

Example

ds = Systems.magnetic_pendulum(γ=1, d=0.2, α=0.2, ω=0.8, N=3)
xg=range(-4, 4, length=150)
yg=range(-4, 4, length=150)
@time basin, attractors = basins_general(xg, yg, ds; idxs = 1:2, reltol = 1e-9)
attractors

Dict{Int16, Dataset{2, Float64}} with 3 entries:
  4 => 2-dimensional Dataset{Float64} with 2 points
  6 => 2-dimensional Dataset{Float64} with 2 points
  2 => 2-dimensional Dataset{Float64} with 3 points

Alright, so far so good, we found 3 attractors (the 3 points of the magnetic pendulum). Let's visualize this beauty now

fig = figure()
pcolormesh(xg, yg, basin'; cmap = "Accent")
fig

Yorke1997H. E. Nusse and J. A. Yorke, Dynamics: numerical explorations Ch. 7, Springer, New York, 1997