Normal modes for piecewise linear vibratory systems

. A method to construct the normal modes for a class of piecewise linear vibratory systems is developed in this study. The approach utilizes the concepts of Poincare maps and invariant manifolds from the theory of dynamical systems. In contrast to conventional methods for smooth systems, which expand normal modes in a series form around an equilibrium point of interest, the present method expands the normal modes in a series form of polar coordinates in a neighborhood of an in variant disk of the system. It is found thal the normal modes, modal dynamics and frequency-amplitude dependence relationship are all of piecewise type. A two degree of freedom example is used to demonstrate the method.


Introduction
In linear vibration theory, modal analysis is a powerful technique that provides a means of reducing the system dynamics to uncoupled individual oscillators. Through the mode shapes and corresponding modal dynamics, the behavior of a linear vibratory system can be completely understood. For this reason, researchers have tried to extend the concept of normal modes to nonlinear systems. Rosenberg took the first steps in this direction [22][23][24]. Since then, many studies have appeared in the literature discussing the existence [5,15,19,40], stability [1,7,10,17,21], and construction [8,[18][19][20]39] of nonlinear normal modes, and the possibility of some type of nonlinear modal analysis [1,32,[35][36][37][38]. Other studies have included the investigation of nonlinear nonnal modes for continuous [2,3,12,13,16], nonconservative [6,28], and internally resonant systems [14]. Shaw and Pierre introduced a constructive method for nonlinear normal modes using invariant manifolds [26,[28][29][30]. This approach is applied here to the case of piecewise linear systems.
While the idea of nonlinear normal modes has drawn much attention, virtually all of the previous works deal with smooth systems. In contrast, nonsmooth systems, in particular piecewise linear (PWL) systems (see, e.g., [4,9,11,25,31]), received little attention in regards to normal modes. However, PWL vibratory systems are useful models for many practical vibratory systems. For example, conservative systems with clearance, backlash, or piecewise linear springs and dissipative systems with simple coulomb friction are all of PWL type.
Recently, Zuo and Cumier reported a study on the modal motions of a class of PWL systems, where the switching hyperplane passes through the origin [41]. As one will see in the sequel, their class is a special case of the present class of PWL systems to be studied, where the switching hyperplane does not necessarily pass through the origin. In particular, the present systems may not satisfy the property of positive homogeneity defined in [41]. As a result, the modal frequency of the present systems depends on the amplitude, in contrast with the systems in [41] whose modal frequency is independent of the amplitude. Also, the construction of the normal mode manifolds and their attendant dynamics were not pursued in [41]. However, they did include systems with gyroscopic terms, which are not considered here.
The objective of this paper is to demonstrate a method for the construction of nonlinear normal modes and their attendant dynamics for a class of general conservative PWL systems with N (finite) degrees of freedom. To this end, we utilize the approach of invariant manifolds and asymptotic expansions. Here the invariant manifold, and hence the normal mode, is a one-parameter family of periodic orbits in the state space. The general procedure involves building a Poincare map with the switching hyperplane as the Poincare section in state space, determining the fixed points of the map, and finally constructing the normal mode by asymptotic expansions. Unlike the case for a smooth system, where the normal mode is expanded in a series form near an equilibrium point of the system, the normal mode for a PWL system is expanded in a series form using polar coordinates in a neighborhood of an invariant disk of the system. This invariant disk is precisely the corresponding normal mode manifold of the linear system for energies below those at which the switching takes place.
The paper is organized as follows. We begin in Section 2 with a description of the class of PWL systems in which we are interested. The general system will then be transformed into a canonical form for easy formulation. Next, in Section 3, the Poincare map is obtained implicitly. The procedure of constructing the normal mode manifolds for PWL systems, including the determination of the fixed points of the Poincare map, using asymptotic expansion is described in Section 4. Their modal dynamics are established in Section 5. In Section 6, a two degrees of freedom example problem with clearance is given to show some of the calculations involved with using this method and a comparison is made with simulation. Some conclusions are drawn in Section 7. never leave the first region, the system is nothing but the first linear subsystem and is well understood. As the motion amplitudes become large, solutions pass through the switching hyperplane and go into the second region, in which case both linear subsystems will be involved and the system behavior is no longer simple. Many complicated behaviors involved in PWL systems have been reported [4,31,34]. In this paper, we are not going to discuss the general behavior of such PWL systems, but only their normal mode motions.

The Piecewise Linear System and Its Canonical
System (1) represents a large class of physical vibratory systems such as systems with clearance, impact, or piecewise linear springs. However, we must point out that it can only represent systems with one clearance, but not those with two or more clearances, which should be modeled by PWL systems with more domains. Normal modes for PWL systems with more than two domains are not considered in the present paper. Also, systems with dry friction do not belong to the class of systems in equation (1) because they are nonconservative oscillators.
Before proceeding to the next step, we make some assumptions about system (1). By premultiply,ing equation (1) with M-1 , the system can be rewritten as (ii) N modes in both linear subsystems are distinct. (2) Assumption (i) says that the right-hand side of (2) is continuous, so that when written in state equation form, i.e., as a system of first order ODEs, the system will possess a continuous vector field. Although the vector field is not C 1 due to its PWL nature, it will satisfy a local Lipschitz condition. Hence, the existence and uniqueness of solutions of (2) are guaranteed. Assumption (ii) is for simplicity and can be removed easily. It is needed to insure that each mode of the linear subsystems is two dimensional in the state space. Without this assumption, we might have a four or more dimensional surface for a certain mode:.. Also, it can be shown that assumption (i) together with the negative definiteness of K and K imply that the origin is the only equilibrium point for system (6) (this can be proved but is not obvious).
With these assumptions in hand, we transform system (2)  denotes the Euclidean 2-norm. Then it is easy to verify that Q = P D 1 1 2 U is as required.
After the transformation x = QT z, system (2) takes the form X -{ ~: + QTbt ::: :: : ~ where K = QT At(QT)-l and K = QT A2(QT)-1 . By our choice ofQ, both K and K are real symmetric negative definite matrices. From the continuity of the vector field, we have for (5) or in a more compact form (6) where Xe = dk-Ik. Moreover, by the symmetry of k, k must be of the form Either equation ( 5) or equation ( 6) can be considered as the canonical form of the PWL system (3) with x; as a parameter, and when x; = 0, they reduce to a purely linear system. Also, for d = 0, the system reduces to that discussed i!l [ 41].
It will be clear in Section 4 that as far as K is negative definite, x; is not restricted in size in our analysis. However, since the method of asymptotic expansion is used, the magnitude of x; might affect the region of validity of the expansions. This will be discussed in Section 6.

The Poincare Map
In this section, we derive the Poincare map for the PWL system in canonical form (6), which will be used to construct the nonlinear normal mode in the next section. The switching hyperplane provides a natural Poincare section. To this aim, we make use of the symmetry of K and k, and the linearity of the subsystems. Since both subsystems are linear, we can obtain their solutions in closed form which enables one to obtain an analytical, although implicit, expression for the Poincare map.
COSwzt, ... , COS WNt), S(t) = diag(sinw1t, sin W2t, ... , sin WNt), n = diag(wl , w2, ... , WN ), and a = [at az . . . aNJT, b = [bt b2 . . . bn]T E i'RN are coefficients determined by the initial conditions. Setting t ::::: 0 in equations (7) and (8) Inserting (9) back into (7) and (8) one obtains (9) or in an integrated form where is a 2N x 2N matrix composed of four N x N matrices. Note that we have used the commutative property of two diagonal matrices to simplify the expression. Furthermore, we can put T(t) in the clearer form It is not difficult to see that T( t) is invertible and that its inverse is given by 0 vr (12) For the solution of the second linear subsystem, a similar procedure gives where we= [xe o]T and T(t) has a definition similar to T(t) and can be obtained from T (t) by replacing wi, Vi with Wi, fJi . Note that forK,= 0, T(t) = T(t) andxe = 0, and equation (13) reduces to equation (10).
We are now in a position to construct the Poincare map. Define the surfaces of section (here, hyperplanes) 2: = { w E ~2N : x 1 = d and Y1 > 0} and 2:* = { w E ~2N : Xt = d and Yt < 0} · The Poincare map F : ~ --+ ~ is thus defined as (14) Vw E :E, where thi is the time required to take w to a point w * on~* and th2 is that required to return from w* back to :E. In general, both tht and th2 satisfy transcendental algebraic equations and only approximate solutions are possible (see [11,27]). Note that although the system is piecewise linear, the corresponding Poincare map is smooth. expansions is provided in this section. Since the systems considered here are conservative, it is well known that each normal mode is a collection of periodic motions [28]. In other words, we are seeking invariant manifolds in state space that are composed of one-parameter families of periodic orbits. Therefore, when viewed in the Poincare map, the invariant manifold is a one-parameter family of fixed points.

Normal Modes for the PWL System
In the linear case, that is, when K = 0 in system (6), a normal mode is a two-dimensional plane in state space, which intersects the switching hyperplane along a straight line. This line constitutes the one-parameter family of fixed points for the Poincare map ( 14) for the corresponding normal mode. To fix the ideas, let us take the first mode as an example. Let the corresponding eigenvector be written as Then the 2-D invariant plane for this mode is: which is composed of a family of ellipses that are periodic solutions with different amplitudes but the same frequency Wt, as shown in Figure 1. Figure 2 shows the corresponding fixed point family for the Poincare map, given by  Figure 2). Correspondingly, the invariant plane will become an invariant manifold consisting of two pieces. One piece is an elliptic disk near the origin which represents the small amplitude vibration and is tangent to the switching hyperplane at Wt. This piece is a portion of the linear normal mode of the first linear subsystem since near the origin the system is linear. The other piece represents large amplitude vibration and involves both linear subsystems. This piece again contains two parts, one in the region xI < d and the other in XJ > d. The two parts must match along the Poincare section in the sense that they meet at the curve of fixed points. Figure 3 is a qualitative sketch of such an invariant manifold.
In order to obtain the normal mode, we only need to construct the large amplitude piece since the elliptic disk is already known. Because it is generally impossible to obtain the entire invariant manifold in closed form, we shall employ the method of asymptotic expansions which allows us to obtain an approximate invariant manifold in the neighborhood of the elliptic disk. For this purpose, it is natural to use polar coordinates, in contrast with the Cartesian coordinates used by Shaw and Pierre [28].
There are several steps involved in the procedure for constructing a normal mode for a PWL systems, and the procedure is exactly the same for each mode. First, an approximation to the curve of fixed points is determined. Next, using the curve obtained as a matching condition, the large amplitude piece of the invariant manifold, which is in the neighborhood of the elliptic disk and belongs to the region Xt < d, is obtained in a series form of polar coordinates. Finally, by matching the curve from the first step again, and by making use of another set of polar coordinates, the other portion of the invariant manifold belonging to the region x 1 > d is expanded in the neighborhood of Wt. Then, by gluing together these two pieces and the elliptic disk from the linear normal mode, we are able to construct an approximation of the entire nonlinear normal mode for the PWL system. In what follows, we will explain these steps in some detail. Again, we take the first mode as the illustrative example, where the extension to the other modes is straightforward.
The fixed points of the Poincare map are those on ~ satisfying subject to (16) where {·}i represents the i-th element of{·}. In general, equation (15)  In order to determine the one-parameter family of fixed points corresponding to the first mode, we can assume a solution in a series in terms of the velocity of the first degree of freedom at the switching point, Yt, up to m-th order as the following . (17) i =2, ... ,N; (18) i = 2, ... ,N; where 0( m + 1) denotes the higher order terms and the f.Li/s, Vi/s, and Ti/s are unknown coefficients, depending on the parameter "'' which are to be determined. Recall that ui = dvit/vu and [d u2 . . . uN]T is an eigenvector associated with the first mode of the first linear subsystem. Note here that we have assumed vn -1-0. In general, if vli = 0 for some 1 < i < N, then the i-th nonnal mode of the PWL system is the same as that of the first linear subsystem, i.e., a 2-D plane. This is because on this mode, Xt -0 and hence it will never exceed d. Note also that as a 8--+ 0 in equations (17) through (21), the assumed solution reduces to Wt, which is a fixed point for all values of K. Substituting equations (17)- (21) into equations (15) and (16) Before moving to the next step, we transform the coordinates ( x1, y 1 ) to a polar form ( r, e) defined by As a result, a periodic solution of the first mode of the first linear subsystem can be represented simply as r = constant when projected onto (x 1 , yt) plane. Then when r < d, the normal mode can be expressed as which is the elliptic disk mentioned above. Generally, in terms of (r, 0), the normal mode for Basically, our aim in this section is to find approximate solutions for fij(r, B) and gij(r, B). Clearly, these functions must satisfythe continuity condition on the boundaries.
Next, consider the region r > d and r cos() < d, which is outside the elliptic disk but still governed by the first linear subsystem. Let P( 8) be a fixed point on 2] with y 1 coordinate being 8 and others being given from the first step in equations (18) and (19). Then there must be a periodic trajectory passing through P. 1 Let Po ( 8) be a point on the trajectory in this region whose coordinates are where the rj's are also functions of() to be determined. In other words, Po is a point on the normal mode. Reference to Figure 4, which is the state space projected onto the ( x 1 , y 1 ) plane, may clarify matters for the reader at this point. If 8 = 0 in equations (25) where the tj 's are again functions of e to be determined.
With Po as the initial point w(O) and substituting expressions (25)-(28) into equation (10), we obtain w(thB) = P or (29) By expanding T(the) with respect to 8 = 0 up to order m and matching coefficients of like order in 8 in equation (29), one can solve for the a ij 's, f3i j 's, r j 's and t j 's. Then, by inverting  the power series of equation (27), 8 can be obtained as an approximate function of r and e. Entering the resulting 8(r, (}) into equations (25) and (26), and comparing the result with equations (23) and (24) (35) As in the procedure described above, the expansion coefficients Cxij 's, ~ij 's, f j 's and ij 's are unknown functions of ¢ determined by the matching condition (36) Transforming (f, ¢) back to (r, B), and from equation (34), we solve for 8 as a function of (r, B). Then, comparing equations (32) and (33) with equations (23) and (24) where one should note that 8(r, B) here is a different function from that given in equations (30) and (31) since they are in different regions of (r, B). This is the final step and it completes the procedure for generating the normal modes of the PWL system (6).
This section is ended with some remarks on the order of approximation. Since we are expanding the invariant manifold around an invariant disk, solutions will be accurate close to the disk. Here the measurement of closeness is 8, the YI coordinate of a fixed point on the Poincare section ~. Throughout the paper, the order of approximation is in terms of 8.
One may relate 8 to the energy level of the modal motion above and near the switching hyperplane. A larger 8 represents a higher energy level on the invariant manifold. However, their relationship cannot be obtained in closed form since a closed form expression for the invariant manifold is not available.

THE EXISTENCE AND UNIQUENESS OF THE NORMAL MODES
The problem of determining the existence and uniqueness of nonlinear normal modes has been the subject of many efforts since the concept of nonlinear normal modes was introduced [5,15,19,40]. For smooth nonlinear systems, the existence problem is well understood. As for uniqueness, the Lyapunov center theorem states that the number of nonlinear normal modes is the same as that of linear modes under nonresonant conditions [30]. In the case of internal resonance, Nayfeh and Chin reported in [14] that the number of nonlinear normal modes may exceed the number of linear ones. For PWL systems, both the existence and uniqueness problems remain open. In [41], Zuo and Cumier numerically obtained the normal modes for a special class of PWL systems and discussed the bifurcation of these modes. However, they did not address the existence problem, nor did they investigate internally resonant cases.
Although we do not have a definite answer to the problem of existence and uniqueness of normal modes for PWL systems either, some evidence (a necessary condition) to the existence of the normal modes has been determined. Specifically, the problem is related to the solution of the fixed point equations (15) and (16). If, for a given linear mode, there is a unique one parameter family of fixed points satisfying equations (15) and (16), then the nonlinear normal mode exists and is unique for the given mode. Suppose this is indeed the case. Then there must be exactly 2N independent equations in the set of 2N + 1 simultaneous equations (equations (15) and (16)), i.e., essentially, only one of the equations in (15) and (16) depends on the others. This is verified by the illustrative example following in Section 6. In other words, the Jacobian matrix of equations (15) and (16) at the given linear mode with K = 0 and 8 = 0 should have rank equal to 2N.It is shown in Appendix A that this is generally true for nonresonant cases. In lieu of equations (15) and (16), a more convenient but equivalent form of fixed point equations is used in Appendix A. Under resonant conditions, the rank of the Jacobian matrix will be less than 2N, indicating the possibility of non unique normal modes for a given linear mode. Interestingly, one should note that the PWL system is linear locally and hence the normal mode is unique near the origin, whether resonant or not. Over the switching hyperplane, however, the normal mode can split into several pieces.

Modal Dynamics for the PWL System
With the normal modes in hand, the next step is to establish the dynamics of the normal modes. As expected, the modal dynamics are also of piecewise type.
Again, the first mode is considered here. Its local coordinate system is taken to be the polar form (r, B) defined by equation (22). Hence · · e 1 · · e r = Xt cos + -Yt sm ; W} (39) e . 1 1 rcosB+h(r,8))}t +K(rcosB-d) forrcosB > d (41) where VI = [1 ( uz/ d) . . . (UN/ d) ]Tis the eigenvector of K associated with eigenvalue -wi and fi(r, B) = [0 fzi(r, B) . .. !Ni(r, B)]T fori = 1 and 2 are as defined in equation (23) and their approximate solutions are given in equations (30) and (37). Also, the reader should recall that { ·} 1 stands for the first element of { ·}. Therefore, -wyr cos e + kf f2(r, B)+ K,(r cos B-d) for r cos() > d where ki is the first column of K. Inserting expressions (41) and (42) into equations (39) and ( 40), one obtains the dynamics of the first mode as forr>dandrcosB<d. (44) As expected, the modal dynamics are of piecewise type and can be regarded as a one degree of freedom nonlinear oscillator. Also, it is easy to see that the right-hand sides of ( 42)-( 44) are continuous. Note that the first piece (r < d) is for motions of the first linear subsystem, which govern for energy levels below that required to exceed the switching hyperplane, while the second and third pieces describe the motion in the two linear domains for higher energy levels. Due to the piecewise nature of the equations, however, conventional methods in determining the frequency-amplitude relationship of nonlinear oscillators fail to apply here. Conventional methods make use of perturbation techniques which require the system to be smooth, which is not the case here. Despite these difficulties, we can still determine the frequency-amplitude relationship of the modal dynamics from the analysis presented in Section 4, rather than deal with equations ( 43) and ( 44) directly.
Consider the motion on the normal mode. When the amplitude is sufficiently small that the vibration stays in the disk region r < d, the system is linear with constant frequency w 1 .
If the vibration is large enough, it will hit the switching hyperplane at a fixed point on the Poincare section ~. Each fixed point is parameterized by the y 1 coordinate and is given by equations (17)- (19) in Section 4. For each Yt = 8, we have one such oscillation whose period is equal to tp = thi + thz, which is also parameterized by 8, as given in equations (20) and (21). Specifically, where tpj = Ttj + 72j. Thus, the frequency can be easily calculated from Equation ( 45) is the frequency-amplitude relationship of the first modal dynamics in terms of 6. We can also convert equation ( 45) to the usual form of relationship between w 1 and A 1 , the peak displacement of XJ. To do this, let iiJ = [x y]T be a fixed point on~ with y 1 = 6.
The peak displacement occurs when the velocity is zero, i.e., Yl = 0. With iJJ as the initial point, we thus seek the value of x1 at Yl = 0. First, the time tho bringing Yl = 8 to Yt = 0 is obtained in terms of a power series in 6 by solving approximately up to the m-th order of 6. Then the peak displacement of x1 is simply given by (46) which is also a series form of 6. Using equation (46), one can obtain 6 as a function of At in a series in terms of 8. The frequency-amplitude dependence in terms of A 1 is then expressed as These ideas are now demonstrated by an example.

Example: A Two Degrees of Freedom System with a Clearance
(47) The system to be considered is shown in Figure 6, which is an undamped, two degrees of freedom system with a clearance. Except for the clearance, no other nonlinearities appear in the system. Following the procedure provided above, the normal modes for the system, the associated modal dynamics and the frequency-amplitude relationship will be given for both normal modes. Numerical simulations are also carried out to compare with the theoretical results.
It is important to point out that most of the manipulations encountered here, such as power series expansions, inversion of power series, and solving for the unknown coefficients of the normal modes, can hardly be done without the aid of computer symbolic manipulators. Indeed, much of the work for this example was done by Mathematica ™ running on a Macintosh II ex computer. All of the numerical simulations presented were carried out using MatLab on a Sun station.
Assume that the springs are unstressed at the equilibrium position, where the displacements XI and x 2 are zero. In addition, all the values for inertia, stiffnesses, and clearance magnitude are taken to be unity except the free spring whose stiffness is assumed to be /1,. The equations of motion are, therefore, In terms of the canonical form given in equation (6), the equations of motion are expressed with the following definitions: The eigenpairs of K are and the hitting times as thl = 1"110 + Tt2D 2 + T13D 3 + 0(4); th2 = 27r + T2t0 + T220 2 + 1"330 3 + 0(4), where the approximations are done up to third order. The coefficients in equations (50)-( 53) have to be chosen to satisfy the fixed point equations (15) and (16). Instead of solving equations (15) and (16) simultaneously, which requires some cumbersome calculations, we divide the job into several steps. First, th 1 is solved (in terms of other unknowns) from equation (16). Next, the hitting point w* of w on :E* is obtained in terms of the unknowns.
First, let P denote a fixed point on :E with the coordinates described in equations (72)-(74). Since w 1 == 1, the coordinate transformation in equation (22) becomes Suppose that a point Po on the periodic orbit passing through P has coordinates in addition to x 1 and y 1 given in equation (75). Here r is not arbitrary. Depending on 8 and e, it must be chosen so that Po lies on the periodic orbit. Hence we take The hitting time from Pe to P is assumed to be of the form th() = B + t1b + t2b 2 + t3b 3 + 0(4).
Substitution of equations (75)-(79) into equation (29) yields a set of 12 algebraic equations from which one can determine the 12 unknown coefficients in equations (76)-(79). The solution procedure is quite similar to that presented in Section 6.1 and is not repeated here. The resulting solution is For the region r cos B > 1, we set as the coordinates of a point P<P on the periodic trajectory through P (see Figure 5). Again, f is not arbitrary. Let f and the backward hitting time th<P of P<P toP be written as ( 1 _ c') (r cos()-1) , Recall that the parameter K, which is the stiffness of the free spring, accounts for the only nonlinearity in our system. In fact, it can be verified that all of the functions fij, 9ij, ffj, and gij are identically zero when r;, = 0. In this case, all of the modal subspaces are simply two-dimensional planes.  The corresponding relationship for the second mode is for At < 1 for At > 1 Again, observe that as "" = 0, the frequencies will be constant, corresponding to the limiting linear system.
It is also interesting to note that the modal motion for this system is synchronous in the sense that both masses simultaneously reach their peaks. However, the displacements do not necessarily vanish simultaneously. Such a phenomenon has been verified in [ 41] through simulations for a special case of PWL systems. To see this, observe that !323 ( 1r) and (3~3 ( 1r) are zeros, which yield 921 (r, e) ~ 0 and g~l (r, e) ~ 0 when e = 7L Moreover, we already have g 2 z(r, ()) ~ 0 and g~2 (r, e) ~ 0. These imply that both YI and Y2 are necessarily zero when sine = 0 by equation (24), the equation for the normal mode. Hence, both modes are synchronous as described above. Thus, we can compute the ratio of peak amplitudes A2/ At for each mode, which are: Although conventional perturbation techniques fail to apply to PWL systems as mentioned previously, some energy-based methods, like the Ritz method, are applicable here. As a comparison, the Ritz method will be employed to the example PWL system ( 48) and ( 49). Since the Ritz method can be found in many books (for example, [33]), only final results are presented here. The detailed calculations are shown in Appendix B.
In what follows, we shall compare the present method with the Ritz method and the exact solutions obtained through numerical simulations. To guarantee that they belong to the same      In simulations it is observed that for small8, both the present method and the Ritz method have good estimates to the modal solution. However, as 8 increases, the Ritz method gives a better approximation than the present method for the periodic solutions. This is because the present method is based on asymptotic expansion which holds only for small 8, i.e., in the neighborhood of the invariant disk.
The frequency and peak amplitudes relationships are given in Figures 9 and 10 for different combinations of r;, and modes. It is easy to see that good approximations are achieved for both methods for small 8 ( 8 < 0.4 ), while for large 8, the Ritz method is more accurate than the present method. Also, the effect of r;, can be clearly seen.

Conclusions
Based on the construction of a convenient Poincare map and making use of invariant manifold theory, a general procedure for constructing the normal modes of a class of piecewise linear (PWL) vibratory systems has been developed. The class of systems considered contains two linear subsystems separated by a switching hyperplane in state space. It is shown that such a PWL system possesses a canonical form upon which the proposed method applies.
1 .1 ·..-.:-: -R1tz· method: · · · · · · · >: · · · · · · · ·  There are several steps in this constructive procedure. First, using the switching hyperplane as a Poincare section ~. a one parameter family of fixed points on :E is obtained for each mode. Then, using a form of polar coordinates, the state space is decomposed into three regions.
The first one (r < d) is the region of small amplitude vibrations in which the usual linear modes hold. While both the second (r > d and r cos B < d) and third (r cos B > d) regions correspond to large amplitude motions, they are governed by two different linear subsystems. In the first region, the normal modes are invariant elliptic disks and are easily obtained. The normal modes in the other two regions are constructed by matching the fixed point family on the Poincare section. Each mode in the second region is expanded in a series form around a neighborhood of the invariant disk in the first region. In the last region, the normal mode is expanded in a series form near the point where the invariant disk in the first region and the switching hyperplane are tangent.
From the expressions for the normal modes, the modal dynamics and corresponding frequency-amplitude dependence relationship on each mode are obtained. As expected, the normal modes, the modal dynamics, and the frequency-amplitude relationships are all of piecewise type (although not piecewise linear).
An undamped, two degrees of freedom system with a clearance is demonstrated as an example to illustrate the procedure. The numerical results demonstrate the effectiveness of the procedure. It must be admitted that the Ritz method is far easier to apply and is more accurate in terms of estimating frequencies. The advantage of the present method is that it allows one to construct the differential equations that give the modal dynamics. Appendix A. A Necessary Condition for the Existence of Normal Modes for PWL Systems The fixed point equations (15) and (16)  where a = dwiJvu, ON stands for theN-dimensional zero vector, and Im is the identity matrix of order m. If J is nonsingular, by the implicit function theorem, there exists only one fixed point for each "' near 0. Then it is impossible to have a one-parameter family of fixed passing through ( ~0 , 0). Thus, J must be singular in order for the normal mode to exist. It will be shown below that J is indeed not of full rank and exactly has rank equal to 4N-1 if the system has no internal resonance. Therefore, without internal resonance, we will have w· 1 -cos _ t 2n =I 0 and . Wi which has 4N -2 columns. It is obvious that j~ is independent of the columns of J". Also, it is not difficult to check thatj~N+l is dependent on j~ and the columns of J". Thus, we deduce that the rank of J" is 4N-3.
Finally, one can prove that adding j 1 and hN+l back to J" will increase the matrix rank by 2, which results in our claim that J has a rank equal to 4N -1.
If wdwt is an integer for some i =f 1, then hN -T(2n jwr) will have a null space of dimension more than 2, which implies the rank of J will be less than 4N-1.

Appendix B. Detailed Calculations for the Ritz Method
The Ritz method assumes a modal solution of the form: where w is the modal frequency and b2 is the ratio of peak amplitudes to be determined. The system equation is rewritten as: Then by minimizing the solution error in the sense 271" I Ci(B) dB= 0, where B = wt, w and b2 can be obtained in terms of At . For the example system given in equations (48) and (49) where + corresponds to the first mode and -to the second mode.