Approximating piecewise nonlinearities in dynamic systems with sigmoid functions: advantages and limitations

In the industry field, the increasingly stringent requirements of lightweight structures are exposing the ultimately nonlinear nature of mechanical systems. This is extremely true for systems with moving parts and loose fixtures which show piecewise stiffness behaviours. Nevertheless, the numerical solution of systems with ideal piecewise mathematical characteristics is associated with time-consuming procedures and a high computational burden. Smoothing functions can conveniently simplify the mathematical form of such systems, but little research has been carried out to evaluate their effect on the mechanical response of multi-degree-of-freedom systems. To investigate this problem, a slightly damped mechanical two-degree-of-freedom system with soft piecewise constraints is studied via numerical continuation and numerical integration procedures. Sigmoid functions are adopted to approximate the constraints, and the effect of such approximation is explored by comparing the results of the approximate system with the ones of the ideal piecewise counter-part. The numerical results show that the sigmoid functions can correctly catch the very complex dynamics of the proposed system when both the above-mentioned techniques are adopted. Moreover, a reduction in the computational burden, as well as an increase in numerical robustness, is observed in the approximate case.


Introduction
Mechanical systems with free-play gaps and piecewise constraints have aroused the interest of the scientific community in the past decades. They find practical engineering applications in the study of impacting capsule systems [1], impact drilling systems [2], cracked systems [3], mechanical oscillators [4][5][6][7] and aeroelastic systems with free-play gaps [8,9], buildings subjected to earthquakes [10], and impact oscillators with rigid walls [11][12][13]. These systems, generally referred to as non-smooth dynam-ical systems, incorporate discrete-events which guide the associated mathematical solution through regions of the domain where different characteristics are applied. They show peculiar dynamic features, named discontinuity-induced bifurcations which cannot be found in the smooth dynamical systems and include: border-collisions, border-equilibrium, grazing, slidingsticking, and corner angle bifurcations [14]. Normally, non-smooth systems with free-play gaps are modelled as hybrid and piecewise-smooth continuous dynamical systems 1 , for which the most common type of bifurcations are the grazing bifurcations [12,15]: such discontinuity-induced phenomena correspond to a topological change of the system phase-plot [9] due to the sudden change of the characteristics in the governing equations, which can generate very complex dynamic behaviours, like stability change, generation of quasi-periodic responses, period-doubling cascade, and chaos [7,9,12].
Although analytical solutions [4,16,17] and mapping techniques [12,15,[18][19][20] exist, numerical methods have been widely adopted in engineering applications. Direct numerical integration has been particularly used in the literature [5,6,[8][9][10]21,22]. It can handle non-smooth properties in two different ways: either locating the exact discontinuity point 2 and restarting the integration process or directly integrating the function over its discontinuity. Missing this point can affect the numerical solution: thus, Henon's method [23] and the MATLAB built-in event function [21] have been used in recent studies. Wright et al. [24] and Dallas et al. [22] studied the capabilities of the built-in MATLAB event function, and they proposed guidelines for selecting the solver parameters. Dai et al. [8] studied an aeroelastic system with a free-play gap; they proved that the adoption of the exact event location with Henon's method allows identifying correctly chaos and transient chaos responses that could not have been found without the event location or with rational polynomial approximation. Recently, Saunders et al. [6] performed similar studies on a single-degree-of-freedom (SDOF) Duffing system with a free-play gap [25]. They used the MATLAB built-in function and proved that the Basins of Attraction (BOAs), the phase plots, and the transition to the chaos of the system are not strongly affected by errors in the event location, however, poorly converged solutions and basins boundaries modifications have been encountered without the exact event location. Despite the event location methods provide more accurate results, the procedure requires high computational demand, strict tolerances, and more complicated algorithms [26].
Other solutions methods for nonlinear systems with piecewise functions include Harmonic Balance Method [27][28][29], coordinate transformations [30][31][32], and path-following continuation methods [1,2,33]. For the last category, authors generally resort to Numerical Continuation toolboxes, like MatCont [34], COCO [35] and AUTO-07 [36]. However, only some of them can deal with non-smooth dynamical systems. Multisegments continuation procedures [35,37] allow solving system with piecewise characteristics: Chávez et al. [1,2] used the multi-segments toolbox TC-HAT [37] to study piecewise-linear impact oscillators. Liu et al. [38] studied the effect of a position feedback control strategy of a vibro-impact drilling via multisegments continuation: they adopted the hspo COCO toolbox to perform their analyses. In other works [33,[39][40][41][42], authors resorted to numerical continuation toolboxes, proving their effectiveness in solving non-smooth dynamical systems. Nonetheless, the analyses available in the literature are generally limited to one or few simple orbits. Moreover, the adoption of the multi-segment continuation method is timeconsuming, complex, and long, especially when multiple re-segmentation procedures are required [35].
In order to face such a problem, authors have suggested adopting smoothing functions to transform the non-smooth dynamical systems into a simplified version [26,[43][44][45]. Narayanan et al. [46] proposed the implementation of sigmoid functions to approximate piecewise characteristics for the first time. Early studies on a piecewise SDOF were conducted numerically and experimentally by Wiercigroch. [43]. The author found that the sigmoid approximation can represent the dynamics of the system, however, only in certain conditions, the approximation does not produce the same chaotic attractors identified by the ideal piecewise counterpart. However, the author did not evaluate the extension of the error caused by the smoothing function, focusing the numerical observations on chaotic responses. Wolf et al. [47,48] studied the dynamical stability of piecewise systems; in particular, they analysed the effect of smoothing functions on the stability of an SDOF mechanical system with free-play gaps when the Harmonic Balance Method is utilised [44]. They demonstrated that the smoothing functions introduce only limited inaccuracies in the computation of the system stability and that the effect of the number of harmonics is more important than the effect of the smoothing approximation. Vasconcellos et al. [45] studied the effect of free-plays on pitch-plunge aeroelastic systems with a trailing-edge flap. They compared the results obtained with the Henon method's [23] with the ones obtained with smoothed functions. The comparison with experimental measurements confirmed that piecewise functions and hyperbolic tangent functions can correctly represent periodic and chaotic behaviours of the system. Dai et al. [8] modelled the free-play gap of an aeroelastic structure with Rational Polynomials (RP): the approximation provided good results when attractors associated with non-chaotic or strongly chaotic behaviours were studied, however, it was not able to correctly identify the dynamic behaviour from transient chaos to periodic solutions. Wang et al. [49] adopted a smoothed C 2 function to simulate the impact of a cantilevered pipe conveying fluid. They numerically validated the proposed approximation by adopting the Runge-Kutta integration algorithm. Recently, Saunders et al. [26] analysed the effect of smoothing approximation functions on an SDOF Duffing system [25] with soft and hard contacts. Numerical integration analyses were performed with MATLAB ode45: the authors demonstrated that the tangent-based approximations with high convergence could describe the majority of the nonlinear phenomena of the system. In other works, authors employed smoothed functions to compute the response of non-smooth dynamical systems with the Harmonic Balance methods: in all these works, the regularisations [50][51][52] have been used to avoid numerical problems during the computation.
From a more theoretical point of view, smoothing approximations play a very important role in the structural stability of non-smooth dynamical systems. This is particularly true for systems with discontinuous characteristics, such as the classical Filippov systems, for which the uniqueness of the solution cannot be guaranteed at the discontinuity boundary [14]. In the literature, authors have demonstrated that there exist conditions for which the use of smoothing functions is admissible and produce commensurate and acceptable dynamic responses [53,54]; nonetheless, more recent studies have highlighted how such regularisations can create ambiguities [55]. Jeffrey et al. [56] suggested that the nature of such ambiguities is precisely the nonuniqueness of the solution which induces the generation of a set of solutions at the discontinuity border. Each of these solutions is equally valid but represents only a portion of the system dynamic response. Nevertheless, in some cases the Hidden Dynamics [57] of the sliding solutions 3 could be missed: in fact, the adopted regularisation (see Examples 1 and 2 of [56]) and solution method (see Example 1 of [57]) can deeply affect the solution of the system. If the solution method can only miss the hidden attractors, the regularisation can affect their existence by preserving the nonlinear properties of discontinuity boundary [57]. This is particularly clear in systems whose switching boundary is nonlinearly defined, e.g. when the border includes functions like sign(x) 3 . In these cases, the adoption of the Filippov regularisation [58] could prevent the generation of hidden attractors; this happens because the modelling assumptions could exclude hidden terms in the definition of the switching boundary [59]. The ambiguity problem of discontinuous systems cannot be solved even if the exact discontinuity point is exactly located, because the mathematical theory does not provide a clear definition of the dynamics that should be applied in that point [57]. Therefore, the problem of the non-uniqueness of the solution for Filippov systems is still an open question with the background mathematical theory still in its infancy.
Despite the large amount of work on non-smooth dynamical systems, most of the practical implementations focused their attention on SDOF systems and numerical integration schemes, neglecting the implications of utilising smoothing functions with multidegree of freedom systems and path-following continuation procedures. In order to investigate this problem, we focus our attention on non-chaotic attractors generated by a nonlinear two-DOF mechanical system with a soft free-play gap. This clearance introduces a piecewise characteristic in the system which is modelled by a continuous function of class C 0 , as typically happens in engineering applications (e.g. see [6,8,9,26,42,43,45,50,60]). Systems with similar properties fall in the category of the piecewise-smooth continuous systems for which the Picard-Lindelöf Theorem guarantees the uniqueness and existence of the solution in the whole domain. Sliding motions, typical of Filippov systems, are not permitted as it is impossible to have solutions stuck on the discontinuity boundary [14]. This also excludes the possibility of having multiple ambiguous solutions at the switching boundary [59], as the system solution is unique and determined. Finally, the presence of a soft-free play gap avoids the generation of chattering sequences (also called Zeno phenomenon [61]) which could generate sticking conditions in hybrid dynamical systems.
The system is studied with both numerical integration techniques and path-following methods, considering two configurations: the first one accounts for the sigmoids approximation and adopts the methods commonly used for solving the smooth dynamical systems, and the second one, instead, implements the nonsmooth piecewise characteristic and computes the solution with the more sophisticated methods, like integration schemes with event-location function and multisegments continuation procedures. The novelty of the paper is two-fold: first, it aims to study the effect of smoothing functions on the mechanical response of a multi-degree of freedom system when path-following methods are employed, secondly, it intends to provide more insight into the dynamic response of the proposed system, which, to the authors' best knowledge, has not been previously studied in the literature. Additionally, a new mathematical definition of the degree of approximation introduced by the smoothing function is proposed. The formulation, named radius of influence, can help researchers to understand the extension of the error introduced when such approximations are adopted, linking the mathematical background theory with practical application problems.

System description
The graphical representation of the two-mass system is reported in Fig. 1. The slightly damped system presents cubic and free-play nonlinearities where k and μ denote, respectively, the linear and nonlinear stiffness coefficients, and c defines the damping coefficient applied between the masses and the external constraints, while m represents the mass. The subscript Fig. 1 Representation of the two-mass system with a free-play gap • d , instead, refers to mechanical properties applied in between the two masses, which are the damping coefficient c d , and the linear, k d , and nonlinear, μ d , stiffness coefficients. The system is excited by the external loads Q 1 and Q 2 , while k p denotes the bi-linear piecewise stiffness applied on the first mass, whose restoring force can be expressed as: where a represents the free-play gap. The II order nondimensional equations of motion are represented by Eq. 2 and the non-dimensional parameters are mathematically represented as in Table 1. The reference values r [rad/s], t r [s], and l r [m], are selected equal to 1 and, for the sake of simplicity, the rest of the paper will refer to non-dimensional equations without using the symbols• and •.
can be easily converted in its autonomous I order version, also called Initial Value Problem (IVP), which is represented by Eq. A1 in the Appendix A. According to the dynamical systems theory, this equation formally represents a non-smooth dynamical system which, in this particular case, is characterised by a C 0 class function on the right-hand side. In the following sections, this definition will be fundamental Dimensional parameter for demonstrating the validity of the approximation, adopted in this work. Finally, Table 2 describes the parameters adopted for the analysis of the nonlinear system. Thesedimensional parameters have been obtained by starting from the data adopted in previous analyses of slightly damped nonlinear systems [62].

Function approximation
In this section, the smoothing sigmoid function is presented and mathematically analysed. This function is used to continuously approximate the piecewise function F p . Firstly, the admissibility of a such approximation is mathematically proven by using the theory developed by Danca [63][64][65] and secondly, a practical tool for estimating the required degree of approximation (the radius of influence) is introduced.

Mathematical admissibility
The considered non-smooth system of Eq. 2 can be represented in a general form with the following IVP: with initial conditions x(0) = x 0 and time t ∈ [0, T ].
The right-hand side of Eq. 3 is represented by the following functions: g(x) : R n → R n , which represents a nonlinear continuous smooth function, K ∈ R n×n which is a square linear real constant matrix, with K x denoting the linear part of f (x), and A(x) ∈ R n×n a square matrix of real function and s(x) : R n → R n a piecewise function, with s i (x) : R → R, i = 1, 2, ..., n a piecewise real constant function, which identify the non-smooth part of f (x). In the theory developed by Danca [63], the admissibility of smoothing approximation functions for systems like the one indicated by Eq. 3, is proved by using the following assumption, here reported for completeness: Assumption 1 [63][64][65] The function A(x)s(x) is discontinuous in at least one of its elements.
The discontinuity boundary M introduced by the piecewise functions s i can be defined with the following formal notation: Notation 1 [64] M denotes the null discontinuity set of f (x) with Lebesgue measure: ν(M) = 0, and it is generated by the discontinuity points of s i .
M separates the domain R n in many disconnected sub-domains D i ⊂ R n where f (x) is continuous and possible differentiable. Thanks to Notation 1, we can define the piecewise-continuous property of the function f (x), as follows: Under Assumption 1, the functions s i in the righthand side of the IVP 3 is piecewise continuous, hence the function f (x) is piecewise-continuous. It is worth noticing that Definition 1 embeds discontinuous functions such as the sign function sign(x) and the Heaviside step function h(x). Therefore, it is important to introduce the following formal definition: and, ∀ x 1 ∈ R n and ∀ x 2 ∈ R n there exists a real constant C > 0, such that: Remark 1 Given the Definitions 1 and 2, it is possible to state that for a closed and bounded real non-trivial interval we have the following inclusion: piecewisesmooth continuous functions ⊂ piecewise continuous functions Remark 1 is observing that the piecewise-smooth continuous functions are a subset of piecewise continuous functions. Essentially, they represent a degenerate particular condition of Definition 1, i.e. when the null discontinuity set M is collapsing to a single point which coincides with the adjacent sub-domains D i and D i+1 . Definition 2 imposes a more restrictive condition on the function f (x) which, now, has to be continuous in the whole domain R n .
In order to prove that the non-smooth system, outlined in Eq. 3 can be approximated with a smooth continuous dynamical system, the Theorem 1 is exploited.
wheres(x) denotes a smoothing continuous approximation of the piecewise continuous functions s(x) and F(x) is a set-valued function representing f (x) [66,67]. The demonstration of Theorem 1 is based on two steps: firstly the discontinuous function f (x) is converted in a set-valued function F(x) through the Filippov regularisation [58]. This allows finding a solution of the piecewise-continuous dynamical system even at the discontinuity point [65]. Then, the Cellina's Theorem [66,67] is utilised to guarantee the existence of a continuous approximationf (x) of the set-valued function F(x) on its domain. Figure 2 graphically describes the Danca's approximation process for a sign(x) function.
For a more detailed explanation of the proof of Theorem 1, the reader should refer to [63].
With Theorem 1, we can approximate the piecewise continuous function f (x) with a smoothed versionf (x), hence we can substitute the original IVP (3) with the following one: Eq. 6 represents a smooth dynamical system that can be solved without the need to consider discontinuity points. Danca [63] suggests using the sigmoid logistic function for globally approximating a piecewise continuous system, thus, by applying Theorem 1, we can write the following example: where δ represents the degree of approximation of the sigmoid function. According to Assumption 1, Theorem 1 is applicable only to dynamical systems with piecewise continuous functions, hence systems with discontinuous functions. Nonetheless, it is easy to extend the validity of Danca's theorem to continuous systems by considering the following corollary.

Corollary 1 Given a piecewise-smooth continuous function f (x), if g(x) is continuous then for every small neighborhood > 0 of F(x) there exists a smooth continuous approximation of f (x),f (x)
The Corollary 1 is demonstrated by Remark 1 and it guarantees the applicability of Theorem 1 to piecewisesmooth continuous dynamical systems. This extension is valid also for the piecewise-smooth continuous system of Eq. 2, which mathematically describes the twomass systems. The Corollary 1 is also in agreement with the numerical results obtained by Saunders et al. [26] which demonstrated that, provided a sufficiently high degree of approximation δ of the smoothing function, the numerical integration of the smoothly approximated system generates the same results than an equivalent piecewise-smooth continuous system.

Sigmoid function
Section. 3.1 demonstrated that continuous smooth functions, like sigmoids, can be used to approximate piecewise characteristics of non-smooth dynamical systems. Thus, to evaluate the effect of smoothing functions on the mechanical response of the proposed two-mass system, the piecewise bi-linear restoring force of Eq. 1 is approximated with the following sigmoid logistic function: where δ is a non-dimensional parameter that denotes the approximation accuracy. The approximated function F p,s introduces an error with respect to the ideal piecewise function F p , as illustrated in Fig. 3. This effect can be reduced by increasing the parameter δ, but it cannot be eliminated entirely. Although Theorem 1 allows approximating piecewise characteristics of a non-smooth system, it does not provide a numerical value for the degree of approximation δ. Nevertheless, it is possible to identify a suitable value for δ by comparing the dynamic responses of the ideal and the approximate smooth dynamical systems for the values prescribed by Table 2. Firstly, the vector fields of the two versions of the system are analytically computed and compared around the discontinuity boundary M. To facilitate the computation and avoid intersecting lines, the vector fields are computed for the undamped unforced version of the two-mass system with the second mass blocked, hence when: x 2 = 0 anḋ x 2 = 0. This condition represents a 2D section plane of the hyperspace of Eq. A1 which is composed of the following dimensions: x 1 , y 1 , x 2 , y 2 , v, u. The smoothed version of the system is obtained by using the sigmoid function described by (8), while the perfect piecewise counterpart is expressed with a linear combination of sign(x) functions, as follows: The results are reported in Fig. 4 where the blue and black arrows represent the system vector fields for the approximate and non-approximate conditions. Three numerically computed orbits of the first mass are also shown: the smallest one (P O A ) represents the noncontact condition, the largest one (P O C ), instead, denotes the penetration condition, and the middle orbit (P O B ) illustrates the grazing condition. The two fields appear to be entirely overlapped in the whole domain, even nearby the discontinuity boundary M. This suggests that the adopted degree of approximation δ is appropriate for the analysis of the two-mass system. Moreover, the error between the two vector fields is shown as colour maps in terms of magnitude and direction. These maps reveal that the distorted field is restricted to the discontinuity boundary with a maximum error in magnitude of 2.4% (Fig. 4a) and in direction of 1.8% (Fig. 4b) for the investigated portion of the domain. Therefore, the selected approximation δ = 1500 can be considered accurate enough for the analytical analysis of the piecewise system. Furthermore, the accuracy of the sigmoid approximation with δ = 1500 is validated with numerical experiments by comparing two numerical attractors, with and without the smoothing approximation. The MATLAB function ode45 is utilised in the numerical simulations with a maximum time step t = 0.4s for two forcing frequencies, namely = 1.2 and = 1.3. The numerical orbits are reported in Fig. 5: the attractors are very similar and practically overlapped thus, the parameter δ = 1500 is considered sufficiently accurate to describe the dynamics of the system, also, from a numerical point of view.

Radius of influence
In order to better understand the link between the -neighbourhood described by Filippov [58] and the approximating parameter δ, a mathematical formulation, called radius of influence R i is proposed. This formulation represents a practical tool to identify the appropriate smoothing parameter δ, avoiding the analytical and numerical analyses reported in Sect. 3.2. Since the approximation error is generally concentrated around the discontinuity point, the radius R i identifies the extension of the error by starting from that point. The radius can be computed by considering the sum of the areas 1 and 2 of Fig. 6: by starting from the discontinuity point we can calculate the error area of a sufficiently large portion of the domain such that we reach an asymptotic numerical value. This is guaranteed by the sigmoid function which globally approximates the piecewise characteristic [63]. Then, we can find the distance from the discontinuity point that constitutes the 68.0% of the total asymptotic error, i.e., the σ 1 error. The error can be numerically evaluated for the discontinuity point at negative displacement as: where r denotes the radius extension around the discontinuity point. Similarly, it is possible to easily obtain the expression of the error for the discontinuity point at positive displacement. The radius is computed in a domain where the horizontal axis indicates the displacement of the mass and the vertical axis denotes the equivalent displacement of the piecewise restoring function, i.e., F p /k p . Such a condition makes δ the only parameter affecting the radius R i . By considering the soft symmetric piecewise constraint, adopted in this work (k p = 10 and a = 0.05) and the sigmoid approximation with a value of δ = 1500, the radius of influence results to be equal to R i = 0.0017: this value is acceptable as it represents only the 3.4% of the non-contact gap. This confirms that the region affected by approximation error is very limited and it allows us to verify that the parameter δ is sufficiently large to accurately describe the dynamics of the piecewise system.
It is worth mentioning that the scope of the paper is to assess the effect of smoothing functions on soft constraints with large free-play gaps. Hard contacts or small gaps would require a quite large value of δ, which can lead to overflow problems as described by Wolf et al. [44]. In that case, tangential-based approximations could be adopted to mitigate the problem; nonetheless, a complete evaluation of the effect of the type of approximation is out of the scope of this paper.

Dynamic response of the approximated system
First, the nonlinear mechanical system of Eq. 2 is analysed using the sigmoid approximation F p,s as described in Eq. 8: path-following methods are used to compute the amplitude response, and numerical integration techniques are adopted to define the BOAs and bifurcation diagrams. In particular, the toolbox po of COCO [35] is adopted to perform the continuation pro-  Fig. 7. Different colours are used to highlight the effect of the ICs on the system: a similar representation allows discovering hidden attractors that could not be found otherwise. This is particularly evident in Figs. 7(a-b) and 7(d) where the IC denoted by the black colour (x0 = [0; 0; 0; 0; 1; 0]) leads the system to generate the same single-period attractor for most of the investi-gated frequency range . On the contrary, the remaining ICs generate complex multi-periodic and chaotic attractors which result in overlapped solutions in the diagram.
The analysis of the bifurcation diagram allows having a clear overview of the system dynamic behaviour, especially from the chaotic point of view; nevertheless, some solutions, with small BOAs, are difficult to obtain resulting in scattered non-continuous solution branches, as illustrated by Fig. 7b. The bifurcation diagram of Fig. 7 shows that chaos is originated nearby the grazing points at x 1 = 0.05, as generally happens for impact oscillators [15,42]: this behaviour is particularly clear in two regions: ≈ 0.93 − −0.98 and ≈ 1.3 − −1.75. The chaotic behaviour of the system could be triggered by the simultaneous presence of free-play gaps and low damping: indeed, the same behaviour has been found in similar piecewise systems [42,68]. Interestingly, chaos does not appear with the typical period-doubling cascade behaviour [69] but, instead, it seems to be generated by crisis bifurca-tions, e.g., internal and boundary crisis bifurcations, which abruptly transform regular attractors into chaotic ones [42,43]. Figures 7a and 7c-d illustrate a mixture of multi-periodic and chaotic responses which could represent the presence of internal and boundary crisis bifurcations.
The bifurcation diagram provides some important information about the system: it demonstrates the ability of the sigmoid function to identify typical nonlinear phenomena encountered in non-smooth dynamical systems, like chaos, and allows locating frequency ranges where stable period-doubling orbits exist. These ranges can be used in path-following continuation methods to identify additional branches, which may result in closed detached solutions (isolas). Finally, it is worth noticing that the presence of the piecewise restoring force generates a very rich dynamic behaviour.

Numerical continuation
Foremost, the system is analysed without the free-play gap: Fig. 8 reports the frequency response diagram for the amplitude of the first degree of freedom. The stable and unstable amplitude are depicted, respectively, with the continuous and dashed black curves, and the dots define bifurcation points, where FP is indicating foldbifurcations. The grey curves, instead, denote the system backbones: these curves are computed by tracking the Hopf bifurcation points, which can be easily continued by using Eq. 11 [35,62]. Indeed, this equation represents the modal version of Eq. 2, and it is obtained by applying the linear modal transformation. The continuation analysis returns the typical hardening response, with only one simple solution branch originating from an elliptic orbit. The mono-periodicity of the system is, then, confirmed via the bifurcation diagram which, for the sake of brevity, is not reported. It is worth noticing that the two-mass system does not develop backbone bifurcations, internal resonances, or other complex dynamics behaviour because the required conditions are not met [62].
In order to understand the effect of the piecewise contact, the complete system is studied again via pathfollowing continuation by adopting the sigmoid function. The results of the analysis are reported in Fig. 9 and Fig. 12 where BP denotes a branch bifurcation point, PD defines a period-doubling bifurcation point, and TR indicates a Neimark-Sacker (torus) bifurcation point, while PO indicates a periodic orbit. For the sake of simplicity, the fold (FP) and saddle-node (SN) bifurcations are not reported. Figure 9a shows the numerical continuation of the single period branches solutions in the projection: amplitude of the first DOF (x 1 ) versus excitation frequency ( ). B 1 represents the continued solution obtained from the elliptic orbit named PO 1 , which is reported on the left side of the figure. This orbit exists both in the contact and non-contact regions and moves through a grazing bifurcation point. B 1 , moreover, presents branch (BP) and torus (TR) bifurcation points, as shown in detail in Fig. 9b-c: the first one bifurcates in a second branch, named B 2 , whose periodic orbits are described by PO 2 A and PO 2B , while the second one originates a region where the system exhibits chaotic behaviour. It is worth noticing that the branch B 2 is originated by the bifurcation of the first backbone curve: the bifurcation appears when the amplitude x 1 overcomes the free-play gap a, as shown in Fig. 10. In the free-play region, below x 1 = 0.05, there is no coupling effect between the modes, i.e., the first backbone S 1 has zero contribution from the second modal DOF and vice versa. Above the free-play gap limit, instead, the backbones, S 1 and S 2 , have both the modal contributions: this generates a backbone bifurcation when the first mode is grown enough to trigger an internal resonance, which results in the branch S 3 . To the authors' best knowledge, the presence of internal resonances in piecewise systems has not been previously shown in the literature, thus its nature should be better investigated in future works.
Interestingly, the presence of piecewise stiffness characteristic induces a symmetry-breaking behaviour in the system orbits: at the same frequency and amplitude, there exist symmetric orbits, as reported in Fig. 11. These orbits can be obtained by applying opposite sign initial conditions and belong to the same branch in the projection amplitude x 1 vs ; nonetheless, they appear in separated branches when the considered projection is composed of the maximum amplitude x 1 vs the frequency , as shown in the detail of Fig. 9a for branch B 2 . Figure 9d shows a single period solution, named I 1 : this solution is closed and has neither branches nor period-doubling bifurcation points, however, it is overlapped to B 1 , and B 2 in the projection amplitude x 1 vs frequency . By taking a different projection, as shown in Fig. 9e, the solution is completely closed and it does not appear to overlap other branches: this type of branch is referred to as isola. Again, the symmetric behaviour of periodic orbits is found in the isola I 1 , as testified by the symmetric periodic orbits P O 3A and P O 3B of Fig. 9.
After the analysis of the single-period orbits, the period-doubling stable solutions are investigated. Although some authors [50,70] have shown that perioddoubling isolas can be tracked with LP code-2 bifurcations, in this works, we have adopted a simpler approach to identify such isolas: firstly we have computed the bifurcation diagram, as shown in Fig. 7, through direct integration to have an idea of the regions of frequency where period-doubling solutions exist, then, for each isola, the solution has been continued with path-following methods. The initial orbits are reported in the left column of Fig. 12, while the continued solutions are reported in Figs. 12(a-l). Figures 12a,e and (i) show the location of the isolas with respect to the main branches B 1 , and B 2 . Most of them appear to be completely detached from other branches, however, even overlapping solutions like I 3 can be demonstrated to be separated by considering different projections, as previously shown for I 1 . The right column of Fig. 12 shows the detailed views of each isola: all the isolas present either torus (TR) or period-doubling bifurcation points (PD), while only I 6 shows bifurcation points. Some isolas are very small (e.g., I 2 and I 6 ), nonetheless, others are larger and stable, e.g. I 10 , representing an issue in the design and the analysis of similar systems. Interestingly, the symmetry-breaking behaviour is not found in multiperiodic orbits which show an anti-symmetry property along the axisẋ 1 = 0, like PO 9 , PO 10 , PO 12 . It is worth noticing that isolas I 4 , I 6 , and I 7 lay partially in the chaotic regions, thus they could be very difficult to identify without path-following continuation procedures. Finally, it is worth remembering that the perioddoubling isolas found are period-doubling 3,4,6,7,9,11, and 17 which testify to the very complex dynamics of such a system.

Basins of attraction
In order to further understand the dynamic behaviour of the mechanical structure with a free-play gap, the BOAs of the system are studied for the frequencies equal to 1.2, 1.25, 1.3, and 1.35. The basins are computed through a direct numerical integration procedure, by integrating Eq. A1 with the ode45 MATLAB function  Fig. 13 and show the coexistence of many periodic orbits at the same frequency . In particular, Fig. 13a shows the presence of two well-defined regions of attraction R 1 and R 2 : the first one is associated with the singleperiod orbits P O 1 of branch B 1 while the second one is generated by the higher amplitude orbits P O 2 of B 2 . At higher frequencies, more complex regions with the coexistence of single-period and multi-period solutions appear: Fig. 13b shows the presence of a quite large single-period region with low amplitude called R 3 , and a high amplitude region, R 4 , with the coexistence of single-period and multi-period orbits. In particular, R 3 shows the presence of only the periodic orbit P O 1 , while R 4 presents three different periodic orbits: PO 2 , PO 3 , and PO 11 ; the first two have period 1 while the last one has period 4. It should be noted that this region does not have a proper shape, and the coexistence of such different solutions in the same zone of the basin is a symptom of the vicinity of chaotic behaviour. The chaos of this system seems to be generated by the continuous switch between different periodic orbits, coexisting in the same region of the BOA. A similar discussion can be done for regions R 5 and R 6 of Fig. 13c, with the exception that the low amplitude region repeats itself even at very large initial conditions. In this case, the orbits P O 1 , P O 2 , P O 3 , P O 5 , and P O 10 , are identified, while P O 7 is not found. This orbit does exist at = 1.3, as shown in Fig. 12, however, being very near to a change of stability, the size of its BOA is very small, thus it is very difficult to match the exact initial conditions that would have led to obtaining such an orbit. Complete chaotic solutions start appearing at = 1.35 as shown in Fig. 13d. Again, we can identify two regions: R 7 with low amplitude and R 8 with high amplitude. Chaos appears in the regions with lower amplitude response, and it is mixed with regular orbits with many periods, while the high amplitude region R 8 shows again the existence of different periodic solutions, but not chaos. However, also in this case, all the previously identified orbits (P O 2 , P O 3 , and P O 9 ) are found in the analysis of the BOA.
It is worth noticing that even if all the identified periodic orbits have a correspondent isola or branch from the previous continuation analysis, it is possible that periodic orbits, with BOA very small, are missing as happened for the orbit P O 7 .

Non-smooth Dynamical System Comparison
To study the effect of the sigmoid approximation, the dynamics of the approximate system are compared with the dynamics of the non-approximate piecewise nonlinear system. Again, numerical integration and pathfollowing analyses are performed, however, different procedures must be adopted to handle the non-smooth piecewise characteristic.

Effect of the sigmoid function on numerical continuation
Numerical continuation procedures of a non-smooth dynamical system can be performed with the hspo COCO toolbox: this toolbox allows multi-segments numerical continuation, which consists of continuing portions of the orbit called segments. The procedure adopted in this work is similar to the one used by Liu et al. [38]: each segment can be described with a smooth equation that is associated with one condition of the piecewise property. This set of smooth equations is generally called vector field equations. The limit between one segment and another is identified by an event function, which describes the discontinuous event. Finally, a restart function defines the initial condition of the next segment. A more detailed explanation of multisegment continuation procedures is reported in [2,35]. The vector field equations, the event function, and the restart function for the two-mass system are reported in Appendix A. Although this procedure does not require any approximation of the piecewise characteristic, its practical application is often difficult, computationally expensive and time-consuming: this is especially true for systems that show many grazing bifurcation points in the same branch. Indeed, grazing bifurcation generates or eliminates segments during the path-following continuation process: this demands a re-segmentation of the orbits to avoid ending in a non-physical dynamic response [35]. The re-segmentation procedure is graph- ically shown in Fig. 14: when the orbit is small enough, two segments can fully describe the orbit, as depicted in Fig. 14a. By increasing the frequency, the orbit reaches the grazing condition, denoted by GR: here, the maximum displacement of the first DOF is equal to the freeplay gap a and the contact occurs at y 1 = 0. This condition can be described with two segments, the upper and lower segments depicted in Fig. 14b, and two additional points at x 1 = ±a, which represent degenerated segments. Increasing the frequency, the two degenerated segments increase their length, describing the portion of orbit exceeding the free play gap, as reported in Fig. 14c. This process is called re-segmentation, and it must be performed at every grazing bifurcation point to avoid the creation of non-physical orbits. Taking into account the difficulties of this numerical procedure, only some solutions previously obtained with the sigmoid approximation, are now replicated for the comparative analysis. The continuation of the smooth approximated dynamical system is performed with polynomial degrees equal to 4 and, respectively, 340, and 500 orbits subdivision (NSTS) for the single period branches and isolas. The continuation of the non-smooth dynamical system with the package hspo, instead, is performed with 85 and 25 NSTS and polynomial degrees equal to 4, for single period branches and isolas. During the simulations, it has been noted, that the approximated system required numerous subdivisions of the orbits: this may be due to numerical difficulties introduced by the sigmoid function. Figure 15 describes the comparative analysis between the results of the continued solutions with the sigmoid approximation (blue) and perfect piecewise characteristic (black): the major dissimilarities are shown in detail, indicating the exact location of the zoom with circle in the top part of each figure. Overall, the comparison shows that the sigmoid function allows obtaining very good results both in terms of amplitude and stability of the solution, as suggested by previous studies of Wolf et al. [44]; however, consistent differences can arise nearby the grazing bifurcation points GR. At the grazing condition, the smoothed function does not perfectly approximate the piecewise characteristic, generating errors in the continued solution. It is worth noticing that the extension of such error is limited to the radius of influence R i , and, does not seriously affect the overall solution of both single-and multi-period solutions. Figure 15a shows two grazing bifurcation points of the branch B 1 : given its nature, the sigmoid approximation can not detect the grazing condition, however, it is able to correctly identify torus bifurcation even nearby a grazing point, as well as the change of stability and the amplitude of the response. In the figure, it is also possible to see that the extension of the error is very limited to the region nearby the grazing bifurcation, especially in terms of amplitude. Figure 15b, instead, shows the bifurcation point between the branch B 1 and B 2 : in this case, the sigmoid function can correctly locate the branch, flip, and saddle-node bifurcation point, nevertheless, it does not perceive the instability of the first branch before the point BP and the occurrence of a torus bifurcation. It should be noted that smaller continuation steps could help in finding missing bifurcation points; moreover, the inaccuracies are located only in a small area around the BP point. Figure 15c-d shows the compared results for the single period isola I 1 ; again stability and amplitude are overall well identified, however, discrepancies appear nearby the grazing bifurcation points: in particular, the S shape solution of  Fig. 15d shows the presence of a torus bifurcation in approximated solution, which is not found with the multi-segment procedure. Finally, Figs. 15e-f shows the compared result for the periodtripling isola I 3 : in particular, Fig. 15e shows that the sigmoid function can correctly locate period-doubling bifurcation even if there are small differences in the stability regions. Figure 15f, instead, shows a point where multiple grazing bifurcation occurs: the maximum differences arise in the vicinity of the grazing conditions with a quite large difference in the amplitude response. Nonetheless, such discrepancies are very limited even in this case, demonstrating that the sigmoid function can provide accurate results.
Finally, it is worth noticing that the amplitude of the approximated solution is systematically larger than the non-approximate counterpart: the error of the smoothed function in the free-play gap, shown in Fig. 3, affects all the system responses, increasing the overall amplitude. This effect is responsible for all the discrepancies nearby grazing points previously addressed, thus, control parameters like the radius of influence R i , the ratio between the linear restoring force and the maximum negative force introduced by the approximation, should be checked to avoid excessive distortion of the response. Despite these negative effects, the sigmoid approximation demonstrates to be suitable for analyses of systems like the one proposed in this paper, especially when soft free-play constraints are adopted.

Effect of sigmoid function on the basins of attraction
The BOAs of the two-mass system are now analysed with and without the sigmoid approximation by using the MATLAB function ode45: the smoothed system can be directly adopted in the numerical integration procedure, while, the ideal counterpart needs to handle the piecewise characteristic in a different way. To this end, the MATLAB built-in event function [21] is used: this function allows the switch of non-smooth properties during the integration procedure, identifying the exact location of the change. The procedure adopted in this paper is similar to the one adopted in [6,22,24]: for each event occurrence, the time integration stops and the final conditions are used as initial conditions for the next integration step. The numerically integrated ODE is modified accordingly with the position of the first mass, x 1 , as described by Eq. 1. In order to compare the dynamic response of the two system configurations, the autonomous version of the equations of motion, described by Eq. A1, is adopted. The such definition allows us to easily switch between the sigmoid and the piecewise function on the base of the definition of F p . The comparison is carried out with the options described in Table 3, which prescribes the relative and absolute tolerances, the maximum time step, and the time span, adopted for the integration procedure. Figure 16 represents the BOAs of the system obtained at different excitation frequencies, in particular: Fig. 16a-d represent the solutions obtained with the exact event location, while Fig. 16e-h represent the solutions obtained with the sigmoid approximation. The results show that the sigmoid approximation is able  to correctly represent the BOAs of the considered two-DOF system with high accuracy without introducing errors or incorrect topologies to the basins. For all the analysed excitation frequencies, the shape, the boundaries, and the maximum amplitude of the basins are not significantly deformed by the introduction of the approximation. The good performance of the proposed smoothed function can be assessed by looking at some details: Fig. 16b and f show the presence of the same small tongue in the lower part of the region R 3 , while Fig. 16f-h report the same granular shaped regions R 4 , R 6 , R 7 , and R 8 obtained in Fig. 16b-d, confirming the presence of co-existing multi-periodic orbits in both the cases. Only the mixed region of R 1 and R 2 of Fig. 16a and e shows some differences between the two cases, especially in the zone between x 2 = 0 and x 2 = −0.1.
The event location function is very accurate [6], however, it has two important drawbacks: it requires long time-consuming numerical simulations, and it must use strict tolerances to correctly locate the event.
In the case of slightly damped systems, the computational burden problem is particularly evident, because such systems require a long simulation time to overcome the transient phase. Table 4 reports the elapsed time during numerical simulations in the two considered cases: with smoothed sigmoid function and with the event location capability. Again, the function ode45 is utilised to simulate the origin point of BOAs 5 by considering the options of Table 3: the results show that the numerical integration with the event function is about 10 times longer.
The second problem is related to the absolute and relative tolerances: while the maximum step can play a certain role in the duration of the chaotic transient phase, the tolerances adopted in the numerical simulations can generate very different BOAs, if wrongly implemented. In order to investigate this numerical problem, we have computed the BOAs at = 1.3 using the event location function with different maximum 5 The origin point denotes the cross point between the abscissa and ordinate in the BOA representation. time steps and absolute/relative tolerances, respectively equal to: 0.05 and 1e-06, 0.08 and 1e-6, and 0.08 and 1e-9. Figure 17 reports the results of the analyses: when the maximum time step is equal to 0.05 and the tolerances are set at the value of 1e-06, numerical instabilities arise in the computation of the BOA, as reported in Fig. 17a. Most of the problems seem to be located at the extrema of the basin, where large initial conditions are adopted. Instead, the central part of the basin looks very similar to its counterpart previously computed with the sigmoid function, as shown in Fig. 16g. Similar results are obtained when the same tolerances and a larger step, equal to 0.08, are adopted, as reported in Fig. 17b. The only significant change is the return of the rings of regions R 6 : however, at the large initial condition, the numerical instabilities are still present. Finally, when the tolerances are increased to the value 1e-09, the correct BOA appears, as Fig. 17c shows. Indeed, the strict tolerances guarantee to identify of the correct location of the event at each crossing point: for large initial conditions, the event function can miss the event occurrence, leading the system to end in another steady-state condition. Figure 17d-f shows the effect of missing the location of the event: the displacement of the first mass is plotted in red when the algorithm applies Eq. 1 with the condition x 1 >= a, in blue when the condition x 1 <= −a is applied, and black when the mass is considered to oscillate in the free-play gap. Figure 17d-e shows that the adoption of tolerances equal to 1e-6 is not enough to capture the correct event location in the first part of the simulation while using stricter tolerances, equal to 1e-9, the event is correctly located every time, as shown by Fig. 17f. The exact location of the events allows for correctly representing the BOA. This phenomenon constitutes a serious issue in the application of numerical simulation procedures with event function: indeed, the problem does not prevent the attainment of the steady-state condition, instead, it leads to a completely different stable attractor which can potentially change the shape of the BOA. Moreover, it could appear even at smaller scales when small free-play gaps are used.

Conclusion
In this paper, a two-degree of freedom system with soft piecewise constraints was studied via numerical continuation and numerical integration. First, the mathematical admissibility of the proposed approximation was demonstrated and verified with analytical and numerical tests. Then, the effect of the sigmoid functions on the dynamics response of the system was analysed: the numerical results showed that the smoothed system can correctly describe the overall system dynamics, in terms of both BOAs boundaries and quality of the continued solution. The robustness and lower computational burden of the approximation were proven, and the main discrepancies between the ideal and the approximate solutions were found nearby grazing bifurcation points, where the approximation is less accurate. No distorted or non-physical solutions were identified during the analysis, however, a systematic increase in the amplitude response of the approximate solution was encountered: this is due to the approximation error introduced by the sigmoids function and remarks the necessity of tools, like the proposed radius of influence, to control the degree of approximation introduced in such systems.
The numerical results also demonstrated that the presence of the piecewise constraint completely changes the dynamic response of the system: indeed, many isolas were found around the two system natural frequencies, and the occurrence of an internal resonance was proven. Interestingly, the internal resonance results in a backbone bifurcation which generates the second branch of stable solution: this behaviour implies an energy exchange between the modes, and it can be identified only in MDOF systems. Moreover, it represents a critical aspect of the design of mechanical structures: indeed, the additional branch cannot be detected with the classical linear design methods, thus the system could experience unforeseen displacements which could end in system failure.
Finally, the authors want to stress the limits of this analysis as this work focused attention on soft piecewise constraints and non-chaotic attractors. Future works could try to address the effect of the contact stiffness or the type of smoothing functions. system implements Eq. 8 in numerical integration and continuation procedures, while the non-approximated version utilises the definition of Eq. 1. Two additional states, namely u and v, allow removing the dependence on time, making the system autonomous.
x 1 = y 1 (A1a) x 2 = y 2 (A1c) The multi-segments continuation, however, requires the definition of the vector field equations, event function, and reset function. The first one represents a set of smooth equations to be used for continuing each orbit segment: it can be considered as a modified version of Eqs. A1 which account for the piecewise function as follows: No Contact -this case occurs when the first mass is not in contact with the external piecewise stiffness, and the piecewise restoring force of Eq. A1b is null: Contact up -this case occurs when the first mass touches the upper limit of the external piecewise stiffness and the associated restoring force becomes: Contact down -this case, instead, occurs when the first mass touches the lower limit of the external piecewise stiffness, which imposes the following restring force: The event function, instead, identifies three possible conditions: Impact up -the segment ends when the first mass touches the upper limit of the external piecewise stiffness. In this case, the event condition r ev is identified by: Impact down -the segment ends when the first mass touches the lower limit of the external piecewise stiffness and the event conditions are identified by: Velocity change -this segment ends when the first mass changes the velocity sign. This is particularly useful for identifying grazing bifurcation points when noncontact orbits are continued. It can be described by the following condition: Finally, the reset function imposes the restart conditions for the next segment continuation as x new = x, where x new indicates the new set of initial values for the continuation of the next segment.