bSTAB: an open-source software for computing the basin stability of multi-stable dynamical systems

The pervasiveness of multi-stability in nonlinear dynamical systems calls for novel concepts of stability and a consistent quantification of long-term behavior. The basin stability is a global stability metric that builds on estimating the basin of attraction volumes by Monte Carlo sampling. The computation involves extensive numerical time integrations, attractor characterization, and clustering of trajectories. We introduce bSTAB, an open-source software project that aims at enabling researchers to efficiently compute the basin stability of their dynamical systems with minimal efforts and in a highly automated manner. The source code, available at https://github.com/TUHH-DYN/bSTAB/, is available for the programming language Matlab featuring parallelization for distributed computing, automated sensitivity and bifurcation analysis as well as plotting functionalities. We illustrate the versatility and robustness of bSTAB for four canonical dynamical systems from several fields of nonlinear dynamics featuring periodic and chaotic dynamics, complicated multi-stability, non-smooth dynamics, and fractal basins of attraction. The bSTAB projects aims at fostering interdisciplinary scientific collaborations in the field of nonlinear dynamics and is driven by the interaction and contribution of the community to the software package.


Introduction
Stability analysis plays a central role in nonlinear dynamics research, whether for structural dynamics, fluid flows, chemical reactions, or others [1]. Linear and local stability concepts are standard methods for the characterization of point attractors and (quasi-) periodic orbits. However, these stability concepts are limited to small perturbations when assessing the longterm system behavior and stability against perturbations. In a multi-stability situation, i.e. when multiple stable solutions co-exist, global metrics are required for studying the stability of a state against any type of perturbation, whether small or non-small. The basin stability proposed by Menck et al. [2,3] is a global stability concept based on the estimated volume of the basins of attraction in the system's phase space, which can overcome some limitations of classical local stability concepts. Practically, the basin stability indicates the probability of the system to arrive on one of the multiple solutions when perturbing the current state of the system. This stability concept is of relevance especially for designing practical application systems which are always subject to uncertainty of initial conditions or instantaneous perturbations. The concept of basin stability has led to a better understanding of some complex dynamical systems. For example, power grids were studied with a special interest in finding those network nodes that are critical for the network stability [4]. Albeit being a rather simplistic theoretical concept, the practical computation of the basin stability involves extensive Monte Carlo simulations and several technical challenges for researchers who want to study the basin stability of their dynamical systems.
As the major contribution of the work at hand, bSTAB aims at equipping researchers with a readily implemented tool to compute the basin stability of a dynamical system with minimal effort and high efficiency. The bSTAB toolbox aims at making the basin stability analysis easier to perform, hence helping this global stability metric make its way to a stateof-the-art analysis tool for nonlinear dynamical systems. bSTAB is designed to be model-agnostic, i.e. the user is required to have only basic knowledge about the system at hand, typical ODE formulations can be re-used, and only very limited coding skills are needed. Hence, this toolbox can be integrated into existing analysis routines without much effort. The toolbox is licensed under the GNU General Public License v3.0 and is freely available through https://github.com/ TUHH-DYN/bSTAB/, which aims at providing easy access and high visibility. The code is under active development, and the community can add functionalities or report issues within the git framework. This project shall foster communication and collaborations across scientific disciplines, and at the same time grow the functionalities of bSTAB further.
Classically, bifurcation diagrams are computed for nonlinear dynamical systems to study the character of different solutions, and their qualitative change under system parameter variations. Subcritical bifurcation points, also denoted as tipping points are of special interest, as they form the connection of unstable and stable solutions in bi-stable system configurations. Bifurcation diagrams are helpful for studying questions related to which dynamics are possible? Even though the computation of unstable solutions, representing the boundaries of the basins of attraction in phase space, is possible, bifurcation analysis can only barely answer questions related to which dynamics are the most probable?. Owing to the possibly intermingled or fractal geometries of the basins of attraction, a simple basin size estimation can become very challenging even for few-degree-of-freedom systems [22]. Lyapunov functions [30,31] can assess the stability against large perturbations, but they are difficult to find, especially for high-dimensional systems. Many open-source software projects exist for computing bifurcation diagram, such as MatCont [32] and NLvib [33], but yet no basin stability toolbox exists. This work aims at addressing this white spot for scientific communities working on nonlinear dynamical systems.
The basin stability concept relies on measuring the volumes of the basins of attraction B in the D− dimensional state space. The basin stability value S B (A) of the attractor A is given by [2,3,7] where κ (y) is an indicator function stating whether a trajectory converged to the attractor A, and thus whether y belongs to the basin B (A). ρ is a density distribution of states that the system may be perturbed to with R ρ (y) dy = 1. Estimates for the basin stability are obtained through Monte Carlo sampling from ρ, resulting in the absolute standard error of the estimate for the N -times repeated Bernoulli experiment, i.e. the square root of the variance. Hence, the choice of N is independent of the state space dimension, which is especially helpful for high-dimensional systems. The relative error of the estimate becomes relevant for small basin stability values, potentially requiring to increase N into the range N ∼ 1/S B (A) [3]. Each state is numerically integrated, and the fraction of trajectories approaching the respective attractor is calculated to yield the basin stability estimate. In the following, we use S B to indicate the estimate of the actual basin stability value for reasons of readability. This work is structured as follows: the basin stability concept is re-visited briefly in Sect. 2 with a focus on computational challenges for practitioners. Next, the structure of bSTAB is introduced in Sect. 3 along with typical workflows for studying the global stability of multi-stable dynamical systems. Sect. 4 illustrates the use of the toolbox for four canonical time-continuous dynamical systems.

Methods
This section re-visits the concept of basin stability from the more technical viewpoint, i.e. discussing how Eqs. (1), and (2) can be translated into efficient programs with flexible choice of parameters. First, the general nomenclature is introduced.
Nonlinear dynamical systemṡ are considered in their D-dimensional phase space. Model parameters are collected in p = p 1 , . . . , p j , and may be varied for parameter studies such as bifurcation diagrams. The time evolution of the states is referred to as the trajectory y (t). Typically, time marching solvers, such as Runge-Kutta schemes, are employed to solve the initial value problem that is associated with finding the trajectory y (t) for some initial condition y 0 = y (t = 0) and the time t. Given that the trajectory converges to either a fixed point, to a (quasi)-periodic orbit, or to a chaotic attractor, the system settles to some steady-state solution, denoted as y, after an initial transient phase. The basin of attraction B (ȳ) for the steady-state solutionȳ denotes the subset of the phase space from which all initial conditions converge towardsȳ as t → ∞. If the system has only one stable solution, the corresponding basin of attraction fills the complete phase space. However, in a multi-stability scenario, the phase space can be subdivided into k different basins of attraction if the system has k stable steady-state solutionsȳ 1,...,k . The boundaries of the basins, i.e. the separatrices, are given by the unstable solutions of the system. The concept of basin stability is a probabilistic approach to estimating the volumes V 1,...,k of the basins of attractions by sampling a finite number N of states from the distribution ρ and classifying the corresponding steady-state trajectories into one class of the k stable solutions. The fraction of each basin's volume V (B i ) in the phase space covered by ρ is denoted as the basin stability value S B (ȳ i ). Methodically, four main steps are involved in the basin stability computation, which are introduced in the following paragraphs. Particularly, a consistent formalization is required for the implementation of a model-agnostic computing toolbox that can run highly automated analysis with minimal user interaction.

Selection of the region of interest and Monte Carlo sampling
Even though the state space has an infinite size, typically only a small portion of this space represents a set of admissible states and potential perturbations for a dynamical system. Therefore, the distribution ρ is typically confined to a hypercube in R D , which is referred to as region of interest Q (y) ⊂ R D . Q should cover the set of possible initial conditions, but also the set of states that would result from perturbing a steadystate solutionȳ + . Hence, the selection of Q is specific to the dynamical system at hand, and relies on the expertise of the domain expert studying this particular system. Naturally, the selected region of interest crucially affects the final absolute basin stability values. One may interpret the a-priori selection of the region of interest as a weak point of the basin stability concept. However, one may on the other hand also argue that for classical local stability metrics the actual size of small, i.e. admissible, perturbations is unknown, too.
In bSTAB, the user specifies the minimum and maximum state values per state space dimension to set Q. N initial conditions are drawn from the region of interest Q following the distribution ρ. For the distribution, a bounded uniform distribution at random is the classical choice [2] that guarantees basin stability values to be proportional to the true basins' volumes for N → ∞. However, in certain situations different sampling strategies may be considered. For example, if a distribution of initial conditions or admissible perturbations is known from previous studies or experiments, they can be supplied through a custom definition of ρ. In the latter case, the basin stability values will no longer be an estimate for the basins' volumes. Instead, the basin stability values will indicate probabilities for observing a specific steady-state solutionȳ under the distribution ρ. bSTAB comes with the following D-dimensional distributions ρ readily implemented: bounded uniform random distribution, bounded multivariate Gaussian distribution, and multivariate uniformly spaced grid sampling. Furthermore, bSTAB allows the user to sample only from a subset of the D dimensions, i.e. exclude some dimensions from the sampling, and thereby create a hypercube Q that has less dimensions than the state space D. Figure 1 illustrates a schematic for the specification of the region of interest, as well as different distributions ρ.

Numerical time integration
For each state y 0 sampled from Q, the steady-state is obtained by solving the initial value problem by numerical time integration. Particularly, the steady-stateȳ is reached at time t * , i.e. after this time the dynamics are qualitatively constant with respect to either statistical properties, spectral properties, geometric [34] and dynamical invariant measures [35][36][37], or other descriptors of the nonlinear dynamics [38]. The integration time span must be chosen large enough such that also slow or even chaotic transients [39] have died out. The numerical time marching scheme, error tolerance settings, the time step length, and other parameters may affect the time-asymptotic behavior, especially for chaotic systems with exponential divergence of nearby starting trajectories. As investigated by Schulz et al. [40], numerical approximation and rounding errors may infer with the basin stability estimation when trajectories start, or come, close to a basin boundary, particularly for fractal boundaries. bSTAB provides functionalities to perform sensitivity analysis studies of the basin stability values against parameters of the numerical integration scheme. Classical Runge-Kutta schemes from Matlab's ode suite are provided: ode45 for fast computations, ode23t for vanishing numerical damping, and ode113 for numerically stiff systems.

Classification of asymptotic trajectories
The third step involves turning the indicator function κ B (y) in Eq. (2) into a data-driven feature extraction and clustering task. Steady-state trajectories y (t > t * ) are characterized by descriptive features X in order to assign a label L of membership to one of the k solutions. For this purpose, the feature extraction function φ converts dynamic data (trajectories y(t)) into static data (descriptive features X), and can be considered the most crucial ingredient to a fully automated basin stability computation. The formulation of the operator φ directly depends on the complexity and the diversity of the dynamical behavior of the system: if only a fixedpoint solutionȳ 1 needs to be distinguished from a periodic steady-stateȳ 2 , a trivial amplitude threshold can serve as discriminative feature X . In situations involving multiple solutions with higher degrees of similarity, the feature extraction function φ may involve advanced signal processing techniques (linear time series analysis, nonlinear time series analysis [36], recurrence plot quantification analysis [37], generalized feature extraction [38]) in order to come up with a set of m characteristic features X 1,...,m that can be used for the classification of different trajectories. The user must specify a feature extraction function in bSTAB taking trajectory data y(t) as input and returning a feature vector X.
Depending on the a-priori knowledge about the multi-stable dynamics, the labeling takes the form of a supervised or unsupervised clustering procedure C.
Assuming complete knowledge about all k solutions (supervised setting), template featuresX 1,...,k can be derived per a-priori known template solutionȳ 1,...,k . Then, a label to a new trajectory y (t) is selected by Fig. 1 Schematic for the definition of the region of interest Q as a subset of the state space (here: R 2 ) for a system that is k = 3 multi-stable in (a). The panel on the right hand side shows N = 500 samples drawn from different distributions ρ: (b) uniform random distribution, (c) multivariate Gaussian distribution and (d) uniformly spaced multivariate grid within the comparing the new feature vector against the template vectors and assigning the class label L according to the minimal distance A k-nearest neighbor classification algorithm (kNN) with a single nearest neighbor is provided with bSTAB, creating Voronoi partitions in the m-dimensional feature space upon which the labels are assigned, see Fig. 2a. If only incomplete knowledge about possible solutions is available, the unsupervised setting aims at finding clusters in the feature space. Particularly, not all solutions are required to be known, hence leaving the number k unknown, and making template feature definitions impossible. The unsupervised clustering algorithm DBSCAN (Density-Based Spatial Clustering of Applications with Noise) [41] is employed to derive optimally separated clusters for the N feature vectors X in the m−dimensional feature space, thereby finding k by the number of identified clusters. The unsupervised setting enables users of the toolbox to blindly put a less-known dynamical system at test, and let bSTAB find out which and how many multi-stable solutions can arise from states in Q.

Computation of the basin stability estimates
Having classified the N trajectories into the k classes of steady-statesȳ 1,...,k , the number of class members N 1,...,k is available. Here, N i denotes the number of initial conditions y 0 that converge to the ith steadystate solutionȳ i with i = 1, . . . , k and k i=1 N i = N . The resulting estimate for the basin stability value, i.e. the estimate for each basin's volume's V (B i ) share in Q, is given by Naturally, the basin stability values add up to unity if all N trajectories were able to be classified belonging to one of the k solutions.

Implementation aspects of bSTAB
This section illustrates the workflow and the readily available analysis routines in bSTAB. Details on the actual implementation, further options and annotated test cases are shown in the user manual (the latest version is accessible through https://github.com/ TUHH-DYN/bSTAB/) for bSTAB. Generally, the toolbox covers three major analysis modes, which cover most of the requirements for the analysis of nonlinear dynamical systems across various scientific disciplines: mode 1: computation of the basin stability values S B at a fixed set of model parameters p using the routine compute_bs. mode 2: computation of the sensitivity of the basin stability values against workflow hyperparameters. Hyperparameters are denoted as settings that are independent of the dynamical system at hand, e.g. number of sampling points and numerical tolerances. Sensitivity studies are performed using the routine compute_bs_ap, which allows to specify a wide range of adaptive parameters. mode 3: computation of the basin stability values along the variation of a model parameter p in analogy to conventional bifurcation analysis using the routine compute_bs_ap.
Mode 1 implements the core methodologies of the toolbox as described in the previous paragraphs, while modes 2 and 3 implement automatic looping through the parameter variations and calling mode 1 iteratively. The pseudo-code displayed in Algorithm 1 shows the implementation of mode 1 in bSTAB. Obtain trajectories y i (t) through time integration from y 0,i : 5: Extract feature vector from steady-state X i = φ (y i ((t > t * )) 6: end for 7: if supervised mode then 8: compute template featuresX 1,...,k 9: end if 10: Cluster feature vectors and assign labels C : To illustrate a typical basin stability computation using mode 1, the necessary commands are briefly discussed. More details can be found in the user manual. Each case study (first line) must be given a unique name in bSTAB for creating respective sub-directories which store the computation results. init_bSTAB initializes the bSTAB toolbox and setup_system is the case definition function. The case definition setup_system serves as anchor for running cases, hence all model parameters and hyperparameters need to be specified here by the user.  Figure 3 displays parallel speedups for one of the case studies shown in Sect. 4, namely the bi-stable chaotic Lorenz system at σ = 0.12, r = 0, b = −0.6. The analysis was performed using N = 5 · 10 5 samples and an extralong time integration span of 10 4 . Even though absolute computing times are affected by multiple factors related to the hardware and software resources, qualitative trends can be read from the analysis. Compute times are compared to the duration of a not-parallelized run on one worker, which in absolute numbers took 0.2938 seconds per sample. The study shows that the speedup linearly scales in a 2 : 1 ratio, such that using 8 workers reduces the computation time by the factor 4 compared to using a single worker. The scaling property makes the toolbox especially valuable for running larger case studies on high-performance machines making use of massive parallelization.

Canonical case studies using bSTAB
A range of four canonical dynamical systems is studied to illustrate the novel insights provided by basin stability analysis and to highlight the main functionalities in bSTAB. From the vast number of multi-stable dynamical systems [18], we select (1) a bi-stable (fixed point and periodic solution) damped driven pendulum, (2) a multi-stable (five periodic solutions) configuration of the Duffing oscillator, (3) a bi-stable (chaotic and chaotic) variant of the Lorenz system, and (4) a bistable (fixed point and periodic solution) non-smooth self-excited friction oscillator. The Appendix 5 lists all parameter settings and user-defined options for those case studies, such that readers can reproduce and verify our results. The different modes of the toolbox, namely basin stability computation, sensitivity analysis and model parameter bifurcation studies, are presented.

Bi-stable damped driven pendulum
Following the original work by Menck et al. [2], first the basin stability of the damped driven pendulum (length l, damping coefficient α, angular acceleration T , gravitational acceleration g, K = g/l) is investigateḋ The global stability picture obtained through this analysis hence adds insight into the dynamical system: Both solutions are stable (as indicated by linear stability analysis), but the periodic response is more than five times more likely to occur than the fixed point solution at the given parameter combination and the chosen distribution ρ, which is the additional piece of information that complements the local stability statements. Even though the system dynamics may seem trivial, manually or analytically estimating the basins' volumes is in fact non-trivial, see Fig. 4b. The basin of attraction of the fixed point is composed of many thin islands in the selected region of interest. The model-agnostic character of bSTAB helps the user to estimate those volumes in a completely automated fashion and enables effortless variations of the region of interest. All bSTAB parameters are given in Appendix A.
In a second study, the torque T is varied to study the local stability and the basin stability of both solutions along the parameter variation, hence making use of mode 3 in bSTAB. The bi-stable regime begins at T = 0.13, after which the stable fixed point FP and the limit cycle LC co-exist, see Fig. 5. The basin stability values are computed for a grid of torque values, and the linear stability analysis for the fixed point is computed along. The Lyapunov exponent λ associated with the fixed pointȳ 1 is constant throughout the complete parameter regime up to T = 0.99999, thereby indicating the stability of the fixed point, but not a transition of the qualitative dynamics picture of the system. On the other hand, the basin stability depicts how the global stability role of the fixed point solution vanishes for increasing driving torques, until eventually the fixed point solution disappears through a fold bifurcation at T = K = 1. The authors in [2] argue that the basin stability analysis can yield a better early warning sign for the critical tipping point of this system than local stability concept would allow.
For all pendulum studies the feature extraction task in Eq. 6 was required to formulate features X that would separate a fixed point solution from a limit cycle solution in the steady-state regimeỹ = y (t > t * ) irrespective of absolute amplitude values. For this purpose, the deviation of the maximum rotational velocity y 2 (t) about the average rotational velocity was considered using a one-hot encoding with respect to a threshold value to generate the feature vector X. Full knowl- edge about all solutions of the system is available apriori, hence the supervised classification using the kNN method was chosen for this system.

Multi-stable Duffing oscillator: five co-existing limit cycles
The Duffing oscillatoṙ is studied in a particular parameter setting (δ = 0.08, k 3 = 1, F = 0.2) that allows this small system to exhibit complicated multi-stability [42]: five periodic solutions co-exist for this configuration. There are two period−1 limit cyclesȳ 1 ,ȳ 2 , two period−2 limit cycles y 3 ,ȳ 4 and one period−3 limit cycleȳ 5 . As all of them are linearly stable, the shape and the volumes of the basins of attraction are of particular interest, but at the same time very challenging to derive analytically. The share of each basin in the state space Q will dictate the probability The largeamplitude period−1 cycleȳ 2 has by far the largest basin stability, and the two period−2 cycles have an almost vanishing probability of less than 3% in the chosen region of interest.ȳ 5 is the second most-likely andȳ 1 is the third most-likely solution of the system in the given configuration. Those results may be used to tune the dynamical system or adapt control strategies for arriving at the desired dynamical behavior. Again, this piece of information is not possible to extract from conventional linear stability analysis, and would not be easy to obtain from bifurcation diagrams or other conventional methods to measure the basins of attraction. Figure 6d depicts the feature space X which is used by the kNN clustering algorithm to assign class labels to the feature vectors X i . For deriving features from the trajectories, the maximum displacement X 1 = max (ỹ 1 (t)) and the displacement's standard deviation X 2 = std (ỹ 1 (t)) in the steady-state regimeỹ (t) = y (900 ≤ t ≤ 1000) are considered. Even though these features are rather simplistic, they turn out to be uniquely separable in the feature space with only minimal inner-class variance. In this case, the supervised clustering mode was used, hence initial conditions for each of the five solutions were supplied to the toolbox, which allowed to compute template feature vectors per class (indicated by cross markers in Fig. 6d). Fig. 6 a displacement trajectories y 1 (t) for the five stable limit cycle solutions y 1,...,5 of the Duffing oscillator in Eq. (11) and b the corresponding state space representation of the steady-state solutions. c N = 5000 samples from the region of interest labeled by their long-term behavior and d the feature space spanned by the maximum displacement (X 1 ) and the displacement standard deviation (X 2 ) being used for clustering the trajectories The bSTAB toolbox allows the user to arrive at these results without additional effort besides setting up the case and running the computation using mode 1 of the software package. All graphical output shown in Fig. 6 is created automatically, and all settings and results are stored in a local folder sub-directory for the current study. Changing the region of interest, the number of samples or other parameters only requires changing the respective values in the case definition file and rerunning the analysis. This design of the toolbox aims at enabling the user to quickly perform several analysis runs, and re-visit the results later on without having to consider version control or extensive logging of case specifications.

Multi-stable Lorenz system: periodic and chaotic states
Li and Sprott [43] study the famous Lorenz system [44] x = σ (y − x) in a parameter setting which exhibits particularly interesting multi-stable dynamics: for r = 0, b = −0.6 and σ = 0.12 the system has two stable co-existing chaotic attractors. Even if the resulting model may no longer represent the dynamics of convection rolls or chemical reactions [45], it is an interesting case study also for the basin stability analysis. Figure 7 shows the bifurcation behavior of the system under a variation of the model parameter σ . As this parameter is decreased, two coexisting limit cycles undergo a period-doubling route to chaos until they eventually are destroyed in a boundary crisis after σ = 0.12. Chaos is born at σ ≈ 0.1225. The two co-existing attractors are symmetric and the chaotic sets were entitled a broken butterfly by Li and Sprott [43]. Even though the attractors approach each other very closely at x = 0, they remain separated throughout the range of σ . Instead of deriving the basin stability at a fixed parameter set, the basin stability of both attractors is studied along a variation of σ using mode 3 (variation of a model parameter) in bSTAB. The specification of all parameters is given in Appendix C.  White areas in the state space relate to states leading into unbounded orbits computed for N = 5 · 10 4 samples from a uniform random distribution each. The basin stability analogue to the bifurcation diagram is depicted in Fig. 8. The analysis shows that more than 80% of the trajectories lead to unbounded orbits. The basin stability of both attractors are very similar owing to the symmetrical structure about the z axis of the basins of attraction shown in the lower panel of Fig. 8. For the four parameter values investigated in detail, the basins of attraction exhibit an intricate and fractal shape [43], which per se is an interesting observation: For σ > 0.01225 regular dynamics arise in a system that has fractal basin boundaries. Moreover, while the dynamics change rather strongly along σ in terms of amplitude and their qualitative characteristics, the basins of attraction remain rather unchanged. As a direct result of this observation, the basin stability values of both attractors are S B ≈ 0.09 throughout the parameter range studied here.
Fractal basin boundaries [46] raise questions with regards to the applicability of the Monte Carlo samplingbased stability concept. Schulz et al. [40] investigate the limitations of the basin stability concept for systems with "'fractal, riddled, or intermingled basin boundaries"'. The authors find that "the numerical basin stability estimation is still meaningful for fractal boundaries"' if large N are chosen. Utilizing bSTAB for hyperparameter studies (mode 2) using the compute_bs_ap routine, Fig. 9a displays the basin stability values of the two coexisting chaotic attractors at σ = 0.12 as a function of the number of sampling point ranging up to N = 20, 000. The basin stability values for this set of parameters saturates quickly at N ≈ 10 3 with only small fluctuations for smaller values of N . The relative tolerance setting of the numerical integration scheme was studied in a second hyperparameter study. The standard fourth-order Runge-Kutta ode45 integrator is used with an absolute error toler-(a) (b) Fig. 9 The basin stability values for the Lorenz system in Eq. 12 at σ = 0.12, r = 0, b = −0.6 as a function of hyperparameters: a variation of the number of sampled N and b variation of the relative error tolerance of the numerical integration scheme ance of 10 −6 . The relative tolerance value is changed in the range of 10 −3 (default setting) to 10 −8 using N = 10 4 points and a bounded uniform random distribution ρ on Q. Figure 9b shows that only minor differences can be observed between the default and very strict tolerance settings. The difference between a tolerance of 10 −3 and a tolerance of 10 −8 is S B = 0.0066 for the first attractor and S B = 0.0004 for the second attractor, which is a small deviation considering a basin stability value S B = 0.0923 (S B = 0.0822) for the first (second) attractor. Larger numbers of samples would even reduce the effect of the numerical integration scheme on the basin stability estimates.

Non-smooth frictional oscillator
A single-mass linear oscillator Mẍ + Cẋ + K x = F fric displayed in Fig. 10a with a non-smooth frictional contact is studied following previous works [20,22]. The friction force is described by velocity-dependent weakening friction slope μ (v rel ) shown in Fig. 10b giving rise to friction-induced self-excited stick-slip vibrations, see Fig. 10c. μ st = μ (0) denotes the static friction coefficient, μ d = μ (v rel → ∞) denotes the dynamic friction coefficient in the asymptotic limit, v d is the driving velocity of the conveyor belt, v 0 is the reference velocity, and the contact is loaded by the normal force F N . The non-dimensional form (·) of the system is obtained by The oscillator exhibits bi-stability for 1.11 ≤ṽ d ≤ 1.84: small initial displacements and velocities die out quickly as the trajectories converge to the stable fixed point of steady sliding. Larger initial states push the system into a stable stick-slip limit cycle which is selfexcited through the negative friction slope Eq. 13. For low velocitiesṽ d < 1.11 only a stable limit cycle exists, while for larger velocitiesṽ d > 1.84 the steady sliding state is globally stable. The correponding bifurcation diagram is shown in Fig. 11a along with the eigenvalue's real part in (b). The real part approaches the positive half plane asṽ d is decreased and hence may indicate a qualitative change of the dynamics in the sense of an early warning sign (in contrast to the pendulum case discussed in Sect. 4.1). Naturally, the coexistence of a stable limit cycle in the range of negative real parts cannot be obtained from the linear stability analysis of the fixed point. The basin stability analysis displayed in Fig. 11c represents a quantitative measure for both stable solutions as a function of the sliding velocity. For the given choice of the region of interest, the basin stability scales linearly with the velocity up toṽ d ≈ 1.83: for 1.11 ≤ṽ d ≤ 1.6 the stick-slip cycles are more likely to occur, and for 1.6 ≤ṽ d ≤ 1.84 the steady sliding state dominates the time-asymptotic behavior under the selected distribution of states. All  sampling from a uniformly spaced grid (c) are studied. The sampling points nicely indicate the unstable periodic orbit that separates both basins of attraction within the area encircled by the stick-slip trajectory.
The maximum difference of basin stability values amounts to 0.097 comparing the Gaussian distribution to the grid sampling. This study illustrates how the user can specify custom sampling strategies and distributions, and thus adapt the analysis to an individual system for which a specific distribution of admissible perturbations may be known. The structure of bSTAB requires changing only a single line (props.roi.samplingPDF = 'uniform') in the case definition for running those different analysis, and hence enables quick interaction of the user with the results, and reduces efforts related to changing programming codes.

Conclusion
This work introduces an open-source computing toolbox for estimating the basin stability of continuous dynamical systems exhibiting multi-stability. Complementing conventional local stability analysis, the basin stability represents a global stability metric for estimating the share of each basin of attraction in the system's state space. Thereby, basin stability translates into probabilities of arriving on one of the competing multi-stable solutions under a distribution of admissible initial states or instantaneous perturbations. Albeit the simplicity of the theoretical concept, practical issues and pitfalls can arise for researchers willing to implement a basin stability analysis of their dynamical systems of interest. Thanks to a plethora of software packages existing for stability and bifurcation analysis of nonlinear dynamical systems, these methods have made their way into the state-of-the-art toolbox in the scientific community. However, an easy-to use toolbox for computing the basin stability has not been proposed yet, potentially hindering the penetration of this rather new stability concept into the toolbox of nonlinear dynamics researchers.
The bSTAB project aims at filling this void. The proposed implementations constitute a highly formalized structure for computing the basin stability with as minimal interaction required by the user as possible. Following a model-agnostic approach, common implementations of dynamical systems can be re-used, such that the basin stability analysis can easily be incorpo-rated into existing analysis procedures. Efficient Monte Carlo simulations are enabled through parallelization strategies for multi-core and high-performance computers. Consistent quantification of multi-stability with error estimates is presented for four canonical nonlinear multi-stable dynamical systems featuring regular, chaotic and non-smooth dynamics, complicated and fractal basin boundary shapes as well as bifurcations along parameter variations.
The concept of open-source software is chosen to enable barrier-free access to the implementation, but also foster interdisciplinary cooperation among nonlinear dynamics researchers. Active code development by the community will help to grow functionalities further, to cover even more efficient implementations and programming languages, and to extent the capabilities of the toolbox to maps and network systems.
Funding Open Access funding enabled and organized by Projekt DEAL.

Declarations
Funding This research was funded by the German Research Foundation (Deutsche Forschungsgemeinschaft DFG) within the Priority Program 1897 'calm, smooth, smart', grant number 314996260.

Conflict of interest
The authors declare that they have no conflict of interest.

Availability of data and material Not applicable.
Code availability All code for reproducing the case studies is available through https://github.com/TUHH-DYN/bSTAB.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

A Definition of the pendulum case
The settings for the driven damped pendulum presented in Sect. 4.1 are given in Table 1 following the imple- mentation in setup_pendulum.m provided with bSTAB.

B Definition of the Duffing oscillator case
The case definition details for the Duffing oscillator case discussed in Sect. 4.2 is given in Table 2 and in setup_duffing.m provided with bSTAB.      Table 4 gives the case definition details employed for the friction oscillator case study.