Sample-Based Methods of Analysis for Multistable Dynamical Systems

In this paper we present how sample based analysis can complement classical methods for analysis of dynamical systems. We describe how sample based algorithms can be utilized to obtain better understanding of complex dynamical phenomena, especially in multistable dynamical systems that are difficult for analytical investigations. Relying on the simple, direct numerical integration algorithms we are able to detect all possible solutions including hidden and rare attractors; investigate the ranges of stability in multiple parameters space; analyse the influence of parameters mismatch or model imperfections; assess the risk of dangerous or unwanted behaviour and reveal the structure of multidimensional phase space. For each mentioned application we present methodology, example on paradigmatic non-linear dynamical system and discuss practical applications. The presented methods of analysis can be applied to solve numerous of scientific problems originating from different disciplines. Moreover, their robustness and efficiency will grow with the upcoming increase of computational power.


Introduction
The stability of dynamical systems have focused the interest of scientist for more than hundred years. One of the oldest example is the solar system. From our point of view, it is stable and we should not care about its destabilization, but looking in the timescale of a universe it becomes an important problem. Copernicus placed the Sun at the centre of the universe [18,19]. Then, Kepller assumed that planets move along the perfect ellipses, hence after full revolution they follow the same trajectory [93]. This beautiful simplicity has been destroyed be the Newtons law of universal gravitation [70] which proofs that interactions occur between the Sun and all objects that orbit it. Thus, there is no reason to assume that the planet orbits are fixed ellipses. In the 18th century Lagrange and Laplace correctly formulated the equations of motion and proofed that Newton's law is the universal explanation for the motion of the celestial bodies [68]. Since then, many scientists have considered this issue. The breakthrough work has been published by Poincaré [4,87] where he showed that it is not possible to integrate the equations of motion of three bodies subjected to mutual interaction, and not possible to find an analytic solution representing the movement of the planets valid over infinite time interval. From this moment, we observe an intensive investigation of this problem, but the real milestone was done by applying numerical methods. Numerical computations show that the Solar system is unstable, but catastrophic phenomena leading to its destruction can take place only in a time comparable with its age (approximately 5 billion years) [47,48,91].
The modern study of stability starts from fundamental works of Poincaré, Lyapunov and Floquet. Henri Poincaré wrote a series of texts under common title ''On curves defined by differential equations'' [76]. He initiated a new branch of mathematics called qualitative theory of differential equations. He showed that even if the differential equation cannot be solved in terms of known functions, one can extract a wealth of information about the properties of solutions. Poincaré main interest was the nature of trajectories of integral curves in the plane space. He provided a classification of singular points (saddle, focus, centre, node) and introduced the concept of a limit cycle. Two years after Poincaré, in 1883 Gaston Floquet published paper [31] in which he define the stability of periodic solutions basing on eigenvalues of monodromy matrix [45]. Nine years after the publication of Floquet, in 1892 Aleksandr Lyapunov defended his PhD thesis entitled ''The general problem of the stability of motion'' [57]. He introduced definition of equilibrium stability, i.e. & P. Brzeski piotr.brzeski@p.lodz.pl 1 Division of Dynamics, Lodz University of Technology, 90-924 Lodz, Poland asymptotically and exponentially stable orbits. The intensive development of numerical techniques let us calculate Lyapunov exponents [6,33,71,86] to characterize the stability of arbitrary type of solution. Basing on the sign of Lyapunov exponents we can also classify the types of quasi-periodic and chaotic solutions [7,39,74].
In aforementioned studies, the stability is calculated only in the close neighbourhood of the considered orbit. Such knowledge is enough if system is monostable, but if the system is multistable some perturbation may cause sudden change in behaviour and result in jump to different coexisting solution. Thus, we should gain better insight on the structure of the phase space in a wider range. It can be achieved via analysis of basins of attraction [5,29,95]. For systems with two generalized coordinates basins of attraction show all solutions as the behaviour of the system depends only on two initial conditions. Thus, we are able to reflect the structure of the phase space in a 2D plot and analyze the boundaries of basins of attraction. Most of practically inspired models have more then two dimensional phase space, which make investigations much more difficult because we can only look at the two-dimensional cuts of multi-dimensional phase space. This shows the importance of the analysis of stability and detection of basin boundaries. Such boundaries can be smooth or fractal [46,60,73]. If the the boundary of basin is a smooth curve or a surface the range of solution stability can be defined precisely. If boundary is more complex -fractal -initial conditions leading to different coexisting attractors are mixed [37] and we cannot predict the effects of even arbitraly small perturbation. To describe such structure and measure its complexity fractal dimension has been introduced. It provides a statistical index comparing how detail in a pattern changes with the scale at which it is measured. We can distinguish several measures like: Hausdorff dimension, box counting dimension, information dimension and ect. All of them can be used to characterize both basin and attractor [82]. If we focus on the structure of basins of attraction, we can distinguish following types: compact, eroded or riddle [40,41]. For compact basins of attraction the attractor is surrounded with its basin. Hence, within this range we can be sure that after perturbation the system will go back to the attractor in finite time. For eroded or riddle basins it is much harder to predict if a perturbation can cause a jump to another solution. The importance of such analysis has been presented in 1989 by Thomson [89]. He discussed the basin erosion and indicated that with varying parameters of the system its basin can change from compact to completely eroded.
The studies of basins of attraction give a detailed knowledge if the dimension of system is low, but if we analyze multidimensional system we can miss important information. To overcome this problem a new methods of investigations have been proposed. Rega and Lenci proposed ''basin integrity measures'' [54,55,80] that describe the erosion of basins, i.e., global integrity measure and integrity factor. These measures enable to asses if the basin has safe, compact structure or it has fractal structure. Still, basin integrity measures are also hard to obtain for highdimensional systems. Nevertheless, for low dimensional systems this method gives precise information about the multistable system dynamics. The advantages of this method has been proved on multiple investigations of various systems including systems with impacts [50][51][52], parametric excitation [53,54], shells [34] and single-mode model of non-contact atomic force microscopy [81]. All these studies show that we can assess the structure of basins of attraction and basing on the obtained results we are able to propose effective control algorithms to reduce the erosion of basins.
The second method called ''basin stability'' has been developed by Menck et al. [62,63]. Basin stability method enables to quantify stability basing on the probability of reaching given attractor from random initial conditions. To calculate basin stability measure one has to perform a significant number of Bernoulli trials and classify the final solutions reached in each trial. Proposed idea has been successfully applied to asses the stability of power grids [63,85], systems with time delay [56], chimera states [79], stabilization of saddle fixed points in coupled oscillators [78] and brain dynamics [61]. In our previous papers we proposed extensions of this method [12]. Additionally to initial conditions we draw the parameters of the system. Such approach let us include in the analysis the uncertainty of parameters or just perform the very efficient investigation of coexisting solutions for systems with varying parameter values. We validate the accuracy of the proposed method with experimental study [15]. The next method ''basin entropy'' has been proposed by Daza et al. [22]. They propose another sample based algorithm to test the structure of basins of attraction to obtain information about its fractality. Basin entropy has been used for investigate the dynamics of propagating matter waves [21,23] and the flux in twist Hamiltonian systems [67].
The three mentioned methods need very efficient algorithms and high computational power. Still, nowadays computational power is constantly increasing and let us use sample-based methods to analyze systems from different disciplines of science. Thanks to that, many complex problems can be solved using numerical algorithms instead of analytical derivations. Advanced computing enables to simplify scientific analysis and gives new possibilities for solving the most important problems.
In the paper we present four possible utilizations of sample-based analysis. In Sect. 2 we briefly present the basics of the basin stability and describe sample applications of this method. In Sect. 3 we show the expansion of the basin stability concept that enables effective analysis of systems with parameters mismatch. It enables to analyze the effects of system imperfections and varying working conditions. Sect. 4 is devoted to the description of how sample based methods can be used to determine ranges of stability and predict conditions that ensures the highest probability of reaching given solution. The presented approach has similar precision to classical methods of analysis but it is straight forward and can be applied for all types of dynamical systems. Moreover, basing on the data obtained with the algorithm we are able to characterize the structure of the phase space. It is described in Sect. 5 where we show how sample based approach supplements basins of attraction and dynamical integrity measures. It is especially important for multi-DoF systems where we cannot project the basins of attraction to present their boundaries and volume. In all sections we start with paradigmatic models to present the underlying methodology and supplement the description with real world example of possible utilization of the method. We conclude the presented material in Sect. 6.

Basin Stability Method
In this section, we show in detail the idea of basin stability method which enables to estimate the stability and number of solutions for given values of system parameters. The idea is strikingly simple but the method has found numerous successful applications and proved to be an efficient tool for dynamical analysis. The growing interest in this method comes from its two main advantages. Firstly, it is suitable for multidimensional and multistable systems where the classical approaches are difficult and time consuming to apply. Secondly, it can be easily applied to all types of systems and reproduces the inherent uncertainty of perturbations. To apply basin stability method we just need to have a reliable direct numerical integration code for the mathematical model of the investigated system. The computational effort does not grow significantly with the increase of the phase space dimensions. This makes this sample-based method even more appealing for the analysis of very high-dimensional systems such as large networks of oscillators, power grids or neural networks. Moreover, in our previous paper we prove that the results form basin stability method match with the results form the path-following analysis and experiment [15].

Methodology
The algorithm is simple and based on the series of N Bernoulli trials. For each trial we randomly select initial conditions and detect the final attractor using direct numerical integration. Based on this one can calculate the chance to reach given solution and determine the distribution of probability for all coexisting solutions. This gives information about the number of stable solutions and the sizes of their basins of attraction. Consider the dynamical system with fixed parameters _ x ¼ f ðx; tÞ, where x 2 R n is the state vector and t 2 R is time. Let B & R n be a set of all possible initial conditions. Let us assume that attractor A exists, is stable and has a basin of attraction bðAÞ. Lets assume that within N overall trials attractor A was reached n A times. Assuming random initial conditions the probability that the system will reach attractor A is given by: This probability is the estimator of the basin stability measure that reflects the relative volume of the basin of attraction of attractor A. If this is equal to one ðp A ð Þ ¼ 1:0Þ this means that considered solution is the only one in the taken range of initial conditions and given values of parameters. Otherwise, other attractors coexist and system is multistable. Due to physical limitation we cannot take the initial conditions form the infinite set but we draw them form set B A & B. We can consider two possible ways to impose the limit to this set and both of them are based on preliminary analysis of the system. The first ensures that set B A includes values of initial conditions leading to all possible solutions. This approach is appropriate if we want to get general overview of the system's dynamics. In the second approach we use narrowed set of initial conditions that corresponds to practically accessible initial states.

Applications
The basin stability method has been already utilized in numerous applications. Originally, it has been created to analyze the synchronizability of Watts-Strogatz networks consisting of paradigmatic Rössler oscillators [62]. Authors vary the topology of the network starting form regular to random structure (which corresponds to diffident applications e.g. visual cortex, power grids and neural networks). The main conclusion is that the synchronous state is much more stable in networks that are more regular. In [85] the relationship between stability against large perturbations and topological properties of a power transmission grid has been analized. Authors introduce a very fast tool to predict the appearance of nodes with poor single-node basin stability. This algorithm can be applied during live operation of a power grid to increase safety of the network. The next study has been focused on details of modelling of power grids [3]. It proves that voltage dynamics should be taken into account in model of synchronous machines. In mountain ranges, big lakes or inland seas we often observe large closed loops in high voltage AC power grids [17]. The basins stability method has been used to identify three mechanisms which create the circulating power flows. Moreover, one can find significant number of studies that use basin stability method to show the way to increase the efficiency of power grids [2,42,43,83] Maslennikov et al. [61] study the basin stability of the burst synchronization regime in small-world networks consisting of chaotic slow-fast oscillators. The results confirm that the coupling density and the coupling strength influence the basin stability similarly and there are threshold values above which the basin stability of the burst synchronization regime significantly increase. The same method has been used to analize the explosive transitions between synchronous and non-synchronous states in networks [96]. Next, in [44] Authors consider synchronized state in time-varying complex networks (the time-varying character is ensured by stochastic rewiring of links with defined frequency). In [64] Mitra et al. propose the general framework of multiple-node basin stability for gauging the global stability and robustness of networks in response to non-local perturbations simultaneously affecting multiple nodes of a system. It helps to estimate the threshold number of simultaneously perturbed nodes that reduce the capacity of the system and triggers the return to the original solution.
In all dynamical systems the final state is always achieved after some transient time. Quite often, it is very hazardous range of the time evolution because the trajectory of the system can move very close to the boundary of stability and even a small perturbation can result in jump to another solution. To study transient dynamics an extension of basin stability method ''constrained basin stability'' has been created [92]. The time of recovery after perturbation has been discussed in [65].
The interesting application of the basin stability method to systems with time delay has been presented in [56]. The time delay is present in population dynamics, epidemiology, neurobiology, control engineering, optoelectronic and many more [1,32,35]. Due to presence of time delay the dimension of system phase space is infinite. Thus, the elaboration of the efficient tool to analyze solutions in this class of equations is extremely important. Authors propose a technique which projects the infinite dimensional initial state space to a finite-dimensional Euclidean space using three different bases: the Bernstein, the trigonometric and the Legendre and show on numerous examples that it is very effective. The basin stability has been also applied to piecewise and discontinuous systems (Amazonian vegetation model [66], externally forced oscillator with impacts [12] and system with dry friction [30]) where it gives information of co-existing solutions and volumes of their basins of attraction in multidimensional phase space.

Investigation of Systems with Parameters Mismatch
Basin stability measure was initially proposed as a tool to characterize the stability of different synchronous states of complex networks of oscillators. Thanks to its paradigmatic simplicity basin stability method can be easily expanded by drawing values of system's parameters. The underlying idea is to take into consideration the fact that values of system's parameters are measured or estimated with some finite accuracy which is often hard to determine.
Moreover, values of parameters can vary during normal operation. Still, in practical applications we usually need certainty that presumed solution has enough surplus of stability to ensure safe operation. Hence, there is a need to describe how small changes of parameters' values influence the behaviour of the system. Our method provides such description and allows to estimate the required accuracy of parameters' values. Moreover, we can also consider the risk of unwanted phenomena under presence of model imperfections and parameters mismatch. Hence, it is time efficient, can be applied in a wide range of systems and does not require high computational power. The method is described in detail in our previous paper [12]. In this Section we briefly recall the proposed methodology, present results from paradigmatic systems and discussed some potential practical application.

Methodology
Consider the dynamical system _ x ¼ f ðx; t; xÞ, where x 2 R n is the state vector, t 2 R is time and x 2 R m are the system parameters, where m is the number of parameters that are taken into account. Let B & R n be a set of all possible initial conditions and C & R m a set of accessible values of system parameters. Let us assume that attractor A exists for x 2 C A & C and has a basin of attraction bðAÞ. In classical approach of Menck et al. [62] the stability of attractor A is assessed basing on the probability of reaching A with random initial conditions but for fixed parameters (see Sect. 2.1).
We propose to not only draw initial conditions but also values of some selected parameters of the system from the set C A & C that consists of all practically accessible values of system's parameters x. This let us ensure that given solution exists in this range (taking into account the mismatch in parameters). Next, we subdivide set C A in to m ¼ 1; 2; . . .M equally spaced subsets. Subsets C m A do not overlap and the relation S m¼1...M C m A ¼ C A is always fulfilled. Then for each subset C m A we randomly pick N sets of initial conditions and value of the considered parameter. For each set we check the final attractor of the system. After sufficient number of Bernoulli trials we can calculate the probability of reaching presumed solution or solutions. As a result we can describe the relation between the value of the system's parameter and volume of basin stability of reachable solutions. Value of N strongly depends on the complexity of analysed system, i.e. mainly on the size of phase space, the number of coexisting solutions and the structure of basins of attraction and size of initial conditions and parameters sets. Moreover, the computation time for single trial and should be adjusted for each system independently to ensure that the final attractor is reached. In general, we recommend that in most cases N should be at least one hundred.

Paradigmatic Examples
In this Section we present the methodology using two paradigmatic nonlinear dynamical systems: pendulum suspended to the Duffing oscillator and bilinear impacting oscillator. Both systems are excited with harmonic external force and have multiple stable attractors.
The first example presented in Fig. 1a is the Duffing oscillator with suspended pendulum that acts as a tuned mass absorber. The main body consists of mass M fixed to the ground with nonlinear spring (hardening characteristic k 1 þ k 2 y 2 ) and viscous damper (damping coefficient c y ). The mass M is forced externally by harmonic excitation with amplitude F and frequency x. The pendulum has length l and mass m. The small viscous damping is present in the pivot of the pendulum that is described by parameter c u .
The derivation of equations of the system's motion can be found in [13]. Here, we present the dimensionless equations of motion: . As control parameters we take the amplitude f and the frequency x of external forcing. The dimensionless parameters have the following values: For simplicity primes in dimensionless equations will henceforth be neglected and we use x and u as dimensionless generalized coordinates of the system. The second considered paradigmatic example is the system with impacts investigated previously in [72]. The system is shown in Fig. 1b and consists of mass M suspended by linear spring with stiffness k 1 and viscous damper with damping coefficient c to harmonically moving frame. The frame oscillates with amplitude A and frequency X. When amplitude of mass M motion reaches value g we observe soft impacts (spring k 2 is significantly stiffer than spring k 1 ). To equation of motion is derived in [72]. Here, we present dimensionless form of equation and to obtain it we use the reference length y 0 ¼ 1 [mm]. Thus, the dimensionless equation of motion is as follow: (a) (b) Fig. 1 The models of the considered paradigmatic systems: externally forced Duffing oscillator with attached pendulum (tuned mass absorber) (a) and bilinear impacting oscillator (b) Sample-Based Methods of Analysis for Multistable Dynamical Systems 1519 n y 0 are the dimensionless vertical displacement, velocity and acceleration of mass M, s ¼ x n t is the dimensionless time, x n ¼ k 1 M is natural frequency of linear system, b ¼ k 2 k 1 is the stiffness ratio, e ¼ g y 0 is the dimensionless gap between equilibrium of mass M and the stop suspended on the spring k 2 , a ¼ A y 0 and x ¼ X x n are dimensionless amplitude and frequency of excitation, n ¼ c 2mx n is the damping ratio and HðÁÞ is the Heaviside function. In our calculations we take the following values of system's parameters: a ¼ 0:7, n ¼ 0:01, b ¼ 29 and e ¼ 1:26. As a controlling parameter we use frequency of excitation x.

Duffing System with Suspended Pendulum
At the beginning, we want to recall the results we have presented in our previous paper [13]. As a summary we show two dimensional bifurcation diagram obtained by path-following method in Fig. 2. It shows bifurcations for varying amplitude f and frequency x of external excitation (see Eq. 2). Lines shown in the plot correspond to different types of bifurcations (period doubling, symmetry breaking, Neimark-Sacker) and resonance tongues. The main difference between solutions is the behaviour of pendulum, hence we introduce the notion which let us distinguish them. The Duffing system follows the excitations, while pendulum oscillates or rotates. The number of full period of pendulum motion differs for one period of excitation. To define it, we introduce the ratio u:z which gives information how many periods of excitation u are performed for z number of pendulum oscillations or rotations (the period of given solution corresponds to u multiply by period of excitation). We present these lines in one style because the structure is too complex to follow bifurcation scenarios and we do not need that data (details are shown in [13]). We mark areas where we observe the coexistence of one (black colour), two (grey colour) and three (hatched area) stable solutions. The remaining part of the diagram (white area) corresponds to situation where there are four or more solutions simultaneously. Additionally, by white colour we also mark areas where only Duffing system is oscillating in 1:1 resonance with frequency of excitation and the pendulum is in stable equilibrium position, i.e., hanging down pendulum (HDP) state. In this case the dynamics of the system is reduced to the oscillations of summary mass (M þ m). As we can see, the range where less than three solutions exist is rather small.
The detailed analysis of the system given by Eq. 2 is time consuming and creation of Fig. 2 was preceded by complex analysis done with large computational effort. Additionally, the obtained results do not give us information about the size of basins of attraction of each solution and their stability margins. Thus, using the path-following analysis we can obtain detailed bifurcation structure of system but cannot assess which solutions are practically accessible and are expected to occur in the real system. For example some solutions may need initial conditions which are not accessible in practice (too large displacement or velocity) or have minimal stability margin such that even small perturbation may cause the jump into another attractor. We will present that sample based analysis serves Þ showing periodic oscillations and rotations of pendulum. Black colour indicates one attractor, grey colour shows two coexisting attractors (the same as for black but with coexisting stable steady state of the pendulum). In the hatched area we observe the coexistence of stable rotations and stable steady state of the pendulum. Figure reprinted from [12]. Detailed analysis can be found in [13] greatly for that purpose. For illustrative purposes we focus on three solutions: 2:1 oscillating resonance, HDP and 1:1 rotating resonance assuming that only they have practical meaning.
Firstly, we show reference results obtained with direct numerical integration. These are two bifurcation diagrams calculated for constant amplitude of excitation f ¼ 0:5 and frequency from the range x 2 ½0:1; 3:0 (see Fig. 3). In Fig. 3a we increase excitation frequency x from 0.1 to 3.0 and in subplot (b) we decrease x from 3.0 to 0.1. As the initial conditions for the first points of the bifurcation diagrams we take equilibrium position ðx 0 ¼ _ In both panels we plot amplitude of pendulum u which reflects the type of attractor. In the considered range there are two dominating solutions: HDP and 2:1 internal resonance. Near x ¼ 1:0 we observe a narrow range of 1:1 and 9:9 resonances and chaotic motion (for details see Fig. 6 in [13]). Ranges where diagrams differ are marked with grey rectangles. Based on previous results we know that we detect all solutions existing in the considered range, however we do not have data about size of their basins of attraction and coexistence. Hence, the analysis with proposed method should give us new important information about system's dynamics. Contrary to bifurcation diagram obtained by path-following presented in Fig. 3 we do not observe a rotating solutions (the other set of initial conditions should be taken). This shows the main disadvantage of classical bifurcation diagrams. Now, we apply the sample based analysis. The initial conditions are random numbers drawn from the following ranges: x 0 2 ½À 2; 2, _ x 0 2 ½À 2; 2, u 0 2 ½À p; p and _ u 0 2 ½À 2:0; 2:0 (ranges there selected basing on the results from [13]). Frequency of excitation is within range x 2 ð0:0; 3:0 (Fig. 4a, c), then to have a zoom in narrower range, we refine it to x 2 ½1:25; 2:75 ( Fig. 4b, d). In both cases we take 15 equally spaced subsets of x and in each subset we calculate the probability of reaching given solution. For each subset we calculate N ¼ 1000 trials each time drawing initial conditions of the system and value of x from the appropriate range. Then we plot the dot in the middle of the subset which indicate the probability of reaching given solution in each considered range. Lines that connect the dots are marked just to show the tendency. For each range we take N ¼ 1000 because we want to estimate the probability of solution with small basin of stability (1:1 rotating periodic solution).
In Fig. 4 we show the probability of reaching three aforementioned solutions obtained using the proposed method. In Fig. 2 we see that the 2:1 resonance solution exists in the area marked by black colour around x ¼ 2:0 and coexists with HDP in neighbouring grey zone. In Fig. 4 we mark the probability of reaching 2:1 resonance (p(2:1) ) using blue dots. As we expect for x\1:4 and x [ 2:2 the solution does not exist. In the range x 2 ½1:4; 2:2 the maximum value of probability pð2 : 1Þ ¼ 0:971 is reached in the subset x 2 ½1:8; 2:0 and outside that range the probability decreases significantly. To check if we can reach pð2:1Þ ¼ 1:0 (so that it is the only stable attractor) we decrease the range of parameter's values to x 2 ½1:25; 2:75 and the size of subset to Dx ¼ 0:1 (to remain with 15 equally spaced subsets). The results are shown in Fig. 4b using blue colour. In the range x 2 ½1:95; 2:05 the probability p(2:1) is equal to unity and in range x 2 ½1:85; 1:95 it is slightly smaller pð2:1Þ ¼ 0:992. Hence, for both subsets we can be nearly sure that system reaches 2:1 solution. This gives us indication of how precise we (a) (b) Fig. 3 Bifurcation diagrams obtained using direct numerical integration. For subplot (a) value of bifurcation parameter x was increased while for subplot (b) we decreased value of x. Gray rectangles mark the ranges of bifurcation parameter x for which different attractors coexist. Figure reprinted from [12] have to set the value of x to be sure that the system will behave in a presumed way. This has practical importance and can be used when designing the system which cannot be multistable. The similar analysis is performed for HDP solution. The values of probability are reflected using red dots on the same plots in Fig. 4. For x 0:8, x 2 ½1:2; 1:4 and x 2 ½2:6; 2:8 the HDP is the only existing solution. The rapid decrease close to x % 1:0 indicate the 1:1 resonance and the presence of other coexisting solutions in this range (see [13]). In the range x 2 ½1:2; 1:4 the basin stability of HDP solution is a unity (pðHDPÞ ¼ 1:0) which corresponds to a border between solutions born from 1:1 and 2:1 resonance. Hence, up to x ¼ 2:0 the probability of HDP solution is a mirror refection of p(2:1). The same tendency is observed in the narrowed range as presented in Fig. 4b. Finally, for x [ 2:0 the third considered solution comes in and we start to observe an increase of probability of rotating solution as shown in Fig. 4c. However, the chance of reaching rotating solution remains small and never exceeds pð1 : 1Þ ¼ 8 Â 10 À3 . We also plot the probability of reaching rotating solution in the narrower range of x in Fig. 4d. The probability is similar to the one presented in Fig. 4c-it is low and do not exceed pð1 : 1Þ ¼ 8 Â 10 À3 . Note, that results presented in Fig. 4a, b and Fig. 4c, d are computed for different sets of random initial conditions and parameter values, hence the obtained probability can be slightly different. Still, it is inherent feature of the sample based methods and does not impede their advantages.

Bilinear Impacting Oscillator
Now, we present the analysis of the bilinear system with impacts. A discontinuity usually increases the number of coexisting solutions. Hence, in the considered system we observe a large number of different stable orbits, but we will consider only periodic solutions of different periods. In Fig. 5 we show two bifurcation diagrams with x as controlling parameter. Both of them are obtained by starting with zero initial conditions x 0 ¼ 0:0 and _ x 0 ¼ 0:0. In panel (a) we increase x from 0.801 to 0.8075; while in panel (b) we decrease x in the same range. We select the range of x basing on the results presented in [72]. As one can see, both diagrams differ in two zones marked by grey colour. Hence, we observe a coexistence of different solutions, i.e., in the range x 2 ½0:8033; 0:8044 solutions with period-3 and -2 are present, while in the range x 2 ½0:8068; 0:8075 we detected solutions with period-2 and -5 (numbers indicate how many times the period of solution is longer then the period of excitation). As presented in [72], some solutions appear from a saddle-node bifurcation and we are not able to detect them with the classical bifurcation diagram. The proposed method solves this problem and shows all existing solutions in the considered range of excitation frequency.
We focus on periodic solutions with periods that are not longer than eight periods of excitation. Solutions with higher periods are observed in the narrow range of x but the probability that they will occur is very small so they can be neglected. All non-periodic solutions are chaotic (quasiperiodic solutions are not present in this system). The results of our calculations are shown in Fig. 6a, b. We take initial conditions from the following ranges x 0 2 ½À 2; 2, _ x 0 2 ½À 2; 2. The controlling parameter x is changed from 0.801 to 0.8075 with step Dx ¼ 0:0005 in Fig. 6a and from 0.806 to 0.8075 with the step Dx ¼ 0:0001 in Fig. 6b. In each subrange of excitation's frequency we pick the exact value of x randomly from this subset, as described earlier.
The probability of reaching periodic solutions is plotted with lines of different colours that are created by joining markers that reflect probability obtained for given range of x. The solutions of different periods were detected and investigated, namely solutions with period-1, -2, -3, -5 (two different attractors with large and small amplitude), -6 and -8. We also calculate the sum of reaching one of all periodic solutions' (also with period higher then eight). It is plotted with yellow markers and lines. When its value is below 1, chaotic solution exists. Dots are drawn for mean value i.e, middle of the subset. For each range we take N ¼ 200 and we increase the calculation time because the transient time is sufficiently larger than in the previous example due to the piecewise smooth characteristic of spring's stiffness.
As we can see, the chance of reaching a given solution strongly depends on x. Hence, in the sense of basin stability we can say that x value strongly affects stability of attractors. In Fig. 6a the probability of a single solution is always smaller than one. Nevertheless, there are two dominant solutions: period-5 with large amplitude in the first half of the considered x range and period-2 in the second half of the range. The maximum registered value of probability is pðperiod À 2Þ ¼ 0:92 and it refers to the period-2 solution for x % 0:80675. To check if we can achieve even higher probability we analyse a narrower range of x and decrease the range from which we draw x value (from Dx ¼ 0:0005 to Dx ¼ 0:0001). In Fig. 6b we see that in range x 2 ½0:8069; 0:807 the probability of reaching the period-2 solution is equal to unity. Hence, in the sense of basin stability it is the only stable solution. In the range x 2 ½0:8065; 0:8072 the probability of reaching this solution is higher then 0.9 and we can say that its basin of attraction is strongly dominant. Hence, the obtained results give the required tolerance of x to omit  [12]. Further analysis can be found in [72] Sample-Based Methods of Analysis for Multistable Dynamical Systems 1523 multistability and be sure that period-2 solution will be observed. This proves the practical importance of the method.
The same information can be also presented using the relative basin volume plots as shown in panels (c, d) in Fig. 6. Such plots are predisposed for presentation when we consider multiple solutions as they enable easy comparison between basins volume and make presentation more readable when many solutions have similar basin probability measure.

Practical Application
Now, we will present an example of practical application of the above technique basing on the hybrid dynamical model of the yoke-bell-clapper system. The model was validated experimentally [11] and thoroughly analysed using classical methods including direct numerical integration [10] and path-following [9,16]. The model possess a plethora of interesting dynamical phenomena including multistability [9,16] and new type of bifurcation called ''impact adding bifurcation'' [9]. This system is highly predisposed to sample analysis because in most cases we cannot determine parameter values with good precision. It is due to the fact that bellfounding (casting of bells) is still based more on traditional methods and local customs than engineering science. Similarly, when designing the yoke and propulsion one does not have the precise knowledge about the shape of the bell and clapper. Thus, we usually observe some imperfections and the real system is just a rough embodiment of the assumed design. This problem can be partly solved with the sample based analysis, in which we can introduce parameters mismatch and find ranges of crucial parameter values which guarantees the highest probability that the system will work properly. The investigated model is based on the existing bell ''The Heart of Lodz'' of the Cathedral Basilica of St Stanislaus Kostka, Lodz, Poland and all parameters values were obtained in a series of dedicated experiments [10,11]. The derivation of the model is described in detail in [11] and now we just briefly recall it. The developed mathematical model is based on the analogy between freely swinging bell and the motion of the equivalent double physical pendulum. The first pendulum has fixed axis of rotation and models the yoke together with the bell that is mounted on it. The second pendulum is attached to the first one and imitates the clapper. Figure 7a, b show schematics indicating the position of the rotation axes of the bell o 1 , the clapper o 2 and presenting parameters involved in the model. For the sake of simplicity, henceforth, the term ''bell'' is used for the bell and its yoke, which is treated as one solid element.
The model has eight physical parameters: L 0 describes the distance between the rotation axis of the bell and its centre of gravity (point C b ), l is the distance between the rotation axis of the clapper and its centre of gravity (point C c ). The distance between the bell's and the clapper's axes of rotation is given by l c0 . The mass of the bell is described by M, while B b0 characterizes the bell's moment of inertia referred to its axis of rotation. Similarly, m describes the mass of the clapper and B c stands for the clapper's moment of inertia referred to its axis of rotation.
We have to remember that we consider a musical instrument, hence we cannot change most of its parameters as it could affect the sound it generates. In real applications we can easily modify the driving motor and the mounting of the bell (by changing the design of the yoke). Therefore, in our investigation we will consider alterations of these two features taking as a reference parameter values that refer to ''The Heart of Lodz''.
Parameter l r is used to describe the modifications of the yoke, as it is presented in Fig. 7b, c. The l r definition is explained in detail in our previous paper [10], where l r ¼ 0 refers to the shape of original yoke used in the Cathedral's Basilica bell. If the centre bell's of gravity is lowered with respect to its axis of rotation, l r \0, otherwise l r [ 0. The change of the yoke design given by the value of l r affects other parameters. Therefore, in the mathematical model, the following parameters that describe the system with the modified yoke are used: The bell is excited with a linear motor. When the deflection of the yoke is smaller than p=15 ½rad ð12 Þ the motor is active and generates the torque. The torque generated by the motor Tðu 1 Þ is given by the piecewise formula: where T max is the maximum torque. Although, the above expression is not fully accurate reflection of the effects generated by the linear motor it is able to reproduce the characteristics of the modern bells' driving mechanisms [11]. We can easily modify the effects generated by the motor by changing the range of bell's deflection in which the motor is active (in our case Àp=15; p=15 ½ ) or by altering the maximum generated torque T max which is much easier to realize in practice [10]. Therefore, we take T max as the second control parameter.
As shown in Fig. 7a we use a planar co-ordinate system where the angle between the bell's axis and the downward vertical is given by u 1 and the angle between the clapper's Fig. 7 Schematics of the physical model of the yoke-bell-clapper system in different planes to show its geometry, kinematics and physical parameters Sample-Based Methods of Analysis for Multistable Dynamical Systems 1525 axis and downward vertical by u 2 . Angular parameter a describes the impact condition as follow: Synonymously, contact between the bell and the clapper occurs when a relative angular displacement between the bell and the clapper is equal to a. The Lagrange equations of the second type are employed to derive two coupled second order ODEs that describe the motion of the system: where g stands for gravity. We use a discreet impact model. If Eq. 6 is fulfilled, the numerical integration process is paused. Then, simulation is restarted with updated initial conditions. The bell's and the clapper's angular velocities are swapped from the values before the impact to the ones after the impact. The angular velocities after the impact are obtained by taking into account the energy dissipation and the conservation of the system's angular momentum that are expressed by the following formulas: where index AI stands for ''after impact'', index BI for ''before impact'' and parameter k is the coefficient of energy restitution and in our simulations we assume k

Results
Depending on local customs and traditions bells can have different working regimes. We can alter them easily by changing the design of yoke and/or propulsion [10]. Still, for some values of parameters the system is multistable and the final attractor depends on initial conditions or starting procedure. Moreover, as mentioned before, the final object is often a rough realization of initial concept, and we cannot predict the actual values of system parameters before it is build. Thus, instead of analyzing the dynamics for some fixed values we can apply the sample based approach. Including some parameters mismatch we can estimate the ranges of parameters that ensure proper working conditions and define the precision required to achieve mandatory level of reliability. As a reference we use the results presented in [9] where we consider alteration of the yoke design given by parameter l r . As a reference we use ''The Heart of Lodz'' which works in ''symmetric falling clapper'' regime (collisions between the bell and the clapper occur when they perform an antiphase motion). We consider the range of parameter l r 2 ½À 0:3; 0:2 [m] in which we observe 13 different types of attractors and multistable regions. The detailed description of all considered solutions is given in [9], where we also present in depth analysis done using the path-following method [20] and direct numerical integration. The obtained results enable to detect the ranges where we presume that there is only one stable attractor. Moreover, for multistable regions we do not have the knowledge about the chances of reaching given solution and structure and size its of basins of attraction. To obtain such information we apply the sample based analysis described above.
We divide the investigated range of l r 2 ½À 0:3; 0:2 [m] into 50 subranges of equal length 0.01 [m] which reflects the maximum expected precision of assembly. For each subset we perform 40,000 trials of direct numerical integration each time drawing the values of l r and initial conditions from accessible ranges which are assumed as follows: additionally condition that ensures that clapper is inside the bell u 1 0 À u 2 0 a must be fulfilled. Then, for each trial we detect the final attractor and estimate probability of reaching given solution in each considered range. The results are shown in Fig. 8. In panel (a) we present the bifurcation diagram obtained using multiple trials of direct numerical integration (it was previously shown and described in [9]). In panel (b) we show how the probability of reaching given solution depends on the value of l r .

Analysis of Multistability
There is a rich variety of mathematical tools to analyze multistable dynamical systems. Still, more sophisticated methods are usually difficult to apply. For example, there are a number of different toolboxes that enable the pathfollowing analysis but their functionality is strictly limited to the type of the investigated system and its dimensionality [25,28,90]. The dynamical analysis is especially challenging for multistable systems, where we have to consider multiple steady states that coexist in the phase space. It is a challenging problem and multistability is widely studied in many disciplines [26,38,58,59,75,84,88,94]. Therefore, new sample based tools for dynamical analysis are being developed such as survivability [36] that includes the analysis of the transient motion and basin entropy [22] that measures the basin compactness. The concept presented in the previous Section is also suitable for investigation of multiple coexisting attractors. Moreover, it enables to perform simultaneous analysis in a multi-parameter space. In this Section we will show analysis in two-parameter space and compare the results from sample-based analysis with detailed two parameter bifurcation diagrams obtained using the path-following method. Then, to validate the method we confront the results with experimental data. By that we are able to critically compare the accuracy of both methods and show their strengths and weaknesses. Finally, we show that the sample-based approach can be applied without a precise knowledge of parameter values and it gives results that can be used in practice.

Investigated System
We consider the specific type of a double pendulum (see Fig. 9) which is a paradigmatic example in nonlinear dynamics. The first pendulum rod is mounted horizontally and connected to the base with a pin joint at one end and via a spring on the second end. Hence, the first pendulum always performs oscillatory motion. The second pendulum is connected to the first one with a rotational pivot at the distance x 1 between both pin joints. The support is mounted on a shaker and excited kinematically in the vertical direction. This is a physical model of the existing experimental rig shown in the photographic picture in Fig. 9b. The dynamics of this system was analysed in the recent article [15].
The system has two degrees of freedom. The angular displacement of the first pendulum is given by generalized coordinate u 1 and the position of the second pendulum by u 2 . The shaker excites the system kinematically with a harmonic function of the amplitude A and frequency x. The upper pendulum has the length l 1 , the mass m 1 , the moment of inertia J 1 and the centre of gravity located at the distance d 1 from its pin joint. The stiffness of the spring that supports the end of the pendulum is given by k. The second pendulum has the mass m 2 , the moment of inertia J 2 and the center of gravity located at d 2 from its pin joint. Due to the relatively small damping in the system we neglect effects of dry friction and other non-linearities of damping characteristics and in both pivots we assume the viscoelastic damping. For the first pendulum the damping coefficient is given by c 1 and for the second one by c 2 .
We have to estimate the ranges of initial conditions that can be applied on our experimental rig. In a number of dedicated experiments we defined the accessible ranges of initial conditions that describe physical restrictions of our rig: u 1 À 0:0015; 0:0015 ½ ½ rad, u 2 À p; p ½ ½rad, _ u 1 ½À 1:2; 1:2 ½rad=s, _ u 2 ½À 60; 60 ½rad=s. Then, we divide the considered ranges of parameter values into a regular two dimensional grid in the following way: for the frequency of excitation (x½0; 60 ½rad=s) we assume 23 subsets with equal width of 2:5 ½rad=s. We intentionally omit extremely small values of x and start from the range 1:25; 3:75 ½ ½ rad=s (for first column of the grid) and finish with 56:25; 58:75 ½ ½ rad=s (for the last column). For the amplitude of forcing A 0; 7:7 ½ ½10 À3 m we take 15 equally spaced subsets with the step of 0:5 ½10 À3 m. Here, we also ignore values near the accessible boundaries and start the lower row of the grid with the range of 0:25; 0:75 ½ ½ 10 À3 m and finish the grid with 7:25; 7:75 ½ ½ 10 À3 m. By that, we receive a lattice of 345 boxes each covering the range of 2:5 ½rad=s and 0:5 ½10 À3 m of the forcing frequencies and amplitudes respectively. For each box in the grid we perform 500 trials of direct numerical integration. Each time we randomly pick the initial conditions from the accessible ranges and draw the values of A and x from the ranges assigned to the grid box. For every trial we recognize the final attractor that is reached by the system. In a procedure described above, we obtain a large data set of 172, 500 Bernoulli trials that we use to characterize the possible behaviors of the system.

Methodology
We investigate the ranges of stability of coexisting solutions in two-parameters space using two different approaches and validate them with experiments. We start with the path-following method to get a precise boundaries of solutions' stability in two-parameters. The path-following analysis is done with the AUTO-07p software [24]. Each time we start with a direct numerical integration to prepare the initial periodic solution for continuation. Then, we perform with one parameter continuation in selected parameter. We detect the bifurcation points in which the stability of solution changes and perform continuation of the bifurcation points in two parameters space. For the solutions detected numerically, we experimentally find the boundaries of the stability similarly in the two parameter space (here we use the amplitude and frequency of excitation A and x, respectively).
Then, for each solution we perform experiments to assess the stability ranges. The starting points for the experimental study are located far from the boundary of stability to ensure that the given solution is reachable experimentally. Depending on the shape and area of the stability range, we select a number of starting points and for each apply the following procedure: • We start the parameter values taken from numerical results and try to reach the given solution by applying proper initial conditions. • When the system reaches the presumed state we mark current parameters values as a starting point (black dots in Fig. 10) in the two parameter space. • Then, we slowly change the value of one parameter with given step. After each step we wait to check if the system remains on the presumed attractor. The range of the parameter value where the solution stays stable is marked as a trace with a black line in Fig. 10d, e. • Eventually, we reach the parameter value for which the solution becomes unstable and the system jumps to a different attractor. We note the parameter value for which the solution changed. • In each case we repeat the above steps four times and compute the average values of parameters for which the solution looses its stability. We take this value as the boundary of stability and mark it with the perpendicular end of the stability trace (see Fig. 10d, e).  Fig. 9 The physical model of the considered double pendulum excited kinematically (Eq. 12) with its parameters

Sample-Based Methods of Analysis for Multistable Dynamical Systems 1529
To estimate the position of a line that indicates the boundary of stability obtained by the path-following method in the two parameter space, we repeat the above procedure for different pairs of parameters values (starting points). Finally, we apply the sample-based approach (as described in Sect. 3) to determine the ranges of parameters for which the given solution can be reached. We divide the considered ranges of parameter values into a regular two dimensional grid. For each box in the grid we perform N trials of direct numerical integration. Each time we randomly pick the initial conditions from the accessible ranges and draw the values of investigated parameters from the ranges assigned to the grid box. For every trial we recognize the final attractor that is reached by the system. The data enable us to estimate the ranges of stability for each considered attractor. For that purpose, we detect all the trajectories that go to the investigated attractor and mark the parameters value (A and x) for which it has been approached. Then, we plot those values as a set of points that indicates the region of stability of the investigated solution. Apart from that, the collected data enable us to calculate the maps of basin stability in the two parameter space, as shown in Fig. 10c and f where we show the probability of occurrence of the considered solution.
In Fig. 10 we show the presentation scheme that we use. It enables to show all of the obtained results and ensures an easy comparison between the methods. In Fig. 10 we present the results yielded for some exemplary solution (1:2 oscillations of system given by equations 12) and explain the presentation layout. Arrows in Fig. 10 show how the data are interchanged between the panels.
We start with the path-following method. We begin with one parameter continuation performed to detect the bifurcation points in which the stability of the solution changes. These results are not presented but we use these points to  To simplify the comparison between the results obtained with different approaches, we combine the results from panels (a), (b) and (d) in a single plot as shown in panel (e) in Fig. 10. However, the maps of basin stability provide an additional knowledge about the attractor stability. Hence, for each solution we present three diagrams corresponding to panels (c), (e) and (f). To underline this fact these three panels are surrounded with a dashed line.

Ranges of Stability
Analyzing the ranges of stability we take as controlling parameters the amplitude A and the frequency x of the external excitation. In the considered ranges of parameters A 2 0; 0:0077 ½ ½ m and x 2 0; 60 ½ ½rad=s we found 12 different solutions. The upper pendulum always performs an oscillatory motion, while for the second pendulum we observe both oscillations and rotations with different locking ratios in respect to the frequency of excitation. Therefore, we name each solution basing on the behaviour of the second pendulum and ratio n:m which means that for m periods of excitation we observe n full oscillations or rotations of the second pendulum. We detected and further consider: 1:1, 1:2, 1:3, 1:4, 1:6, 1:8 oscillations, and 1:1, 1:2, 1:3, 5:5, 6:6, 2:3 rotations. The methodology of the investigation was described in Sect. 4.2 and in Fig. 11 we present the outcome from three considered methods using the described presentation scheme. The first six panels (af) refer to oscillatory solutions while in the last six panels we show the results corresponding to the rotary solutions.
We start the description from the oscillatory case. In panel (a) of Fig. 11 we show the results obtained for 1:1 oscillations of the second pendulum. In this solution the amplitude of oscillations is small, and because of that this solution can be reported as a semi-trivial one. Moreover, if one increase the damping coefficient in the pin joint of the second pendulum its motion would be not or barely observable. This solution is stable in almost the whole analysed range of parameter values except the narrow resonance tongue that occurs for x ¼ 20 ½rad=s. During experimental investigations four runs were performed, as we wanted to find the location of bifurcation lines predicted in the path-following analysis. In each run we slowly changed the value of x until we reach the moment when the system jumps into another solution. The position of the branching bifurcation lines have been detected by pathfollowing and confirmed using both experimental analysis and sample-based approach. The agreement between the results obtained using all three approaches is remarkably good. The second analysed solution is 1:2 oscillations which corresponds to the resonance tongue mentioned earlier (around x ¼ 20 ½rad=s). The results are shown in Fig. 11b. The path following analysis reveals that this solution loses its stability either in a saddle node bifurcation (for A 5 ½10 À3 m) or in a symmetry braking bifurcation (for A [ 5 ½10 À3 m) that is immediately followed by a cascade of period doubling bifurcations. The experimental analysis confirms that the position of the right-hand side border of the resonance tongue has been obtained with excellent precision. However, the detected range of stability differs from the results obtained via the numerical continuation. The difference is especially visible for A ¼ 4 ½10 À3 m (in horizontal direction) and for x ¼ 16 ½rad=s in vertical direction. In these cases the sample-based method enables a better precision, especially when detecting the left-hand side border of stability. We see that the density of points decreases significantly before the bifurcation lines are reached. Thus, using the sample-based approach we can better predict the range in which it is possible to obtain 1:2 oscillations in real life. For the third analysed solution 1:3 oscillations the results from both methods of analysis are in good agreement with the experimental data [see panel (c)]. Similarly, for 1:6 oscillations [panel (e)] and 1:8 oscillations [panel (f)] we observe a concurrence of the results from all the applied methods. The difference between the experimental and numerical results is clearly distinguishable only in two cases. The experimental investigation revealed that the 1:6 oscillations can be achieved also below the line detected using the path-following method. Also, the position of a period doubling bifurcation which destroys the stability of 1:8 oscillations is different from what we observe experimentally. We think that the difference between the numerical results and the experiment mainly comes from the dissipation modelling approach which strongly simplifies the real phenomena [14,49]. In panel (d) we show the results that correspond to 1:4 oscillations. For this solution, with the sample based approach can obtain the boundary of stability that better fits the experimental results then the bifurcation lines obtained via path-following method. It is especially visible for the left and the lower boundaries in the (A, x) plane.
Having considered all oscillatory solutions we can now focus on the rotary ones. The results are presented in the last six panels (g-l) of Fig. 11. It is important to notice that we were unable to perform experimental analysis of last  three rotary solutions, namely 2:3, 5:5, 6:6 rotations. We were only able to imitate 2:3 rotations but this solution is extremely susceptible to perturbations and the system jumps to another solution as soon as we try to change the amplitude or frequency of excitation, so it was impossible to reach the boundary of stability. In panel (g) we show the results obtained for 1:1 rotations of the second pendulum. This solution is stable in a noticeably wide range of excitation parameters and both numerical methods enable to predict the boundary of stability with good precision referring to the experimental result. Similar conclusions can be drawn from the analysis of the results obtained for 1:2 and 1:3 rotations that are shown in panels (h) and (i) respectively. Analysing the last three panels (j, k, l) of Fig. 11 we see that the location of stability boundaries predicted using the sample based procedure coincide with bifurcation lines calculated using the path-following method. To sum up, for all rotary solutions we see a good convergence between the results obtained using both numerical methods which reflects the gathered experimental data with good precision.
The results presented in Fig. 11 prove that for the considered system the path-following method and the samplebased approach offer similar accuracy. Also, when comparing the numerical results with the experimental data we observe good convergence. The advantage of the sample based analysis is that looking at the density of points we can also estimate the basin stability measure and probability density as described in the following sections. Such additional knowledge gives insight on the susceptibility to perturbations and one may expect the difficulties observed during experimental analysis. The important advantage of the sample-based approach is that during a large number of trials we are able to detect hidden attractors [27] or solutions with rather meager basins of attraction. In the investigated case we also found solutions that occurred only once for 172,500 trials, such as for example 2:5 oscillations or 3:15 rotations. Moreover, analysing Fig. 11, we see that with the sample-based approach we find the regions with the maximum density of points. This, in consequence enables to quantify the stability for different values of parameters using the basin stability measure or to prepare 2D probability density plots to find parameter values for which . To further investigate this topic for each solution we obtain maps indicating the changes of basin stability in the two-parameter space.

Basin Stability Maps
In the previous Section we show the results proving that the sample-based approach enables to predict the boundaries of stability with precision comparable to path-following method. Now, we show how we can utilize the results from the sample based method to obtain additional information about the system dynamics. In particular we will focus on the susceptibility to perturbations. The described algorithm is predisposed to such analysis because including random initial conditions we somehow investigate how the system reacts to sudden changes in working conditions. Basin stability measure is based on this concept and proved to be an efficient tool in a number of practical applications. Now we use data obtained from the aforementioned series of Bernoulli trials to analyse the changes of basin stability measure in the two parameter space. In Fig. 12 we demonstrate the results for all 12 considered solutions. It is important to notice that we use the grid that reflects the parameter ranges predefined in the calculation procedure and for each box we perform 500 trials. For all presented diagrams we hold the same colour scale for the basin stability that is given at the bottom of Fig. 12.
Panel (a) refers to the 1:1 oscillations which has the largest range of stability (see Fig. 11a). The plot reveals that also in the large range of parameters this solution has basin stability greater than 0.5 so it has the dominant volume of the basin of attraction. In these ranges, with probability grater than 0.5 we expect that after random perturbation the system will end up in the basin of attraction of 1:1 oscillations. Thus, this solution has relatively low susceptibility to perturbations. Panels (b-f) correspond to the remaining oscillatory solutions. We see that for most of them, basin stability calculated in the assumed grid never exceeds 0.3. Still, for 1:2 and 1:4 oscillations we uncover regions with higher basin stability. Especially for 1:2 oscillations we observe increase of basin stability measure for x 2 ½18:75; 21:25 ½rad=s and A [ 0:004 ½m, i.e., in the middle of the resonance tongue. For rotary solutions [panels (g-l)] we observe the similar structure of basin stability maps with maximum value of basin stability measure never exceeding. The only exception is the region around x 2 ½25; 35 ½rad=s and A [ 0:005 ½m where the probability of reaching 1 : 1 rotations is higher than 0.3.
Using the data obtained during the numerical procedure we can draw the 2D density plots reflecting continuous maps of basin stability (not using the lattice) by using interpolation. For that purpose we use the two-dimensional kernel density estimation in R software (R-project for statistical computing) [77]. With that functions we obtain plots presented in Fig. 13 that enables to detect where the given solution has the biggest relative volume of the basin of attraction. Thus, we are able to find the ranges in parameter space for which the solution is the most likely to appear and has low susceptibility to perturbations. The colour scale reflects the relative probability, so with such plots we cannot asses the value of basin stability measure but we can analyse the changes in the basins volume. Such plots are also convenient for estimating the boundary of stability.
Analysing the results presented in Figs. 12 and 13 we can divide the considered range of parameter values into three regions:    3. The region where, the 1:1 rotations has the dominant basin of attraction and we observe the coexistence of many stable solutions [see panels (g)].
The above results prove that with the sample based procedure we are able to obtain additional knowledge about the systems dynamics that can be presented using basin stability maps (as show in Fig. 12) or probability density plots (as presented in Fig. 13).

Investigation with Parameters Mismatch
The important advantage of sample based methods is the ease to include parameters mismatch as described in [15]. It has huge practical importance, because we often cannot obtain the exact values of the system parameters [69]. The parameter identification procedure is especially difficult when the model contains parameters whose values cannot be measured directly such as, for example, a value of the viscous damping coefficient in a pin joint. Apart from that, when modelling the mechanical and structural objects, we often simplify complex phenomena using simple models. This can decrease the accuracy of simulations and cause the divergence between numerical and experimental results. In this section we show that the extended basin stability method can be applied without the knowledge of actual parameter values and still maintain the high accuracy of the yielded results. Let us assume that we do not know precisely the values of some parameters of the system. In a classical approach, we have to set their values which may lead to wrong results. In the sample-based approach, instead of setting a value of the parameter, we can estimate the range to which this parameter belongs. Then, during a series of trials we draw the values of uncertain parameters. Before we apply that to our problem, we divide the parameters of our system (Eq. 12) into four specific groups basing on the ease of measurements: 1. Parameters that are easy to measure: m 1 , m 2 , l 1 , x 1 .
This group contains the parameters that can be determined easily with good precision, namely masses and lengths of the physical objects. During calculations we assume that we determined values of these parameters precisely and do not draw them. 2. Parameters that can be estimated with good precision: J 1 , J 2 , d 1 , d 2 , k. The methods to measure these parameters values are not straight forward. In this group there are also parameters whose values are given by manufacturer with certain precision -as for example the stiffness of the spring. We assume that the error for the parameters from that group is AE5% and for each trial we draw the values of parameters from a certain range.
3. Parameters that are difficult to estimate: c 1 , c 2 . In this group we put parameters that require complex procedures or approach to infer their values and parameters that refers to phenomena that are difficult to model, such as energy dissipation. Values of parameters from this group are drawn from the range [80%, 120%] of the estimated value. 4. Parameters whose special influence we investigate: A; x. The aim of the study is to analyse the influence of these parameters on the stability of solutions. For these parameters we do not estimate the error. Instead, we draw values of these parameters from a ''grid'' as described in the Methods section.
Apart from the above, we want to maintain the value of the natural frequency of the second pendulum which can be easily measured using a simple stopwatch. For that purpose, we draw the inertia of the second pendulum and recalculate the position of its center of gravity to preserve the natural frequency.
With the above assumptions we repeat all the calculations. In Fig. 14 we present the comparison between the results obtained with fixed values of parameters and when drawing the parameter values. In the first two columns of Fig. 14 we show the results that refer to 1 : 2 oscillations and last two columns refers to 1 : 3 oscillations. In panels (a) and (e) we show the results from the sample based method with fixed parameters and in (b) and (f) the points obtained when drawing the parameter values. Similarly, we compare the maps of the basin stability for fixed parameters (c,g) and obtained with the assumed mismatch (d,h). Results that refer to the remaining 10 solutions are simmilar, hence we do not present them.
The difference in the results from the extended basin stability analysis are barely visible despite the fact, that the positions of the bifurcation lines change when we modify the parameter values within the assumed ranges. This shows that the sample-based approach can be applied even without time consuming detailed measurements of the system's parameters and still ensures sensible results. Instead of precisely measure the parameter values, we only have to asses the precision for each of the determined parameter.

Example of Practical Application
In this Section we apply the above described procedure to analyze the ranges of stability of different working regimes of the yoke-bell-clapper system described in Sect. 3.3. The system is multistable and there are two influencing parameters that can be changed to optimize working conditions while not affecting the generated sound, namely the amplitude of exciting torque T max and geometry of yoke given by parameter l r described in Fig. 7. Altering these two parameters we are able to change the response of the system significantly as described in [10]. Fig. 15a is a fragment of the plot presented and described in detail in [10]. It shows which type of proper working regime will be achieved assuming zero initial conditions (u 1 Fig. 15 Two parameter ringing schemes diagram showing the attractor reached starting from zero initial conditions. In the table we briefly characterize the possible solutions. For more details please see [10] u 2 ¼ 0:0 [rad], _ u 1 ¼ 0:0 [rad/s], _ u 2 ¼ 0:0 [rad/s]). Still, the results presented in Fig. 8 prove that due to multistability, for fixed parameter values we are able to obtain different working regimes by applying non zero initial conditions or using alternative starting procedure. This may be beneficial in practice as given bell can work differently depending on occasion. To enable such property we should know the precise boundaries of stability of solutions with respect to T max and l r as these parameters determine the needed power of linear motor propulsion and geometry of a yoke. This can be achieved using the sample based approach described above.
We consider the following ranges of parameters values: T max 2 ½100; 300 [Nm] and l r 2 ½À 0:3; 0:2 [m] and run 400,000 simulations each time drawing values of T max , l r and initial conditions from the ranges given in Eq. 11. We are especially interested in solutions that correspond to proper working regimes or have significantly large range of stability. Hence, we will only present ranges of stability for solutions given in the table on the right hand side of Fig. 15. The results are presented in Fig. 16. Each panel corresponds to different solution and presents the estimation of probability density of reaching it. Moreover, with white dashed lines we mark the regions in which the solution is reached from zero initial conditions.
In panel (a) we show results for solution 1 (period 1 attractor with no impacts). We see that in the whole stability region the probability of reaching this attractor is similar and stability boundary overlaps with white dashed line which depicts the range of parameters for which solution 1 is reached from zero initial conditions. Such agreement is observed only for solution 1. Panel (b) refers to asymmetric falling clapper regime (solution 6 or 7) and we see that at the right hand side boundary of stability the probability of reaching this solution is gradually decreasing while it drops suddenly on the left boundary. There is a significant range of parameter values for which this solution coexists with symmetric falling clapper (solution 8) as shown in panel (c). In this multistable region solution 8 has dominant volume of basin of attraction. Panels (d) and (e) present the results obtained for solutions 9 and 12 respectively. In both cases the probability of reaching considered solution does not fluctuate strongly in the stability range. The overall probability of reaching a non periodic solution is shown in panel (f). The stability boundary more less overlaps with the white dashed lines which depict the range of parameters for which starting from zero initial conditions we do not reach proper working regime. However, between the white lines the probability of reaching non-periodic solution strongly varies. This means that in that range there are some stable attractors which was also detected in one parameter analysis (see Sect. 3.3) as shown in Fig. 8. We do not prepare additional plots for these solutions as they are not significant from Fig. 16 Two parameters probability density estimations for 6 investigated solutions. White dashed lines mark the region in which given solution is reached from zero initial conditions. Below is the color scale of the relative probability of reaching each solution practical point of view. It is important to notice, that even if there is only one stable attractor we can observe some variability of relative probability density. This is the effect of the sample based approach in which we do not ensure same density of trials in the analyzed parameters region. However, with sufficient number of trails we can neglect it. Summarizing, using the sample based approach we can prepare plots that precisely indicates the ranges of stability in two-parameters space and enables to determine parameter values that ensure biggest relative probability of reaching presumed solution.

Characterization of Phase Space Structure
Another important application of the sample based methods refers to the analysis of phase space organization which is largely unexplored especially for multi-degrees-of-freedom systems due to the computational challenge requested to build basins of attraction. Sample based analysis may help to overcome this problem. We propose to estimate the probability of reaching the attractors using a reasonable number of trials with random initial conditions. Then, we use the obtained results to investigate how this probability depends on particular generalized coordinate or a pair of coordinates. In the analysis we include all initial conditions, and reflect their influence on the final behaviour. The method allows to obtain information about the basins compactness and reveals particular features of the phase space topology. To analyse the trustworthiness of the proposed approach we compare the results with the classical computation of basins of attraction performed in the full range of initial conditions. Finally, we present example of practical application basing on the previously described hybrid dynamical model of the yoke-bell-clapper system.

Methodology and Exemplary Paradigmatic System
In this section we present the methodology basing on the analysis of the phase space structure of externally forced Duffing oscillator with attached pendulum. It is the model presented in Fig. 1a and  We start with explanation of the methodology and description of the presentation scheme for the obtained results. Then, we focus on the advantages of the sample based approach.
In the system we observe two stable solutions. The solution I is period-1 oscillation (it has the same period of the excitation) of mass M with pendulum remaining in the hanging down position. The second possible orbit (solution II) is period-2 solution (its period is double than that of the excitation) in which both mass M and pendulum performs oscillatory motion. At first, we consider the probability of reaching solution I from random initial conditions. We perform 2:5 Â 10 6 trials of direct numerical integration each time drawing the initial conditions from the following ranges x 2 À 2; 2 ½ , u 2 Àp; p ½ , _ x 2 À 2; 2 ½ , _ u 2 À4; 4 ½ . These ranges ensure that all stable solutions can be reached and refer to the practically accessible initial conditions that could be implied in the real world realization of the investigated system. The overall probability of reaching solution I is 54.3% which corresponds to basin stability measure of this attractor. Hence, the second solution has probability of reaching equal to 45.7%. Now, we want to investigate how this value change for different regions of the phase space.

One-Coordinate Histograms
Firstly, we consider how each particular coordinate influences the probability of reaching solution I. For this purpose we calculate histograms for each generalized coordinate separately. Each time we use 42 equal steps receiving on average 600,000 trials for each subset. In Fig. 17a-d we show how the probability of reaching solution I evolves with respect to coordinates x, _ x, u and _ u respectively. Each dot on the histograms corresponds to the series of trials in which all initial conditions are random but for one initial condition given on vertical axis we narrow down the range of drawing its value to 1 42 of the whole accessible range. Hence, we see 42 discrete ranges which show the changes of probability as a function of one selected initial condition. Panel (a) shows how the probability of reaching solution I depends on the initial value of x. We see that there is no dominant trend and the value is scattered around the overall average probability. Hence, we expect that the initial value of x can affect the structure of basins of attraction, but not its overall volume. This is confirmed by the sample 3D projections of the full basins of attraction, obtained for x 0 ¼ À 2:0, x 0 ¼ 0:0 and x 0 ¼ 1:0 presented in panels (a.1, a.2, a.3) respectively. Moreover, in our previous paper [8] we present supplementary material video1 which revels the evolution of 3D projection of the phase space for x 2 À 2; 2 ½ . Panel (b) refers to the influence of initial condition _ x which is much less scattered.Thus, changes in initial value of _ x should have minor influence on the basins structure and compactness. It is also visible in the sample 3D plots (panels (b.1, b.2, b.3)) calculated for _ x 0 ¼ À 2:0, _ x 0 ¼ 0:0 and _ x 0 ¼ 1:0 respectively. Panels (c) and (d) correspond to histograms calculated for u and _ u respectively. In both plots we see similar shape of histograms. The probability is scattered on both ends of the considered ranges of initial conditions while in the middle there is clearly visible trend with the maximum around zero. Hence, we expect that around zero value the basins have large volume (maximum along the histogram) and more compact structure (clearly visible trend). Contrary, on both ends of the considered range we expect basins with smaller volume and more complex structure (lower probability and large dispersion in the histogram). These predictions are confirmed both on the sample 3D projections-panels (c.1, c.2, c.3, d.1, d.2, d.3). Animations showing full evolutions of the 3D basins of attraction can be found in the supplementary material to our recent paper [8]. The presented results show that one coordinate histograms may indicate the changes in the structure and volume of basins of attraction. We see that large dispersion along the histograms indicate fractal structure of basins while clearly visible trends refer to more compact shapes of basins. Apart from that, the proposed histogram plots can be used to detect the ranges of initial conditions for which we observe minimum or maximum volume of the basin of attraction.

Two-Coordinate Density Plots
In this section we consider 2D density plots showing how the probability of reaching solution I changes with respect to two particular coordinates. To obtain the density plot we perform large number of trials with random initial conditions. Then, we clusterize trials with respect to two initial conditions and for each group we calculate the probability of reaching solution I. By that, we obtain drawings that are  Fig. 17 Histograms showing how the probability of reaching solution I changes with respect to four generalized coordinates (a-d) supplemented by examples of 3D plots of basins of attraction somehow squeezed reflections of phase space that can be interpreted similarly as 1D histograms but brings more information. They are especially meaningful for mechanical systems because they are able to reflect the influence of the initial conditions of single degree of freedom (position and velocity). If we want to have the information how other pairs of initial conditions influence the density we can easily plot such diagrams. Moreover 2D plot is the limit of clear visualization method that enables convenient interpretation. The above reasons make the proposed density plots an efficient tool for the analysis of phase space structure.
In Fig. 18 we show 2D density plots obtained for the investigated system. In panels (a, b) we consider initial conditions that refer to separated single degree of freedom. In panel (a) we show how the probability is changing with varying x and _ x. We see that for _ x 2 À1:6; 1:4 ð Þ the probability of reaching assumed solution is much larger than outside this range. The same conclusion can be drawn based on 1D diagram. The advantage of 2D plot is visible in the middle of the plot, where we clearly see that there is low high  Fig. 18. The plot is significantly different from the first one. We observe a high probability in the center of plot and rapid decrease of probability outside this range. The highest probability is marked with dashed line ellipse.
To present the advantage of 2D plots in panels (c) and (d) we show basins of attraction calculated for fixes initial conditions of the Duffing oscillator. The green and yellow colours refer to basins of attraction of solutions I and II respectively. Panel (c) corresponds to the point from panel (a) with minimum probability (x ¼ À 0:24, _ x ¼ 0:74) and panel (d) refers to the point with maximum probability (x ¼ À0:18, _ x ¼ 1:70). Comparing panels (c) and (d) we see that indeed the amount of yellow colour that refers to solution II basins differs strongly. Moreover, on both plots we see unchanged range of high probability of reaching solution I marked with the dashed line ellipse [the same area as in panel (b)]. The above consideration shows that 2D density plots can be effectively utilized to investigate the structure of the phase space in multidimensional systems and provides additional knowledge over traditional presentation methods.

Practical Application
In this Section we one more time refer to the experimentally validated model of the yoke-bell-clapper system described in Sect. 3.3. We analyse the structure of the phase space for the following set of parameter values: T max ¼ 230 [Nm], l r ¼ À 0:075 [m]. Analyzing Figs. 8 and 16 we see that two stable attractors coexist for that case. Both of them are falling clapper hence the collisions between the bell and the clapper occur when they perform an anti-phase motion. The first attractor corresponds to asymmetric falling clapper with one impact per period of motion (solution 6 or 7) and the second one to symmetric falling clapper with two impacts per period-one on each side of the bell (solution 8). Assuming the ranges of applicable initial conditions given by Eq. 11 we analyze the structure of the phase space.
Results are presented in Fig. 19 where the first three panels (a, b, c) refer to the relative probability of reaching solution 6 or 7 and panels (d, e, f) placed in the lower row reflect the results for solution 8. In the first column we present the influence of initial angular position and velocity of the bell while in the second column of the clapper. We see that to obtain the symmetric solution with two impacts per period (solution 8) we have to imply some initial deflection of the bell and the clapper. It is also much easier to obtain precise initial deflection than angular velocity. Hence, for practical application, in the last column we plot how the relative probability depends on the initial deflections of the bell and clapper. The white triangles reflect the geometrical boundaries of accessible initial states of the system (as aforementioned the initial angular potions of bell and clapper must fulfil the relation u 1 0 À u 2 0 \a, see Eq. 11). Panels (c, f) are of practical importance as they indicate that it is impossible to reach symmetric solution without implying some specific initial state of the system or sophisticated starting procedure.

Conclusions
Recent rapid increase in computational power and ensuing development of parallel computing create immense opportunities for novel sample based methods of dynamical analysis. At first glance, these methods seem complex and sophisticated as they need powerful numerical apparatus. Still, they enable to reproduce the inherent uncertainty of perturbations and have numerous advantages which make them appropriate especially for analysis of multistable systems. The first is the fact that we are able to detect all possible solutions also hidden and rare attractors. This feature is extremely important when analyzing systems with unknown number of possible stable solutions. With sample based methods we are able to detect all solutions with basin stability greater than some threshold that depends on the number of trials. We can create basin stability maps or twodimensional probability density plots to show ranges of stability and find parameter values that ensure the highest possible basin stability measure. The presented results prove that sample based algorithm ensures precision similar to classical methods (direct numerical integration, path-following) and the obtained results reflect the experimental data.
Secondly, we are able to investigate the ranges of stability in multiple parameters space including parameters mismatch and variation of parameter values. We can take into account the fact that values of parameters are measured or estimated with some finite precision and also that they can slightly change during normal operation. We can extend the analysis by adding model imperfections or different types of noise. The obtained results enable to assess the robustness of solutions and estimate the risk of unwanted behavior. Moreover, with sample based methods we can concurrently investigate dynamics in multiple parameters space.
The last revealed application of the sample based approach is analysis of the phase space structure of mulit-DOF dynamical systems. The proposed method enables to investigate the influence of particular initial condition. The presented results show that the method enables to obtain additional knowledge of the basin compactness and reveal regions in the phase space that have some specific features. Crucial part of the approach is proposed presentation scheme that enables to overcome some limitations of the currently know methods. Probability histograms are useful to analyze single initial condition and detect ranges with maximum or minimum probability of reaching the solution of interest. Two dimensional density plots correspond to squeezed reflections of phase space and reveal the influence of given degree-of-freedom or a pair of initial conditions. These plots are obtained assuming random initial conditions, hence they uncover the overall influence of considered the initial condition or conditions.
The above advantages prove that sample based methods are efficient tool that can be effectively utilized to investigate the dynamics of wide range of dynamical systems and provide additional knowledge over traditional methods.

Compliance with Ethical Standards
Funding This work has been supported by Lodz University of Technology own Scholarship Fund (PB) and by Stipend for Young Outstanding Scientists from Ministry of Science and Higher Education of Poland (PB).