# A new Galerkin-based approach for accurate non-linear normal modes through invariant manifolds

Eric Pescheck, Christophe Pierre, Steven W. Shaw

► To cite this version:

Eric Pescheck, Christophe Pierre, Steven W. Shaw. A new Galerkin-based approach for accurate non-linear normal modes through invariant manifolds. Journal of Sound and Vibration, 2002, 249 (5), pp.971-993. 10.1006/jsvi.2001.3914. hal-01580920

HAL Id: hal-01580920

https://hal.science/hal-01580920

Submitted on 3 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A NEW GALERKIN-BASED APPROACH FOR ACCURATE NON-LINEAR NORMAL MODES THROUGH INVARIANT MANIFOLDS

E. PESHECK AND C. PIERRE

Department of Mechanical Engineering and Applied Mechanics, University of Michigan, Ann Arbor, MI 48109-2125, U.S.A. E-mail: epesh@adams.com; pierre@umich.edu

AND

S. W. SHAW

Department of Mechanical Engineering, Michigan State University, East Lansing, MI 48824-1226, U.S.A. E-mail: shawsw@egr.msu.edu

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.

## 1. INTRODUCTION

The generation of effective reduced-order models for multi-degree-of-freedom non-linear dynamic systems is difficult, 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 difficult and slow, as well as physically non-intuitive. Alternatively, smaller models often sacrifice 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.
Some early work along these lines concentrated on the definition and existence of non-linear normal modes (NNMs), their properties, and illustrative examples [1, 2]. These concepts were further extended by Rand [3, 4], and more recently by King and Vakakis [5, 6], among others. An interesting feature of these non-linear modes is that they necessarily interact, unless one specifies initial conditions for a pure mode. These interactions can become strong, and even essential to the overall dynamics, leading to instabilities and bifurcations of NNMs. On outcome of these bifurcations is that a system can possess more NNMs than d.o.f.s [16, 18]. Another feature of NNMs is that they do not adhere to superposition or orthogonality properties, in direct contrast with the linear case. Model reduction using NNMs can be achieved, however, due to the important property of invariance, which is common to linear and non-linear normal modes.

The present work builds upon the invariant manifold formulation developed by Shaw and Pierre [7, 8] and extended in reference [9]. This non-linear extension of modal analysis allows one to generate reduced-order models by restricting non-essential modes (or d.o.f.s) to an invariant manifold parameterized by a set of retained modes (or d.o.f.s). This invariant manifold allows the corresponding ROM to account for the participation of all d.o.f.s without requiring their explicit simulation. Although some work has been done using periodicity constraints to locate the invariant manifold numerically [10], most approaches use asymptotic methods to generate polynomials which approximate the manifold locally. This approach has been used successfully for the analysis of systems within the local non-linear regime [11–13]. However, as is typical of asymptotic methods, the resultant dynamics are only accurate for small motions, and the upper bound for these motions is not known a priori. In addition, due to the nature of the polynomial approximation, the divergence of the asymptotic approximation and the actual manifold may occur at small amplitudes—yielding entirely inaccurate results at amplitudes only slightly beyond the domain of linear validity. Hence, though the invariant manifold formulation of the NNM is sound, the asymptotic solution method sometimes yields results with limited applicability.

If non-linear modal analysis is to become a practical tool, it must be accurate, reliable, and general. Earlier efforts have not fully met these criteria. The various asymptotic methods are locally accurate, but the unknown extent of this accuracy, restrictions on the types of non-linearities that can be handled, and the effort necessary result in limited reliability and generality. Other analytical methods have typically restricted either the temporal or spatial behavior of the model, ultimately compromising both accuracy and generality. Lastly, numerical simulations may be used to locate the invariant manifold, and generate ROMs, but this requires extensive simulation and relies upon the periodicity of the system response. Hence, this approach is cumbersome and lacks generality.

The approach developed herein eliminates many of these problems. As with the asymptotic approach, the non-linear partial differential equations governing the manifold are solved approximately. However, instead of the local polynomial approximation, a series approximation is employed over a chosen domain. Non-linear equations in the expansion coefficients are obtained through a Galerkin projection and then solved to achieve an approximation of the manifold which minimizes the error for the selected set of basis functions, i.e., the error is minimized in the domain spanned by the chosen basis functions. Typical results are shown in Figure 1, which illustrates the difference between the asymptotic and Galerkin-based manifolds. The curves shown correspond to cross-sections of the second-mode manifold for the 10-mode finite element beam model discussed in section 4. This approach allows one to determine the manifold geometry to a chosen accuracy through a given (and possibly large) amplitude. Furthermore, the chosen formulation and computational nature of this approach allow it to be used for the analysis of a large variety of non-linear dynamic systems.
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Figure 1. Cross-sections of the exact, asymptotic and Galerkin manifolds at $\phi = 0$ for the second non-linear mode of the finite 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, $\phi$.</figcaption></figure>

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 differential equations (PDEs) developed herein, constitute unique contributions to the field. 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 finite element based beam model with a quadratic/cubic non-linear torsional spring at one end is examined. Results which correspond to the second and fifth 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 effects to remain weak. This accuracy enables one to use the reduced-order model with confidence in its validity, even for amplitudes corresponding to large non-linear effects.

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.

## 2. FORMULATION

This approach assumes that a reduced-order model is sought for a set of non-linear coupled oscillators. Whether from a system of discrete masses, finite element model, or discretized continuous system, these may be expressed in the form

$$ \ddot{\eta} + [Z]\dot{\eta} + [\Omega]\eta = f(\eta, \dot{\eta}), \quad (1) $$

where $\eta$ is a vector of normalized modal co-ordinates, $[Z]$ and $[\Omega]$ are diagonal damping and stiffness matrices, respectively, $f(\eta, \dot{\eta})$ contains any additional linear or non-linear terms, and each overdot indicates a time derivative. The individual elements of $[Z]$ and $[\Omega]$ are written as $2\zeta_i\omega_i$, and $\omega_i^2$, respectively, for the $i$th equation. This form is selected for analytical convenience, and does not reflect 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:

$$ \eta_i = X_i(\eta_k, \dot{\eta}_k), \quad \dot{\eta}_i = Y_i(\eta_k, \dot{\eta}_k) \quad \text{for } i = 1 \cdots N, \quad i \neq k \quad (2) $$

for $N$ original modal co-ordinates. Here, an alternative set of co-ordinates is defined using the transformation

$$ \eta_k = a \cos \phi, \quad \dot{\eta}_k = -a\omega_k \sin \phi \quad (3) $$

and, the original constraints, equation (2) become

$$ \eta_i = P_i(a, \phi), \quad \dot{\eta}_i = Q_i(a, \phi) \quad \text{for } i = 1 \cdots N, \quad i \neq k. \quad (4) $$

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 $k$th mode. These constraint relations are dictated for the $k$th mode (by equation (3)), and must be determined for all others. Determining these unknown relations, that is, the $P_i$ and $Q_i$, will ultimately yield reduced dynamic equations in terms of $a(t)$ and $\phi(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 offer the following advantages:
(1) The dynamics in the transformed co-ordinates are elementary for $f = 0$.

(2) The new constraint equations, ($P_i$, $Q_i$) must be periodic in $\phi$, and thus may be expressed with harmonic functions in $\phi$.

(3) The boundary values for the manifold are reduced from four values ($\pm u_k$, $\pm v_k$) to one—the upper bound of $a$.

The original equation of motion governing $\eta_k$ may be transformed into two first order equations in $a$ and $\phi$ [15]:

$$ \dot{a} = \left( \frac{-f_k}{\omega_k} - 2\zeta_k \omega_k a \sin \phi \right) \sin \phi, \quad \dot{\phi} = \omega_k - \left( \frac{f_k}{a\omega_k} + 2\zeta_k \omega_k \sin \phi \right) \cos \phi \quad (5) $$

and the remaining differential equations may be written in first order form as

$$ \dot{P}_i(a, \phi) = Q_i(a, \phi), \quad \dot{Q}_i(a, \phi) = -2\zeta_i \omega_i Q_i(a, \phi) - \omega_i^2 P_i(a, \phi) + f_i \quad \text{for } i = 1 \cdots N, i \neq k. \quad (6) $$

The time derivatives on the left-hand side of equation (6) may be expanded using the chain rule to obtain

$$ \dot{P}_i(a, \phi) = \frac{\partial P_i}{\partial a} \dot{a} + \frac{\partial P_i}{\partial \phi} \dot{\phi}, \quad \dot{Q}_i(a, \phi) = \frac{\partial Q_i}{\partial a} \dot{a} + \frac{\partial Q_i}{\partial \phi} \dot{\phi} \quad \text{for } i = 1 \cdots N, i \neq k. \quad (7) $$

Equations (5–7) are combined to yield a set of partial differential equations (PDEs) that are independent of time and govern the invariant manifold geometry:

$$ Q_i = \frac{\partial P_i}{\partial a} \left( \frac{-f_k}{\omega_k} - 2\zeta_k \omega_k a \sin \phi \right) \sin \phi + \frac{\partial P_i}{\partial \phi} \left[ \omega_k - \left( \frac{f_k}{a\omega_k} + 2\zeta_k \omega_k \sin \phi \right) \cos \phi \right] \quad \text{for } i = 1 \cdots N, i \neq k, \quad (8) $$

$$ -2\zeta_i \omega_i Q_i - \omega_i^2 P_i + f_i = \frac{\partial Q_i}{\partial a} \left( \frac{-f_k}{\omega_k} - 2\zeta_k \omega_k a \sin \phi \right) \sin \phi + \frac{\partial Q_i}{\partial \phi} \left[ \omega_k - \left( \frac{f_k}{a\omega_k} + 2\zeta_k \omega_k \sin \phi \right) \cos \phi \right] \quad \text{for } i = 1 \cdots N, i \neq k. \quad (9) $$

In previous works, a similar process is used to yield manifold-governing equations in terms of a master position and velocity, ($\eta_k$, $\dot{\eta}_k$) = ($u_k$, $v_k$). A polynomial expansion in $u_k$ and $v_k$ 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 coefficients, but here the functions of amplitude and phase ($a$ and $\phi$), are defined 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

$$ P_i(a, \phi) = \sum_{l=1}^{N_a} \sum_{m=1}^{N_\phi} C_i^{l,m} T_{l,m}(a, \phi), \quad Q_i(a, \phi) = \sum_{l=1}^{N_a} \sum_{m=1}^{N_\phi} D_i^{l,m} U_{l,m}(a, \phi) \quad \text{for } i = 1 \cdots N, i \neq k, \quad (10) $$
where the C's and D's are the unknown expansion coefficients and the $T_{l,m}$ and $U_{l,m}$ are known shape functions, typically composed of products of functions in $a$ and $\phi$. For example, a given $T_{l,m}$ might be a product of a harmonic function in $\phi$ and a polynomial function in $a$, defined over the domain $a \in [0, a_0]$, $\phi \in [0, 2\pi]$. $N_a$ and $N_\phi$ are used to denote the number of expansion functions used in $a$ and $\phi$ 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 = \int_{a,\phi} U_{p,q} \left[ -a \sum_{l,m} D_i^{l,m} U_{l,m} + \sum_{l,m} C_i^{l,m} \frac{\partial T_{l,m}}{\partial a} \left( \frac{-f_k}{\omega_k} - 2\zeta_k \omega_k a \sin \phi \right) a \sin \phi \right. \\ \left. + \sum_{l,m} C_i^{l,m} \frac{\partial T_{l,m}}{\partial \phi} \left[ a \omega_k - \left( \frac{f_k}{\omega_k} + 2\zeta_k \omega_k a \sin \phi \right) \cos \phi \right] \right] da d\phi, \quad (11)$$

$$0 = \int_{a,\phi} T_{p,q} \left[ 2\zeta_i \omega_i a \sum_{l,m} D_i^{l,m} U_{l,m} + \omega_i^2 a \sum_{l,m} C_i^{l,m} T_{l,m} - af_i \right. \\ + \sum_{l,m} D_i^{l,m} \frac{\partial U_{l,m}}{\partial a} \left( \frac{-f_k}{\omega_k} - 2\zeta_k \omega_k a \sin \phi \right) a \sin \phi \\ \left. + \sum_{l,m} D_i^{l,m} \frac{\partial U_{l,m}}{\partial \phi} \left[ a \omega_k - \left( \frac{f_k}{\omega_k} + 2\zeta_k \omega_k a \sin \phi \right) \cos \phi \right] \right] da d\phi \quad (12)$$

for $i = 1 \cdots N$, $i \neq k$, $p = 1 \cdots N_a$, and $q = 1 \cdots N_\phi$. This represents a set of $2(N-1) N_a N_\phi$ 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 coefficients for each $P_i(a, \phi)$, $Q_i(a, \phi)$ pair over the chosen domain, e.g., $\phi \in [0, 2\pi]$, $a \in [0, a_0]$.

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 effort. Also, the domain of convergence is known *a priori*, and the results obtained may be used with confidence throughout that region, once the desired accuracy is achieved.

Once the $P_i$ and $Q_i$ have reached the desired accuracy, $f_k$, in equation (5), may be evaluated for any given $(a, \phi)$, and the dynamics within the chosen mode may be determined by numerical simulations. The resultant time histories, $a(t)$ and $\phi(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 coefficients. For example, an incremental approach which separates each C and D as
<figure><img src="image_4.png" /><figcaption>Figure 2. Two-degree-of-freedom non-linear system with hardening cubic springs.</figcaption></figure>

$C_i = C_i^o + \delta C_i$, and $D_i = D_i^o + \delta D_i$ could be linearized locally (in the $\delta C$'s and $\delta 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 fidelity throughout a predetermined amplitude range.

In the next two sections, this formulation is applied to generate reduced-order models of two example systems. The first, 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 benefits associated with the non-linear Galerkin approach. A non-linear finite 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.

### 3. 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_1$ and $q_2$ are

$$ 
\begin{aligned}
\ddot{q}_1 + q_1 + 5(q_1 - q_2) &= -2q_1^3 - (q_1 - q_2)^3, \\
\ddot{q}_2 + 5(q_2 - q_1) &= -(q_2 - q_1)^3.
\end{aligned}
\qquad (13) $$

A linear eigenanalysis is carried out, and the natural frequencies are found to be $\omega_1 = 0.689$, $\omega_2 = 3.244$, and a transformation matrix of the form

$$ \begin{Bmatrix} q_1 \\ q_2 \end{Bmatrix} = \begin{bmatrix} 0.671 & 0.741 \\ 0.741 & -0.671 \end{bmatrix} \begin{Bmatrix} \eta_1 \\ \eta_2 \end{Bmatrix} $$

may be used to decouple the two equations to linear order. The resulting modal equations are

$$ \ddot{\eta}_1 + \omega_1^2 \eta_1 = f_1(\eta_1, \eta_2), \qquad \ddot{\eta}_2 + \omega_2^2 \eta_2 = f_2(\eta_1, \eta_2), \qquad (14) $$
where the non-linear terms may be written as

$$ 
\begin{aligned}
f_1 &= -0.405 \eta_1^3 - 1.34 \eta_1^2 \eta_2 - 1.51 \eta_1 \eta_2^2 - 0.349 \eta_2^3, \\
f_2 &= -0.448 \eta_1^3 - 1.51 \eta_1^2 \eta_2 - 1.05 \eta_1 \eta_2^2 - 4.58 \eta_2^3.
\end{aligned} 
\qquad (15) 
$$

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 coefficients. Many different sets of basis functions may be used to describe the manifold geometry. Here we use the following double expansion in $a$ and $\phi$:

$$ T_{l,m}(a, \phi) = L_l(a) \cos((m-1)\phi), \quad U_{l,m}(a, \phi) = L_l(a) \sin((m)\phi). $$

The harmonic functions are a natural choice for expanding the $\phi$ dependence over $[0, 2\pi]$, and the use of cosine functions for $T$ (corresponding to a modal position), and sine functions for $U$ (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 $\phi$ may be necessary. The functions $L_l(a)$ were chosen to be a set of polynomials defined over the domain $[0, a_0]$, with zero slope and $a = 0$, and satisfying the orthogonality property

$$ \int_{0}^{a_0} a L_i(a) L_j(a) = \begin{cases} 1 & \text{for } i = j, \\ 0 & \text{for } i \neq j. \end{cases} $$

Analytical expressions for the first seven $L_i(a)$ are given in Appendix A, and are illustrated in Figure 3. The restriction of the slope at $a = 0$ is a reflection 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 $L_i(a)$ must be chosen to allow for linear contributions to the invariant manifold at $a = 0$.

<figure><img src="image_5.png" /><figcaption>Figure 3. Orthogonal polynomials used for the $a$ component of the manifold expansion, shown for $a_0 = 1$.</figcaption></figure>
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

$$D = F(C),$$

where C and D represent vectors containing all the $C_i^{l,m}$, and $D_i^{p,q}$ 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 coefficients. 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.

<figure><img src="image_6.png" /><figcaption>Figure 4. Convergence of the first NNM invariant manifold. In plot (a): —, asymptotic; ---, (1, 2); ----, (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_a, N_\phi$), used in the expansion. Plots (a) and (b) illustrate cross-sections $P_2(a, 0)$, and $P_2(2 \cdot 22, \phi)$ respectively.</figcaption></figure>
### 3.1. THE FIRST NON-LINEAR NORMAL MODE

The first 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_2(a, \phi)$ and $Q_2(a, \phi)$. Two cross-sections of $P_2$, at $\phi = 0$ and $a = a_0 = 2.22$, are illustrated in Figures 4(a) and 4(b), respectively, for a number of different expansion orders, depicted in the form ($N_a, N_\phi$). Also shown are the asymptotic solution and an “exact” reference solution (described below). The asymptotic solution corresponds to a third order manifold in ($\eta_1, \dot{\eta}_1$), generated according to the analytical formulas presented in reference [9], and then transformed into the $(a, \phi)$ 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 configuration space for initial conditions which yield a periodic response, and consequently must lie on the invariant manifold. This search is carried out within the

<figure><img src="image_7.png" /><figcaption>Figure 5. The first NNM invariant manifold: The position (a), and velocity (b) of linear mode 2, as a function of the amplitude and phase of the first NNM (shown for $N_a = 5, N_\phi = 12$).</figcaption></figure>
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 five polynomials and 12 harmonics, there is still some small error near $a_0$. This error appears to be confined to a few regions of $\phi$, 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 $\phi$ 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_2(a, \phi)$ and $Q_2(a, \phi)$ are shown in Figure 5, for the solution corresponding to $N_a = 5$, and $N_\phi = 12$. These surfaces represent an approximate solution to equations (8) and (9), and allow the accurate extension of the first linear mode into the non-linear realm. Given these constraint relations, $f_1(\eta_1, \eta_2)$ in equation (15) becomes $f_1(a, \phi)$, and the amplitude-phase differential equations-equation (5)-may be integrated to find the correct single-mode motion. These manifolds go well beyond the domain of “small motions” and “first order effects”, allowing one to generate results which may be limited by the fidelity of the mathematical model, rather than by the reduction technique. Here, at $a = 2.22$ and $\phi = 0$ the calculated non-linear forces, $f_1(a, \phi)$, is slightly greater in magnitude than the linear force, $\omega_1^2 a$, corresponding to a large non-linear effect.

The degree of non-linearity is easily observed in Figure 6, wherein the response frequency can be seen to change significantly with amplitude. Although not initially apparent, the known manifold does not allow this curve to be extended to an amplitude of $a_0(\approx 2.22)$, as the peak amplitude for a given motion is not at $\phi = 0$, but at $\phi = \pi/2$. Motions with initial amplitudes greater than 1.5 will eventually exceed $a_0 \approx 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 first mode and ignores the effects 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

<figure><img src="image_8.png" /><figcaption>Figure 6. Response frequency for the first non-linear normal mode versus the modal amplitude, $a$, started at $\phi = 0$: ——, Galerkin manifold; ---, asymptotic manifold; ------, one-mode solution: +, exact. The Galerkin manifold corresponds to the $N_a = 5$, $N_\phi = 12$ solution.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7. Response of the first mass, $q_1$, for periodic motions in the first NNM, using various reduced systems: ——, Galerkin manifold; ---, asymptotic manifold; ----, one-mode solution; ·····, exact. The exact and Galerkin manifold ($N_a = 5, N_\phi = 12$) solutions are indistinguishable.</figcaption></figure>

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 different procedure. The “Galerkin Manifold” results are produced by simulating equation (5), using the ($N_a = 5, N_\phi = 12$) manifold to evaluate $f_1$, then using the manifold equations to compute the modal displacements for each $(a, \phi)$, and finally using the eigenvectors to compute the physical displacement. The “Asymptotic Manifold” uses the analytical solutions from reference [9] to determine $\eta_2 = g(\eta_1, \dot{\eta}_1)$, then simulates equation (14) (retaining terms of the proper order), determines the modal deflections, and, finally, the physical deflections. The “One Mode Solution” simply assumes $\eta_2(t) = 0$, simulates equation (14), and uses only the first 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 $\eta_1 = 1.5$ and $\dot{\eta}_1 = 0$, the different constraints on the second linear mode yield different 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 first NNM is simply an extension of the first linear mode into the non-linear realm, and a corresponding extension exists for the system’s second linear mode as well.

## 3.2. 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 $\phi = 0$ (a), and $a = a_0 = 3.0$ (b). For this mode the asymptotic manifold diverges quite quickly from the exact solution, making even the $N_a = 1, N_\phi = 2$ solution a considerable improvement. Plot (b) illustrates that the harmonic content of the manifold is simpler than that of the first normal mode, and this observation is
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8. Convergence of the second NNM invariant manifold. In plot (a): —, asymptotic; ---, (1, 2); —-, (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_a$, $N_φ$), used in the expansion. Plots (a) and (b) illustrate cross-sections $P_1(a, 0)$, and $P_1(3·0, φ)$ respectively.</figcaption></figure>

consistent with the fact that only half as many harmonics are necessary to reach convergence. As with the first mode, convergence appears to be the slowest near $a_0$. However, as was discussed above, the motions which begin at $φ = 0$, and near $a_0$, will exceed $a_0$ during the course of their motion (this is characteristic of systems which stiffen with increased amplitude). Consequently, convergence near the points $(a, φ) = (a_0, 0)$, $(a_0, π)$ and $(a_0, 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_1(a, φ)$ and $Q_1(a, φ)$ are shown in plots (a) and (b) of Figure 9, for an expansion with $N_a = 4$ polynomials, and $N_φ = 6$ harmonics. As before the expansions which describe these surfaces may be used to eliminate $η_1$ 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 different 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.
<figure><img src="image_12.png" /><figcaption>Figure 9. The second NNM manifold: The position (a), and velocity (b) of linear mode 1, as a function of the amplitude and phase of the second NNM (shown for $N_a = 4$, $N_\phi = 6$).</figcaption></figure>

That is, for motions occurring on this invariant manifold, the two linear modes are closer to moving in unison than for motions on the first NNM invariant manifold.

As with the first 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 first mode has little effect. However, an accurate response frequency does not guarantee an 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 first linear mode. Figure 8 shows that the asymptotic solution over-predicts the first linear mode participation, while the “One Mode Solution” assumes no participation of the first linear mode. Once again,
<figure><img src="image_13.png" /><figcaption>Figure 10. Response frequency for the second non-linear normal mode versus the modal amplitude, $a$, started at $\phi = 0$: ——, Galerkin manifold; ---, asymptotic manifold; ----, one-mode solution: +, exact. The Galerkin manifold corresponds to the $N_a = 4$, $N_{\phi} = 6$ solution.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 11. Response of the first mass, $q_1$, 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_a = 4$, $N_{\phi} = 6$) solutions are indistinguishable.</figcaption></figure>

the “Exact” solution is indistinguishable from the non-linear Galerkin results, highlighting
the accuracy of the approach.

4. 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 finite element code PATRAN was used to create a linear model of a beam with transverse deflection $u(x, t)$ using 200 two-noded beam finite 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
<figure><img src="image_15.png" /></figure>

Figure 12. Schematic representation of the finite element model: 200 linear beam elements were used to construct a beam with deflection $u(x, t)$, length $L = 1$ m, density $\rho = 7860$ kg/m$^3$, Young's modulus $E = 2 \times 10^{11}$ N/m$^2$, cross-sectional second moment of area $I = 5 \times 10^{-8}$ m$^4$, spring stiffness $k = 10^8$ N/m, and non-linear torsional stiffness $\gamma_t = 5000u'(0, t)^2 + 20000u'(0, t)^3$ 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 finite element system, each non-linear mode is associated with nine ($P, Q$) pairs. Consequently, a manifold with expansions of order $N_a = 5$, and $N_\phi = 12$, now requires a solution of 540 coupled non-linear equations for 540 coefficients. The computational effort necessary for the solution scales approximately with $(N_a N_\phi (N-1))^2$, 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 efficient solution, the domain in $a$ may be discretized into several smaller intervals. These intervals may then be described as linear in $a$ (fixing $N_a$ at two—one coefficient for each segment endpoint) and solved as individual sub-problems. The harmonic dependence is retained, so that each expansion now contains $2N_\phi$ coefficients, and corresponds to a thin annular strip of the manifold ($a \in [a_i, a_{i+1}]$ while $\phi \in [0, 2\pi]$). As with the polynomial approach, the velocity coefficients (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 benefits to this piecewise amplitude approach. The linear dependence on $a$ simplifies 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 difficulties 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 first example simply may not converge at all in such circumstances.

This piecewise strategy was used to examine the second and fifth NNMs of transverse vibration for the finite element beam system shown in Figure 12. Both solutions used $N_\phi = 8$ harmonics and the second NNM solution used 80 piecewise segments in $a$, while the fifth NNM solution used 40. Both manifolds were obtained for modal amplitudes of $a \in [0, 2]$, which corresponds to a maximum deflection of approximately 0.65 m. If these
<figure><img src="image_16.png" /><figcaption>Figure 13. The asymptotic and Galerkin manifolds for the third linear mode displacement contribution to the second NNM. The asymptotic solution is of third order, while the Galerkin solution uses 80 piecewise linear segments in $a$, and $N_{\phi} = 8$ harmonics in $\phi$.</figcaption></figure>

amplitudes seem unreasonably large for the 1·0 m long beam model, it should be realized that entirely equivalent results may be obtained for a stronger non-linearity through a simple rescaling of parameters. That is, if the slope/moment relationship from the non-linear torsional spring, $\gamma_t$, is adjusted from $[5000u'(0)^2 + (2 \times 10^4)u'(0)^3]$ to $[(5 \times 10^4)u'(0)^2 + (2 \times 10^6)u'(0)^3]$, (where $a'$ indicates a spatial derivative with respect to $x$) the results from $a \in [0, 2]$ are mapped onto the domain $a \in [0, 0.2]$.

A sample pair of cross-sections at $\phi = 0$ of the second-mode invariant manifold are displayed in Figure 1; these correspond to the components of the constrained displacements of the first and third linear modes. As with the first example, the asymptotic manifolds used here are calculated to third order according to the analytical solutions presented in reference [9]. Here, the deficiency 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 magnified by the fact that for strong non-linear effects, 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 $\phi$. Figure 13 shows both the calculated Galerkin invariant manifold and the corresponding asymptotic solution. The curve shown in Figure 1(b) corresponds to the $\phi = 0$ edge of Figure 13, although there is some difference 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 deflections 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
<figure><img src="image_17.png" /><figcaption>Figure 14. Time response for various second NNM solutions: —, Galerkin manifold: ---, asymptotic manifold: -.-.-, one mode solution; ----, exact. Again, the Galerkin manifold and exact solutions are indistinguishable.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 15. Phase space depiction of the periodic linear mode 6 contributions to the second NNM: —, Galerkin manifold; ---, exact.</figcaption></figure>

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 effects are still not entirely accounted for. This could be remedied by increasing the number of harmonics and piecewise linear segments.

In Figure 16, the fifth 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 configurations of the beam for motions within the NNM (at $\phi = 0$). They are normalized for plotting purposes by the amplitude of the first 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><img src="image_19.png" /><figcaption>Figure 16 Normalized mode shape of the beam for the fifth NNM at small ($a = 0.01, \phi = 0$) and large ($a = 2.0, \phi = 0$) amplitudes, obtained through the Galerkin-based solution procedure: ——, small amplitude; ---, large amplitude.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 17. Cross-section of the linear mode 1 contribution to the fifth NNM at $\phi = 0$: ——, asymptotic manifold; ---, Galerkin manifold; +, exact. The asymptotic solution is third order, while the Galerkin solution uses $N_\phi = 8$ harmonics, and 40 linear segments in $a$.</figcaption></figure>

location. Of particular interest here is the increased participation of the end spring at larger amplitudes.

Figure 17 depicts the contribution of the first linear mode to the fifth 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 benefits of the asymptotic approach. In addition, this type of amplitude dependence is quite difficult to capture using the orthogonal basis functions from the first 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 fifth NNM is shown in Figure 18. Several features of this figure are worth noting. First, the asymptotic results obviously suffer greatly
<figure><img src="image_21.png" /><figcaption>Figure 18. Response frequency versus amplitude for the fifth NNM, using various reduction techniques: ——, Galerkin manifold; ---, asymptotic manifold; - - - -, one mode solution; +, exact.</figcaption></figure>

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 first example, the one-mode model leads to
frequency predictions which become qualitatively wrong at moderate amplitudes.

5. FURTHER CONSIDERATIONS

The model fidelity 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 original model as a means of
capturing the leading-order non-linear effects, 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 sufficient for many systems, due to the uncertainties
involved in the measurement and modelling of non-linear effects.

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><img src="image_22.png" /><figcaption>Figure 19. Anomaly for the second NNM manifold, indicating a 6:1 interaction with the fifth linear mode: - - - - - , asymptotic manifold; - - - , Galerkin manifold; - - - - , exact.</figcaption></figure>

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 different 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 significant interactions at greater amplitudes. An instance of this effect was observed during the study of the second NNM for the finite element beam system. These results are illustrated in Figure 19, which depicts the fifth linear mode contribution to the second NNM manifold (at $\phi = 0$). The anomaly occurs close to a 6:1 frequency relation between the second and fifth modes, and the Galerkin manifold contains significant contributions from the sixth harmonic in $\phi$. 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 effect 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 fifth mode to those of the second, and conducting and expansion in four variables. Of course, this requires considerably more expansion coefficients, and the computational effort 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 finite-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 significant promise for the development and understanding of reduced-order models for non-linear dynamic systems.

## ACKNOWLEDGMENTS

This work was supported by a grant from the Army Research Office. Dr Gary Anderson is the grant monitor. The authors are grateful to Dongying Jiang for help with preparation of the manuscript.

## REFERENCES

1. R. M. ROSENBERG 1996 *Advances in Applied Mechanics* 9, 155–242. On nonlinear vibrations of systems with many degrees of freedom.

2. C. P. ATKINSON and B. TASKET 1965 *Journal of Vibration and Acoustics, Transactions of the American Society of Mechanical Engineers* 32, 359–364. A study of the nonlinearly related modal solutions of coupled nonlinear systems by superposition techniques.

3. R. H. RAND 1971 *Journal of Applied Mechanics, Transactions of the American Society of Mechanical Engineers* 38, 561. Nonlinear normal modes in two-degree-of-freedom systems.

4. R. H. RAND 1974 *International Journal of Non-linear Mechanics* 9, 363–368. A direct method for non-linear normal modes.

5. M. E. KING and A. F. VAKAKIS 1994 *Journal of Vibration and Acoustics, Transactions of the American Society of Mechanical Engineers* 116, 332–340. An energy based formulation for computing nonlinear normal modes in undamped continuous systems.

6. M. E. KING and A. F. VAKAKIS 1996 *Journal of Applied Mechanics, Transactions of the American Society of Mechanical Engineers* 63, 810–819. An energy based approach to computing resonant nonlinear normal modes.

7. S. W. SHAW and C. PIERRE 1993 *Journal of Sound and Vibration* 164, 85–124. Normal modes for non-linear vibratory systems.

8. S. W. SHAW and C. PIERRE 1994 *Journal of Sound and Vibration* 169, 319–347. Normal modes of vibration for non-linear continuous systems.

9. E. PESHECK, C. PIERRE, S. W. SHAW and N. BOIVIN 2001 *Nonlinear Dynamics*, to appear. Nonlinear modal analysis of structural systems using multi-mode invariant manifolds.

10. J. C. SLATER 1996 *Nonlinear Dynamics* 10, 19–30. A numerical method for determining nonlinear normal modes.

11. A. H. NAYFEH, C. CHIN and S. A. NAYFEH 1996 *Journal of Vibration and Acoustics, Transactions of the American Society of Mechanical Engineers* 118, 340–345. On nonlinear normal modes of systems with internal resonance.

12. S. W. SHAW, C. PIERRE and E. PESHECK 1999 *Shock and Vibration Digest* 31, 3–16. Modal analysis-based reduced-order models for nonlinear structures—an invariant manifold approach.

13. E. PESHECK, C. PIERRE and S. W. SHAW 2001 *Journal of Vibration and Acoustics, Transactions of the American Society of Mechanical Engineers*, to appear. Modal reduction of a nonlinear rotating beam through nonlinear normal modes.

14. A. STIENDL, H. TROGER and J. ZEMANN 1999 Improved galerkin method in the dimension reduction of nonlinear dynamical systems. In Francis C. Moon, editor, *IUTAM Symposium on New Applications of Nonlinear and Chaotic Dynamics in Mechanics. Proceedings of the IUTAM Symposium, Ithaca, NY, July 27–August 1, 1997*.
15. A. H. NAYFEH and D. T. MOOK 1979. *Nonlinear Oscillations*. New York: John Wiley & Sons.

16. A. F. VAKAKIS, L. I. MIKHLIN, Yu. V. M. PILIPCHUCK and A. A. ZEVIN 1996 *Normal Modes and Localization in Nonlinear Systems*. New York: John Wiley & Sons.

17. W. H. PRESS, S. A. TEUKOLSKY, W. T. VETTERLING and B. P. FLANNERY 1992 *Numerical Recipies in C: The Art of Scientific Computing*. Cambridge: Cambridge University Press.

18. A. F. VAKAKIS 1997 *Mechanical Systems and Signal Processing* **11**, 3–22. Non-linear normal modes (nnms) and their applications in vibration theory: an overview.

19. J. P. CUSUMANO, M. T. SHARKADY and B. W. KINBLE 1994 *Philosophical Transactions of the Royal Society of London* **347**, 421–438. Dynamics of a flexible beam impact oscillator.

20. X. MA, M. F. A. AZEEZ and A. F. VAKAKIS 2000 *Mechanical Systems and Signal Processing* **14**, 37–48. Non-linear normal modes and non-parametric system identification of non-linear oscillators.

APPENDIX A: EXPANSION POLYNOMIALS

The orthogonal polynomials used to express the amplitude dependence in the function expansion for the two-mass example are as follows:

$$L_1(a) = \sqrt{6}(a/a_0)^2, \quad L_2(a) = \sqrt{2}[-12(a/a_0)^2 + 14(a/a_0)^3],$$

$$L_3(a) = \sqrt{10}[-21(a/a_0)^2 + 56(a/a_0)^3 - 36(a/a_0)^4],$$

$$L_4(a) = \sqrt{3}[-112(a/a_0)^2 + 504(a/a_0)^3 - 720(a/a_0)^4 + 330(a/a_0)^5],$$

$$L_5(a) = \sqrt{14}[-126(a/a_0)^2 + 840(a/a_0)^3 - 1980(a/a_0)^4 + 1980(a/a_0)^5] - 715(a/a_0)^6],$$

$$L_6(a) = -1008(a/a_0)^2 + 9240(a/a_0)^3 - 31680(a/a_0)^4 + 51480(a/a_0)^5 \\ - 40040(a/a_0)^6 + 12012(a/a_0)^7,$$

$$L_7(a) = \sqrt{2}[-1386(a/a_0)^2 + 16632(a/a_0)^3 - 77220(a/a_0)^4 + 180180(a/a_0)^5] \\ - 225225(a/a_0)^6 + 144144(a/a_0)^7 - 37128(a/a_0)^8].$$