Direct and inverse problems encountered in vibro-impact oscillations of a discrete system

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


INTRODUCTION
Vibro-impact oscillations are of considerable practical importance and occur in many engineering applications, including systems with clearances, loose joints, or motion confining rigid constraints.Vibro-impacts lead to strong and essential nonlinearities and their analytical or even numerical study presents challenging technical difficulties [1][2][3].
Among various studies in the literature, Masri and Caughey [4] analyzed linear systems with rigid constraints by matching linear solutions computed before and after the time instants of impacts.Studies of piecewise linear and vibro-impact oscillations with analytical/numerical Poincare´maps and geometrical techniques were performed in a series of works by Shaw and Holmes [5], Moon and Shaw [2], Shaw [6], Shaw and Rand [3], and Shaw and Shaw [7].Ivanov [8] studied vibro-impact oscillations by introducing auxiliary phase planes, and Zhuravlev [9] analyzed these oscillations by employing special nonsmooth transformations of variables.In a related work, Vedenova et al. [10] studied localized and nonlocalized vibro-impact periodic solutions of discrete oscillators using nonsmooth coordinate transformations (cf.also reference [11]).Van de Vorst et al. [12], studied periodic motions of a beam undergoing vibro-impacts by using finite elements and reducing the problem to numerically solving two-point boundary value problems.Numerical and experimental investigations of flexible elastic and rotordynamic systems with impact nonlinearities were carried out in references [13,14].
Additionally [5,[15][16][17], local and global bifurcations of periodic orbits and chaotic motions in systems with vibro-impacts were examined.Chaotic motions of such systems were studied in these and other works [18][19][20][21].Whiston [22] presented a study of the ''shredding'' of the invariant manifolds of a periodically excited system with impacts due to the non-differentiable nature of the vibro-impact dynamics, and Brogliato [23] discussed the dynamics and controls of systems with impact nonlinearities.
In this work we study direct and inverse problems related to vibro-impact oscillations of discrete linear or nonlinear systems with rigid constraints.The theoretical part of the work centers on the use of functional relations between coordinates to study periodic solutions and local dynamics, as well as, the formulation of two-point boundary value problems (BVPs) to study analytic continuations and bifurcations of solutions.The numerical part of the work focuses on the numerical solution of the formulated BVPs and on the construction of Poincare´maps to study the global dynamics.An inverse problem is also posed, whose solution leads to the determination of linear or nonlinear systems that are capable of producing a specified vibro-impact motion in configuration space.Bifurcating problems associated with the inverse problem are discussed and numerical simulations that validate the theoretical findings are presented.

analytical results
We analyze the dynamics of the system depicted in Figure 1(a).It possesses two degrees of freedom, but the methodology to be developed can be extended to the multi-degree-of-freedom case.The motion of each particle is restricted by two rigid barriers, and takes place in the interval [−e, e].As a result, for sufficiently high energy of motion impacts occur and the equations of motion are expressed in the form x¨+ 1P(x, y) 1x + P(x) = 0, y¨+ 1P(x, y) 1y + P(y) = 0, where P(x, y) denotes the (smooth) potential energy, and P(u) models purely elastic impacts of the two particles with their rigid constraints.The elastic impacts are mathematically modelled by the relationship [cf. Figure 1(b)] where c(n) is a scaling constant depending on n.
In the direct problem, we seek analytic approximations and numerical solutions to vibro-impact, time-periodic orbits of the system, with a specified number of impacts per cycle.Whereas this problem has been addressed in the past by other authors (especially for the case of single degree of freedom oscillators), the technique that we propose is new since it relies on the concept of a 'nonlinear normal mode' (NNM) [11] to formulate a functional relation between the two coordiantes of the system during the vibro-impact motion.Hence, we assume that during the vibro-impact oscillation the following expression is satisfied at all time instants: The function y[•] is determined by substituting equation (3) into equation ( 1) and eliminating the time variable from the problem.This leads us to the following quasi-linear ordinary differential equation governing yˆ [•] in open intervals of x and y between consecutive impacts, where h denotes the (conserved) energy of the system, and primes denote differentiation with respect to x.Assuming that h − P(x, yˆ) q 0 during the motion, i.e. that the vibro-impacts occur before the system reaches its maximum potential energy value, the equation above is nonsingular, and must be complemented by appropriate boundary conditions at the elastic impacts.We suppose at this point that the periodic motion under consideration involves only impacts in the x-coordinate, and is composed of two branches y = yˆ1(x) and y = yˆ2(x), as depicted in Figure 2(a).Between impacts the functions yˆ1 ,2 [•] are governed by equation ( 4), whereas at each instant of impact they satisfy the compatibility conditions That is, at the time instant t = t* when the elastic impact occurs, there is C 0 and C 1 continuity in y(t), but only C 0 continuity in x(t).The C 1 discontinuity in x(t) satisfies the jump condition, x˙(t* − ) = −x˙(t* + ), which in turn leads to the conditions (5).Hence, the direct vibro-impact problem consists of solving two nonlinear ordinary differential equations governing the two branches yˆ1 ,2 [•] [each similar to expression (4)] subject to the compatibility conditions (5).Since no exact analytical solutions exist for the general nonlinear problem, it will be necessary to resort to asymptotic approximations.We note that more complicated vibro-impact oscillations corresponding to an arbitrarily large number of impacts in the xand/or y-coordinates can also exist, but these will not be considered herein.For such motions, a similar asymptotic analysis can be performed by considering vibro-impact motions possessing more than two branches and/or multiple impacts per period in the xand y-coordinates.
Before we consider the asymptotic solution of the direct problem, we remark that the previous discussion regarding the vibro-impact motion is valid only if the condition h − P(x, yˆ) q 0 is satisfied.If the system reaches at any time instant its maximum equipotential surface [corresponding to the condition P(x, yˆ) = h], the coefficient of the highest derivative in equation ( 4) vanishes and the problem becomes singular.Such problems with singularities at the maximum equipotential surface are encountered in calculations of nonlinear normal modes (NNMs) of discrete oscillators [11]; their solution requires imposing the following boundary orthogonality condition that is valid when P(x, yˆ) = h: This condition must be imposed whenever the orbit of the vibro-impact oscillator reaches the maximum equipotential surface.We now consider the analytic approximation of the two-branch time-periodic orbit depicted in Figure 2(a).It should be clear that such a vibro-impact response can only occur for sufficiently large values of the total energy h.Indeed, for sufficiently small h the two particles cannot reach their rigid constraints and the sought time-periodic motions are merely NNMs of the system with no impacts [11].To be able to perform analytic computations we assume the following specific form for the potential energy: In the above expression, the scalars a, b, d are the coefficients of the quadratic, quartic and cubic parts of the potential functions of the oscillators, whereas g is the coefficient of the weak linear coupling stiffness.At low energies and sufficiently small coupling g, this system possesses two localized NNMs, corresponding to periodic oscillations of the two masses with amplitudes of O (1) and O(g), respectively [27].We are interested in examining the behavior of these localized NNMs as the energy increases and impacts start to occur.Hence, in what follows we will assume that all coefficients in equation ( 7) are of O (1) with the exception of the coupling coefficient which will be assumed small, g 1. Considering the vibro-impact orbit of Figure 2(a), the two branches yˆ1 ,2 (x) must satisfy the relations (4) and (5).For a potential function given by equation (7), and taking into account that the problem possesses a small (perturbation) parameter g, we express the two branches in the following series of successive approximations: With equation ( 8) we seek vibro-impact motions that localize to the particle with coordinate x [note that with x = O(1), we seek y = O(g) x].Substituting equation ( 8) into the equations governing yˆ1 ,2 (x), we obtain for j = 1, 2. This expression is vaid for =x= Q e, i.e. between consecutive impacts, and for sufficiently large values of h so that the condition h − P(x, yˆ) q 0 is satisfied pointwise for each of the two branches of the vibro-impact motion.When =x= = e the elastic impact conditions (5) are imposed and matching of the two branches takes place.Considering terms in equation ( 9) proportional to g we obtain the following quasi-linear ordinary differential equation governing yˆj 1 : Since by assumption the coefficient of the leading derivative is a nonzero quantity, we can express the solution in the power series form with the coefficients of the series determined by substituting equation (11) into equation (10) and matching corresponding powers of x.Following this procedure we obtain the approximation where the leading coefficients of the series remain undetermined up to this point.A similar expression for the second (or higher) order approximations yˆj 2 can be derived in the form, yˆj 2 0 d j0 g 1 (x, h) + d j1 g 2 (x, h) + g 3 (x, h) + O(x 6 ), but this higher order correction will not be pursued in the present work.
Imposing the elastic impact conditions (5) and taking into account only terms of O(g) [expression (12)], we obtain the following matrix equation governing the leading order coefficients c jp , j = 1, 2, p = 0, 1: Provided that the matrix of the determinant of coefficients is nonsingular, the above equation yields a unique solution for the coefficients, which is derived symbolically using Mathematica, as follows: 4 (ae 2 + 3h) a(−9abe 6 + 10d 2 e 6 − 72ae 2 h − 36be 4 h − 288h 2 ) , This solution shows that leading order approximations for the two branches coincide, yˆ1 1 0 yˆ2 1 , and, furthermore, that the slope of the vibro-impact trajectory in the configuration plane at the points of impact is equal to zero, yˆ1 ' 1 (2e) 0 yˆ2 ' 1 (2e) = 0.In Figure 3 we depict the asymptotic approximation ( 11)-( 14) for parameter values a = b = d = 1, e = 1, g = 0•1, and h = 2.In the same figure we show the maximum equipotential energy level corresponding to h = 2, which is sufficiently wide to allow double-sided vibro-impact motions to occur.
Numerical simulations were also performed, integrating the equations of motion (1) with potential energy given by equation (7).The numerical results validated the analytical predictions, indicating that the double-sided vibro-impact localized motion is orbitally stable, and, hence, physically realizable.Moreover, using the numerical simulations it was established that the localized motion persists as a stable solution for arbitrarily high levels of the energy h, and that no bifurcations associated with this motion take place in the system for large h.The above theoretical and numerical results prove the existence of localized vibro-impact oscillations, where the energy of the motion is mainly confined to one of the two particles.Moreover, the motion in the configuration plane is represented by a line (and not by a closed loop), which indicates synchronous motion of the two particles (purely in-phase or out-of-phase).We will show later that this is not the case for systems with other types of potential function, which can support periodic motions possessing nontrivial phase differences between particles.
An interesting question concerns the behavior of the previous localized solution as h decreases, and we now briefly discuss this issue.For the system with (asymmetric) potential function (7), we expect that as the energy decreases below a critical value the particle characterized by the coordinate x will make the transition from double-sided to single-sided impacts, and (after further decrease of energy), eventually, to no impacting motions.Indeed, assuming small g and a localized periodic motion of the form yˆ(x) = gyˆ1(x) + O(g 2 ), the function yˆ[•] is still governed by equation ( 4).The first critical value of the energy h c1 separating the regimes of single-and double-sided impacts is determined approximately by the condition which corresponds to the energy level for which the particle characterized by x reaches its maximum amplitude when x = e, and, hence, the relation h c1 − P(e, yˆ) = 0 is satisfied.Since the potential energy of the system is asymmetric, at x = −e the condition h c1 − P(−e, yˆ) q 0 is still satisfied and impact of the particle with its rigid constraint occurs.For a = b = d = 1, e = 1 and g = 0•1 we find that h c1 = 1•1333.As a result, for a sufficiently small neighborhood of h bounded from above by h c1 the system undergoes single-sided localized vibro-impact oscillations, and the motion in the configuration plane yˆ(x) = gyˆ1(x) + O(g 2 ) is determined by solving equation ( 4) subject to the boundary conditions yˆ'(−e) = 0, (elastic impact at x = −e), (16a) where X 1 is the maximum amplitude attained by x at the right boundary of the single-sided vibro-impact motin (cf. Figure 4).Performing a perturbation analysis similar to the double-sided impact case, we express yˆ1(x) as, yˆ1(x) 6 ), where the functions f k (x, h), k = 1, 2, 3, are defined by equation (12).Imposing the conditions ( 16) we then determine the values of the coefficients c 0 and c 1 .
For the above values of the system parameters and h = 0•93333 Q h c1 , we find c 0 = 0•20042 and c 1 = −0•49791.In Figure 4 the motion in the configuration plane is presented, superimposed to the maximum equipotential energy level.The left boundary curve possesses zero slope since it encounters an impact, whereas the right boundary terminates at the maximum equipotential energy level; at that point the system possesses zero kinetic and purely potential energy.The localized motion of Figure 4 is the analytic continuation of the double-sided vibro-impact motion discussed earlier when the energy is decreased a small distance below the critical value h c1 .
With a further decrease of the energy below a second critical value h c2 , the system ceases to impact with its rigid barriers and performs localized NNM oscillations (Vakakis et al. [11]).The value h c2 is determined from the relation which for the above values of the system parameters is equal to 0•46666.For h Q h c2 the localized NNM is computed by solving equation ( 4) subject to the boundary conditions (6).Expressing the mode as yˆ(x) = gyˆ1(x) + O(g 2 ), a perturbation analysis can be used to compute the first order approximation.The resulting analytical approximation for h = 0•26666 Q h c2 is depicted in Figure 5, superimposed to the maximum equipotential energy level.The previous analysis proves the existence of a localized vibro-impact motion.As the energy decreases this localized motion is preserved in the single-sided impact regime, and ultimately degenerates to a localized NNM of the (smooth) dynamical system with no impacts.Hence, the localized vibro-impact periodic motions discussed in this section can be viewed as analytical continuations of strongly localized NNMs that are encountered in many classes of smooth dynamical systems [11].As mentioned earlier, the localized vibro-impact motion is orbitaly stable at all levels of energy.This result, along with additional ones concerning bifurcations of other (nonlocalized) types of vibro-impact motions, was deduced by means of direct numerical simulations and construction of Poincare´maps.

numerical results
The analytical results of the previous section were local in nature, establishing the existence of localized vibro-impact motions and providing asymptotic approximations for sufficiently small values of the coupling parameter g.The global vibro-impact response of the system will now be considered by numerically integrating the equations of motion (1) and constructing numerical Poincare´maps [24][25][26].
To construct the Poincare´maps we note that the motion of equation ( 1) generally takes place in the four dimensional phase space (x, x˙, y, y˙).Since the impacts are assumed to be purely elastic the energy h is conserved throughout the oscillation, and by fixing it to a specific level we can restrict the flow of the dynamical system to an isoenergetic three-dimensional space.If, in addition, we intersect transversely the three-dimensional isoenergetic flow by a two-dimensional cut-plane, the cut-plane defines a two-dimensional Poincare´map which can be used to study the global dynamics of the system.Choosing the cut plane as T: {x = 0}, the Poincare´map P(•) is defined as P: S:S, (y, y˙):P(y, y˙), where the Poincare´section is given by S = {x = 0, x˙q 0} + {energy h}.The additional restriction of positive x˙at the point of intersection was imposed in order that the Poincareḿ ap be orientation preserving [25].A period-one vibro-impact orbit pierces the cut-section S only once and corresponds to a fixed point of the map.The stability of the periodic orbit can be determined by examining near-by orbits of the Poincare´map: if the fixed point appears as a center and is surrounded by closed curves, the corresponding orbit is stable; otherwise it is unstable.Subharmonic vibro-impact orbits of higher periods pierce S in more than one point, whereas quasiperiodic orbits are represented by closed loops.Chaotic motions appear as randomly distributed points that occupy a definite region of the map.The numerical simulations were carried out by integrating equations (1) between impacts using the fifth-order Runge-Kutta SLATEC routine DDERKF.The energy-conserving elastic impacts were modelled by imposing new initial conditions after each impact, to reflect the abrupt sign change in velocity after the xor y-particle collide with their rigid constraint.The exact time instant when the impact occurs is of great importance in this numerical scheme, since it determines the points where the new initial conditions must be imposed.The time instants of impact were determined numerically by applying Newton's method to the numerical results of the integration close to each instance of impact.Special care was taken to choose a sufficiently small time steps in order to ensure accurate convergence of Newton's method and precise calculation of the instant of impact.This is particularly important for grazing (near impact) oscillations where the accuracy of the numerical integrations is important in order to construct correct bifurcation diagrams and study mode instabilities initiated by impact.

Nonlinear potential function
In Figures 6(a, b) we present the neighborhoods of the Poincare´maps close to the localized double-sided vibro-impact motions of equation ( 1) with potential energy (7), parameters a = b = d = 1, e = 1, g = 0•1, and energy h = 18 and 99.The Poincare´maps consist of centers surrounded by closed curves, indicating that the vibro-impact localized motions are stable even at high values of the energy.The global dynamics of the system are more clearly discerned from the Poincare´maps of Figure 7.In Figure 7(a) we depict the global map of the system for energy h = 2.We note the large central region that is dominated by the localized vibro-impact periodic motion.In addition, the system possesses a stable in-phase periodic orbit represented by the center at the top of the map; this motion corresponds to in-phase synchronous oscillations of the two particles satisfying x = y.We also note the large chaotic ring, composed of chaotic vibro-impact motions, that surrounds the central region.This 'stochastic sea' possess islands of regular subharmonic motions.For comparison purposes, in Figure 7(b) we present the Poincare´map of the same system but with no impacts.The system with no impacts possess a stable localized NNM that also dominates the central part of the map, a smaller chaotic region, and a stable in-phase NNM at the top of the map.Islands of subharmonic motion can also be detected in this case.

Linear potential function
To perform a systematic study of the effects of the vibro-impacts on the global dynamics of system (1) we performed a detailed computation of Poincare´maps for the case of the quadratic potential, with parameter values and varying values of the clearance e.In Figure 8(a-j) we present the results of the numerical simulations, ranging from a system with no vibro-impacts [Figure 8(a)] to a system with vibro-impacts and clearance e = 1.We note that as vibro-impacts start to occur, a region of chaotic motions starts to develop at the boundaries of the domains of influence of the stable in-phase and out-of-phase modes [cf.Figures 8(a-e)].This central chaotic region expands with decreasing clearance and surrounds the domains of influence of the two modes, where regular quasiperiodic motions continue to occur.For e 1 1•78 two saddle-node bifurcations occur and two pairs of stable-unstable vibro-impact modes are generated.In the Poincare´map of Figure 8(f) only the one pair of newly created modes can be observed, since the other pair is very close to the lower boundary curve of the map, and is difficult to observe.The small chaotic region surrounding one of the newly generated stable vibro-impact mode close to the origin of the axes is due to transverse homoclinic intersections of the invariant manifolds of one of the now unstable out-of-phase mode, and, hence, is a direct result of the mode bifurcation.A similar bifurcation of nonlinear normal modes was previously observed in a system with smooth nonlinearities [26].
With the clearance further reduced, the newly created mode becomes localized (it approaches the origin of the axes), and its domain of influence expands.At the same time, the regular regions surrounding the in-phase and out-of-phase vibro-impact modes contract [cf. Figure 8(g)], until the in-phase mode becomes unstable close to e = 1•4  [cf. Figure 8(b)] and the out-of-phase mode becomes unstable below e = z5/3.We note that the in-phase mode regains stability for e E 1•4; the bifurcations associated with the small instability region of this mode is discussed in detail below.The out-of-phase mode loses stability in a pitchfork bifurcation and remains unstable for smaller values of e.This behavior of the out-of-phase mode is consistent with the findings of reference [26] for a similar system with smooth nonlinearities.Further decreases of the clearance results in an expansion of the regular domain of influence of the localized modes, with the in-phase mode continuing to be stable [cf. Figure 8(i, j)].As a general observation we note that the region of chaotic motions contracts as the clearance is further reduced and the stable localized vibro-impact modes dominate the global dynamics of the system.In the following discussion we examine more closely the separate sequences of bifurcations associated with the instability of the in-phase mode, and the generation of the localized mode in the vibro-impact system.
We start with the bifurcations giving rise to the instability of the in-plane vibro-impact mode.These are complicated and generate additional periodic vibro-impact motions that have no counterparts in the corresponding smooth dynamical systems with no impacts.To investigate the series of bifurcations associated with the in-phase mode instability we performed detailed Poincare´map calculations in the neighborhood of that mode, as well as, analytical/numerical boundary value problem (BVP) computations to precisely determine the bifurcation points.We start with a discussion of the Poincare´maps in the vicinity of the in-phase mode, as depicted in Figure 9(a-f).
In Figure 9(a) (corresponding to e = 1•418) the in-phase mode (labelled 'A') is orbitaly stable and non-impacting.That is, for the level of energy h = 2 the in-phase synchronous oscillations of the two particles do not possess enough amplitude to reach the rigid barriers, and are, in fact, in-phase NNMs of the system.We also note the large region of chaotic vibro-impact motions surrounding the in-phase mode, and the plethora of co-existing regular quasi-periodic motions.As e is reduced to 1•417 [cf. Figure 9(b)], the in-phase motion 'A' is still stable and non-impacting, and an additional stable periodic motion appears (labelled as 'B') on the right of the in-phase mode.As shown below, this new orbit corresponds to a double-impacting synchronous motion where both xand y-coordinates undergo single-sided impacts; clearly this type of orbit cannot be realized in the corresponding smooth dynamical system.Moreover, it will be shown that this new motion is generated through a Saddle-node bifurcation, and that there exists an additional unstable orbit that cannot be discerned in the Poincare´map; also, due to symmetry, there exists an additional stable-unstable pair of similar orbits lying on the left side of the in-phase mode, which is not depicted in Figure 9(b).At e = z2 [Figure 9(c)] the in-phase mode 'A' starts impacting and loses stability, whereas, the stable motion 'B' expands its surrounding domain of influence which is composed of regular quasi-periodic orbits.It turns out that at least four unstable branches of vibro-impacting periodic motions converge to the in-phase mode at its point of instability, including the two unstable counterparts of motion 'A'.As discussed below, the in-phase impacting mode remains unstable until e 1 1•40099, when it coalesces with two new unstable branches of vibro-impact motions and regains orbital stability.This can be realized from Figures 9(e) for e = 1•398, where both the in-phase impacting 'A' and the periodic motion 'B' are stable.A further reduction of e to 1•390 reduces the domain of influence of motion 'B' [cf. Figure 9(f)], and at e 1 1•38813 the orbit 'B' ceases to exist when it coalesces with another unstable vibro-impact motion in a saddle-node bifurcation.Hence, the small region of instability of the in-phase motion is associated with a complicated series of bifurcations, which we now proceed to examine in detail using BVP formulations.
The BVP formulation is based on the realization that, since the motions of the two particles between instants of impact is linear, they can be explicitly determined as The undetermined amplitudes and phases can be evaluated by imposing appropriate boundary conditions at instants of impact for each of the vibro-impact orbits.Then, by construction the solution ( 19) will be periodic in t.The first time-periodic solution to be examined is a double-sided impacting synchronous motion, where both xand y-coordinates undergo single-sided impacts.For this type of motion the appropriate boundary conditions that must be imposed on equations (19) are where T is the period of the motion, and the origin of the t-axis is chosen at the time instant when the x-coordinate impacts.Combining equations ( 19) and ( 20), we obtain the set of algebraic equations where v 1 = za and v 2 = za + 2g.Provided that the matrix of coefficients is nonsingular, the above relation provides a unique solution for the coefficients C 1 , C 2 , f 1 , and f 2 and determines the solution (19).Once the solution is determined, the period T of the motion is computed by imposing the additional condition that the total energy of the motion is conserved and equal to h.This completes the calculation.Singularities of the matrix of coefficients indicates points where analytic continuation of this particular branch of solutions cannot take place, i.e. bifurcation points of the solution.Similar BVPs can be formulated for determining the branches of solutions and bifurcation points of other types of vibro-impact modes, such as, in-phase synchronous motions [solution (19) with boundary conditions x(0) = −e, x(T/2) = e, y(0) = −e, y(T/2) = e], out-of-phase synchronous motions [boundary conditions x(0) = −e, x(T/2) = e, y(0) = e, y(T/2) = −e], and 'loop' vibro-impact modes (see discussion and results below, corresponding to boundary conditions The numerical solutions of the aforementioned BVPs enable us to fully understand the series of bifurcations that generate the instability of the in-phase mode, as depicted in the Poincare´maps of Figure 9.In Figure 10 we provide the complete bifurcation structure that gives rise to this instability for the system whose Poincare´maps were depicted previously.In the plots of this figure we graph the values of y and y˙vs e for each branch of periodic solutions at the time instant when x = 0 and x˙q 0. There are four families of periodic solutions in this bifurcation structure.For e e z2 the in-phase mode 'A' is non-impacting and stable.At e = z2 the in-phase mode becomes where it coalesces with four unstable branches: two of loop' vibro-impact modes 'C', and two of double-sided synchronous modes 'B'.The mode 'B' (which were observed in the Poincare´plots of Figure 9) are generated in two saddle-node bifurcations and exist in stable-unstable pairs up to e = z2.At e = (zh/v 1 ) sin (pv 1 /2v 2 ) 1 1•40099 the in-phase mode regains stability in a pitchfork bifurcation and generates two unstable branches of solutions 'D'.These are in-phase synchronous motions with coinciding endpoints with the in-phase mode, and cease to exist at e = (1/v 1 )zh/{1 + 2 tan 2 [(p/2) + (pv 1 /2v 2 )]} 1 1•38813, when they coalesce with the stable modes 'B' in saddle-node bifurcations.The bifurcation structure of Figure 10 agrees completely with the numerical simulations of Figure 9, and leads to a full understanding of the dynamic phenomena associated with the instability of the in-phase mode.
The bifurcation results depicted in Figure 10 reveal the richness of the dynamics of the system.Indeed, in the region where the symmetric in-phase NNM looses stability it gives rise to stable or unstable motions with differing impacting patterns.Solutions C (cf. Figure 10) are of particular interest due to their loop-like appearance, indicating nontrivial phase differences between the two position coordinates; moreover, these motions correspond to anisochronous and single-impacts for x and y.Similar anisochronous and single-impacts occur in solutions B, but for these motions trivial phase differences occur between the coordinates.Nontrivial phase differences between x and y also occur in solutions D, but these represent unstable free oscillations.We point out that it would be of considerable interest to analytically study the degenerate bifurcations associated with the loss of stability of the in-phase mode at e = z2.These bifurcations indicate that the initiation of impacts (grazing) for the in-phase mode is destabilizing, although the mode regains stability when the non-linear impacting effects become stronger for e E 1•40099.Additional bifurcations are expected to occur involving motions in the configuration plane of increased complexity.Such bifurcations can be studied using the techniques outlined herein, by imposing different boundary conditions and impacting restrictions for the motions to be studied.
Similar BVP computations help us understand the bifurcations giving rise to the localized vibro-impact mode that appears in the Poincare´maps of Figure 8.As mentioned previously, due to symmetry there exists two such localized modes, one localizing in the x-coordinate, and the other in the y-coordinate.Here we will analyze only the mode localizing in the x-coordinate, but a similar calculation can be performed for the other mode.To this end, we seek localized solutions governed by equation ( 19) between impacts, with only the x-coordinate impacting, and boundary conditions In Figure 11 we depict the numerical solutions of the BVP.The solutions are depicted as graphs of y vs e for each periodic solution at the time instant when x = −e.We note that the localized mode is generated at e 1 1•79 in a saddle-node bifurcation.The unstable mode of the bifurcation has the same impacting properties with the localized mode and ceases to exist at e = z5/3 1 1•291 when it coalesces with the out-of-phase vibro-impact mode in a pitchfork bifurcation.After this bifurcation the out-of-phase modes loses stability, and remains unstable for lower values of e.It is the transverse intersections of the invariant manifolds of this unstable mode that give rise to the small chaotic regime in the central region of the Poincare´plot of Figure 8(f).
In the previous sections we analyzed some of the dynamic phenomena occurring in the direct vibro-impact problem of a two-degree-of-freedom linear or nonlinear oscillator undergoing purely elastic impacts.We considered both analytical and numerical techniques, which enabled us to study local and global dynamics, and bifurcations generating vibro-impact mode localization and mode instabilities.We now turn our attention to the inverse problem associated with the two degree-of-freedom system.This problem seeks the systems that under vibro-impact conditions generate an a priori specified trajectory in the configuration plane.The solutions of the inverse problem provide interesting insight on the class of dynamical systems that are capable of producing required pre-determined vibro-imapct orbits.

THE INVERSE PROBLEM
The inverse problem associated with the vibro-impact oscillator of Figure 1 is posed as follows.Consider the equations of motion (1) and assume that the trajectory of the system in the configuration plane follows a specified orbit: What is the required potential function P(x, y) required for such an orbit to exist?Moreover, under what conditions does such a potential function exist and when is it unique?In what follows we will discuss in detail a method of solution to the inverse problem and will provide some numerical applications.Although the inverse problem considered herein addresses the two-degree-of-freedom vibro-impact oscillator, it can be easily generalized to higher dimensions.In addition, in the following exposition the inverse problem is formulated for a motion composed of two two-sided impacting branches.Generalization to other types of vibro-impact motion can be similarly performed.
Considering again the equations of motion (1), we assume that the orbit of the system in the configuration space is composed of two double-sided impacting branches.We further assume that each branch is sufficiently smooth to be expanded in power series in x for the entire range of the vibro-impact motion: Assuming that the coefficients a k and b k in the above series are specified, we seek the potential energy function P(x, y) that is required for such a motion.Expressing the potential function in the general form c ij x i y j for i, j = 0, 1, 2, . . ., (24) the solution of the inverse problem reduces to the determination of the coefficients c ij in the above expression.Since the branches (23) represent double-sided vibro-impact orbits, they must satisfy the differential equation ( 4) and the set of compatibility conditions (5).Substituting equations ( 23) and ( 24) into equations ( 4) and ( 5) we obtain the algebraic equations 2 complemented by the compatibility conditions These equations govern the sought coefficients c ij , and their solution provides the answer to the inverse problem.We note that conditions (25b) do not involve the unknowns c ij , and represent mere geometric restrictions on the assumed vibro-impact branches (23).
Considering relations (25a), we make the initial observation that, if we truncate the infinite summations to only the first N terms, by matching coefficients of respective powers x k we obtain an underdetermined set of 2N algebraic equations for N 2 unknowns.Hence, in general, the inverse problem appears to possess non-unique solutions.
Due to the difficulty of analyzing the general equations ( 25), it will be necessary to consider specific geometric forms for the branches (23) in order to proceed with concrete solutions of the inverse problem.Hence, as an application, we consider a vibro-impact time periodic solution composed of the two branches y = yˆ1(x) = a − (a/e 2 )x 2 and y = yˆ2(x) = −yˆ1(x) = −a + (a/e 2 )x 2 for −e E x E e. (26) In terms of the previous notation we have, a 0 = a, a 1 = 0, a 2 = −(a/e) 2 , a k = 0, k e 3, and b 0 = −a, b 1 = 0, b 2 = (a/e) 2 , b k = 0, k e 3. Substituting these values into relations (25a), and considering terms only up to O(x 3 ) we obtain the following matrix relation governing the leading coefficients of the potential function: where the various coefficients are defined in Appendix A. Assuming that all coefficients c i j , i e 4, are set equal to zero, the above matrix relation governs exactly the inverse problem and involves no approximation.The vector on the right-hand-side depends on an infinite number coefficients c ij (cf.Appendix A), and, hence, the inverse problem possesses non-unique solutions.Furthermore, it can be realized that problem (27) can be decomposed into two separate subproblems, one governing the coefficients c 1j and c 3p with j e 1, p e 0, and the other governing the coefficients c 0k and c 2p with k e 2, p e 0. The solutions of these subproblems are summarized below.
(1) Subproblem I governing c 1j and c 3p with j e 1, p e 0. , If c 1j = 0, j e 1, then, either, c 3p = 0, p e 0 or, c 30 = D + 4 /Q + 0 , and c 3j , j e 1 satisfy s In this case we cannot have the trivial solution c 0j = 0, j e 2 and these coefficients satisfy the remaining coefficients satisfy the relation All matrices in the solutions above are invertable, provided that the parameter a of the trajectory ( 26) is nonzero.This eliminates the possibility of a trivial trajectory in the configuration plane.
A solution to the inverse problem can be formed by combining one of the two solutions of subproblem 1 with the solution of subproblem 2. We note that since the above solutions involve an infinite number of coefficients, one can obtain an infinity of solutions to the inverse problem by assigning arbitrary values to all the coefficients except the ones required to satisfy the above relations.Perhaps the simplest of these solutions is obtained by setting c 1j = 0, j e 1 and c 3p = 0, p e 0 (solution of subproblem 1), and c 02 = 2h/(2a 2 + e 2 ), c 0j = 0, j e 3, c 20 = h/2(2a 2 + e 2 ), c 2j = 0, j e 1 (solution of subproblem 2).Hence, we obtain a linear system with clearance nonlinearities, which is predicted to possess the vibro-impact periodic solution (26).In Figure 12 we superimpose the theoretical orbit (26) and the numerical orbit obtained by direct numerical integration of the vibro-impact equations of motion.The agreement between theory and numerical integration validates the previous inverse analysis.
Summarizing, in this section we formulated an inverse vibro-impact problem by computing the required potential function of a system in order to produce a specified trajectory in the configuration plane.Although the analysis herein was performed for a double-branch orbit of parabolic shape, the technique can be applied to a more general class of orbits with single-or double-sided impacts.An interesting aspect of the inverse problem analyzed is the degeneracy of its solutions: indeed, it was found that an infinity of vibro-impact systems can produce the specified trajectory, a result which is due, perhaps, to its relatively simple parabolic shape.It is conjectured that the degeneracy of solutions of the inverse problem decreases as the shape of the trajectory becomes more complicated.This is an interesting topic for some future research.

DISCUSSION
We studied direct and inverse problems associated with vibro-impact oscillations of discrete systems.The usefulness of the presented analytical/numerical techniques stems from the fact that the majority of previous vibro-impact works were direct problems and centered on dynamical systems with a single impacting coordinate.By contrast, in this work we examined systems with two impacting coordinates, undergoing single-or double-sided impacts.The analytical methods employed were a NNM-type analysis and a BVP formulation.In the NNM-type analysis the displacement of one of the particles was expressed as a function of the displacement of the other particle, and the functional relation of the motions of the two particles was asymptotically approximated.The BVP enabled the computation of various branches of bifurcating time-periodic solutions with different impacting properties.Numerical results were obtained by direct integrations of the equations of motion and then used to construct Poincare´maps.
The class of vibro-impact systems considered possesses rich nonlinear dynamical behavior, including vibro-impact localized and nonlocalized time-periodic motions (which can be considered as analytic continuations in the vibro-impact regime of NNMs of corresponding smooth dynamical systems), as well as, complicated bifurcation sequences giving rise to new types of single-or double-sided impacting motions, instabilities, and chaotic responses.These bifurcations can be studied using appropriate BVPs and numerical Poincare´maps.
We also formulated inverse vibro-impact problems, whereby we sought systems that would produce specified orbits in the configuration plane.As expected from experience with similar types of inverse problems for other engineering areas, the solution of the vibro-impact inverse problem is generally non-unique, since it can be reduced to an underdetermined set of algebraic equations with multiple infinities of unknowns.The inverse problems discussed in this work can find applications to tracking problems were the design task calls for the responses of the coordinates of a system to track specified paths.Such problems arise often in kinematics and dynamics of mechanisms, and in other engineering applications.
The presented techniques can be extended to vibro-impact oscillators with more than two degrees-of-freedom, although in such cases the computational and numerical effort is expected to increase.Moreover, the outlined analysis can be used to study bifurcations associated with the inverse problem, i.e., to investigate degeneracies in the class of dynamical vibro-impact systems that produce a specified orbit in configuration space.That would be an interesting problem in inverse nonlinear dynamics.In addition, the mentioned analytical and numerical techniques can be applied to study the bifurcation structure of the periodic solutions of multi-degree-of-freedom vibro-impact systems, and to establish regions of regular and chaotic responses.

Figure 1 .
Figure 1.The vibro-impact oscillator: (a) system configuration, (b) the limiting process defining the elastic impact forces as n:a.

Figure 2 .
Figure 2. Vibro-impact motion with imapcts only in the x-coordinate: (a) motion in the configuration plane, (b) joining of branches at the point of elastic impact.

Figure 6 .
Figure 6.Poincare´map of the neighborhood of the localized vibro-impact motion of a system with a = b = d = 1, e = 1, g = 0•1, and energy (a) h = 18, and (b) h = 99.

Figure 10 .
Figure 10.The bifurcation structure associated with the instability of the in-phase mode.Graphs of (a) y and (b) y˙vs e for each periodic solution at the time instant when x = 0 and x˙q 0: --, stable; -----, unstable motions.

Figure 11 .
Figure 11.The bifurcation structure associated with the generation of the localized vibro-impact mode.Graphs of y vs e for each periodic solution at the time instant when x = −e: --, stable; ---, unstable motions.

,
If not all c 1j , j e 3 are zero, then the coefficients are computed by the relations, Subproblem II governing c 0k and c 2p with k e 2, p e 0.

Figure 12 .
Figure 12.Comparison between theoretical and numerical solutions of the inverse problem.