S W Sha

Pierre

NORMAL MODES FOR NO:N-LINEAR VIBRATORY SYSTE1\1S

A methodology is presented which extends to non-linear systems the cDncept of normal modes of motion which is well developed for linear systems. The method is constructive for weakly non-linear systems and provides the physical natUlre of the normal modes along with the non-linear differential equations which govern their dynamics. It also provides the non-linear co-ordinate transformation which relates the original system co-ordinates to the modal co-ordinates. Using this transformation, we demonstrate how an approximate nunlinear version of superposition can be employed to reconstruct the overall motion from the individual non-linear modal dynamics. The results presented herein for non-linear systems reduce to modal analysis for the linearized system when non-linearities are neglected, even though the approach is entirely different from the traditional one. The tools employed are from the theory of invariant manifolds for dynamical systc~ms and were inspired by the center manifold reduction technique. In this paper the basic ideas are outlined, a few examples are presented and some natural extensions and applications of the method are briefly described in the conclusions.

INTRODUCTION

The idea of developing some type of normal mode represc:ntation for non-linear systems is not new. The work of Rosenberg, the main results of which are gathered together in a review article (1], was ground-breaking and laid the foundations for most of what followed; for example, the works of Atkinson and Taskett [2], Szemplinska [START_REF] Szemplinska | The Behaviour of Non/inear Vibrating Systems[END_REF], Rand [START_REF] Rand | [END_REF]5], Greenberg and Yang [6], van der Varst [7], Yen [8] and Caughey and Vakakis [9]. The thesis of Vakakis [START_REF] Vakakis | Analysis and identification of linear and nonlinear normal modes in vibrating systems[END_REF] contains a nice account of this history and provides the most recent contribution.

Nearly all previous work on non-linear normal modes deals exclusively with conservative systems, and much of it with systems having two degree!; of freedom. In this paper we present a formulation which extends the concept of nonlin.ear normal modes to a class of general systems with N (finite) degrees of freedom which caJO include damping. The method is constructive in a local sense, i.e., near an equilibrium point of the system, so that one can actually obtain the physical nature of the modes, generate the attendant dynamics, and reconstruct a general motion from the modal dynamics. This is demonstrated explicitly in several examples.

A key feature of the present work is that we formulat1e the problem in terms of first order differential equations and include velocities as well as displacements as dependent variables. Not only does this methodology allow for the inclusion of damping in a systematic manner, but we show that our formulation also results in the appearance of important velocity dependent terms for the normal modes of non-linear conservative systems. In the traditional approach for non-linear normal modes of conservative systems, these terms are expressed as functions of the system energy, thereby preventing a straight~ forward extension to non-conservative systems in that approach.

The general approach is inspired by the theory of invariant manifolds for dynamical systems. Throughout the paper we motivate the results as generalizations of the linear theory and relate the two whenever possible. Also, since our primary interest lies in discrete mechanical systems with osciJlatory behavior, we present results for sets of coupled, nonlinear, second order, ordinary differential equations (ODE's), the linearization of which has complex eigenvalues. As described in the closing remarks, the ideas presented in the paper can be easily extended to more general cases.

The paper is arranged as follows. In section 2, we present a review of normal modes for linear systems in order to set the stage for the extension to the non-linear case. This formulation is novel in some respects and provides an alternative approach to linear system analysis. In section 3, the non-linear normal modes are formulated in terms of invariant manifold theory. The method of approximate solution for the normal modes is presented in section 4. Section 5 describes the non-linear transformation from physical to modal coordinates and shows how it can be used to represent a general motion of the full system in terms of the non-linear modal co-ordinates. Three example problems, including extensive simulation results, are given in section 6. The paper closes with a discussion, along with some conclusions and directions for future work, in section 7.

NORMAL MODES FOR LINEAR SYSTEMS

We begin by presenting the definition of normal modes for linear oscillatory systems in a manner best suited to allow for extensions which include non-linear effects. We introduce an approach which is based upon the property of invariance of the modal subspaces and which allows us to represent the system's modal dynamics by a set of uncoupled linear oscillator equations. We show in Appendix A that this representation is strictly equivalent to that obtained from the complex eigensolution of the standard first order eigenvalue problem. An interesting feature of this approach is that the eigenvectors are obtained prior to the evaluation of the eigenvalues.

Consider the free motion of a system with N degrees of freedom which has been linearized about an equilibrium at the origin: Mi+Cx+Kx=O, (I) where x = [x 1 , x 2 , •.• , xN ]T is the displacement vector, a dot denotes a time derivative, and M, C and K are the inertia, damping (or gyroscopic) and stiffness matrices, respectively, obtained by linearizing the dynamics of a mechanical system about an equilibrium point. We require that M be non-singular and allow C and K to be arbitrary, which allows for the inclusion of dissipative, gyroscopic and non-conservative terms. We also assume that the system has independent eigenvectors and, although it is not strictly required, that the system eigenvalues are all complex. The usual procedure is to assume a solution of the form (2) and substitute this into equation (I). This leads to a second order N x N eigenvalue problem which can be solved for the eigenvalues, A., and the eigenvectors, q. When C is non-zero, the 2N eigenvectors are generally complex and the N equations of motion do not uncouple in this second order form except in special circumstances [START_REF] Caughey | [END_REF]. A first order formulation which includes the velocities as separate dependent variables provides the setting needed to uncouple the dynamics [START_REF] Meirovitch | Analytical Methods in Vibrations[END_REF].

Premultiplying equation (1) by the inverse of M and introducing the velocity vector, y, equation ( 1) can be written in first order form as z=Dz, z=(;} D=(O I) AB' [START_REF] Szemplinska | The Behaviour of Non/inear Vibrating Systems[END_REF] where A= -M-1 K and B = -M-1 C. We denote their elements by aij and f3ij (i,j= 1, ... , N), respectively. A solution of the form z( t) = z e..l.t [START_REF] Rand | [END_REF] is assumed which, when substituted into equation [START_REF] Szemplinska | The Behaviour of Non/inear Vibrating Systems[END_REF], leads to a standard, first order, 2N x 2N eigenvalue problem which can be solved for the 2N eigenvectors, z, and the 2N eigenvalues, A. The eigenvectors are not unique in that they provide only a direction in the phase space, but not a magnitude: each contains 2N-1 specified ratios and a free (complex) multiplicative constant. Thus, with this first order representation, a motion purely in a single mode is one for which all co-ordinates have exactly the same time dependence up to constant (complex) scaling factors, i.e., the amplitude ratios. An important property of such a modal motion is that if one knows the motion of a single generalized co-ordinate, then the motions of all other co-ordinates are specified by that mode's eigenvector. Furthermore, because of the independence of the eigenvectors, an initial condition which is aligned with a single mode leads to a purely single-mode motion in which all other modes start with zero amplitude and remain quiescent for all time. This is the important invariance property of the modal subspaces.

A limitation of equations ( 2) and ( 4) is that they implicitly assume that motion in a mode is synchronous, i.e., there are certain relationships between only displacements which must hold for a normal mode motion to take place. This appears to restrict the associated modal subspaces to be one-dimensional. However, normal mode motions for a damped or gyroscopic system are not generally standing waves, but traveling waves such that the system's degrees of freedom do not oscillate in phase. This means that relationships between both the displacements and the velocities must hold for a normal mode to occur, and that the associated modal subspace is two-dimensional. The consequence is that the corresponding eigenvectors in equation ( 4) must be complex and, therefore, in order to obtain a real normal mode, one requires a pair of eigensolutions.

To circumvent these difficulties, we seek the normal modes of the system by a formulation that exploits the invariance property of the modal subspaces. We define a motion in a normal mode as one where all displacements and velocities are related to a single pair of displacement and velocity, say x1 and Y1 (these are selected for simplicity; any pair can be chosen). Because the system is linear, this functional dependence of the x/s and y/s on Xt and Yt is linear:

i=1,2, ... ,N, (5a)
where a11 = 1, a21 = 0, b11 = 0 and b 21 = 1. In vector form, this is written as

1 0 xi ai2 a22 XN a iN a2N [ ;:} (5b) Y1 - 0 1 . b12 b22

YN biN b2N

This equation is that of a plane in the 2N-dimensional phase space. Any motion that takes place in such a linear invariant subspace and satisfies the equations of motion is defined as a linear normal mode. (If x 1 turns out to be a nodal point for a given mode, i.e., x 1 = y 1 = 0, the equations will be singular. This situation is easily remedied by selecting a different displacement/ velocity pair.) Note that this definition allows for phase differences between components, as required for "non-classical, normal modes.

Recalling that the equations of motion (3) can be written as

N x,= y;, y,= L (a0 1 + f3ijYJ ), j=l i= I, ... , N,
we substitute equation ( 5) into equation ( 6) fori= I, 2, ... , Nand make use of

N Yt = 2.: (atj(a tjXt +a2 1 y,) + {3ij(bt]Xt +b2 1 Yt)). j=l (6) (7)
Gathering terms in x 1 and Y1 in the resulting 2N-2 equations and requiring that coefficients of x, and Yt match gives 4N -4 equations that are quadratic in the unknowns a!i, a2,, bli and b21, i = 2, ... , N. These can be written in matrix form as 1 0 1 -0 a,2 an al2 an . . .

(~ ~) 0 1 a,N a2N aJN a2N N N (8) - 0 L (atia1 1 + {31jblj) L (alia2 1 + {31jb2J) . 0 I 1 b,2 b22 bl2 b22 j= l j • l . . . . . . . biN b2N biN b2N
Note that the time dependence of the problem has been eliminated by this procedure and that equation ( 8) can be solved for the geometry of the modal subspaces without solving the standard eigenvalue-eigenvector problem. For an oscillatory system, equation (8) has N real solutions for the linear coefficients (a 11 , a2;, ht 1 , b 21 ) 1 = 2 , ... ,N-one for each normal mode (see Appendix A for a justification; while there is a unique solution for each underdamped mode, there may be multiple solutions for overdamped modes, as their motions take place in one-dimensional spaces). To retrieve the modal dynamics of the system on each of the N invariant planar subspaces, we introduce the co-ordinate transformation

I 0 a\2 a22 ("<tl)-I a,N a2N (u.(t) )= Uw(l), (9)
y(t)

k= I 0 I vk(t) b,2 b22 . . biN b2N k
where U is defined as a 2N x 2N matrix theN pairs of columns of which are theN solutions of equation (8), such that the k th normal mode shape consists of the (2k -1 )th and (2k )th columns ofU. The vector w = [u 1 , v 1 , • • • , uN, vN]T contains the modal co•ordinates, where u1c and vk are the displacement and velocity, respectively, for the kth normal mode. We apply this co•ordinate transformation to the equations of motion (3) and use the property satisfied by the pairs of columns of U (equation ( 8)). This leads to a set of 2N, pairwise coupled, first order equations:

ul UJ 0 vi VI =U-1 DU - (10) UN VN
where in equation [START_REF] Vakakis | Analysis and identification of linear and nonlinear normal modes in vibrating systems[END_REF] we have defined

N N gk=an + L aiJ(aiJ)k+ L fJliblJ)k, j=2 }""2 N N
hk= L a1ia21)k+ fJ11 + L f31;{b2Jh. k= I, ... , N.

(ll)

j=2 j=2
Therefore, the dynamics on each invariant modal subspace is governed by the singledegree-of-freedom oscillator equation

( uk) ( 0 1 )(uk) vk = -roi -2~kf1h vk ' k= 1, ... , N, ( 12
)
where we have let gk = -ror and hk = -2~kmk, where rok and ~k are the undamped natural frequency and damping ratio for the kth normal mode, and uk and Vk are the modal displacement and velocity, respectively. We refer to the above as the oscillator form of the modal dynamics.

It is interesting to note that the selection of ( uk, vk) as the contribution of the kth mode to the displacement/velocity of the first degree of freedom has essentially forced a particular normalization on the eigenvectors. This normalization manifests itself in the following unusual manner: the eigenvalues, as represented in equation [START_REF] Caughey | [END_REF], depend on the eigenvectors. This feature may be initially bothersome, but it is a natural consequence of the formulation which specifies the eigenvectors uniquely.

Equation [START_REF] Vakakis | Analysis and identification of linear and nonlinear normal modes in vibrating systems[END_REF] shows that the dynamics of theN-degree-of-freedom system reduces to that of N uncoupled oscillators. For given initial conditions z(O) = Uw(O), the N second order initial value problems [START_REF] Meirovitch | Analytical Methods in Vibrations[END_REF] can be solved for the modal dynamics and the general motion obtained as a superposition of uncoupled oscillator motions via equation (9).

In this section we have used invariant subspaces to define normal modes and perform modal analysis for linear systems. Although this approach is novel and entirely different from the conventional methodology, we show in Appendix A that it is strictly equivalent to the standard formulation that is based upon the solution of a first order eigenvalue problem. All that is required to retrieve the (real) oscillator form of the modal dynamics is to perform a similarity transformation of the complex eigensolution of the first order representation. However, the power of the invariant subspace approach lies in its natural extension to invariant manifolds which represent normal modes for non-linear systems. This is tackled next.

NORMAL MODES FOR NON-LINEAR SYSTEMS

An invariant set for a dynamical system is defined as a subset S of the phase space such that if the system is given an initial condition in S, the solution of the governing equations of motion remains in S for all time. The non-linear normal modes defined below are invariant subspaces for the non-linear equations of motion. These subspaces are, in general, non-planar manifolds which are tangent to their linear counterparts, the planar eigenspaces, at the equilibrium point.

The equations of motion are assumed to be of the form x;= y;, y; j;(x; y), i= 1, 2, ... , N, [START_REF] Carr | Applications of Centre Manifold Theory[END_REF] where x = [x 1 , x 2 , ... , xN ]T represents generalized co-ordinates (displacements or rota- tions) and y = [y 1 , y 2 , ••• , YN ]T contains the corresponding generalized velocities. The vector f= [Ji, ... ,fN ]T represents the forces and moments acting on the system normalized by the respective inertias. (In some applications one has to invert an inertia matrix which may depend on x and y in order to achieve the form in equation ( 13). Also, it may be that for a given problem it is more natural to use generalized displacements and their conjugate momenta. In this case a more general, first order formulation may be required. The steps below can easily be generalized to handle that case.). We begin by assuming that there exists at least one motion for which all displacements and velocities are functionally related to a single displacement-velocity pair, which we choose arbitrarily here as the first displacement and velocity, x 1 and y 1 • In order to implement this, we write u=x 1 and v= y, and express the other x;'s and y/s functionally in terms of u and v, as follows :

x;=X,{u, v), y;= Y,{u, v), i = 1, 2, 3, . . . , N,

(14a)
where X 1 (u, v) = u and Y 1 (u, v) = v. In vector form this is written as:

XI u Yt V x2 X 2 (u, v) Y2 - Y 2 (u, v) (14b) XN XN(u, v) YN YN(u, v)
(We note that such a representation is not always possible. Cases for which it breaks down include all of those for which the linearized system cannot be uncoupled; for example, if multiplicities of eigenvalues occur with an associated non-zero nilpotent form, and situations in which other types of internal resonances occur, for example when there exists an unremovable non-linear coupling between the interacting modes. However, the method can be extended to handle these situations.) Equation ( 14) is that of a constraint surface of dimension 2 (or, equivalently, of eodimension 2N-2) in the 2N-dimensional phase space. Thus, we define a normal mode of motion for the non-linear, autonomous system as a motion which takes place on a twodimensional invariant manifold in the system's phase space. This manifold has the following properties: it passes through a stable equilibrium point of the system and, at that point, it is tangent to a plane which is an eigenspace of the system linearized about that equilibrium.

A set of equations which can be solved for this constraint surface, that is, for the X/s and Y;'s, can be obtained by requiring that solutions satisfy both the equations of motion and the constraint conditions. This is accomplished by eliminating the time dependence in the problem, which yields a set of equations for the geometry of the manifold. The following procedure for this task is borrowed from center manifold theory [13]. We begin by taking the time derivative of the constraint equations, and using the chain rule for differentiation, to obtain i= 1, 2, 3, ... , N .

(15)

Next, we substitute the equations of motion in for X; and Yi and replace x, and y; everywhere by X; and Y; to obtain 2N-2 equations which can be solved for X, and Y,:

[;(u, X 2 (u, v ), . .. , X N(u, v ); v, Y 2 (u, v ), . . . , Y N(u, v )) oY,(u, v) oY;(u, v) = au V+ av J.(u, X2(u, v), . .. , XN(u, v); V, Yz(U, v), ... ' YN(u. v)) i = 1, 2, 3, ... , N. ( 16
)
Note that fori= 1, the equations are trivially satisfied.

In general, these functional equations are at least as difficult to solve as the original differential equations, but they do allow for an approximate solution in the form of power series expansions. Once they have been solved for the X ;'s and Y/ s, the dynamics on the invariant manifold, which is the normal mode dynamics, can then be generated by simply substituting the X/s and Y;'s in for xi and y; in the first pair of equations of motion; that is, the ones for x 1 and y 1 • This results in the modal dynamic equation u=v, [START_REF] Abraham | Foundations of Mechanics[END_REF] where u and v represent the variables on the invariant manifold and correspond to projec- tions of the modal dynamics onto (x 1 , YJ). In general, at each equilibrium point there are N solutions for the X ;'s and Y/s and N corresponding sets of equations of the form given in equation [START_REF] Abraham | Foundations of Mechanics[END_REF], one set for each mode.

In some cases, for example, for systems which possess certain symmetries, equations ( 16) can be solved exactly, in which case global representations of the modal manifolds and their attendant dynamics can be obtained. These are typically when the manifolds are fiat and identical to the linear normal modes, even though the modal dynamics are nonlinear. These are the so-called "similar normal modes", as defined by Rosenberg, which have received a lot of attention in the past ; see references [1] and [START_REF] Caughey | [END_REF], for examples. Our formulation, on the other hand, generates the equations which must be satisfied for " nonsimilar" as well as " similar" normal modes, and allows for velocity dependence to be included. As is shown in the example in section 6.2, this velocity dependence arises even in conservative systems. For such systems, however, by invoking the conservation of energy, the non-linear velocity terms can be r. eplaced by terms involving the system energy and the modal displacement. This approach, typical of all previous work, allows one to represent a non-linear normal mode as a relationship involving only displacements and the total system energy.

APPROXIMATE SOLUTION FOR THE MODAL DYNAMICS

The approach taken here is local in nature, and yields approximations for the normal mode invariant manifolds and the dynamics on them near the equilibrium point. This restriction permits analysis to be done for systems without symmetry and allows for the inclusion of damping, which covers the most widely studied situations for structural and fluid mechanical systems. The procedure yields exact results for similar normal modes, since they are represented by flat manifolds, and asymptotic results for non-similar normal modes. The approximations are in the form of power series expansions and can, in principle, be generated to any order. For the present, we work to third order in the displacements and velocities. This will be adequate for most problems, although the method can be extended to higher orders in a straightforward manner.

The approach is based on Taylor series expansions about an equilibrium configuration of the system. The development given in section 3 above makes no explicit reference to an equilibrium point. Here we make the assumption that the equilibrium of interest is at x = 0; this is always possible by the proper selection of x. The behavior of the system near the equilibrium point is governed by

.:ij= yj, PJ= UjkXk + PJkYk + OjkmXkXm + BjkmXk)'m + YJkmYkYm + JLJkmqXkXmXq+ VjkmqXkXm)lq+ PJkmqXkYmYq+ ~jkmqYk)'mYq+ • • •, ( 18
)
where the coefficients in the expansions are derived from straightforward differentiation of the forces jj (}= 1, ... , n) with respect to x and y, and where the implicit summation notation is used. Note that, as written, the coefficients of the non-linear terms are not uniquely defined, but can be chosen in any way such that all terms up to the desired order are included. For example, for the x 1 x 2 terms one can take 8 1 1 2 to be the required coefficient, and then 8 1 21 to be zero.

It is now assumed that the normal modes for the non-linear system, as given in equation (14a), can also be expressed in the form of a Taylor series expansion [START_REF] Rand | [END_REF]:

X,{ u, v) = atiu + auv + a 3 ;u 2 + a4;uv + as;V 2 3 2 2 3+ + 06;U + a 7 ;U V + Og;UV + 09;V • • • , Y;( u, v) = bliu + b2;V + b3;U 2 + b4;UV + bs;V 2 + b6;U 3 + b7;U 2 V + bs;UV 2 + b 9 ;V 3 + ' ' ' , i= 1, 2, 3, ... , N. (l9a)
In vector form this series can be expressed as

+ 1 0 0 0 0 1 0 0 0 0 a32u + a42v as2V + h 32 u + b 42 v b 52 v a 3 Nu + a 4 Nv a 5 Nv b 3 NU + b 4 NV b sNV 0 0 a 62 u 2 + a 82 v 2 a 12 u 2 + a 92 v 2 b 62 u 2 + b 82 v 2 b 12 u 2 + b 92 v 2 (19b)
or, more compactly, as

z = {IIlo + m 1 (u, v) + m 2 (u, v)}(:) + • • •, (19c)
where z = [x 1 , y 1 , x 2 , y 2 , ... , xN, YN]T and ~. m 1 and m 2 are 2N x 2 matrices. (Note that this z is a simple reordering of the z in section 2). The matrix IIlo is the linear modal component, and m 1 and m 2 represent the effects of quadratic and cubic non-linear terms, respectively. Also, note that we have factored out a (u, v) vector and that the given representation is not unique; it is, however, convenient for some of the calculations, and all possible representations will yield the same results.

As described in section 2, the coefficients of the linear terms represent the ratios for the usual linear normal modes. For an undamped, non-gyroscopic system it can be shown that the cross-terms between displacement and velocity, a 2 ; and bu, equal zero, while au = b 2 ; represent the usual amplitude ratios (and thus velocity ratios) for a conservative vibratory system. For a damped or a gyroscopic system the linear coupling terms between displacement and velocity are generally non-zero and have no such special relationship with one another (other than being the linear approximation to the modal subspace), except in cases such as proportional damping, where the modes of the damped system are identical to those of its undamped counterpart. These non-zero cross-terms allow for phase differences between the displacements for non-conservative or gyroscopic systems.

The non-linear terms describe the bending of the modal subspace. Their associated coefficients can all be zero for a mode in which the amplitude ratios are fixed constants, as happens in systems which possess similar normal modes, but in general they are nonzero. These terms capture the effects of the non-linear forces, and result in the fact that the displacement and velocity ratios are dependent on the amplitude of motion.

The manipulations given in section 3 are now carried out for the series representations of the equations of motion and the normal modes. The elimination of time, resulting in equation ( 16), yields a set of 2N-2 equations which contain the coefficients of the forces as known quantities and the coefficients of the normal modes as unknowns. They are solved by gathering terms of like powers in u and v and requiring that their coefficients match. For cubic order terms, this leads to 18(N-1) equations for the unknown coefficients of the normal modes. These equations are listed in Appendix B. Note that the use of the chain rule results in coupled algebraic equations for the coefficients at each step. These equations can be solved sequentially, just as is typical in perturbation theory, and this procedure leads to some surprising observations which are stated below.

As can be seen from the equations given in section 2 and in Appendix B, the equations for the linear coefficients are quadratic in the unknowns a 1 j, a 21 , bv and b 21 , and generally have N real solutions, one for each mode. The solution of these quadratic equations presents potential computational difficulties. However, since their solution is equivalent to solving the eigenvalue-eigenvector problem, it is clear that the desired solutions are computationally obtainable by standard methods. The relationships provided in Appendix A can then be used to obtain the a's and b's from the standard eigensolution. The equations to be solved for the non-linear coefficients depend on the solution of the linear coefficients, and we note from Appendix B the very important fact that they are linear in the unknowns and thus generally have a unique solution for each mode. This implies that the methodology should _be easy to incorporate into computational schemes, even for large scale systems.

Solutions of these l8(N -1) equations result in cubic series approximations for theN non-linear normal modes. This procedure provides the geometric structure of the nonlinear normal modes near the equilibrium point. The local approximation of the dynamic equations for each mode can then be constructed via equation [START_REF] Abraham | Foundations of Mechanics[END_REF], by using the series expansion for / 1 given in equation [START_REF] Pierre | [END_REF], with (xh yj) replaced everywhere by the series representations of (Xiu, v), ij(u, v)) for j= I, 2, ... , N. This results in a single-degree-offreedom non-linear oscillator equation which represents the dynamics of the system on an invariant, two-dimensional subspace which is tangent to the linear normal mode eigenspace at the equilibrium point. There will beN of these oscillators, one for each normal mode.

We are now in a position to assemble the complete non-linear modal matrix M and the affiliated modal co-ordinates w. The matrix M is composed by gathering together into a 2N x 2N matrix the N individual modal matrices of size 2N x 2 given in equation (19). Similarly, the modal vector w is constructed by noting that each mode will have its own (u, v) pair associated with it, which we label as (uk! vk) for k= 1, 2, ... , N. These are collected into the vector w = [ u, , v, , u2, v2, ... , uN, v N ]T. The complete transformation between the physical co-ordinates z and the modal co-ordinates w can now be written as

z = M(w)w = {Mo + M 1 (w) + M 2 (w)}w + • • •, ( 20
)
where M(w) = Mo +M, (w) + M 2 (w) +•••,and where the M/s are assembled from the m/s from equation (19c) for p = 0, 1, 2 (note that Mo is simply U 0 of section 2). It is important to note that the coefficients aii and bii will be different for each mode, and if one were to write out the components of equation ( 20), coefficients of the form aiik and bifk should be used where, for the present case, i = 1, 2, ... , 9 corresponds to the number of terms retained in the expansions (here we include all quadratic and cubic terms), j = 2, 3, ... , N corresponds to the number of constraints, and k = 1, 2, ... , N to the number of modes.

There are several important features of this formulation. It is the most natural local description possible for a non-linear normal mode. It is constructive, so that one is able to use it to study specific non-linear systems. It provides information regarding non-linear behavior near an equilibrium, including non-linear frequencies and decay rates. These are extremely useful for determining transient response and set the stage for the behavior which can be expected when external excitation is present. Also, it comes from a more global representation which may allow for other solution strategies to be used to obtain non-local results for non-linear normal modes. One way to do this is by generating numerical solutions of equations ( 16).

TRANSFORMATION FROM PHYSICAL TO NON-LINEAR MODAL CO-ORDINATES

It must first be stated that the approximate superposition described here is simply a non-linear extension of the usual linear concept and is not, of course, strictly the same as that in linear systems. Basically, it is a non-linear co-ordinate transformation which performs the same function as the well-known linear modal transformation in that it permits one to assemble a complete solution from a sum of simpler ones. Furthermore, it reduces to its linear counterpart when non-linearities are not present or are ignored. In effect, linear superposition is the linearized version of the non-linear co-ordinate transformation described here.

The geometry behind this idea is quite simple in concept. In linear systems, one has planar eigenspaces, and a general response can be broken down into modal components via a linear projection onto the eigenspaces, or reassembled by a linear recombination of the modal responses. When non-linearities are present, the modal subspaces are not flat, but curved, and the required projections and recombination must utilize curvilinear coordinates which reflect the nature of these subspaces.

For linear systems the co-ordinat~ transformation between physical and modal coordinates is given by the modal matrix M 0 above. When the equations of motion are written in terms of the modal co-ordinates, they are uncoupled, linear oscillators, as in equation [START_REF] Meirovitch | Analytical Methods in Vibrations[END_REF]. It is of interest to consider what happens if one writes the non-linear equations of motion in terms of the non-linear modal coordinates. In the following, we show how this transformation can be carried out and describe some general features of the modal dynamic equations. Section 6.3 provides an example for which these features are explicitly demonstrated.

Recombination of the dynamics of the original co-ordinates from the modal co-ordinates for the non-linear case is done via the transformation z=M(w)w given in equation (20).

In order to obtain the projection of the physical co-ordinate dynamics onto the modal subspaces, this transformation must be inverted. The inversion is carried out for the general case involving quadratic and cubic non-linearities in Appendix C. The inversion is much simpler when only cubic non-linearities are present, and can be achieved by the following iterative process, which exploits the fact that we are using series expansions.

We begin by premultiplying equation ( 20) (with Mt("w) equal to zero when no quadratic terms are present) by the inverse of the bracketed terms on the right side and expanding the inverse to yield

w = {Mo + M 2 (w)} -Iz + • • • ={I+ M 0 1 M 2 (w)}-1 Mo 1 z + • • • ={I-M; 1 M 2 (w)}M; 1 z + • • •, (21)
where I is the 2N x 2N identity matrix. Note that the right side depends on both w and z.

This can be remedied by substituting the expression for win equation [START_REF] Wiggins | Introduction to Applied Nonlinear Dynamical Systems and Chaos[END_REF] into thew term on the right side and expanding in powers of z. Since M 2 (w) is quadratic in w, the leading order argument of M 2 is simply Mo-1 z, and the dependence on w is pushed out to terms of higher order. This results in the inverse transformation [START_REF] Guckenheimer | Non/inear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF] which is correct up to cubic terms in z, the order to which we are working. This is used in the examples worked out in section 6, where we consider only cubic order non-linearities. We now transform the equations of motion to a modal co-ordinate representation by using a direct co-ordinate transformation. The equations of motion near the equilibrium at z = 0 are written as

z=A(z)z= {Ao+ A1(z) + A2(z)}z+ • • •, ( 23
)
where Aoz represents the linearized dynamics, A 1 (z)z the quadratic terms and A 2 (z)z the cubic terms. The transformation (20) is applied directly to these equations. We write it in a slightly different way for notational convenience here, such that the w vector is not separated out as a factor: z=IVt(w), [START_REF] Shaw | Normal modes of vibration for nonlinear continuous systems[END_REF] where it is implicitly assumed that we use the series expansion for Q. The time derivative of this is required for the transformation; it is given by using the chain rule:

z=(aQ(w)/ow]w. ( 25
)
Substituting of equations ( 24) and ( 25) into the equation of motion [START_REF] Fenichel | [END_REF] yields

[oM(w)jow]w = A(M(w))M(w), (26)
which can be premultiplied by the inverse of the Jacobian of the transformation to give the equations of motion in terms of the modal co-ordinates: w = (oM(w)/ow)-1 A(M(w))M(w).

(27)

By using the series approximation technique given in this paper, these equations can be constructed to the desired order. Each of the N pairs of equations have the form of equation ( 17), plus additional modal coupling terms. These coupling terms are identically zero on each individual modal subspace, thus recovering the modal oscillators on the invariant manifolds.

Even though the modal equations are coupled, we explore in section 6.3 the possible advantages of simply ignoring the coupling terms in the transformed equations in order to recombine the non-linear modal oscillator responses into a general motion. This is done by solving a non-linear initial value problem as follows. Consider the motion which ensues from an initial condition Zo. This initial condition can be transformed to modal co-ordinate initial conditions, w 0 = [u 10 , v 10 , u 20 , V:w, •.• , uNo, vN 0 Y, by using the inverse transformation given in Appendix C (or equation [START_REF] Wiggins | Introduction to Applied Nonlinear Dynamical Systems and Chaos[END_REF] if only cubic non-linearities are present). These initial conditions are for the coupled non-linear modal oscillators given in equation ( 27). These coupled equations could be numerically integrated to obtain the modal dynamics, (uk(t), vk(t)), which can be transformed back into the original co-ordinates via equation ( 20), but since these equations are coupled, the only advantage in doing this is that one can obtain information about modal amplitudes directly. In section 6.3 we compare results from a full simulation of the original equations with simulations and recombinations of the individual modal oscillators in which coupling is ignored, i.e., equation [START_REF] Abraham | Foundations of Mechanics[END_REF], and show that very good agreement is obtained for that example. Such an approach has no sound mathematical basis, but provides a means of constructing the motion of a non-linear system from its non-linear normal modes.

A more important feature of the equations of motion expressed in terms of the nonlinear modal co-ordinates is that very "clean" reduced order models can be obtained by modal truncation. Since the non-linear normal modes were constructed using invariance, the dynamics of a subset of k of theN modes takes place in an approximately invariant subspace of dimension 2k. By ignoring the energy in the truncated modes and considering only k modes to be active, one can generate a reduced order model which accounts for the non-linear coupling between the active modes but is less contaminated by the truncated modes. An explicit demonstration of this is provided in the third example (section 6.3). This fact indicates that these non-linear normal modes are natural candidates for performing modal analysis of non-linear systems, a subject currently under investigation.

EXAMPLES

Three examples are presented. The first is a damped, linear system with two degrees of freedom; it is used to demonstrate how the method handles linear oscillators. The second example is an undamped non-linear system of two degrees of freedom with no special symmetry. The third is a damped non-linear oscillator of two degrees of freedom which combines the damping from the first example with the non-linearity from the second.

Most of the symbolic manipulations and all of the numerical work for the examples were done using Mathematica® [START_REF] Wolfram | Mathematica[END_REF], running on a Macintosh Ilcx computer. Mathematica® uses a variable step size fourth order Runge-Kutta integrator for the numerical solution of ordinary differential equations, which was tricked into using a constant step size of 0•1 in all cases. The results, as seen below, are quite acx:urate for this step size, and the constant step size was useful when manipulating some of the data.

EXAMPLE 1: A LINEAR, NON•CONSERVATIVE SYSTEM

The linear system to be considered is shown in Figure 1. It has no special symmetry. After setting up the general equations which are to be solved for the invariant modal subspaces, we take specific values for the inertia, stiffness and damping elements, since the solution for the coefficients is not easily obtained (if at all possible) in closed form, except in special cases. We assume that the springs are unstressed in the equilibrium position and that the displacement co-ordinates x 1 and x 2 are zero at this point. The equations of motion in first order oscillator form are il=y,, x2= Y2.

Y1 = /1(x1, X2; y,, Y2) = -xi(l + k) +kx2-c(y1-Y2), Y2 = /2(x,, x2; y,, Y2) =kx,-x2(1 + k) +cy,-2cy2, or, equivalently, in matrix form 0 1 0 0 -1-k -c k c :i = Aoz. - 0 1 or 0 0 k c -I-k -2c
A purely modal motion is one which can be expressed in the following form:

y,=v, (28a) (28b) (29)
The procedure for solving for the coefficients is to first take the time derivative of the last two of these equations (the last 2N-2 in general, see equation ( 15)) to obtain (30)

The left side is then equated to (.i 2 , y 2 ) from the equation of motion (28), and in the right side (u, v) is replaced by (.X 1 , YI) from equation (28). This gives (31)

Next, x 1 , y,, X2 and y 2 are replaced everywhere by the constraint (29) and the right side is subtracted from both sides to produce the final equations which are to be solved for a 1 ,

2 (u, v)-a 1 v-a 2 J;(u, X 2 (u, v); v, Y 2 (u, v)) = 0, j 2 (u,X2(u, v); v, Y 2 (u, v))-h 1 v-aJj(u,X2(u, v); v, Y 2 (u, v))=O, (32a)
or, after substituting in for fi , !2, X2 and Y2, u(a2+b1 -a2b,c+a2k-a,a2k) + v(-a, +h2+a2c-a2b2c-~k)=O, u(-a, +b2-2b1c-b,b2c+k -a,k+b2k-a1b2k)

+ v(-a2-bl + c-b2c-bk-a2k-a2b2k) =0. ( 32b
)
There are four unknowns in these equations; these are determined by noting that the equations must hold for all values of u and v.

Gathering the coefficients of u and v in equation (32b) yields a set of four equations to be solved for a1 , a2, b1 and b2 :

(33) Note that these equations are quadratic in the unknowns. Neither the authors nor Mathematica® were able to solve them in dosed form for general values of c and k. Therefore, for the remainder of the example we take the special case c=O• 30 and k = 1•0. For these values, Mathematicaao finds six solutions of equations (33), two of which are real and four of which are complex. The complex roots are meaningless for the problem formulation given. The real roots represent the linear normal modes:

mode 1: Dt =-1•1056, mode2: Ot =0•9471, a 2 = -0•1550, a2 == -0• 1386, b, =0•4599, b, =0•1401, b2 = -0•9891, b2 == 0•9676. ( 34
)
This means simply that a motion in the first mode satisfies the conditions

x2 = -1•1056x,-0•1550y 1 , Y2 = 0•4599x 1 -0• 989ly, , (35)
and a motion in the second mode satisfies

x 2 = 0•947lx 1 -0•1386yl, y2=0•140lx 1 +0•9676y 1 • (36)
As the damping factor c is reduced towards zero, the first mode approaches the mode in which the masses are exactly out of phase, while the second mode approaches the exactly in-phase mode. The modal subspaces are composed of two planes in the four-dimensional phase space. These are depicted in Figure 2 as the planes X2 and Yz versus (x1, y,) for each of the two modes.

The equations of motion which govern the individual modes are given by substituting conditions (35) and (36) into the (x,, y 1 ) component of the equations of motion (28). For the first mode this yields

61 = -2•9676ut-0•7517v,, or u, +0•7517ut +2•9676ut =0,
and for the second mode

v2 = -1•0109u2-0•1483v2, or ih + 0•1483u2 + 1•0109u2 = 0.
At this point the eigenvalues are easily determined to be mode 1: -0•3759 + 1•6812i, where F = -1. Note that by this procedure one obtains the eigenspaces first, and then the eigenvalues associated with them.

The transformation from the physical co-ordinates to modal co-ordinates is given by

X1 1 0 1 0 Y1 0 I 0 1 - or z =Mow, ( 40
) x2 a11 a21 a12 a 22 Y2 b11 b21 h12 b22 v2
where a ij and bij represent a; and b; for the jth mode. The inverse of the transformation is given by u, 0•4602

V1 0•1468 - u2 0•5398 v2 -0•1468 or -0•0714 0•4828 0•0714 0•5172 -0•4865 -0•0795 0•4865 0•0795 w=~-1 z. 0•0041 -0•5104 -0•0041 0•5104 (41)
Here we can explicitly transform the original equations of motion to modal co-ordinates to obtain u• 0

1 0 0 UJ v, -2•9676 -0•7517 0 0 Vt - r.h 0 0 0 I u2 v2 0 0 -1•0109 -0•1483 v2 or w =M;' AoMow, (42)
which is exactly the matrix form of equations ( 37) and (38). These are uncoupled oscillators which are not in the usual Jordan form, but in what we have called the oscillator form (see section 2 and Appendix A).

It is clear that a motion started with initial conditions in one of the modal eigenspaces will result in motion comprised of only that mode. This is demonstrated in detail below for the non-linear examples, but we show one such case here in order to provide a calibration on the precision of our numerical methods, e.g., equation solvers and differential equation integrators. To this end we take an initial condition which is purely in the first mode :

w 0 =(1,0,0,0)T ~ Zo=(1,0,-1•1056,0•4599)T. ( 43
)
A simulation of the equations of motion (28) started at this point results in the response given in Figure 3, which shows the modal phase planes (u 1 , v 1 ) and (u2, v2). If the numerics were exact, the motion should remain exactly at the (u2, v2) origin, and it is observed that the response of the second mode stays bounded within a disk with a radius on the order of 10-17 -----l ------------J -2 x w-17 2 X 10-l? It is worth noting that if the damping coefficient, c, is increased to the point where both modes are overdamped, the quadratic equations (33) have six real solutions. These arise because we are searching for two-dimensional invariant subspaces, while the completely overdamped system possesses four one-dimensional invariant subspaces. In this case the solutions of equations (33) correspond to all possible pairings of the one-dimensional spaces into two-dimensional spaces, of which there are six possibilities. In such a situation, care must be taken in choosing which solution to use, since all four eigenspaces must be accounted for.

The application of superposition for this example is standard and straightforward and is not given. We do so for the non-linear system of Example 3.

EXAMPLE 2: A NON-LINEAR, CONSERVATIVE SYSTEM

A non-linear, conservative system is now considered. For conservative systems such as this one, the existence of invariant modal subspaces with associated non-linear oscillators has been known and is stated in the Lyapunov sub-center manifold theorem [START_REF] Abraham | Foundations of Mechanics[END_REF]. Here we provide a constructive tool for generating them. One important consequence of our formulation is that conservative, velocity dependent terms appear in the normal modes and the modal oscillators.

The physical system is simply that of Example 1, with damping mechanisms removed and the leftmost spring made non-linear with a cubic non-linearity; see Figure 4. Again, this system has no special symmetries to help with the non-linear analysis. The equations of motion are given by x1=y~.

(44a) or, equivalently, in first order form as .XI 0

I 0 0 0 0 0 0 X1 YI -(1 +k) 0 k 0 2 0 0 0 + -gxt Y1 - 0 0 0 i2 0 l 0 0 0 X2 Y2 k 0 -(1 +k) 0 0 0 0 0 Y2 or :i = { Ao + A2(z) }z. (44b)
A normal mode motion is assumed by requiring that XI, Y1, x 2 and Y2 are related as follows:

XI =u, Y1 =v, x2=X2(u, v) =a,u+a2v+a 3 u 2 +a4uv+a 5 v 2 + a6u 3 + a1u 2 v + asuv 2 + a 9 v 3 + • • •, Y2= Y2(u, v) =b1u+b2v+b3~+b4uv+b 5 v 2 + b6u 3 + b 7 u 2 v + bsuv 2 + b 9 v 3 + • • •. ( 45
)
The procedure goes exactly as in Example 1, but involves many more terms. For example, taking the time derivative of X 2 and Y 2 from equation (45) yields, by using the chain rule,

X2 = a1u + a2v + 2a3uu + a4uiJ + a4uv + 2asviJ + • • •, (46)
where we have shown only up to quadratic terms for conciseness of presentation, although cubic terms have been included in the analysis. The left sides are equated to the equations of motion for .i2 and y 2 , respectively, and in the right sides u and v are replaced by .i, and y 1 , respectively, from the equations of motion. This results in y2=aiY1 +a2/J(x,, x2; Y1, Y2) + 2a3XIYI +a4X,/J(xi, x2; y,, y2)

+ 04.YI + 2asy,/J (x,, X2; Y1, Y2) + • • •, h(x, , X2 ; Y1, Y2) = b,yi + b2/J (x, , X2 ; y,, Y2) + 2b3X1Y1 + b4X1/J (x,, x2 ; YI, Y2) + b4.Yi + 2bsy,/J(xi, X2; Y1, Y2) + • • •, (47)
where (x 1 , y 1 ) are used for (u, v) for notational consistency at this point. The substitution of conditions ( 45) into these expressions yields the required conditions for a modal motion.

These are then expanded out to cubic order in u and v, and coefficients of u, v, u 2 , uv, v 2 , u 3

, zlv, uti and v 3 are gathered together to provide the equations for the a, 1 and bj coefficients. The right side is subtracted from the equation in order to give equations which are equal to zero. The results of this calculation are given here in the form where the left column describes the term from which the equation on the right was obtained, such that there are two equations for each term, one each from the two equations in (47):

u term: a2+b1 +a2k-a,a2k=O,

-a1 + b2 + k-a1k + b2k-a1b2k = 0, v term:
-a,+ b2-a~k =0, -a2-b, -a2k-a 2 b2k=O, u 2 term:

a4 + b1-a 2 a 3 k + a.J<_-a,a,J<. = 0, -a3 + b4-a:J<-a3b2k + b,J<.-a,b,Jc = 0, uv term:

-2a 3 + 2a 5 + b4-2a 2 a~ + 2a 5 k-2a 1 ask = 0, -a4-2b3 + 2bs-a~a 4 b2k-a 2 bJc + 2bsk-2a,bsk = 0, -a4 +bs-Ja2ask=O, -as-b4ask-asb2k-2a2bsk = 0, u 3 term:

a,+ b 6 + a,g-a 3 a,J<.-a2adc + a,k-a,a,k = 0, -a6 + b, + bu;-a,)(-aJJ2k-a3bttk + b,k-a.b,k = 0, u 2 v term: -3a6 + 2as + b,-~Jc-2aJask-2a2a1k + 2ask-2a,ask = 0,
-a,-3b6 + 2bsa,k-a,b2k-a4b,Jc-2a3bsk-a2b1k + 2bsk-2a,bsk = 0, uv 2 term: -2a, + 3a 9 + bs-3a~sk-Ja2ask+ 3aJc-3a,aJc= 0, -as-2b7 + 3b9-ask-a 8 b2k-a 5 hJ<-2a4bsk-2a2bsk + 3bJc-3a,bc}c = 0, v 3 term:

-

a 8 + b9-2trsk-4a2ac)< = 0, -a9-hs-a~-a9b2k-2asbsk-3a2b~ = 0. (48)
The first four equations can be solved for the linear mode shape coefficients (again there are six solutions of which we take the two real ones):

a,=b 2 =1, with a 2 =b 1 =0,

Ot = b2 = -1, with 02 = b1 = 0,
for mode 1;

formode2.

These check with the modes of the linea~d, symmetric system.

(49)

The next six equations are uncoupled and are linear in the unknowns a 3 • a 4 , a 5 , b 3 • b 4 and bs. The problem is homogeneous and has a zero solution for these coefficients. This is not unexpected since the physical system has no quadratic terms. The last eight equations cannot be coupled to any higher order coefficients since we have not included any in our analysis. Also) based on the patterns observed in the equations, they will not be in any case. These equations are linear in the coefficients a 6 , a,, a 8 , For mode I:

a6 = g(k-3)/[2k(k-4)], as==-3g /[2k(k-4)],
For mode 2:

a6:::: g(3 + 7k)/[2k( 4 + 9k)], as= 3g/[2k(4+ 9k)],
b7=3g(l +3k)/[2k(4+9k)J, These give, to cubic order, the coefficients which describe the shapes of the invariant manifolds corresponding to the non-linear normal modes. The two normal mode manifolds are therefore defined by, to cubic order:

For mode I:

1 0 0 1 l 0 0 1 0 0 + g[(k-3)ui-3v~]/[2k(k-4)] 0
For mode 2:

1 0 0 1 --1 0 0 -1 0 0 0 0 0 3g[(k-1)uf-vfJ/[2k(k-4)] 0 0 0 + g[(3 + 7k)u~ + 3v~]/[2k(4 + 9k)] 0 3g[( 1 + 3k)u~ + v~]/[2k( 4 + 9k)J ( 51
)
(u 1 , v 1 ) and (u 2 , v 2 ) are the first and second mode displacement and velocity co-ordinates, respectively. These are assembled to form the transformation between physical and modal co-ordinates :

(52)

The inversion of this transformation was obtained in closed form to third order and is given here:

(53a) Xt + x2 [<k-3)(xt + x2) 3 -3(yJ + Y2) 2 (Xt + x2) Ut= 2 g 32k(k-4) (3 + ?k)(x,-x 2 ) 3 + 3(yt-Y2f(x1-x2)] + 32k(4 + 9k) ' v =y1+Y2_ 3 [(k-l)(xi+x2) 2 (yi+Y2)-(y,+y2) 3 I 2 g 32k(k-4) + (1 + 3k)(xt-x2) 2 (YI-Y2) + (Yt-Y2) 3 ] 32k(4+9k) ' Xt-X2 [(k-J)(XJ + X2) 3 -3(yl + Y2) 2 (X1 + X2) u - + g _;__ _ _ _ ____;_ _ ____;_:..__.__._::_.___;__ _ _ _ 2 - 2 32k(k-4) (3 + ?k)(xt-x2) 3 + 3(yl-Y2) 2 (Xt-x2)] + 32k(4 + 9k) , . . v =Y1-Y2+J [(k-l)(xi +x2) 2 (y1 +y2)-(y1 +y2) 3 2 2 g 32k(k-4) (1 + 3k)(xJ-x2) 2 (y,-Y2) + (yJ-Y2) 3 ] + . 32k(4+9k) (53b)
Note that, as required by this formulation, u1 + u2 = x1 and v 1 + v 2 = y 1 . This inverse is needed only if a specific motion is being studied, and that one can obtain the important qualitative features of a general motion by constructing the non-linear modal oscillators, which requires only the forward transformation (52). The non-linear modal oscillators are obtained exactly as was done in the first example, by substituting the constraint equations for x 2 and y 2 into the first pair of equations of motion. For the present example, this gives, in second order form,

mode 1: ii.t +Ut+ gu,[(t + ( 3 -k) )uf + 3 u7] =0. 2(k-4) 2(k-4) d .. ( 2k [( (3 + 7 k) ) 2 3k . 2] mo e 2: u2 + I+ )u2 + gu2 1-2 ( 4 + 9 k) u2-2 ( 4 + 9 k) U2 = 0. ( 54
)
For the present example, we demonstrate the usefulness of the method in determining the frequency-amplitude relationship. This is accomplished using straightforward perturbation methods on the non-linear oscillators. Lindstedt's method, as described in many non-linear oscillations textbooks, is used here. We present only the final results since the technique is well known and rather straightforward (see, for example, the procedure in reference [START_REF] Rand | Perturbation Methods, Bifurcation Theory and Computer Algebra[END_REF]). The results are the following relationships between the frequency of oscillation, {J), and the peak displacement amplitude, A:

mode 1: w = 1 + 3gA 2 /I6 + • • • mode 2: w = [ l + 2k] 1 ' 2 1 + + ... [ 3gA 2 J 16(1 + 2k) . ( 55
)
These have been checked with simulations and are quite accurate out to amplitudes of order unity.

Simulations were performed to check on the accuracy of the invariance properties of the non-linear normal modes, especiaJly in comparison with the linear modes, and for the frequency-amplitude relationships. Four simulations are presented. The first two are results which use the linear modes and the second two use the non-linear modes. For all cases, we took the parameter values of k= 1•0 and g=0•5, a significant non-linearity.

For initial conditions on the first linear mode, we take (x 1 , y 1 , x 2 , y 2 ) = ( 1, 0, I, 0), which are sufficiently large amplitudes to cause non-linear effects. The results of a simulation carried out for 20 time units were stored in a file and then projected onto the linear modal subspaces. They are shown in the two modal phase planes in Figure 5(a). Results for a similar simulation using initial conditions on the second linear mode, (1, 0, -1, 0), are given in Figure 5(b).

Analogous initial conditions on the first non-linear normal mode are given by (1, 0, 1 + 1/6, 0) to cubic accuracy, and for the second non-linear mode the required initial conditions are (1, 0, -1 + 5/26, 0). Results from simulations started from these initial conditions are shown in Figures 5(c) and 5(d), respectively. In these cases, the solutions are projected back onto non-linear modal co-ordinates using equation (53b).

The non-linear normal modes provide an improvement of accuracy of about one order of magnitude in determining the invariant motions for this amplitude. Also, the frequencies of oscillation for the modal motions were well within 1 °/o of the approximations given in equation ( 55), whereas the linearized frequencies were off by 8% and 4%, at this amplitude.

It is also interesting to note that for this (and other) conservative, non-gyroscopic system(s), the modal motions are synchronous in the sense that both (all) masses reach their peak amplitudes simultaneously and pass through zero displacement simultaneously. Such ideas were used by Rosenberg [1] in order to define non-linear normal modes, while in the present formulation they are a special consequence for conservative, non-gyroscopic systems. For this example, this can be seen directly from the equations for the normal modes, equation (51). It is observed that whenever ui is zero, both x 1 and x 2 are necessarily zero, and whenever vi is zero, both y 1 and y 2 are necessarily zero. These facts imply that both modes are synchronous in the sense described above. Such modes allow for the definition of modal amplitude ratios which are amplitude dependent. These can be obtained by noting that x, = ui and then computing x 2 / x 1 for each mode with v; = 0, which yields the ratio of peak amplitudes for a modal motion. The results are:

x 2 g(k-3) 2 mode 1 : -= 1 + x 1 + • • •, (56 a) X1 2k(k-4) mode 2: x 2 = -1 + g( 3 + 7 k) xi+ • • •. X1 2k(4+9k) (56b)
It is interesting to note that the first mode amplitude ratio can be softening or hardening depending on the value of the coupling stiffness:

k < 3, hardening; 3 < k < 4, softening; 4 < k, hardening;
while the second mode is hardening for all k > 0. This behavior is illustrated in Figure 6, which displays the amplitude x 2 versus x 1 for various values of k. - 5. Simulation validation of the invariance of the non-linear normal modal subspaces; in all cases an initial modal amplitude of 1•0 was used: (a) trajectory as projected into the individual linear modal phase planes for an initial condition started in the first linear mode; (b) trajectory as projected into the individual linear modal phase planes for an initial condition started in the second linear mode: (c) trajectory as projected into the individual non-linear modal phase manifolds for an initial condition started in the first non-linear mode; (d) trajectory as projected into the individual non-linear modal phase manifolds for an initial condition started in the second non-linear mode.

1•0 -0•5 0 0•5 1•0 1•5 2•0 -0•3 -0•1 0•1 ().3 0•5 U t Ut i •O 1•5 1•0 0•5 .... .... 0 ::0. ::0. -0•5 -1•5 -2•0 0•5 -2•0 -1•5 -1-0 -0•5 0 0•5 1•0 1•5 2•0 Uz Uz (c) (d) 0•5 I' I I I I - ().J r- - 0•1 r- - ... - ::0. ;a.. -0•1 r- - -0•3 ,... - -2•0L._---I_---l... _ _._ _ _.___.~....____.._----L.._~ -2•0 -1•5 -1•0 -0•5 0 0•5 1•0 1•5 2•0 I J I I I -0•5 -0•5 -0•3 -0•1 0•1 0•3 0•5 Ut 0•5 I I I I 0•3 - - f"4 ::0. 0•1 r- t"- - -0•1 1- \ "' - -0•3 r- - -0•5 I I I I -0•5 -0•3 -0•1 0•1 0•3 0•5 -1•5 -1•0 -0•5 0 0•5 1•0 1•5 2.0 U2 Uz Figure
x, Figure 6. Bending of the first normal mode amplitude ratio for various values of the coupling stiffness. The peak amplitude of x 2 is plotted versus that of x 1 for a modal motion in the first mode. Some closing remarks for this example follow. First, note that when k = 4, the nonlinear terms in the first mode become singular. This is due to the fact that for k = 4, a three-to-one internal resonance exists, in which case the motion of the non-linear system is known to involve important and unremovable coupling between the modes. In the present case, the modal reduction does not help at all since there are only two degrees of freedom. However, in higher order systems, the crucial coupling will be only between the resonant modes and one can reformulate the problem and look for an invariant fourdimensional manifold on which the behavior of the coupled modes can be studied. (The k= -4/9 singularity involves a similar internal resonance for an unstable system).

Another singularity occurs when the coupling spring stiffness approaches zero. This corresponds to mode localization caused by the symmetry-breaking nature of the nonlinearity. This is essentially a point of one-to-one internal resonance. See references [START_REF] Pierre | [END_REF]19] for a treatment of localization in linear systems and reference [20] for results on non-linear localization.

Finally, it is important to note that if one were to assume that x 2 depended on x 1 in some non-linear manner, but was independent of y,, there would be terms missing from the series expansions for the invariant manifolds and the modal oscillators. However, for conservative systems, the velocity dependent terms are accounted for by making use of the conservation of energy. This makes the connection between our approach and the traditional one for conservative systems, where the modes are sought in terms of displacement relationships and the system energy. The invariant manifolds obtained by the two approaches are not the same, and in fact exist in different spaces, since our formulation does not utilize the energy integral to reduce the dimension of the phase space. The present formulation is, however, the one required for qon-conservative systems, and yields correct normal modes for very general systems, including conservative ones.

EXAMPLE 3 ~ A NON~LINEAR. NON~CONSERVA TIVE SYSTEM

In this example, we combine non-linearity and damping to show how the method works in a more general case. We skip many of the details given in the previous examples, since the calculations are very similar, and instead focus on the approximate superposition aspects of the non-linear normal modes.

The physical system, shown in Figure 7, is the same as that for Example 2, except that the dash pots of Example 1 are added in the same places as in Example 1. The equations of motion are given by those in equation ( 28), with the simple addition of a term -gxt in the Y1 equation. The procedure following the development of equations ( 45)-( 47) is repeated, and the coefficients of the various terms involving u and v are collected to produce a set of 18 equations for the 18 ai and b 1 coefficients, just as in equation ( 48).

These equations are more complicated than those in Example 2 and reduce to them when c = 0. We do not give the details here, since a numerical approach is adopted with the parameter values c = 0• 3, k = 1• 0 and g = 0• 5.

For this example we find solutions as shown in Table l. Note that the quadratic terms are all zero but that we now obtain non-zero coefficients for all cubic order terms. The modal subspaces, approximated to cubic order, are presented in Figures 8(a)-(d), just as was done for Example l. Note that the spaces are now non-planar and that they are tangent to the planes in Figure 2 at the origin, as they must be since the system in Example 1 is the linearized version of this one.

A set of simulations, which start in the linear and non-linear modes and are projected back onto the linear and non-linear modes, respectively, is shown in Figure 9 for an initial amplitude of 0• 5. In this case, the non-linear modes again provide an improvement of about one order of magnitude in the accuracy of the normal mode.

In order to demonstrate the power of the approximate non-linear superposition method, a general initial condition which excites both modes is chosen: (xJO, Y10, X2o , Y2o) = (0, 0, 2•0, 0). The simulations are done in three ways. First, the full equations are simulated

as 0 0 a6 0•033 0•155 a1 -0•049 -0•007 Qg -0•011 0•219 a9 -0•019 -0•063 hi 0•460 0•140 b2 -0•989 0•968 b3 0 0 b4 0 0 bs 0 0 b6 0•209 0•053 b, 0•199 0•024 bs 0•082 0•077 b9 0•032 0•245 (b) (a) (c) (d)
Figure 8. The non-linear invariant modal surfaces (approximated to cubic order): (a) X 2 vs. (x 1 , y 1 ) for the first non-linear normal mode; (b) Y 2 vs. (x 1 , y 1 ) for the first non-linear normal mode; (c) X 2 vs. (x 1 , y,) for the second non-linear normal mode; (d) Y 2 vs. (x 1 , y 1 ) for the second non-linear normal mode.

in the usual manner. These results are a very close approximation of the actual solution.

The second set of simulations uses the linear modes in an attempt at superposition. It involves the following steps: projection of the full equations of motion and the initial conditions onto the linear modal co-ordinates via the inverse of the linear transformation; simulation of the resulting individual uncoupled non-linear oscillators which are obtained by ignoring all the coupling terms between the linear modes; recombination using the forward linear transformation; and comparison with the exact results. The third set of simulations repeats these steps using the oscillators obtained from the non-linear normal modes and the associated non-linear co-ordinate transformation, where again all coupling terms between the non-linear modes are ignored. The oscillators and initial conditions are given below.

The results of the simulations are presented in Figure 10, which depicts the four coordinates (x, , y 1 , x 2 , y 2 ) as functions of time. The first trace in each case is labeled "direct" and results from simulations of the original equations of motion. The second and third traces in each case are obtained using the non-linear and linear approximate superposition methods and are labeled "NLNM" and "linear", respectively. All three traces are then plotted together for comparison. The projection of the equations of motion onto the linear modal subspaces is achieved by substituting only the linear parts of the modal expansions for x 2 and Y2 into the first ... ;:.. ----;-'1• . Simulation results for the initial co ndition, (x ,.,, y 1 .,, x 2.,, y 2 u) = (0, 0, 2, 0), used to demonstrate non-linear superposition. In all curves the following line types are used: --, direct simulation of the original system o f equations, labeled "direct" ; ---, simulations o f the oscillators on the no n-linear modal subs paces started with initial conditions obtained from projecting the given initial condition into the non-linear modal subspaces, labeled .. NLNM"; -•-, simulations of the oscillators on the linear modal subspaccs started with initial conditions obtained from projecting the given initial co ndition into the linear modal subspaces, la beled "linea r". (a) Results for x 1 ; (b) results for y 1 ; (c) results for xz; (d) results for y 2 • two equations of motion. The resulting oscillators are : The projection of the equations onto the non-linear normal modes utilizes the cubic approximation of the invariant manifolds and yields the oscillators:

0 -1•0 -0•5 0 •5 Ut 0•01 0.02 (c) l•Or-------,---.....,..---~-----. O•S -1•0'----~-=------~------:--l-::-----:J -1•0 -0•5 0 0•5 l•O "' • 0 - ;:.. (b) 0•02 -O•o:_Ou•-02-----,J.---~0------:0:...L•0:--1---~0•02 "• 1•0 ....----.....----~-------.------, 0•5 -t•O '------1.----'o~---o:'-:----=-"t.o -1•0 -0•5 •5 u2 (d) -0•001 -o~o~.L~------o~.oo~t----~o----~o~.oo~t--~0•002 "• t-o..--------,---------..,.----.------, 0•5 N ;:.. -t•~Lt.-o----o~•..,..5 --~o----:!-o.-=s------:"1•0 "2 "2
mode I: u1 +0•752u1 +2•968ui +0•50ui=0, mode 2 : Uz + 0• 148u 2 + 1•011u 2 + 0• SOu~= 0.
mode 1: ii 1 + 0•752u 1 + 2•968u, + 0•405u~-O•Ollu~u 1 -0•014u 1 u? + 0•009ui = 0, mode 2: ii 2 +0•148u 2 + I•01lu 2 +0•329u~+0•0002u~u 2 -0•242u 2 u~-0•011u~ = 0. (58)
They are significantly different from those•obtained by using the linear projections. Transformation of the initial conditions via the inverse of the non-linear transformation gives (u 10 , v 10 , u 20 , v 20 ) = ( -0• 912, -0• 235, 0•912, 0• 235). Simulations for the oscillators in equation (58) and recombination back into the physical co-ordinates using the non-linear transformation results in the time traces labeled "NLNM" , shown in Figure 10.

We note that the use of non-linear normal modes for the superposition provides an improvement of two orders of magnitude over using the linear modes for this example. It must be pointed out that the improvement will, in general, be very dependent on the problem and on the initial conditions.

For completeness we present the non-linear equations of motion expressed in terms of the linear and the non-linear modal co-ordinates. In terms of the linear modal co-ordinates, the equations of motion are given by

ii1 + 0•752ri 1 + 2•968u 1 + 0•215(u 1 + u 2 ) 3 -0•107(u 1 + u 2 )(u 1 + u 2 ) 2 = 0, u2 + 0•148ri 2 + l•Ollu 2 + 0•264(u 1 + u 2 ) 3 -0•107(u 1 + u 2 )(u 1 + u 2 ) 2 = 0. (59)
In terms of non•linear modal co-ordinates the equation of motion is given by

Ut+ 0• 752u 1 + 2•968u2 + 0•405ui + 0• 724u,u2(uJ + u2) -O•Olluiu1 -0•014ulzii +0•009ui -0•214u 1 u2(riJ + ri2) -0•107(u,u~ + zi2ui) = 0, u2 + 0• 148ri2 + l•Olluz + 0•329u~+ 0•776u,u2(u, + Uz) +0•000lu~u2 -0•242u2u~ -O•Ollu~ +0•214u,u2(ril + u2) + O•lO?(ui u2+ ri1 ~) =0. ( 60
)
Note that these reduce to the modal oscillators on the appropriate subspaces : namely, the first modal oscillator is determined by setting u 2 = v 2 = 0 and the second is obtained by taking u 1 = v 1 :=0. Also, note that equations (60) have special invariance properties not shared by equations (59). This is discussed in more detail in the following section.

CONCLUSIONS AND DIRECTIONS FOR FUTURE WORK

The above examples clearly demonstrate the power of the method for low order systems. Application to larger scale problems is simply a matter of computational power, and it will be possible to go well beyond two degrees of freedom. This is true since the most difficult part of the solution lies in solving for the linear components of the modes (recall that these were coupled quadratic equations), which can be done using standard solution packages for obtaining eigenvalues and eigenvectors. The coefficients required for the oscillator formulation can then be obtained from the standard eigenvectors by using the relationships given in Appendix A. From that point on, only linear problems need to be solved for the coefficients of the non-linear terms in the normal modes, a relatively easy computational task.

An important feature of this work is that it provides quantitative results regarding dynamic behavior, in contrast to much of the current work in dynamical systems theory, which is concerned with qualitative behavior. The non-linear modal oscillators contain useful information about how the system will behave, since they include all the decay rates and natural frequencies, including the manner in which these depend on the amplitude of motion due to non-linear effects. In fact, estimates of the overall long-term transient behavior can be obtained simply by looking at the oscillator(s) with the slowest decay rate.

The method as presented here will not work, or is not convenient, in some circumstances. For example, for some applications it is more natural to use a general first order formulation, in which case there is no special relationship between variables, as there is here between displacements and velocities. In that case, a mixture .of oscillatory and nonoscillatory non-linear modes will occur. Another situation is when an internal resonance exists between two or more modes, in which case the coefficients for the associated nonlinear normal modes will be singular, as in Example 2, indicating that something is amiss. Then one cannot uncouple the resonant modes to first non-linear order, and the problem must be reformulated such that a non-linear modal subspace of dimension 2M, where M is the number of interacting modes, is constructed. In this subspace a set of M coupled oscillators, which cannot be further uncoupled, will govern the dynamics.

This last observation points to a relationship between the present approach and more general invariant manifold theories for dynamical systems [START_REF] Wiggins | Introduction to Applied Nonlinear Dynamical Systems and Chaos[END_REF][START_REF] Guckenheimer | Non/inear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF]. The approximate construction of the invariant manifold is very similar to what is done in center manifold theory [START_REF] Carr | Applications of Centre Manifold Theory[END_REF]. However, we are not constructing center, stable or unstable manifolds, but invariant components in them which are tangent to the linear eigenspaces. These invariant components of the usual center, stable and unstable manifolds span the local phase space near the equilibrium point of interest in a non-linear manner given by the transformations between physical and modal co-ordinates. The generation of the equations of motion on these components is done exactly as it is done in center manifold theory. The existence of these invariant manifolds has been known in the mathematics community for some time, but they have not previously been exploited in the definition and construction of normal modes for non-linear systems. For the mathematically inclined reader, the paper by Fenichel [START_REF] Fenichel | [END_REF] provides existence proofs for these manifolds, including the required conditions on the eigenvalues.

Several extensions, generalizations and applications of the methods described herein are immediately obvious, and the authors are currently engaged in work on many of them. They include: non-linear normal modes and approximate superposition for general, i.e., not necessarily oscillatory, nor even necessarily stable, operating points of dynamical systems; the inclusion of time dependent excitation ; the inclusion of rigid body modes of motion; mode localization in weakly coupled non-linear systems; and the stability and bifurcations of non-linear normal modes [START_REF] Caughey | [END_REF]. Also, a similar method is being developed for continuous systems, by which we can derive non-linear mode shapes for beams, plates, surface waves in fluids, and other distributed parameter systems. This has important implications for the issues of modal convergence and modal truncation for non-linear continuous systems [START_REF] Shaw | Normal modes of vibration for nonlinear continuous systems[END_REF].

Another field in which the method offers substantial benefits is in the area of model reduction. Reduced order models can be obtained from larger scale systems by merely ignoring non-essential modes in the non-linear modal dynamics. The procedure given in this paper yields very "clean" models in that one can, to a given order, systematically eliminate contamination between the modes of interest and the remaining modes, yet retain the modal coupling in the desired model. This can be very important in structural dynamics and in the design and implementation of control systems.

These general ideas are evident in the third example considered in this paper (see section 6.3). In particular, note that by using the equations of motion expressed in terms of the linear modal co-ordinates, equations (59), neither mode can be completely suppressed.

Even if one starts with zero energy in one mode, that is, u, = v, = 0 for i = 1 or 2, the other mode will directly excite it, leading to a fully coupled motion. In contrast, by using the non-linear modal co-ordinates, the equations of motion, equation (60), possess the desired invariance property: if either mode is started with zero energy, it remains at zero energy for all time. These concepts can be generalized to large scale systems in which an arbitrary subset of modes is started with zero energy-they will remain quiescent for all time. This is a direct consequence of invariance.

Finally~ it is important to emphasize that the development of the ideas presented in this paper has been greatly facilitated by the availability of computer-assisted symbolic manipulations. In fact, many of the calculations carried out in the examples are similar to those presented in the chapter on center manifolds in the book by Rand and Armbruster [START_REF] Rand | Perturbation Methods, Bifurcation Theory and Computer Algebra[END_REF].

described by the ith pair of first order equations in the system (equation (A4)). The rate of decay and the pseudofrequency of oscillations are -Rl; = e;w; and lA.;= w;Jl -e~, respectively, where W; and ~;are the mode,s undamped natural frequency and damping ratio, respectively.

A more physical description of the underdamped modal dynamics of the system can be achieved by seeking sets of oscillator equations in the form ( ui) ( 0 1 )(u;) v; --mf -ze,-m; v; ' i= 1, ... , N, (A5) where u; and V; are the "oscillator-form" modal co-ordinates for the ith normal mode. These can be obtained from the first order form modal co-ordinates by the transformation

w(t) = WE(t), W = blockdiag{-• t ~•} • -} (A6)
where w = [u 1 , v 1 , ••• , uN, vNF is the vector of oscillator-form modal co-ordinates (note that for an overdamped mode these modal co-ordinates would not be uniquely defined).

Applying the co-ordinate transformation z(t)=Uw(t), with U=ZW~1 , we transform D into a block diagonal matrix with "oscillator" blocks, which yields the set of pairwise coupled, first order equations (AS). With this oscillator form of the modal dynamics, a mode shape is represented by two (real) columns of U, and a normal mode consists of this mode shape and of the corresponding undamped natural frequency and damping ratio. Note that the oscillator representation features at the outset relationships between both displacements and velocities and thus avoids the use of complex quantities. To obtain the response to initial conditions z(t = 0), these are transformed into w(O) == U~1z(O), which yields a set of N decoupled, second order initial value problems that can be solved for the modal dynamics. The general motion is then obtained as a superposition of N uncoupled oscillator motions: all a2t al2 azz .

(x(l))-r a IN a2N (u.(t) ). The oscillator form of the normal modes can also be obtained directly by formulating a block eigenvalue problem. We simply need to write that U~1DU is block-diagonal. This leads, for the kth normal mode, to which we can regard as a dimension-two block eigenvalue problem, where the "eigenvectors" have two columns and the "eigenvalues" are 2 x 2 matrices with two unknown coefficients. The N eigensolutions are theN oscillator-form normal modes of the system. We note that equation ( 8) in section 2 is simply the block eigenvalue problem, equation (A8), provided that we normalize the kth mode shape such that in equation (A 7) the first displacement is simply taken as uk(t) and the first velocity as vk(t), i.e., a 11 = l , a 21 =0, biJ =0 and b 2 1 = 1. With this normalization, the first and (N+ l)th lines in equation (A8) simply give the eigenvalues of the kth mode as and equations ( 8) and (A8) become identical. Thus, the use of invariant subspace ideas results in the same linear normal modes as those found by transforming the solution of the standard eigenvalue problem to oscillator form and utilizing a special normalization.

APPENDIX B: EQUATIONS FOR THE COEFFICIENTS OF THE NORMAL MODES

In all cases, j = 2, ... , N.

u term, first equation : v term, first equation:

a,kajk-a,ka lkb2i-b11)12Jf3lk + bucf3Jk = 0.

u term, second equation:

v term, second equation:

a2kaJk-b,i-a2ka tkb2i-b2}12kf31k + b2kf3 1 ~c = 0.

u 2 term, first equation:

-a2j13ka lk + b31-a2}13kf3Ik + a4.1{ -a,ka lk -b1~cf3•k) -allc.almtl2.J.Otm~c-al,aqb,kelmk-al}JikbtmYimk = 0.

uv term, first equation:

-2a3J-a2j14katk+ b4i-a2}14k/3lk + asi -2a,kalk-2b,kf3tk) +a4;{ -a2ka lk-b2k/3Jk) + ( -a,m02f12k -alkG2}12nt)8tmk

(Bl)

Figure 1 .

 1 Figure 1. The mechanical model for Example I.

a 2 ,

 2 b 1 and b 2 :

Y

Figure 2 .

 2 Figure 2. The modal invariant planes: (a) X 2 vs. (x~.y 1 ) for the first mode; (b) Y 2 vs. (x.,yd for the first mode; (c) X 2 vs. (x 1 , y 1 ) for the second mode; (d) Y 2 vs. (x 1 , y 1 ) for the second mode.

mode 2 :

 2 -0•0741 + 1•0027i~

• 2 x

 2 w-17 ,....--------r---------. (b) -2 x w-17 L . . . -. -

Figure 3 .

 3 Figure 3. Trajectory in the individual modal phase planes for an initial condition started in the first mode: (a) first mode phase plane; (b) second mode phase plane.

Figure 4 .

 4 Figure 4. The physical model for Example 2.

a 9 ,

 9 b 6 , b 7 , b 8 and b 9 and have two inhomogeneous terms, alK and b2K, in the two u 3 equations. The solution of this linear problem is:

Figure 7 .

 7 Figure 7. The physical model for Example 3.

  ......-----,..-----r----.,.....----, 0•5 -1 ' 0 '-----L....,....---~-----::-10 -=-

Figure 9 .

 9 Figure 9. Simulation validation of the invariance of the non-linear normal modal subspaces; in all cases an initial modal amplitude of 0• 5 was used : (a) trajectory as projected into the individual linear modal phase planes for an initial condition started in the first linear mode; (b) trajectory as projected into the individual linear modal phase planes for an initial condition started in the second linear mode; (c) trajectory as projected into the individual non-linear modal phase manifolds for an initial condition started in the first non-linear mode; (d) trajectory as projected into the individual non-linear modal phase manifolds for an initial condition started in the second non-linear mode.

-t 10

 10 Figure10. Simulation results for the initial co ndition, (x ,.,, y 1 .,, x 2.,, y 2 u) = (0, 0, 2, 0), used to demonstrate

( 57 )

 57 Projection of the initial conditions onto the linear subspaces yields the linear approximation (uw, v10, uzo, v2o)=(-0•973, -0• 159, 0•973, 0•159). Simulations for the oscillators of

  Figure 1 0-continued.

  equation (57) and recombination back into the physical co-ordinates results in the time traces labeled linear in Figure I 0.

  where the kth mode shape consists of the (2k-I )th and (2/c )th columns of V.

  IN a2NalN OzN (-~l -2:.(j)J k= 1, . . . , N,(A8)

  N Nan+ L a,i(alih+ L f3 v(b, 1 )k=-((a2j)k + Pn + I P1 1 (b2i)k = -2ek(J)b k = 1, ... , N,

  a2flonb1k -alma2p2k)Blmk + ( -a2Plmb2k -a2Jblkb2m)Ytmk =0.(B6) v 2 term, first equation:-a4J-a2f1ska lk + bs 1 -a2P5kf3tk +a sA-2a2ka lk-2b2k/3J~c) (B7) u 2 term, second equation:aJk( a Jk-a 1 khzj) + b4i -alka l k -b,kf3Jk) +b3k( -b2j/3lk + f3Jk)-a!kO!mbziO Jmk + alkalmOJmk -almb Jkb2j&Jmk + OJmbik&jmk-blkblmb2JY Jmk + b1kblm YJmk = 0. ( B8)uv term, second equation:a4k( a 1 k-a 1kb 21 ) -2b3 1 + bsi-2alka lk-2bJk/3 lk) +b4j( -a2kalk-bzk/3Jk) + b4k( -bzJPlk+ f3Jk) +( -almazkbz 1 -aikazmbzi )olmk + (a1ma2k+ al~zm)o;mk +( -a2mblkb2J-Olmb2jJ2k)EJmk + (azmblk + O!mb2k)Bjmk + ( -b 1mb2jJ2k-b1kbzp2m) r imk + ( b1mb2k + b,kbZm) YJmk = 0.v 2 term, second equation:

  +< -b,qb2pzJJ2m-btmh2p~2qb,,JJz/12mh2q)g,mkq+ ( b I /J2Jcb2m + b I ,/JlliJ2q + b 1,j}2mb2q) gjmkq = 0.

TABLE 1

 1 Solutions for Example 3

		Model	Mode 2
	a1	-1•106	0•947
	a2	-0•155	-0•139
	a3	0	0
	a4	0	0

  2a1 1 -2a~s 1 a 1 k-a4p5ka lk-azpska Jk + bsJ-2as.P4kf31k -a4Jbskf3lk-az}Jsk/3Ik + agj(-3alka lk-3blk/3Ik) + asi-2a2ka lk-2b2kf31k) si + 2a s.flska1 k + Oz.f19ka 1 k -b91 + 2as}Jskf3 Ik + azpgkp 1 k + a9j( 3a2ka 1 k + 3bzkf3 1 k) +(2a2ka2ma5J + Oz.f12mask + Oz.f12k05m)8Jmk + (2azmOsJb2k + Oz.flsmb2k + Ozp2mbsk)6Jmk +(2asp2kb2m + 02;h2mhsk + OzJb2kbsm)Yimk + Oz.fl2ka2mazq/.llmkq + Ozp2ka2mb2q VJmkq +az.fl2mb2kbzqPimkq + azp2~2mb2q~ Imkq = 0. -at,a2khtrP2i-OJI(fl2mbl/J2J-a1~1mb2p2q) V1mkq +(a,,a2kblq + alkflanb,q + a,lcfl,,b2q) VJmkq +( -aanb")J,qb2 1 -a,mhtrP2P2k-a,mb,kb2p2q)Pimkq +(a2mb1kblq + a,mblqb2k + a,,b,kb2q)PJmkq+< -b.mb,qb2p2k-b!J)11/J2P2m-b,kb1mh2P2q)~,mkq +(b,,JJ,qh2k + b,~,;2m + htkbl,JJ2q)~jmkq = 0.

	+( -a2ka2ma4 1 -Gz;G2mll4k-a2p2ka4m-2almazkas 1 -2alkazmGSJ-almazjJsk-a,/(llz_flsm)81mk
	+( -2azmas./Jik-Ozflsmblk-Gzm04;h2k-az.f14mb2k-2almas}b2k -azflzmb4k-almOzPsk)s,mk
	+(-2as}lmb2k-2as}lkb2m -a4)Jzkb2m-Ozjb2mb4k-Gzjb2kb4m-azplmbsk -OzPtkbsm) y lmk
	+(aJqll2.fl2k02m-a,mazp2k0zq-alkflzjl~zq)P, Imkq
	+( -azjl2ka2mblq-0Jm0z.flzkb2q-GlkfJZ;<I2mhzq) VJmkq
	+( -Ozjl2mbiqb2k-Ozflzmblkb2q-a,,LlJh.kbzq)Pimkq
	u 2 v term, first equation: -3a6J-a4p4ka lk .... 2a3kasJa lk-Ozf17ka lk + b1 1 -2asJ;Jk/3Jk -aJ;4k/3lk-az}hk/3Jk +a si-2alka 1k-2b Jk/3Jk) + a1j( -azka Ik-b2kf3Ik) +( -azp2tna3k-Ozf1zJcll3m-a,,a2ka4J-a,kazm04j-alma2_fl4k-OJk02.f14m-2aJkalmasJ )olmk (B9) (BIO) (Bll) +( -a2ma4J;Jk-Ozf14mblk-2allnflsJblk-az 1 a3mh2k-alm04/Jzk-az_,a2mb3k-a,mOzJb4k) 81mk +( -2as}!khlm -a4/Jlmb2k-aJ;,thzm -az./JzmbJk-az)J2kb3m-azPJmb4k-a2;b1kb4m)y,mk +( -az)Jqb2kb2m-0zPimhzkh2q-Oz}Ikbzmhzq) e I mkq = 0. (B13) v 3 term, first equation: (B 14) u 3 term, second equation: a6k(aJk-a1khz 1 )-a3kalkb4J-b3kb4 1 f3,k + b1j( -alkalk-b,kf3Ik) +b6k( -bzJf3Ik + f3Jk) + ( -alma3kb2J-aikf13mb2J-a1kf11mb4J )8tmk +(alma3k +aJkGJm)8Jmk +( -a3mbl~2j-G!mhz;b'Jk-GJmblkb4J) Glmk + (a3mb!k + OtmbJk) SJmk +( -bimbz}hk-blkbz}3m-b1kblmb4j )Ylmk + (blmb3k + b,kb3m) YJmk -a!/(lllmG!/Jz;Jlimkq + 0Ik0JmllJqJi}mkq-0JI,almblqb2jVImkq +alkillmb!q Vjmkq-almblkblqb2JP1mkq + 0)mbJkbiqPJmkq -bikbinPI/Jlj~lmkq + blkbimb!q~jmkq = 0. u 2 v term, second equation: a7k( a Jk -a Jkbz 1 ) -a4ka Ikb4 1 -2a3ka ikbs 1 -3b6 1 -b4p4k/3Jk-2b3kbsJ/31k + bsi -2alkalk-2b,k/3Jk) +h1i -azka Jk-b2kf31k) + b7k( -b2Jf3lk + f3Jk) (Bl5) +( -azma3kb2}-a2ka3mb2J-GJm04kb2J-atkG4mb2J-GJmG2kb4J-GJkGzmb4J-2allfllmbsj )0 lmk a8k(a 1 k-a1JJ2 1 )-aska,kb 41 -2a4katkhs 1 -2b?J -2h41.hsJ13,k-b4}Js~cf3tk + b9i-3a,ka Jk-3b,k/3Jk) +bsj( -2a2kalk-2b2k/3Jk) + b8k( -b2Jf3lk + f3Jk) + (-02Jnll4kh2J-U2JcU4mb21-a 1 ,a s~21 -a 1 ~ smb21 -a21c02mb4J-2a,,a2kbSJ-2aufl2mhs1 )0 Jmk +(a:zma4k + U2JcU4m + Ot,tlsk + 0Jkflsm)8Jmk +( -asmbu}J2J-a4nihP2k-a~~4J-a2mb2jJ41c -2a2mblkb5J-2a,,JJ2kbSJ-a,mb2Ps~c) Etmk +(asmblk + 04mb2k + a~4k + a,mbs~c)EJmk +( -b2kbu,b4j-b2panh4k-b2puJJ4m-2b,mb2kbsj -2b1JJ2mhsj-b,mb2./JSk-b,kb2jbsm)Ylmk +(b2mb41c + b21cb4m + ht,bsk + blkbSm)YJmk +( -a,qll2JcUanb2J-a,,a21cfll/JzJ-a11,a2nllz/J2J )JlJmkq +(a,qtZ2JcU2m + a1,a21cf12q + Otlc02ma2q)J1 1 m~cq +( -a2Jcll-unb.qh2 1 -a.,a21ch2}2q-a,kll-unb2.P2q) V1mJcq +(a2.k02mhlq + a,,a21cb2q + a,Jc02mb2q) VJmkq +(a2fn03k + uv 2 tenn, second equation: +( -aanbaqh2./J2k-a2mb1kb2p2q-a,,JJ2 1 b2kb2q)p,mkq

ask( aik-a Jkb21) -b4 1 + bsi-2a2ka Jk-2bzkf3Ik) +bsk

( -b21P

lk + PJk) .... 02kG2mhz;D1mk +a].k0.2m0jmk .... 02mbzj;2k&!mk + 02mb2kEjmk -bzp2kb2m r Jmk + b2kb2m YJmk = 0. u 3 term, first equation: 0JJcll4ja lk + azjl6ka tk-b6J+ a4Jb3kf3Jk + OzJh6k/31k +a1;(a1ka tk + blk/3Jk) + (aJmOzflJk

+ 0JJcll2fl3m + alkalma4i )o lmk +(Ozfl3mblk + 01maJ1Ik + OJmOzjb3k)Eimk +(aJY,kblm + a2j;!,b3k + Ozjh!kb3m)Ylmk + alka!ma!q02Jillmkq +al~im02Plq V]mkq + 0)m0z;hlkblqPJmkq + azjb!kblmh!q/glmkq = 0.

+( -aJm0lq02f12k -0Jk0iq02fl2m .... a,kOtm02flzq)JL lmkq +( -alml12.flzkhtq-alk02.flzmhl q-a Jka Jmavbzq) V Jmkq +( -a2_fl2mblkblq-a,ma2jJJqh2k-OJmazPikbzq)Plmkq +( -az.PJmblqb2k-a2p 1kblqb2m-a2jb!kblmb2q); lmkq = 0. (Bl2) uv 2 term, first equation:

-

a a2ka3m + 0Jm04k + a1ka4m)OJmk +( -a4mb!kb2J-GJmb2Jb2k-02mb21b3k-02mblkb4J-GJmb2kb4j-almb2p4k-2almbl~SJ) £Jmk +(a4,/Jik + a3,JJ2k + a2mbJk + 0Jmb4k)EJmk +( -b2.P2mbJk -hzP2kb3m-h1

,JJ2kb 41 -b1kb~4 1 -h1mbz 1 b4k-bikbz.P4m-Zb1kb ,mbs 1 ) r lmk +(bu,b3k + b2kb3m + htmb4k + b,kb4m)Yjmk +( -a,,alq02kb2J-a1kflJqOanb2J-a1kf11ma2/J2J )Jlimkq +(a,mOJq02k + a1~1q02m + alkf11,112q)/J}mkq +( +(aanbtqh2k + aanh,kblq + OtmhulJ2q)pjmkq

ACKNOWLEDGMENTS

The first author is grateful to the Mechanical Engineering and Applied Mechanics Department of the University of Michigan for a Westinghouse Distinguished Lectureship, during which this work was initiated. This work was supported in part by a grant from the National Science Foundation (MSS-9201815). The development of the ideas contained in the paper was aided by discussions with several people including Joe Whitesell, Alan Haddow, Phil Holmes, Galip Ulsoy and Stephen Wiggins. We are especially grateful to Richard Rand for providing us with an historical perspective on non-linear normal modes, and explaining to us how our method relates to the traditional one used for conservative systems.

Because D is real, those eigensolutions that are complex must occur in conjugate pairs and, in general, there are pairs of complex conjugate eigensolutions and otherwise real eigensolutions. Here, we are interested in oscillatory systems with underdamped modes, and for simplicity we assume all eigenvalues of D are complex (our development could easily be extended to handle critically and overdamped modes). Also for simplicity, we assume no repeated eigenvalues, as otherwise complications can arise [START_REF] Hirsch | Differential Equations, Dynamical Systems and Linear Algebra[END_REF].

TheN pairs of complex conjugate eigensolutions are (A. 1 , z 1 ) and (A.f, zi), i= 1, ... , N, where the asterisk denotes a conjugate. Defining a modal matrix Z, the columns of which are the eigenvectors, and changing co-ordinates according to z(t) = Zt(t), yields a set of 2N uncoupled first order equations:

For oscillatory systems a real co-ordinate transformation is often preferred that renders D only block diagonal. This is achieved by introducing the matrix formed by using the real and imaginary parts of the eigenvectors as columns:

where R and I denote real and imaginary parts, respectively, of the eigenvectors z 1 . Letting z(t) == Z'tt(t), equation ( 3) is transformed into first order equations that are pairwise coupled:

These pairs of equations correspond to underdamped modes associated with the complex conjugate pairs of eigenvalues, such that the corresponding dynamics occurs in a twodimensional subspace obtained by the span of the real and imaginary parts of the complex conjugate eigenvectors. (For overdamped modes the modal dynamics occur in one-dimensional spaces.) A motion of the ith underdamped mode consists of decaying oscillations v 3 term, second equation: a9k( a Jk-a 1kb2 1 ) -2aska lkbs 1 -bs 1 -2bsflsk/3 tk + b9 1 ( -3a2ka lk -3b2k/3Jk) +b9k( -b21fJ 1k+ fJJk) + ( -aQnllskb2J-a~smb2J -2a~an/JsJ )Dtmk +(a2mask + 02Jc(Jsm)8Jmk + ( -asmb2./J2k -2aanh~sJ -a2mh2Psk)Elm1c

APPENDIX C: THE INVERSE OF THE MODAL TRANSFORMATION (Bl8)

Index notation and the implicit summation notation are employed for deriving the inverse of the modal transformation. This is helpful since we will encounter a fourth order tensor along the way. The inverse is derived using series expansions.

The forward transformation is written as (Cl) and its inverse is expressed as

The method in section 4 provides the coefficients in transformation (Cl), and the task here is to derive the coefficients for the inverse (C2) in terms of them. Substitution of equation (Cl) into (C2) and expanding out to order three in the w/ s yields Wm ::;; QlmnMinj Wj + QlmnMlnjrWjWr + QlmnMJnjraWjWrWa

Like powers of w are now gathered together on the right side and are equated to the left side. The linear terms yield (C4) the expected result, which says that the linearized inverse must be the inverse of the linarized transformation. The quadratic terms must satisfy (C5) which can be solved for Q 2 by utilizing the transpose of (C4) and post-multiplying twice by Q1 to eliminate the M1 terms on the left side of (C5). The result is (C6)

The cubic tenns in w must satisfy

which can be solved by three post-multiplications of Qa to yield

This entire process will involve the inverse of one matrix and extensive use of that inverse in the higher order tenns. Note that the equations can be solved sequentially, and that after the initial inverse is taken, only tensor multiplications, and no further inverses, are required.