Potentials and Limits to Basin Stability Estimation

Stability assessment methods for dynamical systems have recently been complemented by basin stability and derived measures, i.e. probabilistic statements whether systems remain in a basin of attraction given a distribution of perturbations. This requires numerical estimation via Monte-Carlo sampling and integration of differential equations. Here, we analyze the applicability of basin stability to systems with basin geometries challenging for this numerical method, having fractal basin boundaries and riddled or intermingled basins of attraction. We find that numerical basin stability estimation is still meaningful for fractal boundaries but reaches its limits for riddled basins with holes.

Going back to the path-breaking ideas of Aleksandr M. Lyapunov, dynamical systems are said to be stable if small variations of the initial conditions lead to small reactions of a system, i.e. small perturbations cannot substantially alter the system's state. This is commonly a statement about the asymptotic behaviour, allowing for large transient deviations if only the system eventually returns to the initial state. Multistable systems with several attractors add another subtlety to the problem: perturbations may lead to switching from one attractor to another, substantially altering asymptotic behaviour [1]. While infinitesimal perturbations on an attractor have local effects well-studied in the theory of asymptotic stability, finite (including large) perturbations can be critical by causing non-local effects like the transition to another attractor.
A direct stability method are Lyapunov functions [2][3][4], which decrease along trajectories and have local minima on attractors. Finding Lyapunov functions is, however, difficult in high-dimensional multi-stable systems, although there are recent approaches to determine them from radial basis functions (e.g. [5]).
Here, we pursue an alternative approach to consider non-local perturbations termed basin stability S B . The central idea [6,7] is to use a kind of volume of the basin of attraction to quantify the stability of attractors in multi-stable systems subject to a given distribution of perturbations. An advantage of this measure is that it can be efficiently estimated even in high-dimensional systems and as an intuitive interpretation as a probability to return to an attractor, but it relies on the correct identification of the asymptotic behaviour for a Monte Carlo sample of initial conditions. Basin stability and derived concepts have been successfully applied recently [8], e.g. for power grids [7,[9][10][11], chimera states [12], explosive synchronization [13] delayed dynamics [14] and resilience measures [15].
In numerical simulations, it can be difficult to correctly identify the asymptotic behaviour and determine the attractors. The basin of attraction can practically be defined as the set of all states that enter and stay in some trapping region [16]. Problems may arise if transients are long and chaotic or trajectories stay close to basin boundaries for long, so that numerical errors can move the simulated trajectory across a boundary into a wrong basin and make the simulation converge to a wrong attractor. Principally, three aspects contribute to the overall estimation error: the standard error due to sampling initial conditions, approximation errors in function evaluations or integration of differential equations, and rounding errors due to limited precision. While sampling and approximation errors are controlled by increasing sample size and order of approximating polynomials and by decreasing step size, rounding errors are typically hard to reduce, which is not a problem if they are much smaller.
Our study thus focuses on the critical case of systems where rounding errors cannot be neglected and may even dominate the overall error due to an intricate state space geometry highly sensitive to numerical imprecisions. We put basin stability estimation here to the test by applying it to systems with fractal basin boundaries and riddled or intermingled basins of attraction.
Consider a system of ordinary differential equationṡ that has more than one attractor in its state space X.
Here, we define an attractor as a minimal compact invariant set A ⊆ X whose basin of attraction has positive Lebesgue measure [17]. The basin of attraction of A is the set B(A) ⊆ X of all states from which the system converges to A. arXiv:1603.01844v1 [nlin.CD] 6 Mar 2016 Assume the system moves on an attractor A, yet at t = 0 a random and not necessarily small perturbation pushes the system to a state x(0) outside A. Assume that x(0) is drawn from a probability distribution with measure µ on X that encodes our knowledge about what relevant perturbations are how likely to occur. E.g., µ may be a uniform distribution on some bounded region R ⊃ A.
Will the system converge back to A after the perturbation? To address this, we study the probability mass of B, the probability that the system will return to A. The indicator function 1 B(A) (x) yields 1 if x ∈ B (A) and 0 otherwise. We use S B (A) to quantify just how stable the attractor A is against non-infinitesimal perturbations, and call it the basin stability of A [6,7]. The estimation of volume integrals such as Eq. 2 in high dimensions is a well-known problem, and we assume this is done by simple Monte-Carlo sampling [18,19]. If for each initial state x(0), one can numerically integrate the system x(t) with sufficient precision to decide to which attractor it converges or whether it diverges, the S B estimation procedure is thus: 1. Draw a sample of N > 0 independent initial states from the distribution µ.
2. For each, numerically integrate the system until it is clear whether and where it converges.
3. Count the number M of times the system has converged to A.
4. Use the estimateŜ B = M N . Since this is an N -times repeated Bernoulli experiment with success probability S B , the absolute standard error of the estimateŜ B due to sampling is S B (A) (1 − S B (A))/N , independently of the system's dimension. Thus, the procedure can be applied to highdimensional systems without necessarily increasing the sample size N , although it may take longer to assess convergence. This is of course since we are not interested in the basin of attraction's geometry but only in its volume w.r.t. the measure µ.
Note that when the relative std. err. ofŜ B is more relevant than the absolute std. err., smaller values of S B (A) require larger sample sizes, of the order N ∼ 1/S B (A), since for small S B (A), the rel. std. err. is ∼ 1/ N S B (A). However, even if S B (A) is not small, the geometries of the multiple basins of attraction may still make the estimation of S B difficult for another reason: For some initial conditions x(0), it may be quite difficult to decide where x(t) converges to, since the trajectory may start or come quite close to the boundary between the different basins, so that approximation and rounding errors (rather than sampling errors) in the integration may become relevant and may make the simulated trajectory hop across a basin border, leading to a wrong assessment of where x(t) actually converges to.
Particularly, this is probable if the basins have fractal boundaries, where the nature of the basin boundaries influences the predictability of a system's behaviour in the long run [16,20,21]. Imagine we randomly draw initial states from a box through which the boundary between the basins of two attractors runs. Suppose each initial state is specified up to a certain numerical error ε. Then for an initial state that is closer to the boundary than ε, it is uncertain to which of the two attractors the system will converge. Denote by f (ε) the fraction of initial states for which the outcome is uncertain, i.e. the uncertainty fraction [22,23]. If the boundary is a smooth curve, then f (ε) is just proportional to ε. However, if the boundary is fractal, then f (ε) ∝ ε α . If α < 1, the system exhibits final state sensitivity, i.e., to decrease the uncertainty one needs a substantial improvement in the knowledge of initial conditions. In a way, this power law scaling leads to an obstruction of predictability [22] very similar to the sensitive dependence on initial conditions in chaotic systems. It has been found to occur in Rayleigh-Bénard convection, e.g. numerically [24] and experimentally [25].
Predicting the long-term behaviour -the essence of estimating S B (A) -of systems with fractal basin boundaries may be hard [21] although generally, for most initial conditions, the final state sensitivity is much smaller than the unpredictability of the actual trajectory.
So let us first investigate how fractal basin boundaries impact the accuracy ofŜ B by studying the Wada pendulum [26,27]. Consider a damped, driven pendulum that is subject to a time-dependent forcing: For α = 0.1, K = 1 and X = 7/4, this system has several attractors [28]. The four dominant of them, all limit cycles with period 2π, are shown in Fig. 1a: The black and red attractors correspond to rotations of the pendulum, and the orange and yellow attractors are librations. Their respective basins of attraction at t = 0 are shown in Fig. 1b. Certain regions in this figure appear sprinkled with dots belonging to the different basins, i.e. the boundary between the basins is not easily discernible and remains so when zooming in (Fig. 1c). It is a fractal, resulting from the so-called Wada property of the basins.
Three (or more) subsets of a space are said to have the Wada property if any point on the boundary of one subset is also on the boundary of the two others [16,28]. For the pendulum, the black basin, the red basin and the union of the orange and yellow basins have the Wada property [16,28]. This means that starting within the rounding error ε of the boundary, a trajectory could in principle converge to any of the four attractors.
To verify this empirically, we write ε = 10 −p with p denoting precision, and discard all information after the p-th significant decimal digit in the floating point variables used in all individual operations of the numerical integration. We use 64 bit double precision to allow for a maximum of p = 16, while using untruncated 32 bit single precision would correspond to p ≈ 7. For different values of p, we integrate a fixed set of 50 initial states x(0), drawn uniformly at random from the rectangle R = [−π, π] × [−2, 4]. Fig. 2a, reveals that some initial states, particularly those indicated by arrows, indeed lead to different outcomes for different values of p. To investigate howŜ B depends on p, we let µ be the uniform distribution on R yielding a sample of N = 1, 000 random initial states which are integrated with different precisions p, leading to estimatesŜ B (p). As depicted in Fig. 2b, there seems to be no systematic influence of p onŜ B (p). Indeed, most of the individual values ofŜ B (p) are within one standard error of the most precise valueŜ B (16). This suggests that, in contrast to long-term prediction for individual initial states (cf. Fig. 2a),Ŝ B is robust under variation of p.
Another extreme case are attractors whose basins are not open as for most systems [17] but rather have an empty interior. The complement of such a riddled basin intersects every disk in a set of positive measure [29][30][31][32]. This means that all points in its basin of attraction have pieces of another attractor basin arbitrarily closely nearby [30].
Physical systems exhibiting riddled basins are the damped, periodically-driven particle moving in a special potential [33] or coupled time-delayed systems [34][35][36]. There are also experimental observations for laser-cooled ions in a Paul trap [37] indicating a riddled phase space structure.
In the following, we investigate the impact of riddled basins of attraction onŜ B using a conceptual example [29,38], i.e. the following quadratic map on the complex plane:

Re(z)
Im(z) The black area corresponds to initial conditions for which the dynamics diverge. Below are zoom-ins of two regions, (b)and (c). The locations of the attractors (line segments, see [29]) are highlighted by red/blue/purple bars (not in scale).
This map has three different attractors on the complex plane which are shown in Fig. 3; for simplicity they are referred to as the red/blue/purple attractors with their respective basin of attraction in the following. Interestingly, the three basins of attraction are not just riddled, they are intermingled. A basin of attraction is called intermingled if any open set which intersects one basin in a set of positive measure also intersects each of the other basins in a set of positive measure [39,40].
The fact that there is a positive probability to end up in a different attractor around each initial condition inside a riddled/intermingled basin of attraction renders these systems effectively non-deterministic [33]. As in the case of Wada boundaries, slight variations of initial conditions or numerical imprecisions will affect any forecast of the system's long-term behaviour.
Again, we investigate the effect of limited numerical precision on the significance ofŜ B . In Fig. 4a Fig. 3a. We observe a large variation ofŜ B of up to 50% compared to the most precise estimationŜ B (16) and no systematic dependence on p.
In Fig. 3c we zoomed into the neighbourhood of the red attractor, where the share of the corresponding red basin is increasing in proximity of the attractor. In particular, the measure of this basin of attraction, restricted to an -neighbourhood of the attractor, approaches unit probability for → 0 [29]. This apparent behaviour provides an explanation for Fig. 4c where we determined S B (p) for Fig. 3c. In contrast to our previous observation, the fluctuations ofŜ B (p) almost stay within one standard error and the estimation appears to be more robust. For reference, Fig. 4b depictsŜ B (p) for Fig. 3b not containing any (part of) an attractor. On the one hand, the variation ofŜ B (p) exceeds one standard error, up to about 20% compared toŜ B (16), such that our estimation is more sensitive to numerical imprecisions than in Fig. 4b; on the other hand the variations are smaller than in our first experiment.
In conclusion, we applied the Monte-Carlo estimation procedure of basin stability in two cases, i.e. basins with fractal boundaries and riddled/intermingled basins of attraction. In the former case, we find that while the asymptotic properties of individual trajectories still cannot be determined robustly, the converse is true for the basin stability estimation. It remains an open question for future research, how exactly (in a quantitative sense) the numerical estimation uncertainty might be derived from the actual basin geometry. In the latter case, however, we find that the results can vary drastically with the chosen precision. The effect of rounding errors is comparable or even larger than the standard error of the sampling. Only if the sample region R is chosen in some sense "close enough" to the actual attractor of interest, the foliated structure of the surrounding basins allows for a meaningful numerical estimation.
What are practical implications for the application of basin stability? In general, it is sufficient if the rounding error of an estimation is smaller than its sampling error to get a significant result. However, any numerical procedure is subject to a finite numerical precision and we have to assume that in practice it will not be high enough to reach this goal in dynamical systems with intricate basin geometries. If there is no prior knowledge available, a good starting point is to actually visualize the interesting part of the phase space to get a first idea of the appearance of, e.g., fractal sets. If any are detected, it is necessary to use the highest available numerical precision p h to getŜ B (p h ), potentially avoiding artifacts respectively insignificant estimations. We suggest to repeat the S B estimation at a lower numerical precision p l and take the differenceê p = |Ŝ B (p h ) −Ŝ B (p l )| as a straight-forward (rough) estimator of the variability of S B (p) with p and, by way of extrapolation, as a rough estimate of the remaining standard error ofŜ B (p h ) as an estimate of S B due to finite numerical precision. To assess the influence of rounding errors onŜ B then compareê p with the standard error ofŜ B (p h ) as an estimate of S B (p h ) due to sampling, which can be estimated aŝ s p = Ŝ B (p h )(1 −Ŝ B (p h ))/N . Ifê p <ŝ p , rounding has no significant effect on the estimation quality. For instance, this could be implemented by comparing the results at double and single precision computations.