A NEW GALERKIN-BASED APPROACH FOR ACCURATE NON-LINEAR NORMAL MODES THROUGH INVARIANT MANIFOLDS

A method for producing accurate reduced order models of non-linear vibratory systems is presented based on the invariant manifold description of non-linear normal modes (NNM). This approach makes use of polar co-ordinates to obtain equations which govern the geometry of the invariant manifold. These equations are discretized through a series expansion and Galerkin projection over a chosen amplitude and phase domain, yielding non-linear equations in the expansion coefficients. These equations, when solved numerically, yield an invariant manifold which is accurate to the degree of the expansion, and devoid of the limitations which plague typical asymptotic solutions. Such Galerkin-based solutions may be used to generate accurate reduced-order models for large-amplitude, strongly non-linear motions. This procedure is illustrated using two non-linear examples, a two degree-of-freedom oscillator, and a finite element beam model. The solution convergence and manifold geometry are discussed and the resultant reduced-order models are shown to possess exceptional accuracy over large amplitude ranges. This approach allows the full potential of the invariant manifold formulation to be reached, and is suitably general for application to a wide variety of non-linear systems.


INTRODUCTION
The generation of e!ective reduced-order models for multi-degree-of-freedom non-linear dynamic systems is di$cult, due to the complex interactions between system components. These non-linear interactions typically couple the system dynamics such that accurate results may only be achieved through the inclusion of many modes or degrees of freedom (d.o.f.s). Such large models are cumbersome, making analysis both di$cult and slow, as well as physically non-intuitive. Alternatively, smaller models often sacri"ce accuracy, ultimately yielding questionable results. Ideally a minimal model is sought, which accurately accounts for the non-linear interactions between components, without requiring their explicit simulation. In pursuit of this goal, there has been considerable work on the development of reduced-order models (ROMs) of non-linear dynamic systems. Figure 1. Cross-sections of the exact, asymptotic and Galerkin manifolds at "0 for the second non-linear mode of the "nite element system illustrated in Figure 12: **, asymptotic manifold: } } } , Galerkin manifold: #, exact. The asymptotic manifold is of third order, while the Galerkin manifold uses 80 piecewise linear segments in the amplitude, a, and 8 harmonics in the phase, .
Similar work on the reduction of dynamic systems using Galerkin methods has been conducted by Steindl et al. [14]. Another reduction technique that employs invariant manifolds, but is based on statistical considerations, is the Karhunen}Loeve decomposition. It provides a set of linear bases for the reduced-order model that optimally capture the energy content of a particular system response [19,20].
The present application of the non-linear Galerkin method, as well as the co-ordinate transformation of the manifold-governing partial di!erential equations (PDEs) developed herein, constitute unique contributions to the "eld. Their utility is illustrated through the examination of two example systems. First, a two-degree-of-freedom system with cubic non-linearities is analyzed. The convergence and geometry of the two (numerically obtained) approximate invariant manifolds are discussed, and the Galerkin-based reduced model is shown to accurately predict the system response. These results are compared against those from the asymptotic approach. Second, a "nite element based beam model with a quadratic/cubic non-linear torsional spring at one end is examined. Results which correspond to the second and "fth non-linear normal modes of the beam are presented.
As with the simpler system, these results indicate that the Galerkin-based reduced-order model retains accuracy throughout the chosen amplitude domain, without requiring the non-linear e!ects to remain weak. This accuracy enables one to use the reduced-order model with con"dence in its validity, even for amplitudes corresponding to large non-linear e!ects.
These examples serve to illustrate the potential of this approach. The formulation may be extended to determine manifolds of higher dimension, yielding reduced order, multi-d.o.f., models for motions occurring within multiple, interacting modes. In addition, considerable insight may be obtained through the application of this process to existing problems within the area of non-linear dynamics. Ultimately, this approach promises to be a valuable new tool for the analysis and understanding of non-linear dynamic systems.

FORMULATION
This approach assumes that a reduced-order model is sought for a set of non-linear coupled oscilators. Whether from a system of discrete masses, "nite element model, or discretized continuous system, these may be expressed in the form where is a vector of normalized modal co-ordinates, [Z] and [ ] are diagonal damping and sti!ness matrices, respectively, f( , ) contains any additional linear or non-linear terms, and each overdot indicates a time derivative. The individual elements of [Z] and [ ] are written as 2 G G , and G , respectively, for the ith equation. This form is selected for analytical convenience, and does not re#ect a necessary restriction of the approach.
From this point, previous single-mode invariant manifold formulations [7] have constrained all degrees of freedom to be functions of a chosen co-ordinate pair, as follows, where the subscript k indicates the &&master'' mode of interest: for N original modal co-ordinates. Here, an alternative set of co-ordinates is de"ned using the transformation and, the original constraints, equation (2) become The resulting system uses this hybrid co-ordinate system. That is, all modal positions and velocities are expressed in terms of a &&master'' amplitude and phase corresponding to the motions of the kth mode. These constraint relations are dictated for the kth mode (by equation (3)), and must be determined for all others. Determining these unknown relations, that is, the P G and Q G , will ultimately yield reduced dynamic equations in terms of a(t) and (t) only. As an alternative to using modal positions and velocities for the dependent modes, the entire system may be transformed into an amplitude}phase representation. However, such an approach is unnecessarily complicated.
Unlike earlier methods, this approach is not conducive to solutions through polynomial expansions. However, it does o!er the following advantages: (1) The dynamics in the transformed co-ordinates are elementary for f"0.
(2) The new constraint equations, (P G , Q G ) must be periodic in , and thus may be expressed with harmonic functions in .
(3) The boundary values for the manifold are reduced from four values ($u I , $v I ) to one*the upper bound of a.
The original equation of motion governing I may be transformed into two "rst order equations in a and [15]: #2 I I sin cos (5) and the remaining di!erential equations may be written in "rst order form as The time derivatives on the left-hand side of equation (6) may be expanded using the chain rule to obtain Equations (5}7) are combined to yield a set of partial di!erential equations (PDEs) that are independent of time and govern the invariant manifold geometry: In previous works, a similar process is used to yield manifold-governing equations in terms of a master position and velocity, ( I is then used to asymptotically approximate the local solution to a given order. Herein, the asymptotic approach is replaced by a Galerkin method. As before, the solution is achieved through series expansions with unknown coe$cients, but here the functions of amplitude and phase (a and ), are de"ned over a preselected amplitude domain. The solution of the resulting discretized equations minimizes the error over the chosen domain, yielding manifolds of considerably greater accuracy and range.
The unknown position and velocity constraint relations are expanded as a double series in the amplitude and phase as where the C's and D1s are the unknown expansion coe$cients and the ¹ JK and ; JK are known shape functions, typically composed of products of functions in a and . For example, a given ¹ JK might be a product of a harmonic function in and a polynomial function in a, de"ned over the domain a 3 [0, a ], 3 [0, 2 ]. N ? and N ( are used to denote the number of expansion functions used in a and respectively. This expansion is substituted into the manifold-governing equations (8) and (9), which are multiplied by a to remove the singularity at a"0. Finally, all terms are moved to the right-hand side of the equation, and a Galerkin projection [15] is carried out using the individual shape functions over the chosen domain. This leaves 0" 0" ?( for i"12N, iOk, p"12N ? , and q"12N ( . This represents a set of 2(N!1) N ? N ( non-linear equations in the C's and D's, the solution of which will be optimal (in a least-squares sense) for the chosen basis functions over the given domain.
In principle, this integration could be completed analytically, leaving a set of explicit non-linear equations. However, as the system size and number of expansion functions increase, this quickly becomes impractical. The alternative is to use numerical integration for each function evaluation. Although it is considerably slower than an explicit formulation, this approach requires very little additional analytic work, and the generality of f makes it applicable to a wide variety of non-linear systems. Regardless of the evaluation method, equations (11) and (12) may be solved using an algorithm for the solution of multi-variable non-linear systems of equations. The solution yields the series coe$cients for each P G (a, ), Q G (a, ) pair over the chosen domain, e.g., 3 [0, 2 ], a 3 [0, a ]. Unlike the asymptotic approach, here the expansion order may be increased to produce more accurate results without requiring additional analytical work. Hence, manifold accuracy is reduced to an issue of computational e!ort. Also, the domain of convergence is known a priori, and the results obtained may be used with con"dence throughout that region, once the desired accuracy is achieved.
Once the P G and Q G have reached the desired accuracy, f I , in equation (5), may be evaluated for any given (a, ), and the dynamics within the chosen mode may be determined by numerical simulations. The resultant time histories, a(t) and (t), allow the evaluation of all linear modal positions and velocities (through the P's and Q's*equation (10)), and the global system response may be assembled.
It should be noted that there are other solution methods which may be applied to obtain the coe$cients. For example, an incremental approach which separates each C and D as Figure 2. Two-degree-of-freedom non-linear system with hardening cubic springs.
could be linearized locally (in the C's and D's) and used iteratively until a solution was reached. However, as above, the integration for each step must be done either analytically (using a general formulation) or numerically, and the resulting solution will apply over the chosen domain.
Thus, the non-linear Galerkin method described above allows one to generate accurate reduced-order models for many types of non-linear structures, without restrictions on the form or magnitude of the non-linearity, and these models are expected to maintain "delity throughout a predetermined amplitude range.
In the next two sections, this formulation is applied to generate reduced-order models of two example systems. The "rst, a two-mass system with purely cubic non-linearities, serves to solidify the concepts developed above and uses a single function set to represent the amplitude dependence of the manifold. The second system is intended to showcase the power and bene"ts associated with the non-linear Galerkin approach. A non-linear "nite element beam model is used to generate a set of equations which are then reduced to a single-mode ROM Due to the model size, this reduction is carried out by solving the manifold equations over several local domains. This composite manifold is then used for the ensuing analysis. These examples serve to demonstrate the accuracy and versatility of the general approach.

APPLICATION: A NON-LINEAR TWO-MASS SYSTEM
As an illustrative example, the non-linear undamped two-mass system shown in Figure 2 is considered. The corresponding equations of motion in q and q are A linear eigenanalysis is carried out, and the natural frequencies are found to be "0)689, "3)244, and a transformation matrix of the form may be used to decouple the two equations to linear order. The resulting modal equations are where the non-linear terms may be written as The equations are now of the form assumed in equation (1) and, given a set of expansion functions, equations (11) and (12) may be employed to yield the system of non-linear equations whose solution will produce the manifold coe$cients. Many di!erent sets of basis functions may be used to describe the manifold geometry. Here we use the following double expansion in a and : The harmonic functions are a natural choice for expanding the dependence over [0, 2 ], and the use of cosine functions for ¹ (corresponding to a modal position), and sine functions for ; (corresponding to a modal velocity) is predicated upon the synchronous motions that are expected for this conservative, non-gyroscopic system. That is, in a synchronous periodic motion, all degrees of freedom must simultaneously reach zero displacement, and also simultaneously achieve zero velocity (at a corresponding maximum or minimum displacement). In more general systems this property may not be present, and a more general expansion in may be necessary. The functions¸J(a) were chosen to be a set of polynomials de"ned over the domain [0, a ], with zero slope and a"0, and satisfying the orthogonality property Analytical expressions for the "rst seven¸G(a) are given in Appendix A, and are illustrated in Figure 3. The restriction of the slope at a"0 is a re#ection of the modal form of the original equations of motion. That is, the uncoupling of the governing equations at linear order precludes linear modal contributions to the manifold at a"0. If the equations of motion were in another co-ordinate system, e.g., physical co-ordinates or non-orthogonal component modes, the¸G(a) must be chosen to allow for linear contributions to the invariant manifold at a"0. ; -----, (2, 4); } } } , (3,8); )))))), (4, 10); ------, (5,12); #, exact. In plot (b): **, asymptotic; } ) } ) } ) , (1, 2); } } } , (3,8); )))))) , (4, 10); ------, (5,12). The parenthetical notation refers to the number of basis functions, (N ? , N ( ), used in the expansion. Plots (a) and (b) illustrate cross-sections P (a, 0), and P (2)22, ) respectively. The orthogonality property above, along with the orthogonal properties of the harmonic functions, and the absence of velocity-dependent forces, allows equation (11) to be reduced to where C and D represent vectors containing all the CJK G , and D NO G respectively. Hence, given a guess of the C's, the D's may be determined, and equation (12) may be evaluated. As a consequence, only the C's need to be considered as independent variables.
The resulting set of non-linear equations are polynomials in the unknowns. This follows since the system non-linearity is polynomial, the manifold equations involve products of functions, and the function expansions are linear in the unknown coe$cients. These non-linear equations can be solved computationally. Here, the Galerkin projection was carried out numerically at each iteration, and Powell's Hybrid method [17] was implemented via the Numerical Algorithms Group (NAG) routines to achieve a solution. This numerical algorithm uses successive function evaluations to approximate the Jacobian and progress toward a solution.

THE FIRST NON-LINEAR NORMAL MODE
The "rst non-linear mode of this system may be examined by solving equations (11) and (12) using the above expansions, for k"1, to determine P (a, ) and Q (a, ). Two cross-sections of P , at "0 and a"a "2)22, are illustrated in Figures 4(a) and 4(b), respectively, for a number of di!erent expansion orders, depicted in the form (N ? , N ( ). Also shown are the asymptotic solution and an &&exact'' reference solution (described below). The asymptotic solution corresponds to a third order manifold in ( , R ), generated according to the analytical formulas presented in reference [9], and then transformed into the (a, ) co-ordinates. The conservative nature of this problem requires that the individual non-linear normal modes be periodic, and this can be exploited to obtain the exact solution by searching the con"guration space for initial conditions which yield a periodic response, and consequently must lie on the invariant manifold. This search is carried out within the plane depicted in Figure 4(a) and the exact results are shown. However, these results are not shown in Figure 4(b), as the amplitude does not remain constant through a given motion, and the exact solution is not easily transferred to the plane of constant amplitude. The results indicate that three polynomial functions and eight harmonics are necessary to out-perform the asymptotic results, and that even with "ve polynomials and 12 harmonics, there is still some small error near a . This error appears to be con"ned to a few regions of , as Figure 4(b) indicates convergence over most of the domain. It should be noted that, due to the symmetric nature of the manifold, the expansion yields only contributions from the odd harmonics. Note that further increases in the number of harmonics or polynomials would result in improved convergence, albeit at higher computational cost.
The invariant manifolds for P (a, ) and Q (a, ) are shown in Figure 5, for the solution corresponding to N ? "5, and N ( "12. These surfaces represent an approximate solution to equations (8) and (9), and allow the accurate extension of the "rst linear mode into the non-linear realm. Given these constraint relations, f ( , ) in equation (15) becomes f (a, ), and the amplitude}phase di!erential equations*equation (5)*may be integrated to "nd the correct single-mode motion. These manifolds go well beyond the domain of &&small motions'' and &&"rst order e!ects'', allowing one to generate results which may be limited by the "delity of the mathematical model, rather than by the reduction technique. Here, at a"2)22 and "0 the calculated non-linear forces, f (a, ), is slightly greater in magnitude than the linear force, a, corresponding to a large non-linear e!ect.
The degree of non-linearity is easily observed in Figure 6, wherein the response frequency can be seen to change signi"cantly with amplitude. Although not initially apparent, the known manifold does not allow this curve to be extended to an amplitude of a ("2)22), as the peak amplitude for a given motion is not at "0, but at " /2. Motions with initial amplitudes greater than 1)5 will eventually exceed a "2)22, thereby leaving the known invariant manifold. The plot indicates that both the asymptotic and Galerkin-based manifolds accurately capture the response frequency. The traditional approach, which simply uses a projection onto the "rst mode and ignores the e!ects of the second mode (shown as the &&One Mode Solution''), leads to considerable error at the large amplitudes. The asymptotic and Galerkin-based manifolds yield nearly equivalent response frequencies through an initial amplitude of 1)0, after which some divergence of the asymptotic result can be seen. Time histories for four periodic solutions are shown in Figure 7, where each trajectory is produced by a di!erent procedure. The &&Galerkin Manifold'' results are produced by simulating equation (5), using the (N ? "5, N ( "12) manifold to evaluate f , then using the manifold equations to compute the modal displacements for each (a, ), and "nally using the eigenvectors to compute the physical displacement. The &&Asymptotic Manifold'' uses the analytical solutions from reference [9] to determine "g( , R ), then simulates equation (14) (retaining terms of the proper order), determines the modal de#ections, and, "nally, the physical de#ections. The &&One Mode Solution'' simply assumes (t)"0, simulates equation (14), and uses only the "rst mode shape to determine the physical motion. Lastly, the &&Exact'' solution corresponds to the numerically determined initial conditions which result in periodic behavior of both degrees of freedom. Hence, although each simulation has an initial condition corresponding to "1)5 and R "0, the di!erent constraints on the second linear mode yield di!erent initial conditions. The results indicate that, at this amplitude, the Galerkin-based results are indistinguishable from the &&Exact'', whereas the asymptotic manifold results are nearly correct, and the &&One Mode Solution'' has considerable errors in both amplitude and frequency. These results clearly show that the non-linear Galerkin approach has produced a model which utilizes the full potential of the invariant manifold to reduce the system to a single oscillator whose motions would not incite contamination from the rest of the system. Of course, this "rst NNM is simply an extension of the "rst linear mode into the non-linear realm, and a corresponding extension exists for the system's second linear mode as well.

THE SECOND NON-LINEAR NORMAL MODE
The convergence of the manifold corresponding to the second non-linear normal mode is illustrated in Figure 8, for the cross-sections "0 (a), and a"a "3)0 (b). For this mode the asymptotic manifold diverges quite quickly from the exact solution, making even the N ? "1, N ( "2 solution a considerable improvement. Plot (b) illustrates that the harmonic content of the manifold is simpler than that of the "rst normal mode, and this observation is } } } , (2, 4); -----, (3, 6); )))))) , (4, 6); #, exact. In plot (b): **, asymptotic; } ) } ) } ) , (1, 2); } } } , (2, 4); -----, (3, 6); )))))) , (4, 6). The parenthetical notation refers to the number of basis functions, (N ? , N ( ), used in the expansion. Plots (a) and (b) illustrate cross-sections P (a, 0), and P (3)0, ) respectively. consistent with the fact that only half as many harmonics are necessary to reach convergence. As with the "rst mode, convergence appears to be the slowest near a . However, as was discussed above, the motions which begin at "0, and near a , will exceed a during the course of their motion (this is characteristic of systems which sti!en with increased amplitude). Consequently, convergence near the points (a, )"(a , 0), (a , ) and (a , 2 ) is not necessary for practical use of the manifold. Unfortunately, the question of &&how near?'', must be evaluated for any given problem.
The invariant manifolds corresponding to P (a, ) and Q (a, ) are shown in plots (a) and (b) of Figure 9, for an expansion with N ? "4 polynomials, and N ( "6 harmonics. As before the expansions which describe these surfaces may be used to eliminate from the system, leaving equations of motion in only the amplitude and phase of the second linear mode. As was apparent from the cross-sections shown in Figure 8, the qualitative shape of these surfaces is considerably di!erent from those in Figure 5. The presence of fewer (and lower) harmonics on the surface indicates a closer tie between the phases of the two linear modes. That is, for motions occurring on this invariant manifold, the two linear modes are closer to moving in unison than for motions on the "rst NNM invariant manifold. As with the "rst NNM, simulations may be carried out using this invariant manifold, as well as the other reduction techniques, and the amplitude}frequency relation may be determined. For the second mode, as Figure 10 illustrates, all approaches yield nearly the correct response frequency, indicating that the participation of the "rst mode has little e!ect. However, an accurate response frequency does not guarantee and accurate description of the system motion. Figure 11 illustrates that, though the frequencies are all close, both the &&One Mode Solution'', and the &&Asymptotic Manifold'' contain errors in the response magnitude stemming from their treatment of the "rst linear mode. Figure 8 shows that the asymptotic solution over-predicts the "rst linear mode participation, while the &&One Mode Solution'' assumes no participation of the "rst linear mode. Once again,  Figure 11. Response of the "rst mass, q , for periodic motions in the second NNM, using various reduced systems: **, Galerkin manifold; } } } , asymptotic manifold; ------, one-mode solution: )))))), exact. The exact and Galerkin manifold (N ?
the &&Exact'' solution is indistinguishable from the non-linear Galerkin results, highlighting the accuracy of the approach.

APPLICATION: NON-LINEAR FINITE ELEMENT BEAM
In order to illustrate the general applicability of this approach for the generation of reduced-order models, a more realistic structural model is examined. The "nite element code PATRAN was used to create a linear model of a beam with transverse de#ection u(x, t) using 200 two-noded beam "nite elements. The beam, shown in Figure 12, is pinned at one end, while the other is constrained by a linear spring. A non-linear torsional spring with quadratic and cubic components is located at the pinned end. Given the characteristics of  u(x, t), length¸"1 m, density "7860 kg/m, Young's modulus E"2;10 N/m, cross-sectional second moment of area I"5;10\ m, spring sti!ness k"10 N/m, and non-linear torsional sti!ness R "5000u(0, t)#20 000u(0, t) N (where u indicates the partial derivative of u with respect to x). the non-linear spring and the eigenvectors of the linearized system, one may determine the coupled non-linear forces for each mode due to the torsional spring. For illustration purposes, the lowest 10 (of 400 possible) linear modes are chosen to represent the &&complete'' model, from which an NNM-based reduced-order model will be generated. Note that this 10-d.o.f. model contains non-linear coupling between all linear modes through the non-linear spring.
The increase in example system size from two modes to 10 modes requires a shift in tactics. For the previous example, each non-linear mode was associated with a single (P, Q) pair*representing the contributions of the other linear mode.
For the "nite element system, each non-linear mode is associated with nine (P, Q) pairs. Consequently, a manifold with expansions of order N ? "5, and N ( "12, now requires a solution of 540 coupled non-linear equations for 540 coe$cients. The computational e!ort necessary for the solution scales approximately with (N ? N ( (N!1)), indicating that, for a shift from two to 10 d.o.f., the solution time increases by approximately a factor of 81.
In order to achieve a more e$cient solution, the domain in a may be discretized into several smaller intervals. These intervals may then be described as linear in a ("xing N ? at two*one coe$cient for each segment endpoint) and solved as individual sub-problems. The harmonic dependence is retained, so that each expansion now contains 2N ( coe$cients, and corresponds to a thin annular strip of the manifold (a 3 [a G , a G> ] while 3 [0, 2 ]). As with the polynomial approach, the velocity coe$cients (the D's) may be condensed out, leaving only the C's as variables. This is only slightly more complicated without the orthogonal properties of the polynomial functions used in the previous example.
Besides the reduction of problem size, there are several other bene"ts to this piecewise amplitude approach. The linear dependence on a simpli"es the evaluation of each sub-problem by eliminating the higher order polynomials, and the smaller domain allows fewer points to be used in the numerical integration. In addition, discrete linear segments may capture details or features which would require an impractical number of polynomials. Finally, in the event that numerical di$culties are encountered (e.g. from a bifurcation in the manifold, or a non-removable modal interaction), the anomaly is localized to the sub-problem level, whereas the global solutions used for the "rst example simply may not converge at all in such circumstances.
This piecewise strategy was used to examine the second and "fth NNMs of transverse vibration for the "nite element beam system shown in Figure 12 A sample pair of cross-sections at "0 of the second-mode invariant manifold are displayed in Figure 1; these correspond to the components of the constrained displacements of the "rst and third linear modes. As with the "rst example, the asymptotic manifolds used here are calculated to third order according to the analytical solutions presented in reference [9]. Here, the de"ciency of the asymptotic approach is quite apparent. Although the initial behavior is captured by both methods, the polynomial basis of the asymptotic approach inevitably leads to a rapid departure from the correct solution as the amplitude increases. This tendency is magni"ed by the fact that for strong non-linear e!ects, the exact invariant manifold often approaches a linear dependence on a (the behavior seen in Figure 1(a) is typical, and similar results may be seen in reference [2]). The asymptotic approximation in Figure 1(b) does seem to capture some basic properties of the manifold. However, at the intermediate amplitudes, the relative error is considerable. In addition, the asymptotic solution is considerably worse at other values of . Figure 13 shows both the calculated Galerkin invariant manifold and the corresponding asymptotic solution. The curve shown in Figure 1(b) corresponds to the "0 edge of Figure 13, although there is some di!erence in scale. Clearly, the asymptotic solution is quite limited in its applicability, although it is qualitatively consistent with the Galerkin solution.
In systems such as this, with many d.o.f.s, a poor asymptotic approximation of a single mode's contribution (as in Figure 1(a)) is likely. Even when the remaining modes are accurately approximated, this type of error may quickly dominate the associated reduced-order model. The Galerkin-based approach avoids this error by allowing one to ensure the accuracy of each mode's constraint relations throughout the full domain. Figure 14 illustrates the time response of the various periodic second-mode solutions at the beam's quarter-span, for a given amplitude. The Galerkin Manifold and Exact solutions are indistinguishable, while the One Mode Solution and the Asymptotic Manifold solution yield inaccurate results. The de#ections shown illustrate both the accuracy of this approach, and the inadequacy of the other methods. The ability to accurately account for the dynamic contribution of the other system linear modes allows the corresponding reduced-order model to be nearly exact. Of course, the Galerkin manifold only approaches the exact  manifold as the number of terms in the expansion is increased, and small errors are likely to remain. The degree of these errors (for the present solution) may be seen in Figure 15, which illustrates the periodic contribution of the sixth linear mode to the motions seen in Figure 14. It is apparent that although the exact solution is quite well approximated, some higher order e!ects are still not entirely accounted for. This could be remedied by increasing the number of harmonics and piecewise linear segments.
In Figure 16, the "fth NNM is depicted to highlight the non-similar nature of the physical mode shape of the beam as predicted by the Galerkin approach. The two curves show the small and large amplitude displacement con"gurations of the beam for motions within the NNM (at "0). They are normalized for plotting purposes by the amplitude of the "rst lobe of the mode shape. One should note that this shape not only changes with amplitude, but will also change dynamically throughout the course of a periodic motion. This information may be useful in a number of contexts, such as control algorithms or sensor Figure 16 Normalized mode shape of the beam for the "fth NNM at small (a"0)01, "0) and large (a"2)0, "0) amplitudes, obtained through the Galerkin-based solution procedure: **, small amplitude; ---, large amplitude. location. Of particular interest here is the increased participation of the end spring at larger amplitudes. Figure 17 depicts the contribution of the "rst linear mode to the "fth NNM, and demonstrates a worst-case scenario for the asymptotic approach. Not only is the domain of convergence small, but in addition, slightly outside this domain the error grows extremely rapidly. For this case, these errors quickly overwhelm any bene"ts of the asymptotic approach. In addition, this type of amplitude dependence is quite di$cult to capture using the orthogonal basis functions from the "rst example (as seen in Figure 3). The piecewise linear amplitude approach, however, performs very well, yielding an accurate reduced-order model over a wide amplitude range.
The amplitude}frequency relationship for the "fth NNM is shown in Figure 18. Several features of this "gure are worth noting. First, the asymptotic results obviously su!er greatly from the aforementioned errors. Second, the exact periodic response frequency appears to be approaching a limit with increased amplitude. This is to be expected since, for larger amplitudes, the non-linear torsional spring begins to resemble a clamped boundary condition, and the NNM should document this transition from pinned to clamped-like conditions at the left end for increasing amplitudes. Third, this transition requires interactions among the linear modes and, consequently, cannot be captured through a single-mode truncation. Here, unlike the "rst example, the one-mode model leads to frequency predictions which become qualitatively wrong at moderate amplitudes.

FURTHER CONSIDERATIONS
The model "delity which can be attained using this approach is compelling. However, as with all methods, there are limitations to its application. In particular, further study is necessary to determine what occurs when NNMs become unstable and bifurcate. Some work has been done in this area by Vakakis and others, e.g., reference [18]. Further development is necessary to ensure that the reduced-order models developed herein represent stable motions of the original system. Another limitation of this technique is that it relies on accurate modelling of the non-linearities within the original physical system. That is, if cubic non-linearities were included in the the original model as a means of capturing the leading-order non-linear e!ects, but higher order terms were neglected, the reduced-order model may be accurate only over a limited amplitude range. Even an exact reduction of an approximate system is limited in its applicability. This is in contrast to the asymptotic approaches, where the domain of validity for the reduced model generally does not exceed that of the original model. The consequence of this may be that relatively simple, low-order manifold solutions are su$cient for many systems, due to the uncertainties involved in the measurement and modelling of non-linear e!ects.
The individual NNMs discussed here are useful only when used independently. As with the asymptotic non-linear normal modes developed previously, and unlike linear modes, simultaneous motions within two or more non-linear normal modes are bound to interact. The two primary consequences of this are that (1) larger reduced-order models may not be assembled from the individual NNMs, and (2) sometimes these interactions are not Figure 19. Anomaly for the second NNM manifold, indicating a 6 : 1 interaction with the "fth linear mode: -----, asymptotic manifold; } } }, Galerkin manifold; *z*, exact. removable (typically due to an internal resonance) and certain individual NNMs may not be generated individually.
The ability to approximate the invariant manifold for a range of di!erent amplitudes also raises a host of issues. Internal resonances which are present at very low amplitudes may be avoided at higher amplitudes (due to amplitude-dependent frequency changes). Conversely, systems which have no internal resonances at low amplitudes may have signi"cant interactions at greater amplitudes. An instance of this e!ect was observed during the study of the second NNM for the "nite element beam system. These results are illustrated in Figure 19, which depicts the "fth linear mode contribution to the second NNM manifold (at "0). The anomaly occurs close to a 6 : 1 frequency relation between the second and "fth modes, and the Galerkin manifold contains signi"cant contributions from the sixth harmonic in . Although a complete study of this interaction is beyond the scope of this work, it should be noted that the Galerkin solution indicates its presence, while the asymptotic solution does not, and that the overall e!ect is quite weak. These issues may be addressed through a generalization of the present formulation, by developing multi-mode invariant manifolds using the non-linear Galerkin methodology, wherein several pairs of modal amplitudes and phases are chosen as &&master'' coordinates. For the above case, this would entail adding the amplitude and phase of the "fth mode to those of the second, and conducting and expansion in four variables. Of course, this requires considerably more expansion coe$cients, and the computational e!ort will become correspondingly larger. However, this is the next logical step for this line of work. 6. CONCLUSIONS This work expands upon the invariant manifold formulation of Shaw and Pierre [7], utilizing a co-ordinate transformation and a non-linear Galerkin solution procedure to generate reduced-order models for non-linear systems which are accurate, reliable, general, and based on fundamental principles. The use of the Galerkin projection enables an unprecedented solution accuracy for the manifold-governing equations. This accuracy allows the corresponding reduced-order models to take full advantage of the potential inherent in the invariant manifold formulation. The examples discussed clearly illustrate the accuracy attainable through this approach, with nearly all results showing little or no discrepancy from the numerically determined exact solution. As indicated by the "nite-element beam example, the method is not restricted to simple mass/spring systems. In addition, the computational foundation of the solution procedure allows the system non-linearities to take on a wide variety of forms with little or no additional analytic work. Furthermore, the formulation may be extended to produce reduced-order models which contain several modes or d.o.f.s, allowing the study of more general, multi-mode motions, as well as internal resonances. In summary, the method shows signi"cant promise for the development and understanding of reduced-order models for non-linear dynamic systems.