An invariant manifold approach to nonlinear normal modes of oscillation

Summary. A method for determining the amplitude-dependent mode shapes and the corresponding modal dynamics of weakly nonlinear vibratory systems is de· scribed. The method is a combination of a Galerkin projection and invariant manifold techniques and is applied to a class of distributed parameter vibratory systems. In this paper the general theory for a class of conservative systems is outlined and applied to determine the nonlinear mode shapes and modal dynamics of a linear Euler-Bernoulli team attached to a nonlinear elastic foundation.


Introduction
Traditional methods for determining the effects of nonlinearities on the mode shapes of vibrating systems typically involve an assumed form of harmonic time behavior and series expansions for the frequency of oscillation and the mode shapes in terms of the amplitude of vibration (see, for example, Bennouna and White [1], Benamar et al. [2], or Szemplinska [3]).A more formal approach is found in the classical work of Rosenberg [4], which was followed by many subsequent efforts, for example, those of Rand [5][6][7] and Vakakis [8][9][10] and their co-workers, all of which used definitions of normal modes in terms of periodic motions for finite-dimensional, conservative, nonlinear systems, typically with two degrees-of-freedom.Recent work by Jezequel and Lamarque [11] uses normal form theory to generate nonlinear normal modes, again for systems with two degrees of freedom.Their approach does allow for the inclusion of special types of damping.Vakakis has used nonlinear normal mode methods to analyze nonlinear effects on mode localization in multi degree of freedom, cyclic systems [12].
The method described herein takes a point of view, which is fundamentally different from that found in previous engineering literature on the subject.Rather than view nonlinear normal modes as synchronous motions of conservative systems, they are viewed as motions on invariant manifolds which are tangent to, and of the same dimension as, the linear eigenspaces in the system phase space.These manifolds are the standard foliation of the center-stable manifold of the equilibrium point of interest, and their existence depends crucially on certain nonresonance conditions being satisfied.
These manifolds are familiar to someone trained in dynamical systems theory, but they have only recently been exploited for the construction of nonlinear normal mode shapes and modal dynamics.Shaw and Pierre developed methods of constructing nonlinear normal modes for both lumped-parameter [13] and distributed parameter [14] vibratory systems in this manner.Such an approach allows for the definition and construction of nonlinear normal modes even in the presence of gyroscopic and/ or dissipative terms.These efforts were followed by the work of Nayfeh and co-workers [15,16], who use both perturbation and invariant manifold techniques to generate nonlinear normal modes, including the investigation of a case with internal resonance [17].King and Vakakis [18] have used a combination of Rosenberg's ideas with the method presented in [14] in order to generate normal modes for continuous systems.This work also includes an investigation of the stability properties of nonlinear modes.
Along similar lines, much work has been done on the existence and smoothness of these manifolds for finite-dimensional, conservative, Hamiltonian systems.The main result in the nonresonant case is the Lyapunov center theorem which states that if the Hamiltonian is C 1 and the linearized system has purely imaginary eigenvalues +fij (j = l, ... ,m) such that for j =I= k, nk;nj =I= n (n = 1,2,3, ... ), then there exist m two-dimensional, local C 1 invariant manifolds.(This result is due to Lyapunov, and a simple proof is given by Kelly [19).)These are the nonlinear normal mode manifolds of interest here, and each manifold contains a oneparameter family of periodic solutions.Other work in this area includes that of Weinstein and Moser on bounds for the number of periodic orbits which can exist near an equilibrium point [20,21], and extensive work on the resonant cases; see [22][23][24][25], for example.It is important to point out that in the strong resonance cases (n = 1, 2, 3), the number of families of periodic solutions is not generally the same as the number of pairs of complex conjugate eigenvalues.The nonlinear normal modes defined in the present work do not generally exist in such cases.
The method developed by Shaw and Pierre [14] is a systematic means for generating nonlinear normal mode manifolds for nonresonant partial differential equations which arise in structural vibrations.It involves some rather cumbersome calculation steps and the results depend on the selection of a suitable base point.The method described herein is aimed at providing an alternative approach which avoids some of these difficulties.The first step is to apply a usual Galerkin projection of the partial differential equations of motion onto a set of basis functions which consist of the mode shapes of the linearized system.This leads to a set of ordinary differential equations which describe the dynamics of the nonlinear system in terms of the linear modal amplitudes.This standard first step leads to equations of motion which are uncoupled at linear order.The invariant manifold approach developed in [13] is then applied in order to determine how the nonlinear coupling distorts the mode shapes and dynamics in the weakly nonlinear system.At this step it is important to treat the linear modal amplitudes and velocities as independent variables, as they provide the coordinates which parameterize the manifolds.Overall, this procedure provides a straightforward means for determining the amplitudedependent shapes of natural vibration for weakly nonlinear, distributed parameter systems.The resonance cases are flagged directly by singularities in the equations which are to be solved for the coefficients describing the manifolds.
One important conclusion of the results obtained is that one must be careful when generating reduced order models for the vibration of continua using modal projections.The typical step of projecting the partial differential equation onto a subset of the linear modes must be done with care, especially if one requires more than first-order nonlinear information about the system at hand.This shortcoming was observed by Nayfeh et al. [26] in the analysis of the dynamic response of a valve system, and is systematically explored here.
The paper is arranged as follows.In Section 2 the general method is described for a class of weakly nonlinear partial differential equations which arise in mechanical vibrations of distributed parameter systems.The method is applied to an example problem in Section 3 and the paper closes with some remarks in Section 4. It should be noted that the methods presented herein were motivated by the analysis of nonlinear mechanical structures, but they may have• a wider range of applicability.

Equations of Motion
Consider one-dimensional vibrations of a one-dimensional continuum which occu~ pies the region described by the spatial variable s E (0, 1).The equations of motion are assumed to be of the form a 2 w(s, t) with linear boundary conditions at s = 0, 1 of The equilibrium configuration about which the system vibrates is taken to be w(s, t) = 0, that is, L(O) = N(O) = 0.The operator L is taken to be linear and self-adjoint, while the nonlinear operator N is assumed to possess smoothness properties that allow it to be expanded to any order desired.The linear system, N = 0, is assumed to have nonzero natural frequencies flj with nj;nk::;:: 1, 2, 3, ... and associated normal mode shapes cfJ/s) for j = 1, 2, 3, ... , which are normalized according to the definition given in equation ( 4) below.The system eigenvalues are thus ±i!lr The inner product between real-valued functions on the unit interval is defined in the usual manner as The following general properties then hold for the linear system: where oij is the Kronecker delta.
The solution of the equation of motion is now assumed to be of the form of an infinite series involving the linear normal modes and their amplitudes as functions of time as w(s,t) = L 4>/s)q/t).j ... l (5) Here qi(t) represents the contribution of the jth linear mode to the response.This solution form is substituted into the equation of motion, which is then projected onto the linear modes by using the inner product.Using the properties of the linear system, this results in the following set of differential equations of motion in terms of the linear modal amplitudes: i = 1' 2, 3, ... ' (6) where an overdot denotes a time derivative, q is the infinite-dimensional vector comprised of the q;'s, and the nonlinear coupling terms in Gi are determined from i = 1,2,3, ... .(7) These can, to cubic order in the q/s, be expressed as 00 00 00 00 co where the vilm and J.Lumn coefficients depend on the nonlinearity and the mode number i.
In order to determine the nonlinear normal modes in terms of invariant manifolds by the method developed in [13], these equations of motion are written in first-Jrder form as qi(t) = Pi(t), (9) pi(t) = -OfqJt) -Gt(q(t)), i = 1, 2, 3, . ... We restrict our attention to systems for which OJOi is not equal to an integer for i :f:.j.
In the following, the present approach will be compared with the more commonly used method of "direct projection" onto a single mode.That approach employs a single mode, say c/Jj, in the expansion (5) and projects the equations of motion onto that mode.The i = j equation of motion is the only one of interest in that case, and it has the same form as equation ( 9) with i = j, except that the multiple sums in Gj are reduced to the terms vjjjq](t) + J.Ljjjq}(t) + •• • .The effects of the missing terms are considered in detail below.

The Normal Mode Invariant Manifold Equations
The key observation made in [13} is that for either a linear or a nonlinear system a normal mode motion is one of a family of motions for which the system behaves like a single degree of freedom oscillator; that is, like a second-order nonlinear oscillator.Such a motion will take place on a two-dimensional invariant manifold.For linear oscillatory systems these manifolds are the planar eigenspaces associated with the complex conjugate pairs of eigenvalues.The existence of these eigenplanes for the linear system is guaranteed for this class of systems with the given assumptions on the frequencies.For weakly nonlinear systems of the type being considered here, the curved manifolds are necessarily tangent to the linear eigenspaces at the equilibrium (w(s, t) = 0).Once these manifolds have been determined, the manner in which the linear modes combine to form the nonlinear normal modes will be evident.In addition, the dynamics on these manifolds, which represent the associated dynamics of the nonlinear normal mode, are easily recovered.
The construction of the modal manifolds is carried out in the usual manner.First, one must choose coordinates which describe the manifold; in this case (qk, pk) are a natural choice for the k th normal mode.Next, it is assumed that if the motion of (qk, pk) is known, then the motion of all other (qi, p) pairs can be described in terms of (qk, pk); this reduces the dynamics of the entire system to that of a single degree of freedom system or, equivalently, the dynamics are restricted to a twodimensional manifold.Also, the equations of motion must be satisfied for such a specific type of motion.These ingredients lead to a set of equations for the normal mode manifolds which, although not solvable in general, easily admit solutions in series form in the variables (qk, pk), just as is done in the construction of other types of invariant manifolds (see, for example, [27] or [28]).The solution of the manifolds contains the desired information regarding the composition of the nonlinear normal modes in terms of the linear normal modes.
By defining (qk, pk) = (uk, vk), a normal mode motion is assumed to exist and be (at least locally) expressible in the form where, by definition, (10) (11) In this way, the dynamics of the entire structure is specified by the dynamics of the kth linear modal amplitude and velocity.
The equations for the Qik's and Pik's are obtained by first taking a time derivative of each item in equation (10) to obtain .aQik .aQ;k .
The following substitutions from the equations of motion and equation (10) are then made: where Qk is the vector of Qik's with i = 1, 2, ... .This process eliminates all time derivatives from equation (12), resulting in the following equations for the Qik's and P;k's in terms of uk and vk: i = 1, 2, 3, .. . .(14) Note that the i = k equations are trivially satisfied.
A closed form solution of these equations is not generally attainable (except in the case of so-called similar normal modes [8]), but a local solution near the origin can be obtained in terms of a series expansion.Matters are simplified by noting that for conservative, nongyroscopic systems, only those terms which are consistent with constant amplitude, standing wave normal modes are required.Thus, to cubic nonlinear order, the series is given by (15) i=1,2, 3, ... , where the following identities for i = k are obvious from equation (11): (The missmg afik and b 1 ;k coefficients are identically zero in the conservative, nongyroscopic case; see [13] for examples where they are nonzero.This is easily confirmed by including them in the above and subsequent steps.) Substituting the expansions (15) into equations ( 13) and ( 14), expanding in terms of uk and vk, and gathering the coefficients of the u'f: 1 uf:' 1 terms, where m 1 , m 2 = 0, 1, 2, 3 and m 1 + m 2 < 4 (we are working only to cubic order here), yields a set of equations in terms of the unknowns afik and buk• These equations are as follows.
At linear order: from equation ( 13), the uk coefficient (18) from equation ( 14), the uk coefficient (19) At quadratic order: from equation ( 13), the ukuk coefficient (20) from equation ( 14), the u~ coefficient (21) from equation ( 14), the v'f coefficient (22) At cubic order: from equation ( 13), the u~uk coefficient (23) from equation ( 13), the vt coefficient from equation ( 14), the u~ coefficient = -b7ikO%-bz;k L L Vktm(a3lm + a3ml)-1Lkkkkb2ik-b4ikvkkk; (25) l=lm=l from equation ( 14), the ukv~ coefficient Several features of these equations are worth noting.First, they are linear in the unknowns, and nonhomogeneous terms arise from the nonlinearity.Also, the known solution for the i = k case [see equation (11)] is consistent with these equations.The solution for the coefficients of the linear, quadratic, and cubic terms are now considered in turn.
Regarding the linear terms, it is immediately clear from equations ( 18) and ( 21) that (27) since nk :fo ni and using identity (16) above.This is not unexpected, since it simply says that the only contribution at the linear order to the kth nonlinear normal mode comes only from the kth linear mode.As described in detail below, this implies that the linear modes are sufficient for correctly capturing the leading order nonlinear dynamics, as nonlinear corrections to the mode shapes will riot affect the leading order nonlinear terms in the equations of motion.However, the modal distortions must be accounted for if one considers higher order dynamics.
Regarding the quadratic terms, these represent three linear equations in three unknowns, and are easily solved when the equations are nonsingular.The solution is given by vikkca The coefficients for i = k are given in equation ( 16) above.Note that these equations are singular in two cases: (1) If any ni = 2flk, an internal resonance condition arises in which the i and k modes cannot be uncoupled. ( 2 If any ni = 0, a condition not allowed by the original restrictions on the natural frequencies.
The solution for the cubic terms involves four linear equations in four unknowns, with nonhomogeneous terms which arise from the cubic order nonlinearity in the equation of motion and from the solution at quadratic order.The appearance of terms from lower order at the current order is typical of such procedures and is well known in, for example, perturbation and normal form calculations.The solution for the cubic coefficients, for i -:f:.k, are given by The coefficients for i = k are given in equation ( 17) above.Note that these equations are singular whenever fl; = 3flk or nj = nk, resonance conditions not allowed by the original restrictions on the natural frequencies.
The uij terms can be simplified by noting that, although the a 11 m terms are not symmetric in the l and m subscripts, they are summed over both l and m, so that This is useful for calculations in specific examples . .

Nonlinear Normal Mode Shapes and Dynamics
The nonlinear mode shapes are now reconstituted by using equation ( 5) with the q/s replaced by their modal versions, the Q/s.
where the reader is reminded that u k and u k represent the displacement and velocity of the kth linear modal amplitude during the kth nonlinear normal mode motion.These normal modes represent standing wave motions of the system as can be seen by the following consideration.At any instant, say t = t *, at which the linear modal velocity is zero, that is, vk(t*) = 0, the entire velocity field is zero: (awk(s, t*))/ at = 0 'Vs E [0, 1].Thus every point of the system reaches its maximum displacement at the same instant.This fact allows mode shapes, denoted here by Wk(s ), to be plotted in a manner just as is done for linear systems.They are simply the zero-velocity configurations of the system, which are given by 00 where u\ = uk(t*) is the maximum value that uk(t) achieves during the motion.
Note that at any instant, say t = t 0 , at which the linear modal amplitude is zero, uk(t 0 ) = 0, the entire displacement field is generally not zero, wk(s, t 0 ) =I= 0. This is due to the presence of the nonsymmetric, quadratic nonlinearities, and in particular the term a 5 ;kvf(t 0 ).For systems which do not possess quadratic nonlinearities, wk(s, t 0 ) = 0 to leading order, since all a 5 ik terms are zero.In that case, every point of the system passes through zero at the same instant.
Note that, in contrast to the linear mode shapes, the nonlinear mode shapes depend on the peak amplitude of motion.The distortion from nonlinearities is contained in the quadratic and cubic terms, and in particular is expressed here in terms of the linear modes shapes with the contribution of the ith linear mode shape on the kth nonlinear mode shape being captured in the coefficients a 3 ;k and a 6 ik• In addition, during a given modal motion, the shape of the structure varies in time, as can be seen by considering equations ( 30) and ( 31) and assuming that uk(t) and vk (t) are periodic in time.This aspect of the modal response is not allowed by the methods used in [1][2][3], which essentially assume a fixed shape whose amplitude varies according to a second-order nonlinear equation of motion.
The dynamics of a nonlinear normal mode are obtained by simply restricting the equation of motion to the modal manifold.This is accomplished by substituting the modal expressions for q and p given in equation ( 15) into the equations of motion (6).The result is which, due to the aki terms, differs from what one would obtain by simply projecting the equation of motion onto the kth linear mode, <f>k(s).However, a projection onto c/>k(s) gives exactly the same dynamics as this result if one is •working only to first nonlinear order; this is discussed in more detail below.The amplitude-dependent frequency of nonlinear oscillation can be determined form equation (34) by using a perturbation method (here Lindstedt's method is used).It is found to be where uZ is the peak amplitude of uk during the motion.

Remarks
At this point it is interesting to examine some special cases.First, consider the general system with quadratic and cubic nonlinearities, but for which only the first-order nonlinear modal dynamics are of interest.These are trivially obtained by simply using a projection onto the single linear mode c/>k(s), which directly yields the quadratic coefficient ak = vkkk in equation (34).In this case there is no need to compute the nonlinear mode shapes, even to first nonlinear order.However, note that the frequency information will not be correct to leading nonlinear order; see equation ( 35).If one is interested in second-order nonlinear modal dynamics, which are required to obtain the correct frequency estimate in this case, it is clear from equation (34) that the first-order nonlinear corrections for the mode shapes must be known.These are expressed in terms of a 3 ;k and a 5 ;k , which combine with the vilm to give ak 3 and aks• For a system which possesses symmetric nonlinearities, in which case Pum = a 3 ;k = a 5 ik = b 4 ik = ak 3 = aks = 0, a projection onto the single mode </>k(s) will give the correct nonlinear modal dynamics to cubic order, via the required coefficient {3k = P..kkkk• Again, the correct higher (e.g., fifth-) order dynamics will require the appropriate nonlinear mode shape terms.Continuing this line of reasoning, a generalization of the development above provides the cubic order nonlinear terms for the mode shapes, and this information can be used to generate the correct nonlinear modal dynamics which includes terms up through and including quintic order in uk and ilk.Boivin et al. have carried out such calculations and have demonstrated the benefits of including such higher order terms in simulation studies [29].The fact that shape corrections are required to generate correct higher order dynamics has been used in estimating wave speed-wave amplitude relationships for water waves (see, e.g., Lamb [30], article 250), but has been overlooked in the structural dynamics community, where reduced order models for systems with quadratic and cubic nonlinearities are typically obtained by direct projections onto the linear modes of interest.

Example
The example considered here is a simply supported, linear, Euler-Bernoulli beam attached to a nonlinear elastic foundation.Figure 1 depicts the physical system.The linear problem is well known and can be found in textbooks (e.g., [31]).This problem was selected since the linear mode shapes are simple functions, making closed form solutions possible.The general approach will be made clear from this relatively simple example, and it will be evident where one will be required to carry out numerical computations in less straightforward cases.
The nondimensional equation of motion, expressed out to cubic order, is a 2 w(s,t) a 4 w(s , t) with boundary conditions as as Here K, p, and TJ represent the linear, quadratic, and cubic coefficients of the foundation, respectively.The linear problem ( p = TJ = 0) has natural frequencies and normal mode shapes, respectively, of w(s,t) ( Two cases will be considered.In Case I the system with both quadratic and cubic nonlinearities will be studied.First-order nonlinear terms for the mode shapes will be generated, and from these the second-order dynamics will be obtained.For Case II, the system with symmetric nonlinearities, p = 0, will be studied, and for this system the first -order nonlinear terms will be computed for both the mode shapes and the modal dynamics.

Case I
The coefficients vilm and ILilmn need to be generated; this is carried out as follows for this system.First, the nonlinear operator acting on the assumed form of the solution, equation ( 5), is written as Next, this is projected onto the ith linear mode according to equation (7) in order to obtain Gi: 00 00 G;(q(t)) = L 1 lf>Jt)p L L <J>,(t)lf>m(t)q,(t)qm(t)dt+ L 1 lf>;(~)1J 1=1 m=l n=l The order of summations and integration is switched in order to obtain the coefficients vilm and J.Lamn from equation (8).The result is (41) given in equation ( 38).The nontrivial a 3 ik and a 5 ik simplify to Here the internal resonance condition is obvious: K = l1r 4 (i 4 -4k 4 )/31 is exactly the condition for which ai = 2Gk.At this point it is also seen that the contributions of the ith linear mode shape to the kth nonlinear normal mode decreases at the rate of i-7 for large i.It is also worthwhile to note that a 3 ik and a 5 ik are zero for all even values of i, including i = 2k, and for i = k (even or odd).Thus, only the odd, that is, the symmetric, linear modes contribute to the nonlinear modes at this order.
The summation terms uk 3 and uk 5 are obtained by direct substitution of the above expressions for vilm' a 3 ik, and a 5 ik into the expression for uij given in equation ( 28), and making the appropriate subscript assignments.These equations are rather unenlightening and are not presented explicitly here.It is sufficient to note that they converge since the coefficients of uk 3 decay in l and m like 1-3 m -a for large l, m while for uk 5 the decay is at the very rapid rate of z-12 m -7 .
Convergence is also aided by the fact that many of the terms in the sums are zero.For plotting purposes, it is convenient to define the following function which, for each mode number k, are functions of only the linear foundation coefficient K. Figure 2 depicts these functions for the first three modes.In each case the resulting plots indicate a close approximation to the exact value of the doubly infinite sum.The number of terms was varied and convergence was inspected visually; the typical required number of nonzero terms was in the range 3-6.Note that the sums are less significant for beams with a large linear stiffness, and it is easily shown that for large K, Sk 3 -K-1 and Sks -K-2 • These decay rates are interrupted by values of K at which internal resonances arise, at which points the S kJ are singular.In the K range shown in Figure 2, only the 2:1 resonance between the first and the third linear modes appears [at K = 7T 4 (4 X 3 4 -1) /3 :::: 10, 488].
All terms required for the mode shapes and modal dynamics are now in hand.
In this case the peak-amplitude nonlinear normal mode shapes are easily expressed to first nonlinear order by using equation ( 33), where the a 3 ik are given in equation ( 43) above, the cubic order terms are ignored, and six nonzero terms were used in each summation.Plots of the peak-amplitude mode shapes for the cases k = 1, 2 are depicted in Figure 3 for various amp1itudes, specified here by u't.In each case the zero-velocity beam configurations for both the "up" and the "down" positions are shown, indicating the lack of symmetry which arises since the foundation is softer (relative to the linear foundation) in the "down" direction and stiffer in the "up" direction.(See the Appendix for the amplitudes used for the "down" configuration.)It is worth noting that the effect of the elastic foundation does not have a decreasing effect on higher modes.This is in contrast with the linear case [cf.In order to determine the modal dynamics, the coefficients vkkk' J.Lkkkk' uk 3 , and ak 5 are needed.They are given above for the general case and are used to produce the oscillator (46) which governs the dynamics of the kth nonlinear normal mode.The amplitudedependent frequency of nonlinear oscillation can be written as (47) where ut is the peak amplitude of uk during the motion and the coefficient Ak is obtained directly from equation ( 35) and is given by A.CK, p, 1Jl = (9(K + e"•l( ~ 1J + s.,CK)p 2 ) where the terms S k/ K) are given in Figure 2. (48) For a fixed amplitude of motion, the frequency of oscillation is seen to be linear in the cubic foundation coefficient TJ, quadratic in the quadratic foundation coefficient p, and has a complicated dependence on the linear foundation coefficient K.In order to plot Ak and examine the difference between the usual direct projection result and the present result, it is convenient to express A k as where /{ and ft are the terms one would obtain by a direct projection, and fff 2 involves the S kj terms and arises from the nonlinear distortions to the mode shapes.The net effect of the nonlinear mode shape depends on the relative sizes of iJ and p and the size of fff 2 relative to ff and ft.Plots of these three functions for the first two modes are given in Figure 5.
The importance of including the modal distortions in the modal dynamics can be demonstrated by a single case: that for the second mode, k = 2, with 17 = 0. Since f!P = 0, the only nonzero term in A 2 is proportional to f!J 2 , which is precisely the term that is neglected if one ignores the coupling between the linear modes.The magnitude of A 2 in this case is given by p 2 fi_ 2 , where f!J 2 is shown in Figure 5.It is seen that, while the correct result of zero frequency shift is approached as K ~ co, in the lower range of K the direct projection method can lead to significant discrepancies.
It is thus confirmed that, although the calculations are nontrivial, the first-order nonlinear corrections to the mode shapes provide the correct modal dynamics to second order, and that these are nontrivial in some cases.We now turn our attention to a simpler case, in which the foundation contains only cubic nonlinearities.Note that fi. 1 = 0.

Case II
The calculations here follow the above with p set equal to zero.Since only first-order information regarding mode shapes and dynamics is to be obtained, the calculations are much less involved.The modal dynamics are immediately obtained from the results of Case I with p set equal to zero, and are precisely what one would obtain from a simple projection of the equations of motion onto <f>is ).
For the nonlinear mode shapes, the following quantities are required and are easily obtained using equation ( 41):  It is clear that the above procedure can be generalized in a straightforward manner to other situations.In fact, in [13] the complete theory has been completely worked out for nonresonant systems with a finite number of degrees of freedom possessing all possible quadratic and cubic nonlinearities.This includes cases which are nonconservative and/or gyroscopic in nature.
In [14] an alternative formulation is used for defining nonlinear normal modes of continuous systems in which the dynamics of the entire structure is specified by the dynamics of a single point on the structure, s = s 0 • While the method in [14] is more naturally suited to handle some problems, such as those with nonlinearities in boundary conditions, the solution procedure is rather cumbersome and the results are sensitive to the choice of the base point s 0 • The present formulation is equivalent, yet significantly simpler to implement.The present approach is capable of handling discrete nonlinear elements.For example, a nonlinear spring attached at some point along a beam can be handled by simply using delta functions in the spatial variable (see [29]).If nonlinear effects exist at the boundaries, the method described in [14] may be easier to apply, although a suitable change of variables may render the equation of motion in the form of equation (1).
A possible extension of these ideas is in the area of nonlinear waves.Here, one would seek special solutions of the governing partial differential equation which represent the dynamics of some low dimensional dynamical system.The method described in [14] is well suited for this.For example, for determining the "single mode" dynamics of a nonlinear wave equation, one simply says that if the displacement and velocity of a single point are known, then the entire displacement and velocity field are determined by these.
One issue which has remained unexplored in the current work is a consideration of the stability of nonlinear normal modes.While much work has been done in this area using traditional methods for two degree-of-freedom systems, it has not been considered in the context of the present method.This has important implications with regard to whether or not these nonlinear normal modes will be observable.The usual parametric instabilities between coupled modes are flagged by the present method in the form of singularities.However, other instabilities may occur, especially at large amplitudes (see [9], for example).
The method described herein can be thought of as a means of generating the "best" single-mode model of the given system.A generalization of the procedure has been developed which can be used to generate invariant, reduced-order models of arbitrary order.In order to achieve this, one expresses the response in a form in which the unmodeled modes are expressed as functions of the modeled modes.This idea is of potential use in developing low-order models which may be helpful in computational algorithms or in the development of control schemes.These ideas have some relation to those of the nonlinear Galerkin methods which are used to approximate inertial manifolds [32], although there is an important difference since that method is applicable to systems with dissipation whereas the systems of interest here are conservative.
A final note: While the above method is in some respects similar to standard perturbation techniques, a comparison between the results obtained by this method and those obtained using a combination of harmonic balance and eigenfunction expansions (see, e.g., [1][2][3]) shows that different results are obtained for the nonlinear mode shapes, even for the simple example presented here.The approxi-mations for the frequency of nonlinear modal oscillations agree to first order, but the discrepancies in the mode shape approximations at first order will lead to differences in higher order frequency estimates.Since the present method is based on the fundamental principle of dynamic invariance, the results obtained here are the correct ones, and naive perturbation strategies should be applied with caution.

Case[
The first-order nonlinear solution for the modal oscillator (34), subject to initial conditions (uk(O), vk(O)) = (X, 0), is obtained using Lindstedt's method and is given as follows: Note that while this starts at a displacement of X with zero velocity, one-half cycle later the velocity is again zero but the displacement is not equal to -X, but is given by The first-order nonlinear solution for the modal oscillator (34) with ak = 'Yk = 0, subject to initial conditions (uk(O), vk(O)) = (X, 0), is obtained using Lindstedt's method and is given as follows: where wk = nk + X 2 (3,Bkj80V + ... .This starts at a displacement of X with zero velocity and one-half cycle later the velocity is again zero and the displacement is -X.

1 Fig. 7 .
Fig.7.Nonlinear normal mode shape at peak deflection for k = 2, K = 13,800, p = 0, and "' = 10 6 , which is close to an internal resonance.Here u~ = 0.04 and the linear mode shape for the same ui value is shown as a dashed curve.