# Normal Modes for Non-Linear Vibratory Systems

Steven Shaw, Christophe Pierre

► To cite this version:

Steven Shaw, Christophe Pierre. Normal Modes for Non-Linear Vibratory Systems. Journal of Sound and Vibration, 1993, 164 (1), pp.85-124. 10.1006/jsvi.1993.1198. hal-01544500

HAL Id: hal-01544500

https://hal.science/hal-01544500

Submitted on 8 Jul 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
NORMAL MODES FOR NON-LINEAR
VIBRATORY SYSTEMS

S. W. SHAW AND C. PIERRE

*Department of Mechanical Engineering and Applied Mechanics, University of Michigan,*
*Ann Arbor, Michigan 48109, U.S.A.*

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

# 1. INTRODUCTION

The idea of developing some type of normal mode representation for non-linear systems is not new. The work of Rosenberg, the main results of which are gathered together in a review article [1], was ground-breaking and laid the foundations for most of what followed; for example, the works of Atkinson and Taskett [2], Szemplinska [3], Rand [4, 5], Greenberg and Yang [6], van der Varst [7], Yen [8] and Caughey and Vakakis [9]. The thesis of Vakakis [10] contains a nice account of this history and provides the most recent contribution.

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

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

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

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

## 2. NORMAL MODES FOR LINEAR SYSTEMS

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

Consider the free motion of a system with $N$ degrees of freedom which has been linearized about an equilibrium at the origin:

$$ M\ddot{x} + C\dot{x} + Kx = 0, \quad (1) $$

where $x = [x_1, x_2, \dots, x_N]^T$ is the displacement vector, a dot denotes a time derivative, and $M, C$ and $K$ are the inertia, damping (or gyroscopic) and stiffness matrices, respectively, obtained by linearizing the dynamics of a mechanical system about an equilibrium point. We require that $M$ be non-singular and allow $C$ and $K$ to be arbitrary, which allows for the inclusion of dissipative, gyroscopic and non-conservative terms. We also assume that the system has independent eigenvectors and, although it is not strictly required, that the system eigenvalues are all complex. The usual procedure is to assume a solution of the form

$$ x(t) = q e^{\lambda t}, \quad (2) $$

and substitute this into equation (1). This leads to a second order $N \times N$ eigenvalue problem which can be solved for the eigenvalues, $\lambda$, and the eigenvectors, $q$. When $C$ is non-zero, the $2N$ eigenvectors are generally complex and the $N$ equations of motion do not uncouple in this second order form except in special circumstances [11]. A first order formulation which includes the velocities as separate dependent variables provides the setting needed to uncouple the dynamics [12].
Premultiplying equation (1) by the inverse of M and introducing the velocity vector, y, equation (1) can be written in first order form as

$$ \dot{z} = Dz, \quad z = \begin{pmatrix} x \\ y \end{pmatrix}, \quad D = \begin{pmatrix} 0 & I \\ A & B \end{pmatrix} \qquad (3) $$

where $A = -M^{-1}K$ and $B = -M^{-1}C$. We denote their elements by $\alpha_{ij}$ and $\beta_{ij}$ ($i, j=1, \dots, N$), respectively. A solution of the form

$$ z(t) = \bar{z} e^{\lambda t} \qquad (4) $$

is assumed which, when substituted into equation (3), leads to a standard, first order, $2N \times 2N$ eigenvalue problem which can be solved for the $2N$ eigenvectors, $\bar{z}$, and the $2N$ eigenvalues, $\lambda$. The eigenvectors are not unique in that they provide only a direction in the phase space, but not a magnitude: each contains $2N-1$ specified ratios and a free (complex) multiplicative constant. Thus, with this first order representation, a motion purely in a single mode is one for which all co-ordinates have exactly the same time dependence up to constant (complex) scaling factors, i.e., the amplitude ratios. An important property of such a modal motion is that if one knows the motion of a single generalized co-ordinate, then the motions of all other co-ordinates are specified by that mode's eigenvector. Furthermore, because of the independence of the eigenvectors, an initial condition which is aligned with a single mode leads to a purely single-mode motion in which all other modes start with zero amplitude and remain quiescent for all time. This is the important *invariance property* of the modal subspaces.

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

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

$$ x_i(t) = a_{1i}x_1(t) + a_{2i}y_1(t), \quad y_i(t) = b_{1i}x_1(t) + b_{2i}y_1(t), \quad i=1, 2, \dots, N, \qquad (5a) $$

where $a_{11}=1$, $a_{21}=0$, $b_{11}=0$ and $b_{21}=1$. In vector form, this is written as

$$ \begin{pmatrix} x_1 \\ \vdots \\ x_N \\ y_1 \\ \vdots \\ y_N \end{pmatrix} = \begin{bmatrix} 1 & 0 \\ a_{12} & a_{22} \\ \vdots & \vdots \\ a_{1N} & a_{2N} \\ 0 & 1 \\ b_{12} & b_{22} \\ \vdots & \vdots \\ b_{1N} & b_{2N} \end{bmatrix} \begin{pmatrix} x_1 \\ y_1 \end{pmatrix}. \qquad (5b) $$
This equation is that of a plane in the 2N-dimensional phase space. Any motion that takes place in such a linear invariant subspace and satisfies the equations of motion is defined as a linear normal mode. (If $x_1$ turns out to be a nodal point for a given mode, i.e., $x_1 = y_1 \equiv 0$, the equations will be singular. This situation is easily remedied by selecting a different displacement/velocity pair.) Note that this definition allows for phase differences between components, as required for “non-classical” normal modes.

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

$$ \dot{x}_i = y_i, \quad \dot{y}_i = \sum_{j=1}^{N} (\alpha_{ij}x_j + \beta_{ij}y_j), \quad i=1, \dots, N, \qquad (6) $$

we substitute equation (5) into equation (6) for $i=1, 2, \dots, N$ and make use of

$$ \dot{x}_1 = y_1, \quad \dot{y}_1 = \sum_{j=1}^{N} (\alpha_{1j}(a_{1j}x_1 + a_{2j}y_1) + \beta_{1j}(b_{1j}x_1 + b_{2j}y_1)). \qquad (7) $$

Gathering terms in $x_1$ and $y_1$ in the resulting $2N-2$ equations and requiring that coefficients of $x_1$ and $y_1$ match gives $4N-4$ equations that are quadratic in the unknowns $a_{1i}$, $a_{2i}$, $b_{1i}$ and $b_{2i}$, $i=2, \dots, N$. These can be written in matrix form as

$$ \begin{pmatrix} \mathbf{0} & \mathbf{I} \\ \mathbf{A} & \mathbf{B} \end{pmatrix} \begin{bmatrix} 1 & 0 \\ a_{12} & a_{22} \\ \vdots & \vdots \\ a_{1N} & a_{2N} \\ 0 & 1 \\ b_{12} & b_{22} \\ \vdots & \vdots \\ b_{1N} & b_{2N} \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ a_{12} & a_{22} \\ \vdots & \vdots \\ a_{1N} & a_{2N} \\ 0 & 1 \\ b_{12} & b_{22} \\ \vdots & \vdots \\ b_{1N} & b_{2N} \end{bmatrix} \begin{pmatrix} 0 & 1 \\ \sum_{j=1}^{N} (\alpha_{1j}a_{1j} + \beta_{1j}b_{1j}) & \sum_{j=1}^{N} (\alpha_{1j}a_{2j} + \beta_{1j}b_{2j}) \end{pmatrix}. \quad (8) $$

Note that the time dependence of the problem has been eliminated by this procedure and that equation (8) can be solved for the geometry of the modal subspaces without solving the standard eigenvalue-eigenvector problem. For an oscillatory system, equation (8) has $N$ real solutions for the linear coefficients $(a_{1i}, a_{2i}, b_{1i}, b_{2i})_{i=2, \dots, N}$—one for each normal mode (see Appendix A for a justification; while there is a unique solution for each underdamped mode, there may be multiple solutions for overdamped modes, as their motions take place in one-dimensional spaces). To retrieve the modal dynamics of the system on each of the $N$ invariant planar subspaces, we introduce the co-ordinate transformation

$$ \begin{pmatrix} \mathbf{x}(t) \\ \mathbf{y}(t) \end{pmatrix} = \sum_{k=1}^{N} \begin{bmatrix} 1 & 0 \\ a_{12} & a_{22} \\ \vdots & \vdots \\ a_{1N} & a_{2N} \\ 0 & 1 \\ b_{12} & b_{22} \\ \vdots & \vdots \\ b_{1N} & b_{2N} \end{bmatrix} \begin{pmatrix} u_k(t) \\ v_k(t) \end{pmatrix} = \mathbf{U} \mathbf{w}(t), \qquad (9) $$

where $\mathbf{U}$ is defined as a $2N \times 2N$ matrix the $N$ pairs of columns of which are the $N$ solutions of equation (8), such that the $k$th normal mode shape consists of the $(2k-1)$th and $(2k)$th columns of $\mathbf{U}$. The vector $\mathbf{w} = [u_1, v_1, \dots, u_N, v_N]^T$ contains the modal co-ordinates, where $u_k$ and $v_k$ are the displacement and velocity, respectively, for the $k$th normal mode. We apply this co-ordinate transformation to the equations of motion (3) and use the property
satisfied by the pairs of columns of U (equation (8)). This leads to a set of 2N, pairwise coupled, first order equations:

$$ \begin{pmatrix} \dot{u}_1 \\ \dot{v}_1 \\ \vdots \\ \dot{u}_N \\ \dot{v}_N \end{pmatrix} = U^{-1}DU \begin{pmatrix} u_1 \\ v_1 \\ \vdots \\ u_N \\ v_N \end{pmatrix} = \begin{pmatrix} 0 & 1 & & & \\ g_1 & h_1 & & & \\ & & \ddots & & \\ & & & 0 & \\ & & & h_N & g_N \end{pmatrix} \begin{pmatrix} u_1 \\ v_1 \\ \vdots \\ u_N \\ v_N \end{pmatrix}, \quad (10) $$

where in equation (10) we have defined

$$ \begin{gathered} g_k = \alpha_{11} + \sum_{j=2}^{N} \alpha_{1j}(a_{1j})_k + \sum_{j=2}^{N} \beta_{1j}(b_{1j})_k, \\ h_k = \sum_{j=2}^{N} \alpha_{1j}(a_{2j})_k + \beta_{11} + \sum_{j=2}^{N} \beta_{1j}(b_{2j})_k, \quad k=1, \dots, N. \end{gathered} \quad (11) $$

Therefore, the dynamics on each invariant modal subspace is governed by the single-degree-of-freedom oscillator equation

$$ \begin{pmatrix} \dot{u}_k \\ \dot{v}_k \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -\omega_k^2 & -2\xi_k\omega_k \end{pmatrix} \begin{pmatrix} u_k \\ v_k \end{pmatrix}, \quad k=1, \dots, N, \quad (12) $$

where we have let $g_k = -\omega_k^2$ and $h_k = -2\xi_k\omega_k$, where $\omega_k$ and $\xi_k$ are the undamped natural frequency and damping ratio for the kth normal mode, and $u_k$ and $v_k$ are the modal displacement and velocity, respectively. We refer to the above as the *oscillator form* of the modal dynamics.

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

Equation (10) shows that the dynamics of the N-degree-of-freedom system reduces to that of N uncoupled oscillators. For given initial conditions $z(0) = Uw(0)$, the N second order initial value problems (12) can be solved for the modal dynamics and the general motion obtained as a superposition of uncoupled oscillator motions *via* equation (9).

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

### 3. NORMAL MODES FOR NON-LINEAR SYSTEMS

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

The equations of motion are assumed to be of the form

$$ \dot{x}_i = y_i, \dot{y}_i = f_i(x; y), \quad i=1, 2, \dots, N, \qquad (13) $$

where $x=[x_1, x_2, \dots, x_N]^T$ represents generalized co-ordinates (displacements or rotations) and $y=[y_1, y_2, \dots, y_N]^T$ contains the corresponding generalized velocities. The vector $f=[f_1, \dots, f_N]^T$ represents the forces and moments acting on the system normalized by the respective inertias. (In some applications one has to invert an inertia matrix which may depend on $x$ and $y$ in order to achieve the form in equation (13). Also, it may be that for a given problem it is more natural to use generalized displacements and their conjugate momenta. In this case a more general, first order formulation may be required. The steps below can easily be generalized to handle that case.).

We begin by assuming that there exists at least one motion for which all displacements and velocities are functionally related to a single displacement-velocity pair, which we choose arbitrarily here as the first displacement and velocity, $x_1$ and $y_1$. In order to implement this, we write $u=x_1$ and $v=y_1$ and express the other $x_i$'s and $y_i$'s functionally in terms of $u$ and $v$, as follows:

$$ x_i = X_i(u, v), \quad y_i = Y_i(u, v), \quad i=1, 2, 3, \dots, N, \qquad (14a) $$

where $X_1(u, v)=u$ and $Y_1(u, v)=v$. In vector form this is written as:

$$ \begin{bmatrix} x_1 \\ y_1 \\ x_2 \\ y_2 \\ \vdots \\ x_N \\ y_N \end{bmatrix} = \begin{bmatrix} u \\ v \\ X_2(u, v) \\ Y_2(u, v) \\ \vdots \\ X_N(u, v) \\ Y_N(u, v) \end{bmatrix}. \qquad (14b) $$

(We note that such a representation is not always possible. Cases for which it breaks down include all of those for which the linearized system cannot be uncoupled; for example, if multiplicities of eigenvalues occur with an associated non-zero nilpotent form, and situations in which other types of internal resonances occur, for example when there exists an unremovable non-linear coupling between the interacting modes. However, the method can be extended to handle these situations.)

Equation (14) is that of a constraint surface of dimension 2 (or, equivalently, of co-dimension $2N-2$) in the $2N$-dimensional phase space. Thus, we define a normal mode of motion for the non-linear, autonomous system as a motion which takes place on a two-dimensional invariant manifold in the system's phase space. This manifold has the following properties: it passes through a stable equilibrium point of the system and, at that point, it is tangent to a plane which is an eigenspace of the system linearized about that equilibrium.

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

$$ \dot{x}_i = \frac{\partial X_i}{\partial u} \dot{u} + \frac{\partial X_i}{\partial v} \dot{v}, \quad \dot{y}_i = \frac{\partial Y_i}{\partial u} \dot{u} + \frac{\partial Y_i}{\partial v} \dot{v}, \quad i=1, 2, 3, \dots, N. \quad (15) $$

Next, we substitute the equations of motion in for $\dot{x}_i$ and $\dot{y}_i$ and replace $x_i$ and $y_i$ everywhere by $X_i$ and $Y_i$ to obtain $2N-2$ equations which can be solved for $X_i$ and $Y_i$:

$$ Y_i(u, v) = \frac{\partial X_i(u, v)}{\partial u} v + \frac{\partial X_i(u, v)}{\partial v} f_1(u, X_2(u, v), \dots, X_N(u, v); v, Y_2(u, v), \dots, Y_N(u, v)), \\ f_1(u, X_2(u, v), \dots, X_N(u, v); v, Y_2(u, v), \dots, Y_N(u, v)) \\ = \frac{\partial Y_i(u, v)}{\partial u} v + \frac{\partial Y_i(u, v)}{\partial v} f_1(u, X_2(u, v), \dots, X_N(u, v); v, Y_2(u, v), \dots, Y_N(u, v)) $$

$$ i=1, 2, 3, \dots, N. \quad (16) $$

Note that for $i=1$, the equations are trivially satisfied.

In general, these functional equations are at least as difficult to solve as the original differential equations, but they do allow for an approximate solution in the form of power series expansions. Once they have been solved for the $X_i$'s and $Y_i$'s, the dynamics on the invariant manifold, which is the normal mode dynamics, can then be generated by simply substituting the $X_i$'s and $Y_i$'s in for $x_i$ and $y_i$ in the first pair of equations of motion; that is, the ones for $x_1$ and $y_1$. This results in the modal dynamic equation

$$ \dot{u} = v, \quad \dot{v} = f_1(u, X_2(u, v), \dots, X_N(u, v); v, Y_2(u, v), \dots, Y_N(u, v)), \quad (17) $$

where $u$ and $v$ represent the variables on the invariant manifold and correspond to projections of the modal dynamics onto $(x_1, y_1)$. In general, at each equilibrium point there are $N$ solutions for the $X_i$'s and $Y_i$'s and $N$ corresponding sets of equations of the form given in equation (17), one set for each mode.

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

## 4. APPROXIMATE SOLUTION FOR THE MODAL DYNAMICS

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

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

$$ 
\begin{aligned}
\dot{x}_j &= y_j, \\
\dot{y}_j &= a_{jk}x_k + \beta_{jk}y_k + \delta_{jkm}x_k x_m + \epsilon_{jkm}x_k y_m + \gamma_{jkm}y_k y_m \\
&\quad + \mu_{jkmq}x_k x_m x_q + \nu_{jkmq}x_k x_m y_q + \rho_{jkmq}x_k y_m y_q + \xi_{jkmq}y_k y_m y_q + \dots,
\end{aligned} 
\quad (18) 
$$

where the coefficients in the expansions are derived from straightforward differentiation of the forces $f_j$ ($j=1, \dots, n$) with respect to $x$ and $y$, and where the implicit summation notation is used. Note that, as written, the coefficients of the non-linear terms are not uniquely defined, but can be chosen in any way such that all terms up to the desired order are included. For example, for the $x_1x_2$ terms one can take $\delta_{j12}$ to be the required coefficient, and then $\delta_{j21}$ to be zero.

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

$$ 
\begin{aligned}
X_i(u, v) &= a_{1i}u + a_{2i}v + a_{3i}u^2 + a_{4i}uv + a_{5i}v^2 \\
&\quad + a_{6i}u^3 + a_{7i}u^2v + a_{8i}uv^2 + a_{9i}v^3 + \dots, \\
Y_i(u, v) &= b_{1i}u + b_{2i}v + b_{3i}u^2 + b_{4i}uv \\
&\quad + b_{5i}v^2 + b_{6i}u^3 + b_{7i}u^2v + b_{8i}uv^2 + b_{9i}v^3 + \dots,
\end{aligned} 
\quad i=1, 2, 3, \dots, N. 
\quad (19a) 
$$

In vector form this series can be expressed as

$$ 
\begin{bmatrix} x_1 \\ y_1 \\ x_2 \\ y_2 \\ \vdots \\ x_N \\ y_N \end{bmatrix} = \left\{ \begin{bmatrix} 1 & 0 \\ 0 & 1 \\ a_{12} & a_{22} \\ b_{12} & b_{22} \\ \vdots & \vdots \\ a_{1N} & a_{2N} \\ b_{1N} & b_{2N} \end{bmatrix} + \begin{bmatrix} 0 & 0 \\ 0 & 0 \\ a_{32}u + a_{42}v & a_{52}u \\ b_{32}u + b_{42}v & b_{52}u \\ \vdots & \vdots \\ a_{3N}u + a_{4N}v & a_{5N}u \\ b_{3N}u + b_{4N}v & b_{5N}u \end{bmatrix} \right. \\
\left. + \begin{bmatrix} 0 & 0 \\ 0 & 0 \\ a_{62}u^2 + a_{82}v^2 & a_{72}u^2 + a_{92}v^2 \\ b_{62}u^2 + b_{82}v^2 & b_{72}u^2 + b_{92}v^2 \\ \vdots & \vdots \\ a_{6N}u^2 + a_{8N}v^2 & a_{7N}u^2 + a_{9N}v^2 \\ b_{6N}u^2 + b_{8N}v^2 & b_{7N}u^2 + b_{9N}v^2 \end{bmatrix} \right\} \begin{pmatrix} u \\ v \end{pmatrix} + \dots, 
\quad (19b) 
$$
or, more compactly, as

$$z = \{\mathbf{m}_0 + \mathbf{m}_1(u, v) + \mathbf{m}_2(u, v)\} \begin{pmatrix} u \\ v \end{pmatrix} + \dots, \qquad (19c)$$

where $z = [x_1, y_1, x_2, y_2, \dots, x_N, y_N]^T$ and $\mathbf{m}_0, \mathbf{m}_1$ and $\mathbf{m}_2$ are $2N \times 2$ matrices. (Note that this $z$ is a simple reordering of the $z$ in section 2). The matrix $\mathbf{m}_0$ is the linear modal component, and $\mathbf{m}_1$ and $\mathbf{m}_2$ represent the effects of quadratic and cubic non-linear terms, respectively. Also, note that we have factored out a $(u, v)$ vector and that the given representation is not unique; it is, however, convenient for some of the calculations, and all possible representations will yield the same results.

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

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

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

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

Solutions of these $18(N-1)$ equations result in cubic series approximations for the $N$ non-linear normal modes. This procedure provides the geometric structure of the non-linear normal modes near the equilibrium point. The local approximation of the dynamic equations for each mode can then be constructed via equation (17), by using the series
expansion for $f_1$ given in equation (18), with $(x_j, y_j)$ replaced everywhere by the series representations of $(X_j(u, v), Y_j(u, v))$ for $j=1, 2, \dots, N$. This results in a single-degree-of-freedom non-linear oscillator equation which represents the dynamics of the system on an invariant, two-dimensional subspace which is tangent to the linear normal mode eigenspace at the equilibrium point. There will be $N$ of these oscillators, one for each normal mode.

We are now in a position to assemble the complete non-linear modal matrix $M$ and the affiliated modal co-ordinates $w$. The matrix $M$ is composed by gathering together into a $2N \times 2N$ matrix the $N$ individual modal matrices of size $2N \times 2$ given in equation (19). Similarly, the modal vector $w$ is constructed by noting that each mode will have its own $(u, v)$ pair associated with it, which we label as $(u_k, v_k)$ for $k=1, 2, \dots, N$. These are collected into the vector $w=[u_1, v_1, u_2, v_2, \dots, u_N, v_N]^T$. The complete transformation between the physical co-ordinates $z$ and the modal co-ordinates $w$ can now be written as

$$z = M(w)w = \{M_0 + M_1(w) + M_2(w)\}w + \dots, \quad (20)$$

where $M(w) = M_0 + M_1(w) + M_2(w) + \dots$, and where the $M_p$'s are assembled from the $m_p$'s from equation (19c) for $p=0, 1, 2$ (note that $M_0$ is simply $U_0$ of section 2). It is important to note that the coefficients $a_{ij}$ and $b_{ij}$ will be different for each mode, and if one were to write out the components of equation (20), coefficients of the form $a_{ijk}$ and $b_{ijk}$ should be used where, for the present case, $i=1, 2, \dots, 9$ corresponds to the number of terms retained in the expansions (here we include all quadratic and cubic terms), $j=2, 3, \dots, N$ corresponds to the number of constraints, and $k=1, 2, \dots, N$ to the number of modes.

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

## 5. TRANSFORMATION FROM PHYSICAL TO NON-LINEAR MODAL CO-ORDINATES

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

The geometry behind this idea is quite simple in concept. In linear systems, one has planar eigenspaces, and a general response can be broken down into modal components via a linear projection onto the eigenspaces, or reassembled by a linear recombination of the modal responses. When non-linearities are present, the modal subspaces are not flat, but curved, and the required projections and recombination must utilize curvilinear co-ordinates which reflect the nature of these subspaces.
For linear systems the co-ordinate transformation between physical and modal co-ordinates is given by the modal matrix $M_0$ above. When the equations of motion are written in terms of the modal co-ordinates, they are *uncoupled*, linear oscillators, as in equation (12). It is of interest to consider what happens if one writes the non-linear equations of motion in terms of the non-linear modal coordinates. In the following, we show how this transformation can be carried out and describe some general features of the modal dynamic equations. Section 6.3 provides an example for which these features are explicitly demonstrated.

Recombination of the dynamics of the original co-ordinates from the modal co-ordinates for the non-linear case is done *via* the transformation $z=M(w)w$ given in equation (20). In order to obtain the projection of the physical co-ordinate dynamics onto the modal subspaces, this transformation must be inverted. The inversion is carried out for the general case involving quadratic and cubic non-linearities in Appendix C. The inversion is much simpler when only cubic non-linearities are present, and can be achieved by the following iterative process, which exploits the fact that we are using series expansions.

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

$$ 
\begin{aligned}
w &= \{M_0 + M_2(w)\}^{-1}z + \dots \\
&= \{I + M_0^{-1}M_2(w)\}^{-1}M_0^{-1}z + \dots \\
&= \{I - M_0^{-1}M_2(w)\}M_0^{-1}z + \dots,
\end{aligned} 
\quad (21) $$

where $I$ is the $2N \times 2N$ identity matrix. Note that the right side depends on both $w$ and $z$. This can be remedied by substituting the expression for $w$ in equation (21) into the $w$ term on the right side and expanding in powers of $z$. Since $M_2(w)$ is quadratic in $w$, the leading order argument of $M_2$ is simply $M_0^{-1}z$, and the dependence on $w$ is pushed out to terms of higher order. This results in the inverse transformation

$$ w = \{I - M_0^{-1}M_2(M_0^{-1}z)\}M_0^{-1}z + \dots, \quad (22) $$

which is correct up to cubic terms in $z$, the order to which we are working. This is used in the examples worked out in section 6, where we consider only cubic order non-linearities.

We now transform the equations of motion to a modal co-ordinate representation by using a direct co-ordinate transformation. The equations of motion near the equilibrium at $z=0$ are written as

$$ \dot{z} = A(z)z = \{A_0 + A_1(z) + A_2(z)\}z + \dots, \quad (23) $$

where $A_0z$ represents the linearized dynamics, $A_1(z)z$ the quadratic terms and $A_2(z)z$ the cubic terms. The transformation (20) is applied directly to these equations. We write it in a slightly different way for notational convenience here, such that the $w$ vector is not separated out as a factor:

$$ z = \bar{M}(w), \quad (24) $$

where it is implicitly assumed that we use the series expansion for $\bar{M}$. The time derivative of this is required for the transformation; it is given by using the chain rule:

$$ \dot{z} = [\partial \bar{M}(w)/\partial w] \dot{w}. \quad (25) $$
Substituting of equations (24) and (25) into the equation of motion (23) yields

$$[\partial \bar{M}(\mathbf{w}) / \partial \mathbf{w}] \dot{\mathbf{w}} = A(\bar{M}(\mathbf{w})) \bar{M}(\mathbf{w}), \qquad (26)$$

which can be premultiplied by the inverse of the Jacobian of the transformation to give the equations of motion in terms of the modal co-ordinates:

$$\dot{\mathbf{w}} = (\partial \bar{M}(\mathbf{w}) / \partial \mathbf{w})^{-1} A(\bar{M}(\mathbf{w})) \bar{M}(\mathbf{w}). \qquad (27)$$

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

Even though the modal equations are coupled, we explore in section 6.3 the possible advantages of simply ignoring the coupling terms in the transformed equations in order to recombine the non-linear modal oscillator responses into a general motion. This is done by solving a non-linear initial value problem as follows. Consider the motion which ensues from an initial condition $z_0$. This initial condition can be transformed to modal co-ordinate initial conditions, $w_0 = [u_{10}, v_{10}, u_{20}, v_{20}, \dots, u_{N0}, v_{N0}]^T$, by using the inverse transformation given in Appendix C (or equation (21) if only cubic non-linearities are present). These initial conditions are for the coupled non-linear modal oscillators given in equation (27). These coupled equations could be numerically integrated to obtain the modal dynamics, $(u_k(t), v_k(t))$, which can be transformed back into the original co-ordinates via equation (20), but since these equations are coupled, the only advantage in doing this is that one can obtain information about modal amplitudes directly. In section 6.3 we compare results from a full simulation of the original equations with simulations and recombinations of the individual modal oscillators in which coupling is ignored, i.e., equation (17), and show that very good agreement is obtained for that example. Such an approach has no sound mathematical basis, but provides a means of constructing the motion of a non-linear system from its non-linear normal modes.

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

## 6. EXAMPLES

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

Most of the symbolic manipulations and all of the numerical work for the examples were done using Mathematica® [15], running on a Macintosh IIcx computer. Mathematica® uses a variable step size fourth order Runge-Kutta integrator for the numerical solution of
ordinary differential equations, which was tricked into using a constant step size of 0.1 in all cases. The results, as seen below, are quite accurate for this step size, and the constant step size was useful when manipulating some of the data.

## 6.1. EXAMPLE 1: A LINEAR, NON-CONSERVATIVE SYSTEM

The linear system to be considered is shown in Figure 1. It has no special symmetry. After setting up the general equations which are to be solved for the invariant modal subspaces, we take specific values for the inertia, stiffness and damping elements, since the solution for the coefficients is not easily obtained (if at all possible) in closed form, except in special cases.

<figure><img src="image_2.png" /><figcaption>Figure 1. The mechanical model for Example 1.</figcaption></figure>

We assume that the springs are unstressed in the equilibrium position and that the displacement co-ordinates $x_1$ and $x_2$ are zero at this point. The equations of motion in first order oscillator form are

$$ 
\begin{aligned}
\dot{x}_1 &= y_1, & \dot{y}_1 &= f_1(x_1, x_2; y_1, y_2) = -x_1(1+k) + kx_2 - c(y_1 - y_2), \\
\dot{x}_2 &= y_2, & \dot{y}_2 &= f_2(x_1, x_2; y_1, y_2) = kx_1 - x_2(1+k) + cy_1 - 2cy_2,
\end{aligned}
\quad (28a) $$

or, equivalently, in matrix form

$$ \begin{pmatrix} \dot{x}_1 \\ \dot{y}_1 \\ \dot{x}_2 \\ \dot{y}_2 \end{pmatrix} = \begin{bmatrix} 0 & 1 & 0 & 0 \\ -1-k & -c & k & c \\ 0 & 0 & 0 & 1 \\ k & c & -1-k & -2c \end{bmatrix} \begin{pmatrix} x_1 \\ y_1 \\ x_2 \\ y_2 \end{pmatrix} \quad \text{or} \quad \dot{z} = Az. \quad (28b) $$

A purely modal motion is one which can be expressed in the following form:

$$ 
\begin{aligned}
x_1 &= u, & y_1 &= v, \\
x_2 &= X_2(u, v) = a_1 u + a_2 v, & y_2 &= Y_2(u, v) = b_1 u + b_2 v.
\end{aligned}
\quad (29) $$

The procedure for solving for the coefficients is to first take the time derivative of the last two of these equations (the last $2N-2$ in general, see equation (15)) to obtain

$$ \dot{X}_2 = a_1 \dot{u} + a_2 \dot{v}, \quad \dot{Y}_2 = b_1 \dot{u} + b_2 \dot{v}. \quad (30) $$

The left side is then equated to $(\dot{x}_2, \dot{y}_2)$ from the equation of motion (28), and in the right side $(\dot{u}, \dot{v})$ is replaced by $(\dot{x}_1, \dot{y}_1)$ from equation (28). This gives

$$ 
\begin{aligned}
& y_2 = a_1 y_1 + a_2 f_1(x_1, x_2; y_1, y_2), \\
& f_2(x_1, x_2; y_1, y_2) = b_1 y_1 + b_2 f_1(x_1, x_2; y_1, y_2).
\end{aligned}
\quad (31) $$
Next, $x_1$, $y_1$, $x_2$ and $y_2$ are replaced everywhere by the constraint (29) and the right side is subtracted from both sides to produce the final equations which are to be solved for $a_1$, $a_2$, $b_1$ and $b_2$:

$$Y_2(u, v) - a_1v - a_2f_1(u, X_2(u, v); v, Y_2(u, v)) = 0,$$

$$f_2(u, X_2(u, v); v, Y_2(u, v)) - b_1v - a_2f_1(u, X_2(u, v); v, Y_2(u, v)) = 0, \quad (32a)$$

or, after substituting in for $f_1$, $f_2$, $X_2$ and $Y_2$,

$$u(a_2 + b_1 - a_2b_1c + a_2k - a_1a_2k) + v(-a_1 + b_2 + a_2c - a_2b_2c - a_2^2k) = 0,$$

$$\begin{aligned} & u(-a_1 + b_2 - 2b_1c - b_1b_2c + k - a_1k + b_2k - a_1b_2k) \\ & \quad + v(-a_2 - b_1 + c - b_2c - b_2^2c - a_2k - a_2b_2k) = 0. \end{aligned} \quad (32b)$$

There are four unknowns in these equations; these are determined by noting that the equations must hold for all values of $u$ and $v$.

Gathering the coefficients of $u$ and $v$ in equation (32b) yields a set of four equations to be solved for $a_1$, $a_2$, $b_1$ and $b_2$:

$$a_2 + b_1 - a_2b_1c + a_2k - a_1a_2k = 0,$$

$$-a_1 + b_2 + a_2c - a_2b_2c - a_2^2k = 0,$$

$$-a_1 + b_2 - 2b_1c - b_1b_2c + k - a_1k + b_2k - a_1b_2k = 0,$$

$$-a_2 - b_1 + c - b_2c - b_2^2c - a_2k - a_2b_2k = 0. \quad (33)$$

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

$$\text{mode 1: } a_1 = -1.1056, \quad a_2 = -0.1550, \quad b_1 = 0.4599, \quad b_2 = -0.9891$$

$$\text{mode 2: } a_1 = 0.9471, \quad a_2 = -0.1386, \quad b_1 = 0.1401, \quad b_2 = 0.9676$$

$$\begin{aligned} a_2 &= -0.1550, & a_2 &= -0.1386, \\ b_1 &= 0.4599, & b_1 &= 0.1401, \\ b_2 &= -0.9891, & b_2 &= 0.9676. \end{aligned} \quad (34)$$

This means simply that a motion in the first mode satisfies the conditions

$$x_2 = -1.1056x_1 - 0.1550y_1, \quad y_2 = 0.4599x_1 - 0.9891y_1, \quad (35)$$

and a motion in the second mode satisfies

$$x_2 = 0.9471x_1 - 0.1386y_1, \quad y_2 = 0.1401x_1 + 0.9676y_1. \quad (36)$$

As the damping factor $c$ is reduced towards zero, the first mode approaches the mode in which the masses are exactly out of phase, while the second mode approaches the exactly in-phase mode. The modal subspaces are composed of two planes in the four-dimensional
<figure><img src="image_3.png" /></figure>

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

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

<figure><img src="image_6.png" /><figcaption>Figure 2. The modal invariant planes: (a) $X_2$ vs. $(x_1, y_1)$ for the first mode; (b) $Y_2$ vs. $(x_1, y_1)$ for the first mode; (c) $X_2$ vs. $(x_1, y_1)$ for the second mode; (d) $Y_2$ vs. $(x_1, y_1)$ for the second mode.</figcaption></figure>

phase space. These are depicted in Figure 2 as the planes $X_2$ and $Y_2$ versus $(x_1, y_1)$ for each of the two modes.

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

$$ \dot{u}_1 = v_1, \quad \dot{v}_1 = -2.9676u_1 - 0.7517v_1, $$

or

$$ \ddot{u}_1 + 0.7517\dot{u}_1 + 2.9676u_1 = 0, \quad (37) $$

and for the second mode

$$ \dot{u}_2 = v_2, \quad \dot{v}_2 = -1.0109u_2 - 0.1483v_2, $$

or

$$ \ddot{u}_2 + 0.1483\dot{u}_2 + 1.0109u_2 = 0. \quad (38) $$

At this point the eigenvalues are easily determined to be

$$ \text{mode 1: } -0.3759 \pm 1.6812i, \quad \text{mode 2: } -0.0741 \pm 1.0027i, \quad (39) $$

where $i^2 = -1$.
Note that by this procedure one obtains the eigenspaces first, and then the eigenvalues associated with them.

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

$$ \begin{pmatrix} x_1 \\ y_1 \\ x_2 \\ y_2 \end{pmatrix} = \begin{bmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ a_{11} & a_{21} & a_{12} & a_{22} \\ b_{11} & b_{21} & b_{12} & b_{22} \end{bmatrix} \begin{pmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \end{pmatrix} \quad \text{or} \quad z = M_0 w, \qquad (40) $$

where $a_{ij}$ and $b_{ij}$ represent $a_i$ and $b_i$ for the jth mode. The inverse of the transformation is given by

$$ \begin{pmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \end{pmatrix} = \begin{bmatrix} 0.4602 & -0.0714 & -0.4865 & 0.0041 \\ 0.1468 & 0.4828 & -0.0795 & -0.5104 \\ 0.5398 & 0.0714 & 0.4865 & -0.0041 \\ -0.1468 & 0.5172 & 0.0795 & 0.5104 \end{bmatrix} \begin{pmatrix} x_1 \\ y_1 \\ x_2 \\ y_2 \end{pmatrix}, $$

or

$$ w = M_0^{-1} z. \qquad (41) $$

Here we can explicitly transform the original equations of motion to modal co-ordinates to obtain

$$ \begin{pmatrix} \dot{u}_1 \\ \dot{v}_1 \\ \dot{u}_2 \\ \dot{v}_2 \end{pmatrix} = \begin{bmatrix} 0 & 1 & 0 & 0 \\ -2.9676 & -0.7517 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & -1.0109 & -0.1483 \end{bmatrix} \begin{pmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \end{pmatrix}, $$

or

$$ \dot{w} = M_0^{-1} A_0 M_0 w, \qquad (42) $$

which is exactly the matrix form of equations (37) and (38). These are uncoupled oscillators which are not in the usual Jordan form, but in what we have called the oscillator form (see section 2 and Appendix A).

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

$$ w_0 = (1, 0, 0, 0)^T \quad \Rightarrow \quad z_0 = (1, 0, -1.1056, 0.4599)^T. \qquad (43) $$

A simulation of the equations of motion (28) started at this point results in the response given in Figure 3, which shows the modal phase planes $(u_1, v_1)$ and $(u_2, v_2)$. If the numerics were exact, the motion should remain exactly at the $(u_2, v_2)$ origin, and it is observed that the response of the second mode stays bounded within a disk with a radius on the order of $10^{-17}$.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 3. Trajectory in the individual modal phase planes for an initial condition started in the first mode: (a) first mode phase plane; (b) second mode phase plane.</figcaption></figure>

It is worth noting that if the damping coefficient, $c$, is increased to the point where both modes are overdamped, the quadratic equations (33) have six real solutions. These arise because we are searching for two-dimensional invariant subspaces, while the completely overdamped system possesses four one-dimensional invariant subspaces. In this case the solutions of equations (33) correspond to all possible pairings of the one-dimensional spaces into two-dimensional spaces, of which there are six possibilities. In such a situation, care must be taken in choosing which solution to use, since all four eigenspaces must be accounted for.

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

## 6.2. EXAMPLE 2: A NON-LINEAR, CONSERVATIVE SYSTEM

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

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

$$ 
\begin{aligned}
\dot{x}_1 &= y_1, & \dot{y}_1 &= f_1(x_1, x_2; y_1, y_2) = -x_1(1+k) - gx_1^3 + kx_2, \\
\dot{x}_2 &= y_2, & \dot{y}_2 &= f_2(x_1, x_2; y_1, y_2) = kx_1 - x_2(1+k),
\end{aligned} 
\quad (44a) $$
<figure><img src="image_9.png" /><figcaption>Figure 4. The physical model for Example 2.</figcaption></figure>

or, equivalently, in first order form as

$$ \begin{pmatrix} \dot{x}_1 \\ \dot{y}_1 \\ \dot{x}_2 \\ \dot{y}_2 \end{pmatrix} = \left\{ \begin{bmatrix} 0 & 1 & 0 & 0 \\ -(1+k) & 0 & k & 0 \\ 0 & 0 & 0 & 1 \\ k & 0 & -(1+k) & 0 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 & 0 \\ -gx_1^2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \right\} \begin{pmatrix} x_1 \\ y_1 \\ x_2 \\ y_2 \end{pmatrix}, $$

or

$$ \dot{z} = \{A_0 + A_2(z)\}z. \quad (44b) $$

A normal mode motion is assumed by requiring that $x_1, y_1, x_2$ and $y_2$ are related as follows:

$$ \begin{align} x_1 &= u, & y_1 &= v, \nonumber \\ x_2 &= X_2(u, v) = a_1 u + a_2 v + a_3 u^2 + a_4 uv + a_5 v^2 \nonumber \\ & \phantom{X_2} + a_6 u^3 + a_7 u^2 v + a_8 uv^2 + a_9 v^3 + \dots, \nonumber \\ y_2 &= Y_2(u, v) = b_1 u + b_2 v + b_3 u^2 + b_4 uv + b_5 v^2 \nonumber \\ & \phantom{Y_2} + b_6 u^3 + b_7 u^2 v + b_8 uv^2 + b_9 v^3 + \dots. \tag{45} \end{align} $$

The procedure goes exactly as in Example 1, but involves many more terms. For example, taking the time derivative of $X_2$ and $Y_2$ from equation (45) yields, by using the chain rule,

$$ \begin{align} \dot{X}_2 &= a_1 \dot{u} + a_2 \dot{v} + 2a_3 u \dot{u} + a_4 u \dot{v} + a_5 \dot{u} v + 2a_5 v \dot{v} + \dots, \nonumber \\ \dot{Y}_2 &= b_1 \dot{u} + b_2 \dot{v} + 2b_3 u \dot{u} + b_4 u \dot{v} + b_5 \dot{u} v + 2b_5 v \dot{v} + \dots, \tag{46} \end{align} $$

where we have shown only up to quadratic terms for conciseness of presentation, although cubic terms have been included in the analysis. The left sides are equated to the equations of motion for $\dot{x}_2$ and $\dot{y}_2$, respectively, and in the right sides $\dot{u}$ and $\dot{v}$ are replaced by $\dot{x}_1$ and $\dot{y}_1$, respectively, from the equations of motion. This results in

$$ \begin{align} y_2 &= a_1 y_1 + a_2 f_1(x_1, x_2; y_1, y_2) + 2a_3 x_1 y_1 + a_4 x_1 f_1(x_1, x_2; y_1, y_2) \nonumber \\ & \phantom{= a_1 y_1 + a_2 f_1(x_1, x_2; y_1, y_2) + } + a_4 y_1^2 + 2a_5 y_1 f_1(x_1, x_2; y_1, y_2) + \dots, \nonumber \\ f_2(x_1, x_2; y_1, y_2) &= b_1 y_1 + b_2 f_1(x_1, x_2; y_1, y_2) + 2b_3 x_1 y_1 \nonumber \\ & \phantom{=} + b_4 x_1 f_1(x_1, x_2; y_1, y_2) + b_4 y_1^2 \nonumber \\ & \phantom{=} + 2b_5 y_1 f_1(x_1, x_2; y_1, y_2) + \dots, \tag{47} \end{align} $$

where $(x_1, y_1)$ are used for $(u, v)$ for notational consistency at this point. The substitution of conditions (45) into these expressions yields the required conditions for a modal motion. These are then expanded out to cubic order in $u$ and $v$, and coefficients of $u, v, u^2, uv, v^2, u^3, u^2v, uv^2$ and $v^3$ are gathered together to provide the equations for the $a_j$ and $b_j$
coefficients. The right side is subtracted from the equation in order to give equations which are equal to zero. The results of this calculation are given here in the form where the left column describes the term from which the equation on the right was obtained, such that there are two equations for each term, one each from the two equations in (47):

$$ u \text{ term: } \begin{aligned} & a_2 + b_1 + a_2 k - a_1 a_2 k = 0, \\ & -a_1 + b_2 + k - a_1 k + b_2 k - a_1 b_2 k = 0, \end{aligned} $$

$$ v \text{ term: } \begin{aligned} & -a_1 + b_2 - a_2^2 k = 0, \\ & -a_2 - b_1 - a_2 k - a_2 b_2 k = 0, \end{aligned} $$

$$ u^2 \text{ term: } \begin{aligned} & a_4 + b_3 - a_2 a_3 k + a_4 k - a_1 a_4 k = 0, \\ & -a_3 + b_4 - a_3 k - a_3 b_2 k + b_4 k - a_1 b_4 k = 0, \end{aligned} $$

$$ uv \text{ term: } \begin{aligned} & -2a_3 + 2a_5 + b_4 - 2a_2 a_4 k + 2a_5 k - 2a_1 a_5 k = 0, \\ & -a_4 - 2b_3 + 2b_5 - a_4 k - a_4 b_2 k - a_2 b_4 k + 2b_5 k - 2a_1 b_5 k = 0, \end{aligned} $$

$$ v^2 \text{ term: } \begin{aligned} & -a_4 + b_5 - 3a_2 a_5 k = 0, \\ & -a_5 - b_4 - a_5 k - a_5 b_2 k - 2a_2 b_5 k = 0, \end{aligned} $$

$$ u^3 \text{ term: } \begin{aligned} & a_7 + b_6 + a_2 g - a_3 a_4 k - a_2 a_6 k + a_7 k - a_1 a_7 k = 0, \\ & -a_6 + b_7 + b_2 g - a_6 k - a_6 b_2 k - a_3 b_4 k + b_7 k - a_1 b_7 k = 0, \end{aligned} $$

$$ u^2 v \text{ term: } \begin{aligned} & -3a_6 + 2a_8 + b_7 - a_4^2 k - 2a_3 a_5 k - 2a_2 a_7 k + 2a_8 k - 2a_1 a_8 k = 0, \\ & -a_7 - 3b_6 + 2b_8 - a_7 k - a_7 b_2 k - a_4 b_4 k - 2a_3 b_5 k - a_2 b_7 k + 2b_8 k - 2a_1 b_8 k = 0, \end{aligned} $$

$$ uv^2 \text{ term: } \begin{aligned} & -2a_7 + 3a_9 + b_8 - 3a_4 a_5 k - 3a_2 a_8 k + 3a_9 k - 3a_1 a_9 k = 0, \\ & -a_8 - 2b_7 + 3b_9 - a_8 k - a_8 b_2 k - a_5 b_4 k - 2a_4 b_5 k - 2a_2 b_8 k + 3b_9 k - 3a_1 b_9 k = 0, \end{aligned} $$

$$ v^3 \text{ term: } \begin{aligned} & -a_8 + b_9 - 2a_5^2 k - 4a_2 a_9 k = 0, \\ & -a_9 - b_8 - a_9 k - a_9 b_2 k - 2a_5 b_5 k - 3a_2 b_9 k = 0. \end{aligned} \quad (48) $$

The first four equations can be solved for the linear mode shape coefficients (again there are six solutions of which we take the two real ones):

$$ \begin{aligned} a_1 &= b_2 = 1, & \text{with } a_2 &= b_1 = 0, & \text{for mode 1;} \\ a_1 &= b_2 = -1, & \text{with } a_2 &= b_1 = 0, & \text{for mode 2.} \end{aligned} \quad (49) $$

These check with the modes of the linearized, symmetric system.

The next six equations are uncoupled and are linear in the unknowns $a_3, a_4, a_5, b_3, b_4$ and $b_5$. The problem is homogeneous and has a zero solution for these coefficients. This is not unexpected since the physical system has no quadratic terms.

The last eight equations cannot be coupled to any higher order coefficients since we have not included any in our analysis. Also, based on the patterns observed in the equations, they will not be in any case. These equations are linear in the coefficients $a_6, a_7, a_8$,
$a_9, b_6, b_7, b_8$ and $b_9$ and have two inhomogeneous terms, $a_2g$ and $b_2g$, in the two $u^3$ equations. The solution of this linear problem is:

For mode 1:

$$ 
\begin{aligned}
a_6 &= g(k-3)/[2k(k-4)], & a_8 &= -3g/[2k(k-4)], \\
b_7 &= 3g(k-1)/[2k(k-4)], & b_9 &= a_8,
\end{aligned} 
\qquad a_7 = a_9 = b_6 = b_8 = 0. \tag{50}
$$

For mode 2:

$$ 
\begin{aligned}
a_6 &= g(3+7k)/[2k(4+9k)], & a_8 &= 3g/[2k(4+9k)], \\
b_7 &= 3g(1+3k)/[2k(4+9k)], & b_9 &= a_8,
\end{aligned} 
\qquad a_7 = a_9 = b_6 = b_8 = 0.
$$

These give, to cubic order, the coefficients which describe the shapes of the invariant manifolds corresponding to the non-linear normal modes. The two normal mode manifolds are therefore defined by, to cubic order:

For mode 1:

$$ 
\begin{pmatrix} x_1 \\ y_1 \\ x_2 \\ y_2 \end{pmatrix} = \left\{ \begin{bmatrix} 1 & 0 \\ 0 & 1 \\ 1 & 0 \\ 0 & 1 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ g[(k-3)u_1^2 - 3v_1^2]/[2k(k-4)] & 0 & 0 & 0 \\ 0 & 3g[(k-1)u_1^2 - v_1^2]/[2k(k-4)] & 0 & 0 \end{bmatrix} \right\} \begin{pmatrix} u_1 \\ v_1 \end{pmatrix}.
$$

For mode 2:

$$ 
\begin{pmatrix} x_1 \\ y_1 \\ x_2 \\ y_2 \end{pmatrix} = \left\{ \begin{bmatrix} 1 & 0 \\ 0 & 1 \\ -1 & 0 \\ 0 & -1 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ g[(3+7k)u_2^2 + 3v_2^2]/[2k(4+9k)] & 0 & 0 & 0 \\ 0 & 3g[(1+3k)u_2^2 + v_2^2]/[2k(4+9k)] & 0 & 0 \end{bmatrix} \right\} \begin{pmatrix} u_2 \\ v_2 \end{pmatrix}. \tag{51}
$$

$(u_1, v_1)$ and $(u_2, v_2)$ are the first and second mode displacement and velocity co-ordinates, respectively. These are assembled to form the transformation between physical and modal co-ordinates:

$$ z = \{M_0 + M_2(w)\}w + \dots \tag{52} $$
The inversion of this transformation was obtained in closed form to third order and is given here:

$$ w = \{Q_0 + Q_2(z)\}z + \dots, \qquad (53a) $$

$$ u_1 = \frac{x_1 + x_2}{2} - g \left[ \frac{(k-3)(x_1 + x_2)^3 - 3(y_1 + y_2)^2(x_1 + x_2)}{32k(k-4)} + \frac{(3+7k)(x_1 - x_2)^3 + 3(y_1 - y_2)^2(x_1 - x_2)}{32k(4+9k)} \right], $$

$$ v_1 = \frac{y_1 + y_2}{2} - 3g \left[ \frac{(k-1)(x_1 + x_2)^2(y_1 + y_2) - (y_1 + y_2)^3}{32k(k-4)} + \frac{(1+3k)(x_1 - x_2)^2(y_1 - y_2) + (y_1 - y_2)^3}{32k(4+9k)} \right], $$

$$ u_2 = \frac{x_1 - x_2}{2} + g \left[ \frac{(k-3)(x_1 + x_2)^3 - 3(y_1 + y_2)^2(x_1 + x_2)}{32k(k-4)} + \frac{(3+7k)(x_1 - x_2)^3 + 3(y_1 - y_2)^2(x_1 - x_2)}{32k(4+9k)} \right], $$

$$ v_2 = \frac{y_1 - y_2}{2} + 3g \left[ \frac{(k-1)(x_1 + x_2)^2(y_1 + y_2) - (y_1 + y_2)^3}{32k(k-4)} + \frac{(1+3k)(x_1 - x_2)^2(y_1 - y_2) + (y_1 - y_2)^3}{32k(4+9k)} \right]. \qquad (53b) $$

Note that, as required by this formulation, $u_1 + u_2 = x_1$ and $v_1 + v_2 = y_1$. This inverse is needed only if a specific motion is being studied, and that one can obtain the important qualitative features of a general motion by constructing the non-linear modal oscillators, which requires only the forward transformation (52).

The non-linear modal oscillators are obtained exactly as was done in the first example, by substituting the constraint equations for $x_2$ and $y_2$ into the first pair of equations of motion. For the present example, this gives, in second order form,

$$ \text{mode 1: } \ddot{u}_1 + u_1 + gu_1 \left[ \left( 1 + \frac{(3-k)}{2(k-4)} \right) u_1^2 + \frac{3}{2(k-4)} \dot{u}_1^2 \right] = 0. $$

$$ \text{mode 2: } \ddot{u}_2 + (1 + 2k)u_2 + gu_2 \left[ \left( 1 - \frac{(3+7k)}{2(4+9k)} \right) u_2^2 - \frac{3k}{2(4+9k)} \dot{u}_2^2 \right] = 0. \qquad (54) $$

For the present example, we demonstrate the usefulness of the method in determining the frequency-amplitude relationship. This is accomplished using straightforward perturbation methods on the non-linear oscillators. Lindstedt's method, as described in many non-linear oscillations textbooks, is used here. We present only the final results since the technique is well known and rather straightforward (see, for example, the procedure in reference [16]). The results are the following relationships between the frequency of
oscillation, $\omega$, and the peak displacement amplitude, $A$:

$$ \text{mode 1: } \omega = 1 + \frac{3gA^2}{16} + \dots $$

$$ \text{mode 2: } \omega = [1 + 2k]^{1/2} \left[ 1 + \frac{3gA^2}{16(1 + 2k)} \right] + \dots \quad (55) $$

These have been checked with simulations and are quite accurate out to amplitudes of order unity.

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

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

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

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

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

$$ \text{mode 1: } \frac{x_2}{x_1} = 1 + \frac{g(k-3)}{2k(k-4)} x_1^2 + \dots, \quad (56a) $$

$$ \text{mode 2: } \frac{x_2}{x_1} = -1 + \frac{g(3+7k)}{2k(4+9k)} x_1^2 + \dots \quad (56b) $$

It is interesting to note that the first mode amplitude ratio can be softening or hardening depending on the value of the coupling stiffness:

$$ k < 3, \text{ hardening}; \quad 3 < k < 4, \text{ softening}; \quad 4 < k, \text{ hardening}; $$

while the second mode is hardening for all $k>0$. This behavior is illustrated in Figure 6, which displays the amplitude $x_2$ versus $x_1$ for various values of $k$.
<figure><img src="image_10.png" /><figcaption>Figure 5. Simulation validation of the invariance of the non-linear normal modal subspaces; in all cases an initial modal amplitude of 1.0 was used: (a) trajectory as projected into the individual linear modal phase planes for an initial condition started in the first linear mode; (b) trajectory as projected into the individual linear modal phase planes for an initial condition started in the second linear mode; (c) trajectory as projected into the individual non-linear modal phase manifolds for an initial condition started in the first non-linear mode; (d) trajectory as projected into the individual non-linear modal phase manifolds for an initial condition started in the second non-linear mode.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 6. Bending of the first normal mode amplitude ratio for various values of the coupling stiffness. The peak amplitude of $x_2$ is plotted versus that of $x_1$ for a modal motion in the first mode.</figcaption></figure>

Some closing remarks for this example follow. First, note that when $k=4$, the non-linear terms in the first mode become singular. This is due to the fact that for $k=4$, a three-to-one internal resonance exists, in which case the motion of the non-linear system is known to involve important and unremovable coupling between the modes. In the present case, the modal reduction does not help at all since there are only two degrees of freedom. However, in higher order systems, the crucial coupling will be only between the resonant modes and one can reformulate the problem and look for an invariant four-dimensional manifold on which the behavior of the coupled modes can be studied. (The $k=-4/9$ singularity involves a similar internal resonance for an unstable system).

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

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

### 6.3. EXAMPLE 3: A NON-LINEAR, NON-CONSERVATIVE SYSTEM

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

The physical system, shown in Figure 7, is the same as that for Example 2, except that the dashpots of Example 1 are added in the same places as in Example 1. The equations
<figure><img src="image_12.png" /><figcaption>Figure 7. The physical model for Example 3.</figcaption></figure>

of motion are given by those in equation (28), with the simple addition of a term $-gx_1^3$ in the $\dot{y}_1$ equation. The procedure following the development of equations (45)-(47) is repeated, and the coefficients of the various terms involving $u$ and $v$ are collected to produce a set of 18 equations for the 18 $a_j$ and $b_j$ coefficients, just as in equation (48). These equations are more complicated than those in Example 2 and reduce to them when $c=0$. We do not give the details here, since a numerical approach is adopted with the parameter values $c=0.3$, $k=1.0$ and $g=0.5$.

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

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

In order to demonstrate the power of the approximate non-linear superposition method, a general initial condition which excites both modes is chosen: $(x_{10}, y_{10}, x_{20}, y_{20}) = (0, 0, 2.0, 0)$. The simulations are done in three ways. First, the full equations are simulated

<table><caption>TABLE 1
Solutions for Example 3</caption><thead><tr><th></th><th>Mode 1</th><th>Mode 2</th></tr></thead><tbody><tr><td>$a_1$</td><td>-1.106</td><td>0.947</td></tr><tr><td>$a_2$</td><td>-0.155</td><td>-0.139</td></tr><tr><td>$a_3$</td><td>0</td><td>0</td></tr><tr><td>$a_4$</td><td>0</td><td>0</td></tr><tr><td>$a_5$</td><td>0</td><td>0</td></tr><tr><td>$a_6$</td><td>0.033</td><td>0.155</td></tr><tr><td>$a_7$</td><td>-0.049</td><td>-0.007</td></tr><tr><td>$a_8$</td><td>-0.011</td><td>0.219</td></tr><tr><td>$a_9$</td><td>-0.019</td><td>-0.063</td></tr><tr><td>$b_1$</td><td>0.460</td><td>0.140</td></tr><tr><td>$b_2$</td><td>-0.989</td><td>0.968</td></tr><tr><td>$b_3$</td><td>0</td><td>0</td></tr><tr><td>$b_4$</td><td>0</td><td>0</td></tr><tr><td>$b_5$</td><td>0</td><td>0</td></tr><tr><td>$b_6$</td><td>0.209</td><td>0.053</td></tr><tr><td>$b_7$</td><td>0.199</td><td>0.024</td></tr><tr><td>$b_8$</td><td>0.082</td><td>0.077</td></tr><tr><td>$b_9$</td><td>0.032</td><td>0.245</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>Figure 8. The non-linear invariant modal surfaces (approximated to cubic order): (a) $X_2$ vs. $(x_1, y_1)$ for the first non-linear normal mode; (b) $Y_2$ vs. $(x_1, y_1)$ for the first non-linear normal mode; (c) $X_2$ vs. $(x_1, y_1)$ for the second non-linear normal mode; (d) $Y_2$ vs. $(x_1, y_1)$ for the second non-linear normal mode.</figcaption></figure>

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

The results of the simulations are presented in Figure 10, which depicts the four co-ordinates $(x_1, y_1, x_2, y_2)$ as functions of time. The first trace in each case is labeled "direct" and results from simulations of the original equations of motion. The second and third traces in each case are obtained using the non-linear and linear approximate superposition methods and are labeled "NLNM" and "linear", respectively. All three traces are then plotted together for comparison.

The projection of the equations of motion onto the linear modal subspaces is achieved by substituting only the linear parts of the modal expansions for $x_2$ and $y_2$ into the first
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>Figure 9. Simulation validation of the invariance of the non-linear normal modal subspaces; in all cases an initial modal amplitude of 0.5 was used: (a) trajectory as projected into the individual linear modal phase planes for an initial condition started in the first linear mode; (b) trajectory as projected into the individual linear modal phase planes for an initial condition started in the second linear mode; (c) trajectory as projected into the individual non-linear modal phase manifolds for an initial condition started in the first non-linear mode; (d) trajectory as projected into the individual non-linear modal phase manifolds for an initial condition started in the second non-linear mode.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 10. Simulation results for the initial condition, $(x_{10}, y_{10}, x_{20}, y_{20}) = (0, 0, 2, 0)$, used to demonstrate non-linear superposition. In all curves the following line types are used: ——, direct simulation of the original system of equations, labeled "direct"; ---, simulations of the oscillators on the non-linear modal subspaces started with initial conditions obtained from projecting the given initial condition into the non-linear modal subspaces, labeled "NLNM"; ···, simulations of the oscillators on the linear modal subspaces started with initial conditions obtained from projecting the given initial condition into the linear modal subspaces, labeled "linear". (a) Results for $x_1$; (b) results for $y_1$; (c) results for $x_2$; (d) results for $y_2$.</figcaption></figure>

two equations of motion. The resulting oscillators are:

$$ \text{mode 1: } \ddot{u}_1 + 0.752\dot{u}_1 + 2.968u_1 + 0.50u_1^3 = 0, $$

$$ \text{mode 2: } \ddot{u}_2 + 0.148\dot{u}_2 + 1.011u_2 + 0.50u_2^3 = 0. $$

(57)

Projection of the initial conditions onto the linear subspaces yields the linear approximation $(u_{10}, v_{10}, u_{20}, v_{20}) = (-0.973, -0.159, 0.973, 0.159)$. Simulations for the oscillators of
<figure><img src="image_17.png" /><figcaption>Figure 10—continued.</figcaption></figure>

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

The projection of the equations onto the non-linear normal modes utilizes the cubic approximation of the invariant manifolds and yields the oscillators:

$$ \text{mode 1: } \ddot{u}_1 + 0.752\dot{u}_1 + 2.968u_1 + 0.405u_1^3 - 0.011u_1^2\dot{u}_1 - 0.014u_1\dot{u}_1^2 + 0.009\dot{u}_1^3 = 0, $$

$$ \text{mode 2: } \ddot{u}_2 + 0.148\dot{u}_2 + 1.011u_2 + 0.329u_2^3 + 0.0002u_2^2\dot{u}_2 - 0.242u_2\dot{u}_2^2 - 0.011\dot{u}_2^3 = 0. \quad (58) $$

They are significantly different from those obtained by using the linear projections. Transformation of the initial conditions via the inverse of the non-linear transformation gives
$(u_{10}, v_{10}, u_{20}, v_{20}) = (-0.912, -0.235, 0.912, 0.235)$. Simulations for the oscillators in equation (58) and recombination back into the physical co-ordinates using the non-linear transformation results in the time traces labeled "NLNM", shown in Figure 10.

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

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

$$ 
\begin{aligned}
\ddot{u}_1 + 0.752\dot{u}_1 + 2.968u_1 + 0.215(u_1 + u_2)^3 - 0.107(\dot{u}_1 + \dot{u}_2)(u_1 + u_2)^2 &= 0, \\
\ddot{u}_2 + 0.148\dot{u}_2 + 1.011u_2 + 0.264(u_1 + u_2)^3 - 0.107(\dot{u}_1 + \dot{u}_2)(u_1 + u_2)^2 &= 0.
\end{aligned} 
\quad (59) 
$$

In terms of non-linear modal co-ordinates the equation of motion is given by

$$ 
\begin{aligned}
\ddot{u}_1 + 0.752\dot{u}_1 + 2.968u_2 + 0.405u_1^3 + 0.724u_1u_2(u_1 + u_2) - 0.011u_1^2\dot{u}_1 \\
- 0.014u_1\dot{u}_1^2 + 0.009\dot{u}_1^3 - 0.214u_1u_2(\dot{u}_1 + \dot{u}_2) - 0.107(\dot{u}_1u_2^2 + \dot{u}_2u_1^2) &= 0, \\
\ddot{u}_2 + 0.148\dot{u}_2 + 1.011u_2 + 0.329u_2^3 + 0.776u_1u_2(u_1 + u_2) + 0.0001u_2^2\dot{u}_2 \\
- 0.242u_2\dot{u}_2^2 - 0.011\dot{u}_2^3 + 0.214u_1u_2(\dot{u}_1 + \dot{u}_2) + 0.107(u_1^2\dot{u}_2 + \dot{u}_1u_2^2) &= 0.
\end{aligned} 
\quad (60) 
$$

Note that these reduce to the modal oscillators on the appropriate subspaces: namely, the first modal oscillator is determined by setting $u_2 = v_2 = 0$ and the second is obtained by taking $u_1 = v_1 = 0$. Also, note that equations (60) have special invariance properties not shared by equations (59). This is discussed in more detail in the following section.

## 7. CONCLUSIONS AND DIRECTIONS FOR FUTURE WORK

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

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

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

This last observation points to a relationship between the present approach and more
general invariant manifold theories for dynamical systems [21, 22]. The approximate con-
struction of the invariant manifold is very similar to what is done in center manifold theory
[13]. However, we are not constructing center, stable or unstable manifolds, but invariant
components in them which are tangent to the linear eigenspaces. These invariant compo-
nents of the usual center, stable and unstable manifolds span the local phase space near
the equilibrium point of interest in a non-linear manner given by the transformations
between physical and modal co-ordinates. The generation of the equations of motion on
these components is done exactly as it is done in center manifold theory. The existence of
these invariant manifolds has been known in the mathematics community for some time,
but they have not previously been exploited in the definition and construction of normal
modes for non-linear systems. For the mathematically inclined reader, the paper by
Fenichel [23] provides existence proofs for these manifolds, including the required condi-
tions on the eigenvalues.

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

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

These general ideas are evident in the third example considered in this paper (see section
6.3). In particular, note that by using the equations of motion expressed in terms of the
linear modal co-ordinates, equations (59), neither mode can be completely suppressed.
Even if one starts with zero energy in one mode, that is, $u_i=v_i=0$ for i=1 or 2, the other
mode will directly excite it, leading to a fully coupled motion. In contrast, by using the
non-linear modal co-ordinates, the equations of motion, equation (60), possess the desired
invariance property: if either mode is started with zero energy, it remains at zero energy
for all time. These concepts can be generalized to large scale systems in which an arbitrary
subset of modes is started with zero energy—they will remain quiescent for all time. This
is a direct consequence of invariance.
Finally, it is important to emphasize that the development of the ideas presented in this paper has been greatly facilitated by the availability of computer-assisted symbolic manipulations. In fact, many of the calculations carried out in the examples are similar to those presented in the chapter on center manifolds in the book by Rand and Armbruster [16].

ACKNOWLEDGMENTS

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

REFERENCES

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

2. C. P. ATKINSON and B. TASKETT 1965 *Journal of Applied Mechanics* **32**, 359-364. A study of the nonlinearly related modal solutions of coupled nonlinear systems by superposition techniques.

3. W. SZEMPLINSKA 1990 *The Behaviour of Nonlinear Vibrating Systems* (Vol. 2). Dordrecht: Kluwer.

4. R. H. RAND 1971 *Journal of Applied Mechanics* **38**, 561. Nonlinear normal modes in two degree-of-freedom systems.

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

6. H. J. GREENBERG and T.-L. YANG 1971 *International Journal of Nonlinear Mechanics* **6**, 311-326. Modal subspaces and normal mode vibrations.

7. P. VAN DER VARST 1982 *Doctoral Thesis, Technical University of Eindhoven, Eindhoven, Germany*. On normal mode vibrations of nonlinear conservative systems.

8. D. YEN 1974 *International Journal of Non-linear Mechanics* 9, 45-53. On the normal modes of nonlinear dual-mass systems.

9. T. K. CAUGHEY and A. F. VAKAKIS 1991 *International Journal of Non-linear Mechanics* **26**(1), 89-103. A method for examining steady-state solutions of forced discrete systems with strong nonlinearities.

10. A. F. VAKAKIS 1990 Ph.D. Dissertation, California Institute of Technology. Analysis and identification of linear and nonlinear normal modes in vibrating systems.

11. T. K. CAUGHEY 1960 *Journal of Applied Mechanics* **27**, 269-271. Classical normal modes in damped linear dynamics systems.

12. L. MEIROVITCH 1967 *Analytical Methods in Vibrations*. New York: Macmillan.

13. J. CARR 1981 *Applications of Centre Manifold Theory*. New York: Springer-Verlag.

14. T. K. CAUGHEY, A. F. VAKAKIS and J. SIVO 1990 *International Journal of Non-linear Mechanics* **25**, 521-533. Analytical study of similar normal modes and their bifurcations in a class of strongly nonlinear systems.

15. S. WOLFRAM 1988 *Mathematica*. Redwood City, California: Addison Wesley.

16. R. H. RAND and D. ARMBRUSTER 1987 *Perturbation Methods, Bifurcation Theory and Computer Algebra*. New York: Springer-Verlag.

17. R. ABRAHAM and J. MARSDEN 1978 *Foundations of Mechanics*. Reading, Massachusetts: Benjamin/Cummings.

18. C. PIERRE 1988 *Journal of Sound and Vibration* **126**, 485-502. Mode localization and eigenvalue loci veering phenomena in disordered structures.
19. C. PIERRE 1990 *Journal of Sound and Vibration* **139**, 111-132. Weak and strong vibration localization in disordered structures: A statistical investigation.

20. C. PIERRE and S. W. SHAW 1991 *International Journal of Bifurcation and Chaos* **1**(2), 471-475. Mode localization due to symmetry-breaking nonlinearities.

21. S. WIGGINS 1990 *Introduction to Applied Nonlinear Dynamical Systems and Chaos*. New York: Springer-Verlag.

22. J. GUCKENHEIMER and P. HOLMES 1986 *Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields*. New York: Springer-Verlag.

23. N. FENICHEL 1971 *Indiana University Mathematical Journal* **21**, 193-225. Persistence and smoothness of invariant manifolds for flows.

24. S. W. SHAW and C. PIERRE *Journal of Sound and Vibration* (to appear). Normal modes of vibration for nonlinear continuous systems.

25. M. HIRSCH and S. SMALE 1974 *Differential Equations, Dynamical Systems and Linear Algebra*. New York: Academic Press.

## APPENDIX A: PROOF OF THE EQUIVALENCE IN THE LINEAR CASE

Consider the $2N \times 2N$ eigenvalue problem obtained from equations (3) and (4):

$$D\bar{z} = \lambda\bar{z}. \quad (A1)$$

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

The N pairs of complex conjugate eigensolutions are $(\lambda_i, \bar{z}_i)$ and $(\lambda_i^*, \bar{z}_i^*)$, $i=1, \dots, N$, where the asterisk denotes a conjugate. Defining a modal matrix Z, the columns of which are the eigenvectors, and changing co-ordinates according to $z(t) = Z\epsilon(t)$, yields a set of $2N$ uncoupled first order equations:

$$\dot{\epsilon} = \Lambda\epsilon, \quad \Lambda = \text{diag}\{\lambda_1, \lambda_1^*, \dots, \lambda_i, \lambda_i^*, \dots, \lambda_N, \lambda_N^*\}. \quad (A2)$$

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

$$Z' = [R\bar{z}_1 \quad \bar{z}_1 \quad \dots \quad R\bar{z}_i \quad \bar{z}_i \quad \dots \quad R\bar{z}_N \quad \bar{z}_N], \quad (A3)$$

where R and I denote real and imaginary parts, respectively, of the eigenvectors $\bar{z}_j$. Letting $z(t) = Z'\eta(t)$, equation (3) is transformed into first order equations that are pairwise coupled:

$$\dot{\eta}(t) = \Lambda'\eta(t), \quad \Lambda' = \begin{bmatrix} R\lambda_1 & I\lambda_1 & & 0 \\ -I\lambda_1 & R\lambda_1 & & \\ & & \ddots & \\ 0 & & & R\lambda_N & I\lambda_N \\ & & & -I\lambda_N & R\lambda_N \end{bmatrix}. \quad (A4)$$

These pairs of equations correspond to underdamped modes associated with the complex conjugate pairs of eigenvalues, such that the corresponding dynamics occurs in a *two-dimensional* subspace obtained by the span of the real and imaginary parts of the complex conjugate eigenvectors. (For overdamped modes the modal dynamics occur in one-dimensional spaces.) A motion of the ith underdamped mode consists of decaying oscillations
described by the *i*th pair of first order equations in the system (equation (A4)). The rate of decay and the pseudofrequency of oscillations are $-R\lambda_i = \xi_i\omega_i$ and $I\lambda_i = \omega_i\sqrt{1-\xi_i^2}$, respectively, where $\omega_i$ and $\xi_i$ are the mode's undamped natural frequency and damping ratio, respectively.

A more physical description of the underdamped modal dynamics of the system can be achieved by seeking sets of oscillator equations in the form

$$ \begin{pmatrix} \dot{u}_i \\ \dot{v}_i \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -\omega_i^2 & -2\xi_i\omega_i \end{pmatrix} \begin{pmatrix} u_i \\ v_i \end{pmatrix}, \quad i=1, \dots, N, \qquad (A5) $$

where $u_i$ and $v_i$ are the “oscillator-form” modal co-ordinates for the *i*th normal mode. These can be obtained from the first order form modal co-ordinates by the transformation

$$ w(t) = W\epsilon(t), \quad W = \text{blockdiag}\left\{ \cdots \begin{pmatrix} 1 & 1 \\ \lambda_i & \lambda_i^* \end{pmatrix} \cdots \right\}, \qquad (A6) $$

where $w = [u_1, v_1, \dots, u_N, v_N]^T$ is the vector of oscillator-form modal co-ordinates (note that for an overdamped mode these modal co-ordinates would not be uniquely defined).

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

$$ \begin{pmatrix} x(t) \\ y(t) \end{pmatrix} = \sum_{k=1}^{N} \begin{bmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} \\ \vdots & \vdots \\ a_{1N} & a_{2N} \\ b_{11} & b_{21} \\ b_{12} & b_{22} \\ \vdots & \vdots \\ b_{1N} & b_{2N} \end{bmatrix}_k \begin{pmatrix} u_k(t) \\ v_k(t) \end{pmatrix}, \qquad (A7) $$

where the *k*th mode shape consists of the $(2k-1)$th and $(2k)$th columns of U.

The oscillator form of the normal modes can also be obtained directly by formulating a block eigenvalue problem. We simply need to write that $U^{-1}DU$ is block-diagonal. This leads, for the *k*th normal mode, to

$$ D \begin{bmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} \\ \vdots & \vdots \\ a_{1N} & a_{2N} \\ b_{11} & b_{21} \\ b_{12} & b_{22} \\ \vdots & \vdots \\ b_{1N} & b_{2N} \end{bmatrix}_k = \begin{bmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} \\ \vdots & \vdots \\ a_{1N} & a_{2N} \\ b_{11} & b_{21} \\ b_{12} & b_{22} \\ \vdots & \vdots \\ b_{1N} & b_{2N} \end{bmatrix}_k \begin{pmatrix} 0 & 1 \\ -\omega_k^2 & -2\xi_k\omega_k \end{pmatrix}, \quad k=1, \dots, N, \qquad (A8) $$
which we can regard as a dimension-two block eigenvalue problem, where the "eigenvec-
tors" have two columns and the "eigenvalues" are 2 × 2 matrices with two unknown
coefficients. The N eigensolutions are the N oscillator-form normal modes of the system.

We note that equation (8) in section 2 is simply the block eigenvalue problem, equation
(A8), provided that we normalize the kth mode shape such that in equation (A7) the first
displacement is simply taken as u_k(t) and the first velocity as v_k(t), i.e., a_11 = 1, a_21 = 0,
b_11 = 0 and b_21 = 1. With this normalization, the first and (N+1)th lines in equation (A8)
simply give the eigenvalues of the kth mode as

$$
\begin{gather}
\alpha_{11} + \sum_{j=2}^{N} \alpha_{1j}(a_{1j})_k + \sum_{j=2}^{N} \beta_{1j}(b_{1j})_k = -\omega_k^2, \nonumber \\
\sum_{j=2}^{N} \alpha_{1j}(a_{2j})_k + \beta_{11} + \sum_{j=2}^{N} \beta_{1j}(b_{2j})_k = -2\xi_k\omega_k, \quad k = 1, \dots, N, \tag{A9}
\end{gather}
$$

and equations (8) and (A8) become identical. Thus, the use of invariant subspace ideas
results in the same linear normal modes as those found by transforming the solution of
the standard eigenvalue problem to oscillator form and utilizing a special normalization.

APPENDIX B: EQUATIONS FOR THE COEFFICIENTS OF THE NORMAL MODES

In all cases, $j=2, \dots, N$.

u term, first equation:

$$
-a_1 k a_2 j a_{1k} + b_{1j} - a_2 j b_{1k} \beta_{1k} = 0. \quad (B1)
$$

v term, first equation:

$$
a_{1k}a_{jk} - a_{1k}a_{1k}b_{2j} - b_{1k}b_{2j}\beta_{1k} + b_{1k}\beta_{jk} = 0. \quad (B2)
$$

u term, second equation:

$$
-a_{1j} - a_{2j}a_{2k}a_{1k} + b_{2j} - a_{2j}b_{2k}\beta_{1k} = 0. \quad (B3)
$$

v term, second equation:

$$
a_{2k}a_{jk} - b_{1j} - a_{2k}a_{1k}b_{2j} - b_{2j}b_{2k}\beta_{1k} + b_{2k}\beta_{jk} = 0. \quad (B4)
$$

$u^2$ term, first equation:

$$
\begin{align}
& -a_{2j}a_{3k}a_{1k} + b_{3j} - a_{2j}b_{3k}\beta_{1k} + a_{4j}(-a_{1k}a_{1k} - b_{1k}\beta_{1k}) \nonumber \\
& \qquad -a_{1k}a_{1m}a_{2j}\delta_{1mk} - a_{1m}a_{2j}b_{1k}\varepsilon_{1mk} - a_{2j}b_{1k}b_{1m}\gamma_{1mk} = 0. \tag{B5}
\end{align}
$$

uv term, first equation:

$$
\begin{align}
& -2a_{3j} - a_{2j}a_{4k}a_{1k} + b_{4j} - a_{2j}b_{4k}\beta_{1k} + a_{5j}(-2a_{1k}a_{1k} - 2b_{1k}\beta_{1k}) \nonumber \\
& \quad + a_{4j}(-a_{2k}a_{1k} - b_{2k}\beta_{1k}) + (-a_{1m}a_{2j}a_{2k} - a_{1k}a_{2j}a_{2m})\delta_{1mk} \nonumber \\
& \quad + (-a_{2j}a_{2m}b_{1k} - a_{1m}a_{2j}b_{2k})\varepsilon_{1mk} + (-a_{2j}b_{1m}b_{2k} - a_{2j}b_{1k}b_{2m})\gamma_{1mk} = 0. \tag{B6}
\end{align}
$$

$v^2$ term, first equation:

$$
\begin{align}
& -a_{4j} - a_{2j}a_{5k}a_{1k} + b_{5j} - a_{2j}b_{5k}\beta_{1k} + a_{5j}(-2a_{2k}a_{1k} - 2b_{2k}\beta_{1k}) \nonumber \\
& \quad -a_{2j}a_{2k}a_{2m}\delta_{1mk} - a_{2j}a_{2m}b_{2k}\varepsilon_{1mk} - a_{2j}b_{2k}b_{2m}\gamma_{1mk} = 0. \tag{B7}
\end{align}
$$
u² term, second equation:

$$
\begin{align}
& a_{3k}(\alpha_{jk} - \alpha_{1k}b_{2j}) + b_{4j}(-a_{1k}\alpha_{1k} - b_{1k}\beta_{1k}) \nonumber \\
& \qquad + b_{3k}(-b_{2j}\beta_{1k} + \beta_{jk}) - a_{1k}a_{1m}b_{2j}\delta_{1mk} + a_{1k}a_{1m}\delta_{jmk} \nonumber \\
& \qquad - a_{1m}b_{1k}b_{2j}\varepsilon_{1mk} + a_{1m}b_{1k}\varepsilon_{jmk} - b_{1k}b_{1m}b_{2j}\gamma_{1mk} + b_{1k}b_{1m}\gamma_{jmk} = 0. \tag{B8}
\end{align}
$$

uv term, second equation:

$$
\begin{align}
& a_{4k}(\alpha_{jk} - \alpha_{1k}b_{2j}) - 2b_{3j} + b_{5j}(-2a_{1k}\alpha_{1k} - 2b_{1k}\beta_{1k}) \nonumber \\
& \quad + b_{4j}(-a_{2k}\alpha_{1k} - b_{2k}\beta_{1k}) + b_{4k}(-b_{2j}\beta_{1k} + \beta_{jk}) \nonumber \\
& \quad + (-a_{1m}a_{2k}b_{2j} - a_{1k}a_{2m}b_{2j})\delta_{1mk} + (a_{1m}a_{2k} + a_{1k}a_{2m})\delta_{jmk} \nonumber \\
& \quad + (-a_{2m}b_{1k}b_{2j} - a_{1m}b_{2j}b_{2k})\varepsilon_{1mk} + (a_{2m}b_{1k} + a_{1m}b_{2k})\varepsilon_{jmk} \nonumber \\
& \quad + (-b_{1m}b_{2j}b_{2k} - b_{1k}b_{2j}b_{2m})\gamma_{1mk} + (b_{1m}b_{2k} + b_{1k}b_{2m})\gamma_{jmk} = 0. \tag{B9}
\end{align}
$$

v² term, second equation:

$$
\begin{align}
& a_{5k}(\alpha_{jk} - \alpha_{1k}b_{2j}) - b_{4j} + b_{5j}(-2a_{2k}\alpha_{1k} - 2b_{2k}\beta_{1k}) \nonumber \\
& \quad + b_{5k}(-b_{2j}\beta_{1k} + \beta_{jk}) - a_{2k}a_{2m}b_{2j}\delta_{1mk} \nonumber \\
& \quad + a_{2k}a_{2m}\delta_{jmk} - a_{2m}b_{2j}b_{2k}\varepsilon_{1mk} + a_{2m}b_{2k}\varepsilon_{jmk} \nonumber \\
& \quad - b_{2j}b_{2k}b_{2m}\gamma_{1mk} + b_{2k}b_{2m}\gamma_{jmk} = 0. \tag{B10}
\end{align}
$$

u³ term, first equation:

$$
\begin{align}
& a_{3k}a_{4j}\alpha_{1k} + a_{2j}a_{6k}\alpha_{1k} - b_{6j} + a_{4j}b_{3k}\beta_{1k} + a_{2j}b_{6k}\beta_{1k} \nonumber \\
& \quad + a_{7j}(a_{1k}\alpha_{1k} + b_{1k}\beta_{1k}) + (a_{1m}a_{2j}a_{3k} + a_{1k}a_{2j}a_{3m} + a_{1k}a_{1m}a_{4j})\delta_{1mk} \nonumber \\
& \quad + (a_{2j}a_{3m}b_{1k} + a_{1m}a_{4j}b_{1k} + a_{1m}a_{2j}b_{3k})\varepsilon_{1mk} \nonumber \\
& \quad + (a_{4j}b_{1k}b_{1m} + a_{2j}b_{1m}b_{3k} + a_{2j}b_{1k}b_{3m})\gamma_{1mk} + a_{1k}a_{1m}a_{1q}a_{2j}\mu_{1mkq} \nonumber \\
& \quad + a_{1k}a_{1m}a_{2j}b_{1q}\nu_{1mkq} + a_{1m}a_{2j}b_{1k}b_{1q}\rho_{1mkq} + a_{2j}b_{1k}b_{1m}b_{1q}\xi_{1mkq} = 0. \tag{B11}
\end{align}
$$

u²v term, first equation:

$$
\begin{align}
& -3a_{6j} - a_{4j}a_{4k}\alpha_{1k} - 2a_{3k}a_{5j}\alpha_{1k} - a_{2j}a_{7k}\alpha_{1k} + b_{7j} - 2a_{5j}b_{3k}\beta_{1k} \nonumber \\
& \quad - a_{4j}b_{4k}\beta_{1k} - a_{2j}b_{7k}\beta_{1k} + a_{8j}(-2a_{1k}\alpha_{1k} - 2b_{1k}\beta_{1k}) + a_{7j}(-a_{2k}\alpha_{1k} - b_{2k}\beta_{1k}) \nonumber \\
& \quad + (-a_{2j}a_{2m}a_{3k} - a_{2j}a_{2k}a_{3m} - a_{1m}a_{2k}a_{4j} - a_{1k}a_{2m}a_{4j} - a_{1m}a_{2j}a_{4k} - a_{1k}a_{2j}a_{4m} - 2a_{1k}a_{1m}a_{5j})\delta_{1mk} \nonumber \\
& \quad + (-a_{2m}a_{4j}b_{1k} - a_{2j}a_{4m}b_{1k} - 2a_{1m}a_{5j}b_{1k} - a_{2j}a_{3m}b_{2k} - a_{1m}a_{4j}b_{2k} - a_{2j}a_{2m}b_{3k} - a_{1m}a_{2j}b_{4k})\varepsilon_{1mk} \nonumber \\
& \quad + (-2a_{5j}b_{1k}b_{1m} - a_{4j}b_{1m}b_{2k} - a_{4j}b_{1k}b_{2m} - a_{2j}b_{2m}b_{3k} - a_{2j}b_{2k}b_{3m} - a_{2j}b_{1m}b_{4k} - a_{2j}b_{1k}b_{4m})\gamma_{1mk} \nonumber \\
& \quad + (-a_{1m}a_{1q}a_{2j}a_{2k} - a_{1k}a_{1q}a_{2j}a_{2m} - a_{1k}a_{1m}a_{2j}a_{2q})\mu_{1mkq} \nonumber \\
& \quad + (-a_{1m}a_{2j}a_{2k}b_{1q} - a_{1k}a_{2j}a_{2m}b_{1q} - a_{1k}a_{1m}a_{2j}b_{2q})\nu_{1mkq} \nonumber \\
& \quad + (-a_{2j}a_{2m}b_{1k}b_{1q} - a_{1m}a_{2j}b_{1q}b_{2k} - a_{1m}a_{2j}b_{1k}b_{2q})\rho_{1mkq} \nonumber \\
& \quad + (-a_{2j}b_{1m}b_{1q}b_{2k} - a_{2j}b_{1k}b_{1q}b_{2m} - a_{2j}b_{1k}b_{1m}b_{2q})\xi_{1mkq} = 0. \tag{B12}
\end{align}
$$
uv² term, first equation:

$$
\begin{aligned}
& -2a_{7j} - 2a_{4k}a_{5j}\alpha_{1k} - a_{4j}a_{5k}\alpha_{1k} - a_{2j}a_{8k}\alpha_{1k} + b_{8j} - 2a_{5j}b_{4k}\beta_{1k} \\
& \quad - a_{4j}b_{5k}\beta_{1k} - a_{2j}b_{8k}\beta_{1k} + a_{9j}(-3a_{1k}\alpha_{1k} - 3b_{1k}\beta_{1k}) + a_{8j}(-2a_{2k}\alpha_{1k} - 2b_{2k}\beta_{1k}) \\
& \quad + (-a_{2k}a_{2m}a_{4j} - a_{2j}a_{2m}a_{4k} - a_{2j}a_{2k}a_{4m} - 2a_{1m}a_{2k}a_{5j} - 2a_{1k}a_{2m}a_{5j} - a_{1m}a_{2j}a_{5k} - a_{1k}a_{2j}a_{5m})\delta_{1mk} \\
& \quad + (-2a_{2m}a_{5j}b_{1k} - a_{2j}a_{5m}b_{1k} - a_{2m}a_{4j}b_{2k} - a_{2j}a_{4m}b_{2k} - 2a_{1m}a_{5j}b_{2k} - a_{2j}a_{2m}b_{4k} - a_{1m}a_{2j}b_{5k})\epsilon_{1mk} \\
& \quad + (-2a_{5j}b_{1m}b_{2k} - 2a_{5j}b_{1k}b_{2m} - a_{4j}b_{2k}b_{2m} - a_{2j}b_{2m}b_{4k} - a_{2j}b_{2k}b_{4m} - a_{2j}b_{1m}b_{5k} - a_{2j}b_{1k}b_{5m})\gamma_{1mk} \\
& \quad + (a_{1q}a_{2j}a_{2k}a_{2m} - a_{1m}a_{2j}a_{2k}a_{2q} - a_{1k}a_{2j}a_{2m}a_{2q})\mu_{1mkq} \\
& \quad + (-a_{2j}a_{2k}a_{2m}b_{1q} - a_{1m}a_{2j}a_{2k}b_{2q} - a_{1k}a_{2j}a_{2m}b_{2q})\nu_{1mkq} \\
& \quad + (-a_{2j}a_{2m}b_{1q}b_{2k} - a_{2j}a_{2m}b_{1k}b_{2q} - a_{1m}a_{2j}b_{2k}b_{2q})\rho_{1mkq} \\
& \quad + (-a_{2j}b_{1q}b_{2k}b_{2m} - a_{2j}b_{1m}b_{2k}b_{2q} - a_{2j}b_{1k}b_{2m}b_{2q})\xi_{1mkq} = 0.
\end{aligned}
\quad (B13)
$$

v³ term, first equation:

$$
\begin{aligned}
& a_{8j} + 2a_{5j}a_{5k}\alpha_{1k} + a_{2j}a_{9k}\alpha_{1k} - b_{9j} + 2a_{5j}b_{5k}\beta_{1k} + a_{2j}b_{9k}\beta_{1k} + a_{9j}(3a_{2k}\alpha_{1k} + 3b_{2k}\beta_{1k}) \\
& \quad + (2a_{2k}a_{2m}a_{5j} + a_{2j}a_{2m}a_{5k} + a_{2j}a_{2k}a_{5m})\delta_{1mk} + (2a_{2m}a_{5j}b_{2k} + a_{2j}a_{5m}b_{2k} + a_{2j}a_{2m}b_{5k})\epsilon_{1mk} \\
& \quad + (2a_{5j}b_{2k}b_{2m} + a_{2j}b_{2m}b_{5k} + a_{2j}b_{2k}b_{5m})\gamma_{1mk} + a_{2j}a_{2k}a_{2m}a_{2q}\mu_{1mkq} + a_{2j}a_{2k}a_{2m}b_{2q}\nu_{1mkq} \\
& \quad + a_{2j}a_{2m}b_{2k}b_{2q}\rho_{1mkq} + a_{2j}b_{2k}b_{2m}b_{2q}\xi_{1mkq} = 0.
\end{aligned}
\quad (B14)
$$

u³ term, second equation:

$$
\begin{aligned}
& a_{6k}(\alpha_{jk} - \alpha_{1k}b_{2j}) - a_{3k}\alpha_{1k}b_{4j} - b_{3k}b_{4j}\beta_{1k} + b_{7j}(-a_{1k}\alpha_{1k} - b_{1k}\beta_{1k}) \\
& \quad + b_{6k}(-b_{2j}\beta_{1k} + \beta_{jk}) + (-a_{1m}a_{3k}b_{2j} - a_{1k}a_{3m}b_{2j} - a_{1k}a_{1m}b_{4j})\delta_{1mk} \\
& \quad + (a_{1m}a_{3k} + a_{1k}a_{3m})\delta_{jmk} \\
& \quad + (-a_{3m}b_{1k}b_{2j} - a_{1m}b_{2j}b_{3k} - a_{1m}b_{1k}b_{4j})\epsilon_{1mk} + (a_{3m}b_{1k} + a_{1m}b_{3k})\epsilon_{jmk} \\
& \quad + (-b_{1m}b_{2j}b_{3k} - b_{1k}b_{2j}b_{3m} - b_{1k}b_{1m}b_{4j})\gamma_{1mk} + (b_{1m}b_{3k} + b_{1k}b_{3m})\gamma_{jmk} \\
& \quad - a_{1k}a_{1m}a_{1q}b_{2j}\mu_{1mkq} + a_{1k}a_{1m}a_{1q}\mu_{jmkq} - a_{1k}a_{1m}b_{1q}b_{2j}\nu_{1mkq} \\
& \quad + a_{1k}a_{1m}b_{1q}\nu_{jmkq} - a_{1m}b_{1k}b_{1q}b_{2j}\rho_{1mkq} + a_{1m}b_{1k}b_{1q}\rho_{jmkq} \\
& \quad - b_{1k}b_{1m}b_{1q}b_{2j}\xi_{1mkq} + b_{1k}b_{1m}b_{1q}\xi_{jmkq} = 0.
\end{aligned}
\quad (B15)
$$

u²v term, second equation:

$$
\begin{aligned}
& a_{7k}(\alpha_{jk} - \alpha_{1k}b_{2j}) - a_{4k}\alpha_{1k}b_{4j} - 2a_{3k}\alpha_{1k}b_{5j} - 3b_{6j} \\
& \quad - b_{4j}b_{4k}\beta_{1k} - 2b_{3k}b_{5j}\beta_{1k} + b_{8j}(-2a_{1k}\alpha_{1k} - 2b_{1k}\beta_{1k}) \\
& \quad + b_{7j}(-a_{2k}\alpha_{1k} - b_{2k}\beta_{1k}) + b_{7k}(-b_{2j}\beta_{1k} + \beta_{jk}) \\
& \quad + (-a_{2m}a_{3k}b_{2j} - a_{2k}a_{3m}b_{2j} - a_{1m}a_{4k}b_{2j} - a_{1k}a_{4m}b_{2j} - a_{1m}a_{2k}b_{4j} - a_{1k}a_{2m}b_{4j} - 2a_{1k}a_{1m}b_{5j})\delta_{1mk} \\
& \quad + (a_{2m}a_{3k} + a_{2k}a_{3m} + a_{1m}a_{4k} + a_{1k}a_{4m})\delta_{jmk} \\
& \quad + (-a_{4m}b_{1k}b_{2j} - a_{3m}b_{2j}b_{2k} - a_{2m}b_{2j}b_{3k} - a_{2m}b_{1k}b_{4j} - a_{1m}b_{2k}b_{4j} - a_{1m}b_{2j}b_{4k} - 2a_{1m}b_{1k}b_{5j})\epsilon_{1mk} \\
& \quad + (a_{4m}b_{1k} + a_{3m}b_{2k} + a_{2m}b_{3k} + a_{1m}b_{4k})\epsilon_{jmk} \\
& \quad + (-b_{2j}b_{2m}b_{3k} - b_{2j}b_{2k}b_{3m} - b_{1m}b_{2k}b_{4j} - b_{1k}b_{2m}b_{4j} - b_{1m}b_{2j}b_{4k} - b_{1k}b_{2j}b_{4m} - 2b_{1k}b_{1m}b_{5j})\gamma_{1mk}
\end{aligned}
$$
$$
\begin{align*}
& + (b_2mb_{3k} + b_{2k}b_{3m} + b_{1m}b_{4k} + b_{1k}b_{4m})\gamma_{jmk} \\
& + (-a_{1m}a_1qa_2kb_{2j} - a_{1k}a_1qa_2mb_{2j} - a_{1k}a_{1m}a_2qb_{2j})\mu_{1mkq} \\
& + (a_{1m}a_1qa_{2k} + a_{1k}a_1qa_{2m} + a_{1k}a_{1m}a_{2q})\mu_{jmkq} \\
& + (-a_{1m}a_2kb_{1qb_{2j}} - a_{1k}a_2mb_{1qb_{2j}} - a_{1k}a_{1mb_{2jb_{2q}}})\nu_{1mkq} \\
& + (a_{1m}a_2kb_{1qb} + a_{1k}a_2mb_{1qb} + a_{1k}a_{1mb_{2q}})\nu_{jmkq} \\
& + (-a_2mb_{1kb_{1qb_{2j}}} - a_{1m}b_{1qb_{2jb_{2k}}} - a_{1m}b_{1kb_{2jb_{2q}}})\rho_{1mkq} \\
& + (a_2mb_{1kb_{1qb}} + a_{1m}b_{1qb_{2k}} + a_{1m}b_{1kb_{2q}})\rho_{jmkq} \\
& + (-b_{1mb_{1qb_{2jb_{2k}}}} - b_{1kb_{1qb_{2jb_{2m}}}} - b_{1kb_{1mb_{2jb_{2q}}}})\xi_{1mkq} \\
& + (b_{1mb_{1qb_{2k}}} + b_{1kb_{1qb_{2m}}} + b_{1kb_{1mb_{2q}}})\xi_{jmkq} = 0. \tag{B16}
\end{align*}
$$

$uv^2$ term, second equation:

$$
\begin{align*}
& a_{8k}(a_{jk} - a_{1k}b_{2j}) - a_{5k}a_{1k}b_{4j} - 2a_{4k}a_{1k}b_{5j} - 2b_{7j} \\
& \quad - 2b_{4k}b_{5j}\beta_{1k} - b_{4k}b_{5k}\beta_{1k} + b_{9j}(-3a_{1k}a_{1k} - 3b_{1k}\beta_{1k}) \\
& \quad + b_{8j}(-2a_{2k}a_{1k} - 2b_{2k}\beta_{1k}) + b_{8k}(-b_{2j}\beta_{1k} + \beta_{jk}) \\
& \quad + (-a_{2m}a_{4k}b_{2j} - a_{2k}a_{4m}b_{2j} - a_{1m}a_{5k}b_{2j} - a_{1k}a_{5m}b_{2j} \\
& \quad - a_{2k}a_{2m}b_{4j} - 2a_{1m}a_{2k}b_{5j} - 2a_{1k}a_{2m}b_{5j})\delta_{1mk} \\
& \quad + (a_{2m}a_{4k} + a_{2k}a_{4m} + a_{1m}a_{5k} + a_{1k}a_{5m})\delta_{jmk} \\
& \quad + (-a_{5m}b_{1kb_{2j}} - a_{4m}b_{2jb_{2k}} - a_{2m}b_{2kb_{4j}} - a_{2m}b_{2jb_{4k}} \\
& \quad - 2a_{2m}b_{1kb_{5j}} - 2a_{1m}b_{2kb_{5j}} - a_{1m}b_{2jb_{5k}})\epsilon_{1mk} \\
& \quad + (a_{5m}b_{1k} + a_{4m}b_{2k} + a_{2m}b_{4k} + a_{1m}b_{5k})\epsilon_{jmk} \\
& \quad + (-b_{2kb_{2mb_{4j}}} - b_{2jb_{2mb_{4k}}} - b_{2jb_{2kb_{4m}}} - 2b_{1mb_{2kb_{5j}}} \\
& \quad - 2b_{1kb_{2mb_{5j}}} - b_{1mb_{2jb_{5k}}} - b_{1kb_{2jb_{5m}}})\gamma_{1mk} \\
& \quad + (b_{2mb_{4k}} + b_{2kb_{4m}} + b_{1mb_{5k}} + b_{1kb_{5m}})\gamma_{jmk} \\
& \quad + (-a_{1qa_2ka_2mb_{2j}} - a_{1ma_2ka_2qb_{2j}} - a_{1ka_2ma_2qb_{2j}})\mu_{1mkq} \\
& \quad + (a_{1qa_2ka_2m} + a_{1ma_2ka_2q} + a_{1ka_2ma_2q})\mu_{jmkq} \\
& \quad + (-a_{2ka_2mb_{1qb_{2j}}} - a_{1ma_2kb_{2jb_{2q}}} - a_{1ka_2mb_{2jb_{2q}}})\nu_{1mkq} \\
& \quad + (a_{2ka_2mb_{1qb}} + a_{1ma_2kb_{2q}} + a_{1ka_2mb_{2q}})\nu_{jmkq} \\
& \quad + (-a_{2mb_{1qb_{2jb_{2k}}}} - a_{2mb_{1kb_{2jb_{2q}}}} - a_{1mb_{2jb_{2kb_{2q}}}})\rho_{1mkq} \\
& \quad + (a_{2mb_{1qb_{2k}}} + a_{2mb_{1kb_{2q}}} + a_{1mb_{2kb_{2q}}})\rho_{jmkq} \\
& \quad + (-b_{1qb_{2jb_{2kb_{2m}}}} - b_{1mb_{2jb_{2kb_{2q}}}} - b_{1kb_{2jb_{2mb_{2q}}}})\xi_{1mkq} \\
& \quad + (b_{1qb_{2kb_{2m}}} + b_{1mb_{2kb_{2q}}} + b_{1kb_{2mb_{2q}}})\xi_{jmkq} = 0. \tag{B17}
\end{align*}
$$
v³ term, second equation:

$$
\begin{align}
& a_{9k}(\alpha_{jk} - \alpha_{1k}b_{2j}) - 2a_{5k}\alpha_{1k}b_{5j} - b_{8j} - 2b_{5j}b_{5k}\beta_{1k} + b_{9j}(-3a_{2k}\alpha_{1k} - 3b_{2k}\beta_{1k}) \nonumber \\
& \quad + b_{9k}(-b_{2j}\beta_{1k} + \beta_{jk}) + (-a_{2m}a_{5k}b_{2j} - a_{2k}a_{5m}b_{2j} - 2a_{2k}a_{2m}b_{5j})\delta_{1mk} \nonumber \\
& \quad + (a_{2m}a_{5k} + a_{2k}a_{5m})\delta_{jmk} + (-a_{5m}b_{2j}b_{2k} - 2a_{2m}b_{2k}b_{5j} - a_{2m}b_{2j}b_{5k})\epsilon_{1mk} \nonumber \\
& \quad + (a_{5m}b_{2k} + a_{2m}b_{5k})\epsilon_{jmk} + (-2b_{2k}b_{2m}b_{5j} - b_{2j}b_{2m}b_{5k} - b_{2j}b_{2k}b_{5m})\gamma_{1mk} \nonumber \\
& \quad + (b_{2m}b_{5k} + b_{2k}b_{5m})\gamma_{jmk} - a_{2k}a_{2m}a_{2q}b_{2j}\mu_{1mkq} + a_{2k}a_{2m}a_{2q}\mu_{jmkq} \nonumber \\
& \quad - a_{2k}a_{2m}b_{2j}b_{2q}\nu_{1mkq} + a_{2k}a_{2m}b_{2q}\nu_{jmkq} - a_{2m}b_{2j}b_{2k}b_{2q}\rho_{1mkq} \nonumber \\
& \quad + a_{2m}b_{2k}b_{2q}\rho_{jmkq} - b_{2j}b_{2k}b_{2m}b_{2q}\xi_{1mkq} + b_{2k}b_{2m}b_{2q}\xi_{jmkq} = 0. \tag{B18}
\end{align}
$$

APPENDIX C: THE INVERSE OF THE MODAL TRANSFORMATION

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

The forward transformation is written as

$$
z_k = M_{1kj}w_j + M_{2kjr}w_j w_r + M_{3kjrm}w_j w_r w_m + \dots, \quad (C1)
$$

and its inverse is expressed as

$$
w_m = Q_{1mn}z_n + Q_{2mno}z_n z_o + Q_{3mnoq}z_n z_o z_q + \dots \quad (C2)
$$

The method in section 4 provides the coefficients in transformation (C1), and the task here
is to derive the coefficients for the inverse (C2) in terms of them. Substitution of equation
(C1) into (C2) and expanding out to order three in the $w_j$'s yields

$$
\begin{equation}
\begin{split}
w_m = {}& Q_{1mn} M_{1nj} w_j + Q_{1mn} M_{2njr} w_j w_r + Q_{1mn} M_{3njra} w_j w_r w_a \\
& + Q_{2mno} M_{1nj} M_{1or} w_j w_r + Q_{2mno} M_{1nj} M_{2ora} w_j w_r w_a \\
& + Q_{2mno} M_{1oj} M_{2nra} w_j w_r w_a + Q_{3mnoq} M_{1nj} M_{1or} M_{1qa} w_j w_r w_a + \dots
\end{split}
\tag{C3}
\end{equation}
$$

Like powers of w are now gathered together on the right side and are equated to the left side. The linear terms yield

$$
Q_{1mn} M_{1nj} = \delta_{mj}, \tag{C4}
$$

the expected result, which says that the linearized inverse must be the inverse of the
linearized transformation. The quadratic terms must satisfy

$$
Q_{2mno} M_{1nj} M_{1or} = -Q_{1mn} M_{2nrj}, \quad (C5)
$$

which can be solved for Q₂ by utilizing the transpose of (C4) and post-multiplying twice
by Q₁ to eliminate the M₁ terms on the left side of (C5). The result is

$$
Q_{2miq} = -Q_{1mn} M_{2njr} Q_{1rq} Q_{1ji}. \quad (C6)
$$
The cubic terms in $w$ must satisfy

$$Q_{3mnoq}M_{1nj}M_{1or}M_{1qa} = -[Q_{1mn}M_{3njra} + Q_{2mno}M_{1nj}M_{2ora} + Q_{2mno}M_{1oj}M_{2nra}], \quad (C7)$$

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

$$Q_{3mich} = -[Q_{1mn}M_{3njra} + Q_{2mno}M_{1nj}M_{2ora} + Q_{2mno}M_{1oj}M_{2nra}][Q_{1ah}Q_{1rc}Q_{1ji}]. \quad (C8)$$

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