Framework for global stability analysis of dynamical systems

Multistable dynamical systems exhibit two or more co-existing stable states, formally called attractors. Multistability is ubiquitous in nature and in mathematical models,^1,2 with examples ranging from power grids,^3–5 the climate,^6,7 ecosystems like the Amazon rain forest,^8,9 the brain and neuronal circuits therein,^10–12 or metabolic systems.^13–15 Some attractors of these systems can be desirable, such as synchronized oscillations in power grids, crucial for their proper functioning.¹⁶ However, they can also be undesirable, as, for example, the extinction of a certain species in ecological models or the collapse of circulation in climate models.¹⁷ In a multistable system, the attractor at which the system ends up depends on the initial conditions (ICs), but perturbations of the state may enforce switching between attractors, a phenomenon called “tipping.”^1,9,18 Alterations in the parameters of a dynamical system can trigger tipping. Hence, it becomes important to evaluate how “resilient” attractors are to perturbations, either to parameters or to the system’s variables. This is a crucial problem of practical importance in several areas of research.^9,19

A traditional solution to this problem is the continuation-based bifurcation analysis (CBA). It identifies fixed points and (under some requirements) limit cycles and describes their linear (also called exponential) stability dependence on a system parameter via the eigenvalues of the Jacobian.²⁰ One major downside is that, by definition, it cannot be applied to chaotic attractors.²⁰ In the majority of cases, this analysis must be done numerically via one of several software, e.g., AUTO,²¹ MATCONT,²² CoCo,²³ or Bifurcationkit.jl.²⁴ Continuation tools can be combined with other numerical techniques to enrich the exploration of the parameter space as shown in Ref. 25. The information provided by these frameworks is useful but incomplete: rigorously speaking, local stability only conveys information about the system’s response to infinitesimally small perturbations. It cannot yield insight into the response to finite perturbations in the state space, which are predominant in practice.

For such responses, it is necessary to study the global stability²⁶ of the system’s attractors, which involves the nonlinear dynamics over the full state space of the system.²⁷ A proxy for global stability of an attractor is the portion of all possible initial conditions ending up at this attractor, i.e., the fraction of the state space that is in the basin of said attractor. When the state space is infinite, the concept of the state space fraction becomes a pragmatic one: we need to define a finite-volume box of physically plausible initial conditions for the system under study, and we are concerned about the fractions of these plausible initial conditions. In this analysis, attractors with larger basins fractions are globally more stable because stronger perturbations are typically needed to switch the system to another attractor.²⁶ Frequently, this measure is also a much better indicator of the loss of stability as a system parameter is varied when compared to the local stability analysis of the system (see, e.g., Ref. 26 or Chap. 12 of Ref. 20).

Analyzing global stability as a function of a parameter demands extensive effort from researchers, as it requires the creation of algorithms that can find system attractors and their global stability, “continue” them across a parameter and also perform the expensive numerical simulations required for such algorithms. In the literature, the only framework so far that can aid this analysis is the featurizing and grouping approach, proposed first by Gelbrecht et al.²⁸ as MCBB (Monte Carlo Basin Bifurcation Analysis) and then later very similarly by Stender and Hoffmann²⁹ as bSTAB (basin stability). The method integrates randomly sampled initial conditions (ICs) of a dynamical system for a preset time span. The trajectories of these ICs are then mapped onto features, numbers describing the trajectories, such as the mean or standard deviation of some of the system variables. All the feature vectors are clustered using the DBSCAN algorithm³⁰ so that ideally each cluster corresponds to an attractor of the system. The fractions of ICs in each cluster approximate the basin fractions and, hence, the global stability of each attractor. More details on the method are given in Sec. IV.

This method works well in a variety of circumstances and can also be applied across a parameter range. However, it comes with significant downsides. One is that it is not clear a priori which features should be chosen to correctly separate the attractors into clusters, requiring a lot of trial and error. Another downside is that the method cannot guarantee that the clusters of features really correspond to unique attractors and that it is not mixing two or more attractors together. An alternative method for finding attractors and their basins of attraction is the recurrence-based approach recently proposed.³¹ The method locates attractors by finding recurrences in the system’s state space, assuming the Poincaré recurrence theorem holds for the system attractors. The input to this method is a state space box, and its tessellation, defining a grid to search for recurrences. Hence, the method will only find attractors within the given box, although the box can initially be arbitrarily large. We describe the method in more detail in Sec. IV B and provide a comparison between the two techniques in Sec. III B. The main advantage of the recurrences method is that it locates the actual system attractors and only requires as an input a state space box that may contain the attractors. So far, however, it has not been clear how to “continue” attractors across a parameter range with this method.

In this work, in Sec. III A, we present a novel global stability analysis and continuation algorithm that utilizes the recurrences-based method for finding attractors.³¹ In Sec. III C, we apply it to exemplary models of climate, ecosystem dynamics, and more. As detailed in Sec. III, this novel continuation algorithm is the most accurate in finding the actual attractors of a dynamical system, the most transparent in matching attractors across parameter values, and requires the least amount of guesswork from the researcher. We believe that this novel continuation of global stability, much like global stability analysis itself,²⁶ is a crucial new tool for the analysis of dynamical systems. In some cases, it supersedes CBA, and in others, it can be complemented by CBA, as we discuss in Sec. III A.

This continuation method is part of a novel automated framework that performs global stability analysis and continuation, which we present in Sec. III B. Our framework significantly advances existing methodology, including the featurizing methods, thereby including all upsides of current literature while addressing most downsides (Secs. IV C and IV D describe the improvements in detail). Its design is based on modular components that can be configured or extended independently. This allows researchers to simply compose the methodology that is best suited to their problem and then let an automated algorithm execute the process. The framework is accompanied by an extensively tested, well documented, and highly optimized open-source software implementation, part of the DynamicalSystems.jl³² general purpose library for nonlinear dynamics (see Sec. IV for code example and documentation).

A major contribution of our framework is the novel algorithm for global stability analysis and continuation that we name recurrences-based attractor find-and-match continuation, RAFM for short. This algorithm can be applied to any system whose attractors satisfy the Poincaré recurrence theorem. As illustrated in Fig. 1(a), it works as follows.

1. the starting point of the algorithm. Attractors and their basins, or basins fractions, are already known at a parameter $p= p 1$ using the recurrences-based algorithm.³¹ 

2. new initial conditions are seeded from the existing attractors. Then, we set the system parameter to $p= p 2$⁠.

3. evolve the seeded initial conditions according to the dynamic rule of the system. The seeds are evolved until they converge to an attractor using the recurrences-based method (the grid reflects the tessellation of the state space that is decided by the user; the finer the grid, the more accurate the results³¹). The main performance bottleneck of the recurrences-based method is finding the attractors. Once found, convergence of other initial conditions is generally much faster.³¹ To address this in the continuation, we use the observation that, unless a bifurcation is occurring, attractor size, shape, and position, typically depend smoothly on $p$⁠. Hence, the seeded initial conditions at each new parameter will most likely converge the fastest to the new attractors.

4. with the main bottleneck of the algorithm (finding the attractors) being taken care of, now compute the basins fractions by formally applying the recurrence-based algorithm³¹ to randomly sampled initial conditions. It can be arbitrarily configured how to randomly sample initial conditions, but typically the samples are drawn uniformly from a state space box. Importantly, the algorithm may still find new attractors during this step (here the “dark green” one) that did not exist before.

5. match attractors in current parameter $p 2$ to those in parameter $p 1$⁠. Matching is arguably the most sophisticated part of the algorithm. Attractors are matched by their “distance” in state space, with distance any arbitrary metric on the space of state space sets, see Sec. III D for more details. In this illustration, the “red” attractor is matched to the “purple” one of Step 0, while neither the “yellow” or “dark green” attractors match to the previous “teal” one, because their state space distance is beyond a pre-defined threshold.

The end result is the (matched) system attractors, and their basins fractions (or full basins if computationally feasible), as functions of the parameter. The attractors and basins are labeled with the positive integers (enumerating the different attractors), and the basins always sum to 1. Note that because RAFM works on a parameter-by-parameter basis, it can be used to perform continuation across any number of parameters, not just one (we present one here as the simplest conceptual example).

To perform global stability analysis, several tasks need to be taken in sequence [see Fig. 1(b) for an overview]. We have abstracted and generalized the tasks to allow researchers different possibilities of how to achieve them.

The first task is the creation of a dynamical system for the global stability analysis. For our framework, this is achieved for free simply by making the implementation part of the DynamicalSystems.jl library³² (see Sec. IV E).

The second task is the creation of a mechanism to find attractors and map initial conditions to them. Possibilities for this mechanism are: (1) featurizing and grouping initial conditions into attractors (as discussed in the introduction), (2) finding attractors using the recurrences algorithm,³¹ or (3) mapping initial conditions to previously known attractors by proximity: once the evolution of an initial condition comes close enough to a pre-determined attractor, the initial condition is mapped to that attractor. Note that in (1) or (2), attractors are found via random sampling in the state space. For the typically used uniform sampling, the probability to find an attractor with basin fraction $f$ after $n$ samples is $1− ( 1 − f ) n$⁠. However, users may use system-specific knowledge to adjust the sampling into something that may find attractors with better scaling with respect to sample number $n$⁠.

Mechanism (1) is paired with instructions on how to group features. Currently, the possibilities are: (a) clustering features into groups using DBSCAN (as done in MCBB or bSTAB), (b) grouping features by histograms in the feature space, so that features that end up in the same histogram bin belong to the same group (novel grouping approach), or (c) mapping features to their nearest feature in feature space, from a set of pre-defined features (as done in bSTAB). Having more grouping possibilities than only clustering can be useful, as discussed in Sec. III D.

At this point, one can analyze global stability at a given parameter and also analyze the basin boundaries of attractors for fractal properties. From here, the third task is to “continue” the attractors and basins across a parameter range. Our framework currently has two continuation algorithms; however, due to an extendable design, more may be added in the future. The first continuation algorithm is what has been employed so far by the MCBB or bSTAB algorithms, but with significantly increased accuracy (see MatMeth), and extended to allow any kind of instructions for how to group features. We will call this “Featurize and Group Across Parameter” (FGAP). In this approach, trajectories from the dynamical system are generated by sampling random ICs across all parameter values of interest. All these trajectories are mapped to features, and all these feature vectors are then grouped using one of currently three grouping instructions (clustering, histogram, nearest feature). Each group is representing an attractor. The grouped ICs are then re-distributed into the parameter slices they came from, providing the fractions of each group at each parameter value. The second continuation algorithm is the RAFM algorithm that we described in Sec. II A. These two approaches are compared in detail in Sec. III B.

In Fig. 2, we apply RAFM on some exemplary systems. We stress that we could characterize the different attractors in accurate detail because RAFM finds the actual system attractors, not some incomplete representation of them (i.e., features used in the featurizing-and-grouping approach).

To illustrate a concrete application of the method, we discuss how the example of the Lorenz 84 dynamical system in Fig. 2(a) was generated. Along with the following step by step “tutorial,” we provide the computer code for this example in Listing 1. First, a continuous time dynamical system object ds (whose equations are provided in Sec. IV F) is created, and we choose the ODE integrator according to the system specifics (here we use the Verner 9th order solver³⁸). Next, we choose the state space box, and its tessellation, that will be used both for the recurrences-algorithm³¹ as well as for sampling random initial conditions when estimating the basin fractions. For this example, the box ranges from $−$3 to $+3$ divided into 600 points along each dimension.

With the grid and the dynamical system, we construct a mapper object that maps initial conditions to attractors using the recurrences method.³¹ To decide the meta-parameters of this algorithm, we consult our previous publication specifically about the algorithm.³¹ For this example, the presence of a chaotic attractor involves longer recurrence times, so we set high values for the parameters mx_chk_fnd_att and mx_chk_loc_att to improve the precision of the characterization of the attractors. Sometimes, transient trajectories may spend some time outside the defined grid. In this case, we can either provide a larger state space box or set the parameter mx_chk_lost to a higher value. Very small time steps may lead to the identification of false attractors if the orbit stays too long in a grid cell. To prevent these artifacts, we specify that during time integration, a non-adaptive method with a fixed time step of Dt should be used (note that the value of Dt will depend on the system’s internal timescale.

We chose 100 values of the parameter $G$ in the range $[1.34;1.37]$ and construct a continuation object rsc, which instructs how to perform the continuation (in this case, this corresponds to using the RAFM algorithm). During the construction of rsc, we could also specify how to match attractors (which is something we only discuss in Sec. II D and in this example, we used the default matching). We also specify how to sample random initial conditions via the sampler object, which in this case is randomly drawn initial conditions within the limits of the provided state space box. The continuation output is obtained after the call of the continuation function. This workflow is identical for the other examples used, in addition to the obvious changes of the dynamical system, state space box, parameter range, and possible tuning of meta-parameters. We also stress that the best place to learn how to choose optimal meta-parameters of all these algorithms is the documentation of the associated software, which is constantly updated with best practices.³⁹ 

Traditional CBA has a rigid “matching” procedure: it always matches the next point found along a “continuation curve” to the previous point. This is often correct for infinitesimal perturbations of fixed points but becomes problematic for global stability analysis, which attempts to find all attractors in a state space box and then continue them. In this case, matching attractors from one parameter to the next becomes a crucial part of the algorithm. For instance, the analysis presented in Fig. 2 is only coherent because of the powerful matching procedure implemented in our framework. Without it, the colors would alternate arbitrarily at each parameter value.

In the featurize-and-group algorithm, matching and grouping are the same process. In RAFM, matching is rather sophisticated and operates on a parameter-by-parameter basis. Each time the parameter is incremented and the new attractors are found, a matching sub-routine is launched. The distance between attractors before and after the parameter change is estimated, with “distance” being any symmetric positive-definite function defined on the space of state space sets. By default, the Euclidean distance of the attractor centroids is used because it is extremely fast and in the majority of cases, it works very well. A more rigorous metric is the Hausdorff distance,⁴⁰ which is also provided out of the box. Additionally, “distance” is not limited to state space distances. It could be the distance across dynamic invariants. For example, one can track the Lyapunov spectrum or the fractal dimension²⁰ of each attractor and define a distance in terms of their absolute difference. This is easily possible in our code implementation because it is part of DynamicalSystems.jl, which offers algorithms for computing Lyapunov spectra or fractal dimensions out of the box.

After the distance is computed between all new-old attractor pairs, the new attractor labels are matched to the previous attractor labels that have the smallest distance to them, prioritizing pairs with the smallest distance. The matching respects uniqueness, so that once an attractor from the previous parameter has been matched, this attractor is removed from the matching pool and cannot be matched to an additional new attractor. Additionally, a distance threshold value can be provided, so that old-new pairs of attractors whose distance is larger than this threshold are guaranteed to get assigned different IDs. Note that in principle finding the attractors and matching them are two completely independent processes. If after the continuation process is finished the user decides that the chosen matching procedure was unsuitable, they can launch a “re-matching” algorithm with different matching “distance” function, without having to re-do any computations for finding the attractors or their fractions (i.e., matching only renames attractor labels but leaves the attractors themselves untouched).

The last thing to highlight in this section is the desirable post-processing of grouping similar enough attractors. This happens automatically if one uses the featurize-and-group continuation method. However, taking as an example Fig. 2(e), the RAFM method finds countless individual attractors. For the researcher, the individual attractors may be useful for careful analysis, but it is sometimes desirable to group similar enough attractors. In our framework, it is possible to use exactly the same grouping infrastructure utilized by the featurizing-and-grouping continuation, but now applied to the outcome of RAFM as a post-processing step. In Fig. 3, we highlight examples that utilize the powerful matching and/or grouping components offered by our framework.

In Table I and Fig. 4, we provide a careful comparison between CBA and RAFM. A direct comparison of the two approaches is difficult, since their operational basis, main output, and even the way stability is quantified, are fundamentally different. On one hand, CBA finds and continues the curves of individual fixed points or limit cycles across the joint state-parameter space on the basis of Newton’s method. The stability is quantified in terms of local stability (also called exponential stability via the Jacobian eigenvalues). On the other hand, RAFM first finds attractors at all parameter values using the original system equations and then matches appropriately similar attractors across different parameters, giving the illusion of continuing them individually. Additionally, the curves of stable fixed points in the joint parameter space (Fig. 4) are only a small part of the information our framework provides. Important, provided information is the basin fractions and how they change, which is completely absent in CBA.

Based on this comparison, we argue that the RAFM algorithm and the global stability framework we provide is an essential tool for stability analysis of dynamical systems. We further believe that in some application scenarios, RAFM will supersede CBA, especially given the difference in required user expertise, required interventions, and steepness of the learning curve that CBA has over RAFM. In other scenarios, we envision that RAFM can be used as the default analysis method, providing the majority of information, and CBA then becomes a more in-depth analysis of fixed points, limit cycles, bifurcations, and even stable/unstable manifolds, if such analysis is required.

Our framework provides two radically different methods for finding attractors: the recurrence-based and the featurize-and-group approach. Generally speaking, the recurrence-based method should be preferred when possible, because of its accuracy (finding the actual attractors), and the possibility for follow-up analysis of found attractors. However, the featurize-and-group method should be preferred when the recurrence-based method fails, because, e.g., the provided state space tessellation is ill-defined, or because computational demands exceed what is available, such as in higher-dimensional systems with chaotic attractors. In Table II, we provide a comprehensive comparison between the two methods.

We group together two similar methods that have been recently proposed in the literature for finding attractors. One is called Monte Carlo basin bifurcation analysis (MCBB),²⁸ the other is basin stability analysis (bSTAB).²⁹ Both methods work by identifying attractors as clusters of user-defined-features of trajectories. Their first step is to integrate $N$ randomly chosen initial conditions inside a certain box in state space. The integration is done for some time $T$⁠, after a transient $T tr$⁠. Both $T$ need to be sufficiently large so that the trajectories correspond to the systems’ long-term behavior and also $T tr$ needs to be large to avoid the transient regime. Each trajectory $x →(t)$ is then transformed into a vector of $K$ features $F →$⁠, specified by some featurizer function $f →$ such that $f → ( x →(t))= F →$ that has to be defined by the user. Each vector of features $F →$ describes a point in the $K$-dimensional space of features $f 1× f 2×⋯× f K$⁠. The key idea is that features belonging to the same attractor cluster together in state space, so that each attractor forms a distinct cluster (a cloud of points) in feature space. The final step in the method is to therefore cluster the features. The clustering algorithm chosen for this is the Density-based Spatial Clustering of Applications with Noise (DBSCAN).³⁰ It first classifies two points as neighbors if their distance if smaller than a radius $ε$⁠. Then, it clusters together points with many neighbors (equal or more than a parameter minPts), and leave as outliers points with too few neighbors (less than minPts). The radius $ε$ is a crucial parameter for the algorithm and often needs fine tuning for proper clustering. The methods by Refs. 28 and 29 use two different ways to identify a value for $ε$⁠. Authors in Ref. 28 use a method that looks at the ordered distance of the $k$ nearest neighbors to each point in the dataset and finds the first knee (high derivative point).^30,45 Authors in Ref. 29 iteratively search for the $ε$ that maximizes a criterion of clustering quality. To calculate this criterion, they evaluate the silhouette of each cluster, which measures how similar each point is to the cluster it currently belongs to, compared to the other clusters, and ranges from $−1$ (worst matching) to $+1$ (ideal matching). This leads to one silhouette value per feature; the authors then take the minimum value as the representative for the clustering quality for each radius. The chosen radius is, thus, the value that leads to the highest minimum silhouette. In both methods, the clusters found by DBSCAN are considered then as attractors.

The inputs to this method are a dynamical system, a state space box that may contain the attractors (although initially it may be arbitrarily large) and a tessellation of the given box into cells. An initial condition of the system is evolved step-by-step with time step $Δt$⁠. At each step, the location of the trajectory in state space is mapped to its cell and that the cell is labeled as visited. If the dynamical system has attractors and they satisfy the Poincaré recurrence theorem (Chap. 9 of 20), the trajectory is guaranteed to revisit cells it has visited before. Once a pre-decided number of recurrences $n f$ have been accumulated consecutively (i.e., enough previously visited cells are visited again), the method claims to have found an attractor. It then proceeds to locate the attractor accurately, by collecting a pre-decided number $n l$ of points on the attractor [Fig. 1 of Datseris and Wagemakers³¹ and panel (3) of Fig. 1]. A finite state machine formulation keeps track of coexisting attractors, so that each attractor is unique. It also keeps track of divergence to infinity by counting steps $n d$ outside the box, ensures algorithm termination by setting a total $n m$ of max amount of $Δt$ iterations, and makes convergence faster by utilizing information already encoded in the grid: if the trajectory visits consecutively a relatively small number $n r$ of cells already labeled as an attractor, convergence is already eagerly decided, i.e., converging to an already found attractor is much faster than finding that attractor for the first time. Hence, $Δt, n f, n l, n d, n m, and n r$ are the meta-parameters of the algorithm and have sensible default values that work in most cases. More information on the method can be found in Ref. 31. Notice that the recurrence method is different at a fundamental level from the Global Analysis of Invariant Objects (GAIO)⁴⁶ and other cell mapping techniques.⁴⁷ We expand more on this in Sec. IV. Also, note that the method is not perfect; it may identify two attractors when only one exists due to, e.g., choosing too low convergence criteria $n l, n f$⁠, or due to a commensurate period of attractor and integrator time step. However, once again the importance of finding the “actual” attractors becomes apparent: further analysis by, e.g., plotting the attractors, immediately highlights such a failure and how to deal with it, and we also provide several tips in the documentation of our method in the code implementation.³⁹ 

A large drawback of the recurrences method was that it scaled poorly with the dimension $D$ of the dynamical system. If an $ε$-sized tessellation of the state space is chosen, then memory allocated scaled as $1/ ε D$⁠. We now use sparse arrays to store accessed grid locations. This changes the memory scaling to $1/ ε Δ$⁠, with $Δ$ the capacity dimension²⁰ of the attractor, which is typically much smaller than $D$ (and is only $1$ for limit cycles).

First, we have changed the criterion used for finding the optimal radius in the clustering method. We have found the knee method consistently more unreliable than the iterative search. We have also found that the mean, instead of the minimum, silhouette value as the measure of clustering quality leads to better clustering. For instance, this leads to correct clustering in the Lorenz86 system, whereas the minimum value criterion did not. Furthermore, our method searches for the optimal radius with a bisection method, instead of the linear method used by authors in Ref. 29. This significantly speeds up the code. Another simple modification we introduced is to rescale the features in each dimension (⁠ $f 1, f 2,…, f K$⁠) into the same interval, for instance, $[0,1]$⁠. We noticed that the clustering method performs poorly if the features span different ranges of values, and this simple modification proved to be a very powerful solution. Third, we allow the integration of all initial conditions to be done in parallel, using several computer cores, which speeds up the solution. Last, grouping in our framework can also happen based on a histogram in feature space.

The code implementation of our framework is part of the DynamicalSystems.jl library³² as the Attractors.jl package.³⁹ The code is open-source code for the Julia language, has been developed following best practices in scientific code,⁴⁸ is tested extensively, and is accompanied by a high quality documentation. An example code snippet is shown in Listing 1.

In addition to the quality of the implementation, three more features of the code are noteworthy. First, that it is part of DynamicalSystems.jl instead of an isolated piece of code. This integration makes the simplicity and high-levelness of Listing 1 possible and makes the input for the code easy to set up. Moreover, the direct output of the code can be used with the rest of the library to further analyze attractors in terms of, e.g., Lyapunov exponents or fractal dimensions. Indeed, in the provided code example Listing 1, we compute the Lyapunov spectra of all found attractors, across all parameter values, in only two additional lines of code. Second, utilizing the Julia language’s multiple dispatch system,⁴⁹ the code is extendable. It establishes one interface for how to map initial conditions to attractors and one for how to group features, both of which can be extended, and yet readily be usable by the rest of the library such as the continuation methods. Third, a lot of attention has been put into user experience, by establishing a short learning curve via a minimal user interface, by carefully considering how to provide the output in an intuitive format, as well as providing easy-to-use plotting functions that utilize the code output. More overview and information on the code or its design can be found in its online documentation or source code.³⁹ 

In this section, we state the dynamic rules (equations of motion), parameter values, and state space boxes used for all systems in the main text.

Global parameters are $L x=1.75π$⁠, $L z=1.2π$⁠, $α=2π/ L x$⁠; and $β=π/2$⁠, $γ=2π/ L z$⁠.

Listing 1: Julia code snippet showcasing the usage of the DynamicalSystems.jl implementation of our framework. The code produces panel (a) of Fig. 2. The main output of the code are two vectors, containing the basins fractions and attractors at each parameter value, respectively. The fractions and attractors are formulated as dictionaries, mapping attractor labels (the integers) to basin fractions and sets of points on the attractor, respectively. At its end, the code snippet computes the Lyapunov spectra of all found attractors by using the first point on each attractor as initial condition for the computation of the Lyapunov spectrum.

The benchmarks presented in Fig. 5 provide a comparison between techniques for finding attractors for a discrete and continuous dynamical system.

The numerical approximations of the attractors and their basins of attraction can be achieved with other well known numerical tools. The cell mapping and the GAIO algorithms rely on a subdivision of the state space. In its simplest form,⁴⁷ the cell mapping technique transforms the dynamical system into a discrete mapping. Each cell of the state space is mapped to another cell following the dynamics during a fixed time $T$⁠. The new representation of the dynamical system is a directed graph where each node represents an initial condition on the tessellated phase space and a single edge starts from this node to another one in the graph. Once the full mapping has been obtained, efficient search algorithms inspired from graph theory approximate the basins and the attractors. The computational effort is centered around the construction of the mapping and depends directly on the discretization of the state space. The computational complexity explodes with the system dimension and, hence, limits its practical use to low-dimensional systems.

The Global Analysis of Invariant Objects (GAIO) algorithm⁵⁰ is also based on the discretization on a region but only a subset of cells is subdivided following the dynamics of the system. An iterative process allows us to approximate accurately the attractor manifold and also other invariant sets embedded in the state space. The complexity only depends on the dimension of the manifold not on the dimension of the state space, which is a considerable gain over the cell mapping. Since the method is focused on global attracting sets, it is not usable for multistable systems, which are the main focus of our work.

These two techniques are hard to apply in estimating the basins fractions across a parameter value. The cell mapping technique requires a full description of the state space for each parameter, making the random sampling of the phase space impossible. The GAIO technique can continue a global attractor but will fail if two or more exist at any parameter value.⁴⁶ 

We should mention another algorithm using interval arithmetic for continuation of the dynamics.⁵¹ The continuation is performed in a database space using Conley indices and Morse decomposition. Although very effective for low-dimensional maps, it lacks the flexibility and ease of use of our proposed method.

The authors would like to thank Ulrike Feudel, Peter Ashwin, Ulrich Parlitz, and Harry Dankowicz for helpful discussions. G.D. was supported by the Royal Society International Newton Fellowship. K.L.R. was supported by the German Academic Exchange Service (DAAD). A.W. was supported by the Spanish State Research Agency (AEI) and the European Regional Development Fund (ERDF) under Project No. PID2019-105554GB-I00 (MCIN/AEI/10.13039/501100011033).

The authors have no conflicts to disclose.

G.D. coordinated the project. All authors contributed in compiling the results, creating the code, writing, and revising the manuscript.

George Datseris: Conceptualization (lead); Formal analysis (equal); Software (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). Kalel Luiz Rossi: Conceptualization (equal); Formal analysis (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Alexandre Wagemakers: Conceptualization (equal); Formal analysis (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal).

The code we used to create the figures of this article is fully reproducible and available online, Ref. 52.