Basin stability approach for quantifying responses of multistable systems with parameters mismatch

In this paper we propose a new method to detect and classify coexisting solutions in nonlinear systems. We focus on mechanical and structural systems where we usually avoid multistability for safety and reliability. We want to be sure that in the given range of parameters and initial conditions the expected solution is the only possible or at least has dominant basin of attraction. We propose an algorithm to estimate the probability of reaching the solution in given (accessible) ranges of initial conditions and parameters. We use a modified method of basin stability (Menck et. al., Nature Physics, 9(2) 2013). In our investigation we examine three different systems: a Duffing oscillator with a tuned mass absorber, a bilinear impacting oscillator and a beam with attached rotating pendula. We present the results that prove the usefulness of the proposed algorithm and highlight its strengths in comparison with classical analysis of nonlinear systems (analytical solutions, path-following, basin of attraction ect.). We show that with relatively small computational effort (comparing to classical analysis) we can predict the behaviour of the system and select the ranges in parameter's space where the system behaves in a presumed way. The method can be used in all types of nonlinear complex systems.


Introduction
In mechanical and structural systems the knowledge of all possible solutions is crucial for safety and reliability. In devices modelled by linear ordinary differential equations we can predict the existing solutions using analytical methods [29,25]. However, in case of complex, nonlinear systems analytical methods do not give the full view of system's dynamics [31,20,3,24]. Due to nonlinearity, for the same set of parameters more then one stable solution may exist [16,11,32,18,7,23,22]. This phenomenon is called multistability and has been widely investigated in all types of dynamical systems (mechanical, electrical, biological, neurobiological, climate and many more). The number of coexisting solutions strongly depends on the type of nonlinearity, the number of degrees of freedom and the type of coupling between the subsystems. Hence, usually the number of solutions vary strongly when values of system's parameters changes.
As an example, we point out the classical tuned mass absorber [2,9,1,17,30,21,12,6]. This device is well known and widely used to absorb energy and mitigate unwanted vibrations. However, the best damping ability is achieved in the neighbourhood of the multistability zone [7]. Among all coexisting solutions only one mitigates oscillations effectively. Other solutions may even amplify an amplitude of the base system. So, it is clear that only by analyzing all possible solutions we can make the device robust.
Similarly, in systems with impacts one solution can ensure correct operation of a machine, while others may lead to damage or destruction [5,4,14,28,15]. The same phenomena is present in multi-degree of freedom systems where interactions between modes and internal resonances play an important role [8,19,10,26].
Practically, in nonlinear dynamical systems with more then one degree of freedom it is impossible to find all existing solutions without huge effort and using classical methods of analytical and numerical investigation (path-following, numerical integration, basins of attractions), especially in cases when we analyse a wider range of system's parameters and we cannot precisely predict the initial conditions. Moreover, solutions obtained by integration may have meager basins of attraction and it could be hard or even impossible to achieve them in reality. That is why we propose here a new method basing on the idea of basin stability [23]. The classical basin stability method is based on the idea of Bernoulli trials, i.e., equations of system's motion are integrated N times for randomly chosen initial conditions (in each trial they are different). Analyzing the results we asses the stability of each solution. If there exist only one solution the result of all trials is the same. But, if more attractors coexist we can estimate the probability of their occurrence for a chosen set of initial conditions. In mechanical and structural systems we want to be sure that a presumed solution is stable and has the dominant basin of attraction in a given range of system's parameters. Therefore, we build up a basin stability method by drawing values of system's parameters. We take into account the fact that values of parameters are measured or estimated with some finite precision and also that they can slightly vary during normal operation.
The paper is organized as follows. In Section 2 we introduce simple models which we use to demonstrate the main idea of our approach. In the next section we present and describe the proposed method. Section 4 includes numerical examples for systems described in Section 2. Finally, in Section 5 our conclusions are given.

Model of systems
In this section we present systems that we use to present our method. Two models are taken from our previous papers [7,13] and the third one was described by Pavlovskaia et. al. [27]. We deliberately picked models whose dynamics is well described because we can easily evaluate the correctness and efficiency of the method we propose.

Tuned mass absorber coupled to a Duffing oscillator
The first example is a system with a Duffing oscillator and a tuned mass absorber. It was investigated in [7] and is shown in Figure 1. The main body consists of mass M fixed to the ground with nonlinear spring (hardening characteristic k 1 + k 2 y 2 ) and a viscous damper (damping coefficient c 1 ). The main mass is forced externally by a harmonic excitation with amplitude F and frequency ω. The absorber is modelled as a mathematical pendulum with length l and mass m. A small viscous damping is present in the pivot of the pendulum. The equations of the system's motion are derived in [7], hence we do not present their dimension form. Based on the following transformation of coordinates and parameters we reach the dimensionless form: . The dimensionless equations are as follows: where µ is the frequency of the external forcing and we consider it as controlling parameter. The dimensionless parameters have the following values: f = 0.5, a = 0.091, b = 3.33, α = 0.031, d 1 = 0.132 and d 2 = 0.02. Both subsystems (Duffing oscillator and the pendulum) have a linear resonance for µ = 1.0.

System with impacts
As the next example we analyse a system with impacts [27]. It is shown in Figure 2 and consists of mass M suspended by a linear spring with stiffness k 1 and a viscous damper with the damping coefficient c to harmonically moving frame. The frame oscillates with amplitude A and frequency Ω. When amplitude of mass M motion reaches the value g, we observe soft impacts (spring k 2 is much stiffer than spring k 1 ).
The dimensionless equation of motion is as follow (for derivation see [27]) : where x = y y0 is the dimensionless vertical displacement of mass M , τ = ω n t is the dimensionless time, ω n = k1 M , β = k2 k1 the stiffness ratio, e = g y0 the dimensionless gap between equilibrium of mass M and the stop suspended on the spring k 2 , a = A y0 and ω = Ω ωn are dimensionless amplitude and frequency of excitation, ξ = c 2mωn is the damping ratio, y 0 = 1.0[mm] and H(·) the Heaviside function. In our calculations we take the following values of system's parameters: a = 0.7, ξ = 0.01, β = 29, e = 1.26. As a controlling parameter we use the frequency of excitation ω.

Beam with suspended rotating pendula
The last considered system consists of a beam which can move in the horizontal direction and n rotating pendula. The beam has the mass M and supports n rotating, excited pendula. Each pendulum has the   same length l and masses m i (i = 1, 2, . . . , n). We show the system in Figure 3 [13]. The rotation of the i-th pendulum is given by the variable ϕ i and its motion is damped by the viscous friction described by the damping coefficient c ϕ . The forces of inertia of each pendulum acts on the beam causing its motion in the horizontal direction (described by the coordinate x). The beam is considered as a rigid body, so we do not consider the elastic waves along it. We describe the phenomena which take place far below the resonances for longitudinal oscillations of the beam. The beam is connected to a stationary base by a light spring with the stiffness coefficient k x and viscous damper with a damping coefficient c x . The pendula are excited by external torques proportional to their velocities: N 0 −φ i N 1 , where N 0 and N 1 are constants. If no other external forces act on the pendulum, it rotates with the constant velocity ω = N 0 /N 1 . If the system is in a gravitational field (where g = 9.81 [m/s 2 ] is the acceleration due to gravity), the weight of the pendulum causes the unevenness of its rotation velocity, i.e., the pendulum slows down when the centre of its mass goes up and accelerates when the centre of its mass goes down. The system is described by the following set of dimensionless equations: In our investigation we analyze two cases: a system with two pendula (where n = 2 and i = 1, 2) and with 20 pendula (n = 20 i = 1, 2, ..., n).. The values of the parameters are as follows: m i = 2.00 n , l = 0.25, c ϕ = 0.02 n , N 0 = 5.00, N 1 = 0.50, M = 6.00, g = 9.81, c x = ln(1.5) parameter. The derivation of the system's equations can be found in [13]. We present the transformation to a dimensionless form in Appendix A.

Methodology
In [23] Authors present a "basin stability" method which let us estimate the stability and number of solutions for given values of system parameters. The idea behind basin stability is simple, but it is a powerful tool to assess the size of complex basins of attraction in multidimensional systems. For fixed values of system's parameters, N sets of random initial conditions are taken. For each set we check the type of final attractor. Based on this we calculate the chance to reach a given solution and determine the distribution of the probability for all coexisting solutions. This gives us information about the number of stable solutions and the sizes of their basins of attraction.
We consider the dynamical systemẋ = f (x, ω), where x ∈ R n and ω ∈ R is the system's parameter. Let B ⊂ R n be a set of all possible initial conditions and C ⊂ R a set of accessible values of system's parameter. Let us assume that an attractor A exists for ω ∈ C A ⊂ C and has a basin of attraction β(A). Assuming random initial conditions the probability that the system will reach attractor A is given by p (A). If this probability is equal to p (A) = 1.0 this means that the considered solution is the only one in the taken range of initial conditions and given values of parameters. Otherwise other attractors coexist. The initial conditions of the system are random from set B A ⊂ B. We can consider two possible ways to select this set.
I The first ensures that set the B A includes values of initial conditions leading to all possible solutions. This approach is appropriate if we want to get a general overview of the system's dynamics.
II In the second approach we use a narrowed set of initial conditions that corresponds to practically accessible initial states.. . .
In our method we chose the second approach because it let us take into account constrains imposed on the system and because in engineering we usually know or expect the initial state of the system with some finite precision.
In the classical approach of Menck et. al. [23] the values of system's parameters are fixed and do not change during calculations. The novelty of our method is that we not only draw initial conditions but also values of some selected parameters of the system. We assume that the initial conditions and some of the system's parameters are chosen randomly. Then using N trials of numerical simulations we estimate the probability that the system will reach a given attractor A (p (A)). The idea is to take into consideration the fact that the 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. Therefore drawing valus of parameters we can describe how a mismatch in their values influences the dynamics of the system and estimate the risk of failure. In many practical applications one is interested in reaching only one presumed solution A, and the precise description of other coexisting attractors is not necessary. We usually want to know the probability of reaching the expected solution p (A) and the chance that the system behave differently. If p (A) is sufficiently large, we can treat the other attractors as an element of failure risk.
In our approach we perform the following steps: I We pick values of system's parameters from the set C A ⊂ C.
II We select the set C A so that it consists of all practically accessible values of system's parameters ω . This let us ensure that a given solution indeed exists in a practically accessible range (taking into account the mismatch in parameters).
III We subdivide the set C A in to m = 1, 2, . . . M equally spaced subsets. The subsets C m A do not overlap and the relation m=1...M C m A = C A is always fulfilled.
IV 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.
V After a suficient number of trials we calculate the probability of reaching a presumed solution or solutions.
VI Finally we describe the relation between the value of the system's parameter and the "basin stability" of reachable solutions.. . .
In our calculations for each range of parameter values (subset C m A ) we draw from N = 100 up to N = 1000 sets of initial conditions and parameter. The value of N strongly depends on the complexity of the analysed system. Also the computation time for a single trial should be adjusted for each system independently such that it can reach the final attractor. In general, we recommend that in most cases N should be at least 100.

Tuned mass absorber coupled to a Duffing oscillator
At the beginning we want to recall the results we present in our previous paper [7]. As a a summary we show Figure 4 with a two dimensional bifurcation diagram obtained by the path-following method. It gives bifurcations for varying amplitude f and frequency µ of the external excitation (see Eq. 1). Lines shown in the plot correspond to different types of bifurcations (period doubling, symmetry breaking, Neimark-Sacker and resonance tongues). 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 [7]). We mark areas where we observe the existence of one solution (black colour), or the coexistence of two (grey) and three (hatched area) stable solutions. The remaining part of the diagram (white area) corresponds to situations where there are four or more solutions. Additionally, by white colour we also mark areas where only the Duffing system is oscillating in 1:1 resonance with the frequency of excitation and the pendulum is in a stable equilibrium position, i.e., HDP (hanging down pendulum) state. In this case the dynamics of the system is reduced to the oscillations of summary mass (M + m).
The detailed analysis of system 1 is time consuming and creation of Figure 4 was preceded by complex analysis done with large computational effort. Additionally, the obtained results give us no information about the size of the basins of attraction of each solution -which practically means that some of the solutions may occur only very rarely in the real system (i.e. due to not accessible initial conditions). Nevertheless, such analysis gives us an in-depth knowledge about the bifurcation structure of the system. As we can see, the range where less then three solutions exist is rather small, especially for µ < 2.0. To illustrate our method of analysis, we focus on three solutions: 2 : 1 oscillating resonance, HDP and 1 : 1 rotating resonance assuming that only they have some practical meaning.   (1) in the plane (f, µ) showing periodic oscillations and rotations of the pendulum. Black colour indicates one attractor, grey colour shows two coexisting attractors (the same as for black but with a coexisting stable steady state of the pendulum). In the hatched area we observe the coexistence of stable rotations and a stable steady state of the pendulum. A detailed analysis is presented in [7].
To show our results obtained with integration, we compute bifurcation diagrams for f = 0.5 in the range µ ∈ [0.1, 3.0] (see Figure 5). In Figure 5(a) we increase µ from 0.1 to 3.0 and in Figure 5(b) we decrease µ from 3.0 to 0.1. As the initial conditions we take the equilibrium position (x 0 =ẋ 0 = 0.0 and γ 0 =γ 0 = 0.0). In both panels we plot the amplitude of the pendulum γ. Ranges where the diagrams differ we mark by grey rectangles. It is easy to see that there are two dominating solutions: HDP and 2 : 1 internal resonance. Near µ = 1.0 we observe a narrow range of 1 : 1 and 9 : 9 resonances and chaotic motion (for details see Figure 6 in [7]). Based on previous results we know that we detected all solutions existing in the considered range, however we do not have information about the size of their basins of attraction and coexistence. Hence the analysis with the proposed method should give us new important information about the system's dynamics. Contrary to the bifurcation diagram obtained by path-following in Figure 5, we do not observe rotating solutions (the other set of initial conditions should be taken). In Figure 6 we show the probability of reaching the three aforementioned solutions obtained using the proposed method. The initial conditions are random numbers drawn from the following ranges: x 0 ∈ [−2, 2], x 0 ∈ [−2, 2], γ 0 ∈ [−π, π] andγ 0 ∈ [−2.0, 2.0] (ranges there selected basing on the results from [7]). The frequency of excitation is within a range µ ∈ [0, 3.0] (Figure 6(a,c) ), then we refine it to µ ∈ [1.25, 2.75] (Figure 6(b,d) ). In both cases we take 15 equally spaced subsets of µ and in each subset we calculate the probability of reaching a given solution. For each subset we calculate 100 trials each time drawing initial conditions of the system and a value of µ from the appropriate range. Then we plot the dot in the middle of the subset which indicate the probability of reaching a given solution in each considered range. Lines that connect the dots are shown just to ephasize the tendency. For each range we take N = 1000 because we want to estimate the probability of a solution with small a basin of stability (1:1 rotating periodic solution).
As we can see in Figure 4, the 2 : 1 resonance solution exists in the area marked by black colour around µ = 2.0 and coexists with HDP in the neighbouring grey zone. In Figure 6 we mark the probability of reaching the 2 : 1 resonance using blue dots. As we expected, for µ < 1.4 and µ > 2.2 the solution does not exist. In the range µ ∈ [1.4, 2.2] the maximum value of probability p(2 : 1) = 0.971 is reached in the subset µ ∈ [1.8, 2.0] and outside that range the probability decreases. To check if we can reach p(2 : 1) = 1.0, we decrease the range of parameter's values to µ ∈ [1.25, 2.75] and the size of subset to ∆µ = 0.1 (we still have 15 equally spaced subsets). The results are shown in Figure 6(b) similarly in blue colour. In the range µ ∈ [1.95, 2.05] the probability p(2 : 1) is equal to unity and in the range µ ∈ [1.85, 1.95] it is slightly smaller p(2 : 1) = 0.992. Hence, for both subsets we can be nearly sure that the system reaches the 2 : 1 solution. This gives us indication of how precise we have to set the value of µ to be sure that the system will behave in a presumed way.
A similar analysis is performed for HDP. The values of probability is indicated by the red dots. As one can see for µ < 0.8, µ ∈ [1.2, 1.4] and µ ∈ [2.6, 2.8], the HDP is the only existing solution. The rapid decrease close to µ ≈ 1.0 indicates the 1 : 1 resonance and the presence of other coexisting solutions in this range (see [7]). In the range µ ∈ [1.2, 1.4] the probability p(HDP) = 1.0 which corresponds to a border between solutions born from 1 : 1 and 2 : 1 resonance. Hence, up to µ = 2.0 the probability of the HDP solution is a mirror refection of p(2 : 1). The same tendency is observed in the narrowed range as presented in Figure 6(b). Finally, for µ > 2.0 the third considered solution comes in and we start to observe an increase of probability of the rotating solution S(µ, HDP) as shown in Figure 6(c). However, the chance of reaching the rotating solution remains small and never exceeds p(1 : 1) = 8 × 10 −3 . We also plot the probability of reaching the rotating solution in the narrower range of µ in Figure 6(d). The probability is similar to the one presented in Figure 6(c) -it is low and does not exceed p(1 : 1) = 8 × 10 −3 . Note that the results presented in Figure 6(a,b) and Figure 6(c,d)

System with impacts
In this subsection we present our analysis of different periodic solutions in the 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 and their classification is necessary. In Figure 7 we show two bifurcation diagrams with ω as controlling parameter. Both of them start with initial conditions x 0 = 0.0 andẋ 0 = 0.0. In panel (a) we increase ω from 0.801 to 0.8075; while in panel (b) we decrease ω in the same range. We select the range of ω basing on the results presented in [27]. 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  (2). For subplot (a) the value of the bifurcation parameter ω was increased while for subplot (b) we decreased the value of ω. Grey rectangles mark the range of the bifurcation parameter ω for which different attractors coexist. Further analysis can be found in [27].
range ω ∈ [0.8033, 0.8044] solutions with period-3 and -2 are present, while in the range ω ∈ [0.8068, 0.8075] we detected solutions with period-2 and -5. As presented in [27] 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. We observe periodic solutions with higher periods in the narrow range of ω but the probability that they will occur is very small and we can neglect them. All non-periodic solutions are chaotic (quasiperiodic solutions are not present in this system). The results of our calculations are shown in Figure 8(a,b). We take initial conditions from the following ranges x 0 ∈ [−2, 2],ẋ 0 ∈ [−2, 2]. The controlling parameter ω is changed from 0.801 to 0.8075 with step ∆ω = 0.0005 in Figure 8(a) and from 0.806 to 0.8075 with the step ∆ω = 0.0001 in Figure 8(b) (in each subrange of excitation's frequency we pick the exact value of ω randomly from this subset). The probability of periodic solutions is plotted by lines with different colours and markers. We detect the following solutions: period-1, -2, -3, -5 (two different attractors with large and small amplitude), -6 and -8. The dot lines indicate the sum of all periodic solutions' probability (also with period higher then eight). Hence, when its value is below 1, chaotic solution exist. 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 ω. Hence, in the sense of basin stability we can say that stability of solutions rely upon the ω value. In Figure 8(a) the probability of a single solution is always smaller than one. Nevertheless, we observe two dominant solutions: period-5 with large amplitude in the first half of the considered ω 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 ω ≈ 0.80675. To check if we can achieve even higher probability we analyse a narrower range of ω and decrease the step (from ∆ω = 0.0005 to ∆ω = 0.0001). In Figure 8(b) we see that in range ω ∈ [0.8069, 0.807] the probability of reaching the period-2 solution is equal to 1. Hence, in the sense of basin stability it is the only stable solution. Also in the range ω ∈ [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.
Other periodic solutions presented in Figure 8

Beam with suspended rotating pendula
The third considered system consists of a beam that can move horizontally with two (n = 2) or twenty (n = 20) pendula suspended on it. As a control parameter we use k x which describes the stiffness of the beam's support. For the considered range of k x ∈ [100, 5000] two stable periodic attractors exist in that system. One corresponds to complete synchronization of the rotating pendula. The second one is called anti-phase synchronization and refers to the state when the pendula rotate in the same direction but are shifted in phase by π.
In Figure 9 we show four bifurcation diagrams with k x as the controlling parameter and a Pioncare map of rotational speed of the pendula. The subplots (a,b) refer to the system with two pendula (n = 2). We start with zero initial conditions: x 0 = 0.0,ẋ 0 = 0.0, ϕ 10 = 0.0,φ 10 = 0.0, ϕ 20 = 0.0,φ 20 = 0.0 and take k x ∈ [100, 5000]. The parameter k x is increasing in subplot (a) and decreasing in (b). We see that in the range marked by grey rectangle both complete and anti-phase synchronization coexist. In subplots (c,d) we present results for twenty pendula (n = 20). We start the integration from initial conditions that refer to anti-phase synchronization (two clusters of 10 pendula shifted by π) i.e. x 0 = 0.1,ẋ 0 = 0.00057, ϕ k0 = 0.0,φ k0 = 9.81, ϕ j0 = 3.09,φ j0 = 9.784 where: k = 1, 2, . . . 10 and j = 11, 12, . . . 20. The value of k x is increasing in subplot (c) and decreasing in (d). Similarly as in the two pendula case, we observe the region (k x ∈ [100, 750]) where two solutions coexist: anti-phase synchronization and non-synchronous state. To further analyse multistability in that system we use proposed method.
In Figure 10 we present how the probability of reaching a given solution depends on the parameter k x . In subplot (a) we show the results for the system with 2 pendula, while in subplot (b) results obtained for the system with 20 pendula suspended on the beam are given. In both cases we consider k x ∈ [0, 5000] and assume the following ranges of initial conditions:  For subplots (a,c) the value of the bifurcation parameter kx was increased, while for subplots (b,d) we decreased the value of kx. Grey rectangles mark the ranges of the bifurcation parameter kx for which different attractors coexist. Further analysis of number of solutions can be found in [13]. ϕ i0 ∈ [−π, π],φ i0 ∈ [−π, π] where i = 1 . . . 20 in Figure 10(b). We take 20 subsets of parameter k x values with the step equal to ∆k x = 250 and mark their borders with vertical lines. For each set we run N = 100 simulations; each one with random initial conditions and k x value drawn from the respective subset. Then, we estimate the probability of reaching given solution. The dots in Figure 10 indicate the probability of reaching a given solution in the considered range (dots are drawn for mean value, i.e, middle of subset). Contrary to both already presented systems, this one has a much larger dimension of phase space (six and forty two), hence we decide to decrease number of the trials to N = 100 in order to minimise the time of calculations.
In Figure 10(a) we show the results for 2 pendula. When k x ∈ [0, 250] only anti-phase synchronization is possible. Then, with the increase of k x we observe a sudden change in the probability and for k x ∈ [750, 1750] only complete synchronization exists. For k x > 2000 a probability of reaching both solutions fluctuates around p(complete) = 0.7 for complete and p(anti − phase) = 0.3 for anti-phase synchronization. Further increase of k x does not introduce any significant changes.
In Figure 10(b) we show the results for twenty pendula. For k x ∈ [0, 250] the system reaches solutions different from the two analysed (usually chaotic). Then, the probability of reaching complete synchronization drastically increases and for k x ∈ [750, 5000] it is equal to p(complete) = 1.0 which means that the pendula always synchronize completely. We also present the magnification of the plot where we see that in fact for k x ∈ [715, 5000] we will always observe complete synchronization of the pendula. Please note that for calculating both plots we use random initial conditions and k x value hence, the results for a narrower range may differ. Anti-phase synchronization was never achieved with randomly chosen initial conditions. This means that even though this solution is stable for k x ∈ [100, 750] (see Figure 9(c)) it has a much smaller basin of attraction and is extremely hard to obtain in reality. The results presented in Figure10 prove that by proper tuning of the parameter k x we can control the systems behaviour even if we can only fix the k x value with finite precision.

Conclusions
In this paper we propose a new method of detection of solutions' in non-linear mechanical or structural systems. The method allows to get a general view of the system's dynamics and estimate the risk that the system will behave behave differently than assumed. To achieve this goal we extend the method of basin stability [23]. We build up the classical algorithm and draw not only initial conditions but also values of system's parameters. We take this into account because the identification of parameters' values is quite often not very precise. Moreover values of parameters often slowly vary during operation. Whereas in practical applications we usually need certainty that the presumed solution is stable and its basin of stability is large enough to ensure its robustness. Hence, there is a need to describe how small changes of parameters' values influence the behaviour of the system. Our method provides such a description and allows us to estimate the required accuracy of parameters values and the risk of unwanted phenomena. Moreover it is relatively time efficient and does not require high computational power.
We show three examples, each for a different class of systems: a tuned mass absorber, a piecewise smooth oscillator and a multi-degree of freedom system. Using the proposed method we can estimate the number of existing solutions, classify them and predict their probability of appearance. Nevertheless, in many cases it is not necessary to distinct all solutions existing in a system but it is enough to focus on an expected solution, while usually other periodic, quasi-periodic and chaotic solutions are classified as undesirable. Such a strategy simplifies the analysis and reduces the computational effort. We can focus only on probable solutions and reduce the number of trials omitting a precise description of solutions with low probability.
The proposed method is robust and can be used not only for mechanical and structural systems but also for any system given by differential equations where the knowledge about existing solutions is crucial.