# Simultaneous normal form transformation and model-order reduction for systems of coupled nonlinear oscillators

## Abstract

In this paper, we describe a direct normal form decomposition for systems of coupled nonlinear oscillators. We demonstrate how the order of the system can be reduced during this type of normal form transformation process. Two specific examples are considered to demonstrate particular challenges that can occur in this type of analysis. The first is a 2 d.f. system with both quadratic and cubic nonlinearities, where there is no internal resonance, but the nonlinear terms are not necessarily *ε*^{1}-order small. To obtain an accurate solution, the direct normal form expansion is extended to *ε*^{2}-order to capture the nonlinear dynamic behaviour, while simultaneously reducing the order of the system from 2 to 1 d.f. The second example is a thin plate with nonlinearities that are *ε*^{1}-order small, but with an internal resonance in the set of ordinary differential equations used to model the low-frequency vibration response of the system. In this case, we show how a direct normal form transformation can be applied to further reduce the order of the system while simultaneously obtaining the normal form, which is used as a model for the internal resonance. The results are verified by comparison with numerically computed results using a continuation software.

### 1. Introduction

The idea of a normal form transformation was first introduced by Poincaré as a method for reducing a system of equations to their simplest or ‘normal’ form [1]. This method was subsequently developed by other researchers throughout the early twentieth century, and the interested reader can find details of the historical background and analytical development of the theory of normal form transformations in a number of comprehensive texts including [2–7]. A crucial part of the transform process is the identification of ‘resonant’ terms in the sense of terms that can grow without bound when the system has eigenvalues at (or very close to) integer, 1, 2, 3… or fractional (1/2), (1/3)… multiples. These resonant terms cannot be removed during the transformation, and therefore represent an irreducible component of the system. However, using the Poincaré–Dulac theorem it is typically possible to remove all other non-resonant terms during the normal form transformation [3].

Normal form transformations have been used extensively to treat a range of practical applications in physics and engineering, particularly single-degree-of-freedom nonlinear oscillators, and systems of two or more coupled oscillators [4,6]. Relevant to the type of applications considered here, Jezequel & Lamarque [8] proposed a normal form decomposition for a system of two coupled oscillators with cubic nonlinearities and both forcing and damping. Following this, the relationship between the normal form transformation and nonlinear normal modes was established by Touzé and co-workers [9,10], using examples of coupled oscillator systems that included both quadratic and cubic nonlinear terms. In addition to establishing this connection, Touzé's method has the important advantage that it maintains the final result in the form of second-order oscillator equations, and in the process keeps the system eigenvalues in a real form—as a result, this is sometimes called the ‘real normal form’.

Retaining the form of second-order oscillator equations is very useful for vibration-based applications where there are systems of multiple oscillator equations, as for example, in modal analysis [11]. Also following this approach, Neild & Wagg proposed a normal form transform applied directly to second-order oscillator equations [12] which is sometimes called the ‘direct normal form’ in that it is applied directly to second-order oscillator equations without having to transform the equations into an equivalent set of first-order differential equations. Neild and co-workers subsequently used this approach to analyse the resonant behaviour of mechanical systems [13–17] using the established idea of backbone curves [4,18,19] that also provided a link to nonlinear normal modes [9,10,20,21]. This enabled phenomena such as out-of-unison nonlinear normal modes [22] to be captured and recast as bifurcations of the backbone curves [23].

More recently, numerical approximations to the spectral sub-manifolds associated with backbone curves have been developed by Haller and co-workers—see [24] and references therein—which can be used to model the dynamics in the non-internally resonant case. Of particular interest, in a recent study Breunung & Haller [24] carried out a comparison between their spectral sub-manifold method and the methods of Touzé & Amabili [10] and Neild & Wagg [12]. This comparison considered the system we will discuss in §4, and it showed that the *ε*^{1}-order direct normal form gave the incorrect approximations for this example. In fact, as we will show in §4, the *ε*^{2} terms are required to give the correct solutions for systems of this type—a fact that has been previously pointed out in [10]. In addition, a useful study of the upper limit of normal form validity was carried out by Lamarque *et al.* [25], which clarifies the parameter region within which any normal form decomposition can be considered valid.

In terms of reduction of the system order (i.e. size in terms of numbers of degrees-of-freedom), there is a large body of the literature spanning several different domains. From a physical applications perspective, there has been much work in the area of reducing high-dimensional systems such as models of solids and fluids, a few examples from the very large literature on this topic can be found in [26–33], and in some cases, explicit mention of reduction as part of a normal form procedure is discussed [28]. The concept of nonlinear normal modes has been strongly motivated by reducing the size of the system, once again the literature on this topic is large and a selection of examples are given in [20,21,34–36]. For example, Touzé and co-workers demonstrated how the application of the real normal form in this context can be used to reduce the system dimension of shell structures [37,38].

In this paper, we present a modification to the direct normal form transformation for systems of nonlinearly coupled oscillators, based on the method of [12], which also includes systematic model order reduction as part of the process. The paper is structured as follows. In §2, we first describe the normal form transform directly applied to the general forced and damped systems governed by *N* nonlinearly coupled second-order ordinary differential equations. Following this, a variant of the method is introduced which includes the ability to reduce the order of the model. Then in §4, the proposed technique is used to transform and reduce a 2 d.f. discrete example system of quadratic and cubic nonlinearities to an equivalent single-degree-of-freedom normal form model. The resonant curves and backbone curves of the example system are then found and compared with numerically computed results, in order to validate the proposed technique up to order *ε*^{2} and also resolve the issues raised by Breunung & Haller [24]. Then in §5, the application of the new method to a thin plate system is considered in order to demonstrate the ability of this method to work for (i) systems with many degrees-of-freedom, and (ii) systems that have internal resonance. Conclusions are drawn in §6.

### 2. Direct normal form method

The direct normal form method uses the same philosophy as the classical Poincaré–Dulac theorem [3], but without initially transforming the system of *N* second-order differential equations into a system of 2*N* equivalent first-order differential equations, (thereby avoiding doubling the size of the system being considered). Instead, the formulation is carried out directly on the second-order differential equations, and incorporates the benefits of (i) not increasing the system dimension to 2N, and (ii) keeping the eigenvalues analysis real instead of complex, as for the method of Touzé & Amabili [10].

There are several other differences with these existing methods that are notable. First to enable a compact formulation, an *ε* parameter, which can be assumed to be small, is used so that order *ε*^{γ}, for *γ* > 1, terms can be neglected as required. However, the method is entirely consistent with the classical method derived without *ε* and using a power series approach [3], and the interested reader can find this formulation of the direct normal form in [15]. Note also that, as a result of the compact formulation the *ε* order does not correspond directly to the order of the nonlinear term. This is because for a wide class of coupled oscillator systems including just the *ε*^{1} terms will provide sufficient levels of accuracy, although a counter example to this is given in §4.

Furthermore, the type of detuning used in the direct normal form is different from that typically used in either classical normal form methods or multiple scales methods—see for example [4] for the traditional application of detuning. In fact, the detuning used in the direct normal form can be interpreted as giving an increased accuracy at a lower *ε* order than the traditional detuning, because the *ε* orders are different, although we note that the same accuracy can be achieved in all methods if sufficient terms in the expansion are used. A detailed explanation of how this mechanism works is given in [39], for both the direct normal form and multiple scales method. Note also that, rather than using this type of detuning approach, it is possible to use the direct normal form with the detuning parameter treated as a (small) *unfolding* parameter—as shown in [15]—where additional unfolding parameters representing damping and forcing were also included.

Here, we consider the nonlinear oscillating system, whose dynamic behaviour could be described using the equations of motion given by

**x**is the {

*N*× 1} vector of physical displacements, i.e. $\mathbf{x}={[{x}_{1},\hspace{0.17em}\dots ,\hspace{0.17em}{x}_{N}]}^{\u22ba}$,

**M**and

**K**are the {

*N*×

*N*} mass and stiffness matrices respectively,

**N**

_{x}is an {

*N*× 1} vector of linear damping terms and nonlinear terms related to stiffness, damping and external forcing,

*ε*is a dimensionless book-keeping parameter that indicates which terms are small,

*N*is the number of degrees-of-freedom of the system, and an over-dot represents the differentiation with respect to time

*t*. The external force is expressed as the product of an {

*N*× 2} forcing amplitude matrix

**P**

_{x}and a {2 × 1} forcing frequency related vector

**r**, written as

*Ω*is the frequency of the external excitation and i is the imaginary unit. Here, the analysis is restricted to the case where the nonlinear terms in

**N**

_{x}are polynomial functions of

**x**, $\dot{\mathbf{x}}$ and

**r**only, other cases are detailed in [40].

The first step in dealing with the problem is to carry out a transform which can decouple the linear part of equation (2.1). To do this, the vibration modes obtained from the underlying unforced and undamped linear system defined by setting the damping, forcing and nonlinear terms in equation (2.1) to zero, i.e. $\mathbf{M}\ddot{\mathbf{x}}+\mathbf{K}\mathbf{x}=\mathbf{0}$, are used. Then the linear transform can be achieved by transforming equation (2.1) using the substitution **x** = *Φ***q**, where **q** is an {*N* × 1} vector of modal coordinates, i.e. $\mathbf{q}={[{q}_{1},\dots ,{q}_{N}]}^{\u22ba}$, and ** Φ** is the {

*N*×

*N*} matrix of linear modeshapes found from the eigenvalue problem

*Φ***=**

*Λ***M**

^{−1}

**K**

**, and**

*Φ***is a {**

*Λ**N*×

*N*} diagonal matrix containing the square of the natural frequencies of the vibration modes, i.e.

*ω*

^{2}

_{ni}for

*i*= 1, 2, ….,

*N*. The linear modes are scaled so that the largest element in each mode vector is ||1||.

Making the substitution **x** = *Φ***q** into equation (2.1) and pre-multiplying by ${\mathit{\Phi}}^{\u22ba}$ gives

**M**

^{−1}

**K**

**=**

*Φ*

*Φ***, leading to**

*Λ***N**

_{q}(

**q**) is the {

*N*× 1} vector of nonlinear terms in modal coordinates, derived from

The next step is to remove non-resonant forcing terms via a transform written as

**e**is an {

*N*× 2} force transform matrix. Full details of this forcing transformation process can be found in [41], but for completeness, the main steps are outlined here. Substituting equation (2.6) into equation (2.4) gives

**W**is a {2 × 2} diagonal matrix with leading elements +i

*Ω*and −i

*Ω*. After the forcing transform, equation (2.7) is expected to be arranged as

where **N**_{v} is an {*N* × 1} vector of nonlinear and damping terms in **v** and **P**_{v} is an {*N* × 2} matrix of near-resonant forcing amplitudes in which the near-resonant one may be retained and non-resonant one has been removed. Combining equation (2.7) and equation (2.8) leads to

For a nonlinear system, when the forcing frequency is commensurably close to one of the system natural frequencies, a significant response would occur, related to the corresponding mode, regarded as near-resonant. Therefore, the matrix of the near-resonant forcing amplitudes may be set as

**P**

_{v,i}and

**P**

_{q,i}denote the

*i*th row of

**P**

_{v}and

**P**

_{q}, respectively. Now making use of equation (2.11) and equation (2.10), the terms in the force transform matrix,

**e**, can be determined using,

**e**

_{i}is the

*i*th row of

**e**. The detailed derivation process is given in [41]. Once the matrix

**e**is determined, the transformed linear damping and nonlinear terms vector

**N**

_{v}can be computed using equation (2.9).

Now a nonlinear transform is applied to simplify equation (2.8) by cancelling the non-resonant terms, which then leads to the normal form of the system. This nonlinear transform is often referred to as a near-identity transform [41] and may be expressed as

**u**and $\mathbf{H}(\mathbf{u},\dot{\mathbf{u}},\mathbf{r})$ are the {

*N*× 1} vectors of fundamental (resonant) and non-resonant components of

**v**respectively. Here, the non-resonant components are assumed to be small when compared with the resonant responses, thus $\mathbf{H}(\mathbf{u},\dot{\mathbf{u}},\mathbf{r})$ is denoted to be of order

*ε*

^{1}. The resulting normal form after the near-identity transform has been applied may be expressed as

*i*th state,

*u*

_{i}where

*i*= 1, …,

*N*, responds at its corresponding resonant response frequency,

*ω*

_{ri}, where

*ω*

_{ri}is close to, but not equal to

*ω*

_{ni}in general.

Here, the $\mathbf{H}(\mathbf{u},\dot{\mathbf{u}},\mathbf{r})$ vector in equation (2.13) and the ${\mathbf{N}}_{u}(\mathbf{u},\dot{\mathbf{u}},\mathbf{r})$ vector in equation (2.14) are approximated by a series of functions of reducing levels of significance, i.e.

*a*

*b*

*N*× 1} vectors in the same form as $\mathbf{H}(\mathbf{u},\dot{\mathbf{u}},\mathbf{r})$ and ${\mathbf{N}}_{u}(\mathbf{u},\dot{\mathbf{u}},\mathbf{r})$, respectively. Additionally, a Taylor series expansion is applied to the nonlinear term vector in equation (2.8), such that after the substitution of equation (2.13), about the equilibrium point $[\mathbf{v},\dot{\mathbf{v}}]=[\mathbf{u},\dot{\mathbf{u}}]$, we obtain

*ω*

_{ri}values are close to, but not equal to

*ω*

_{ni}in general, a frequency detuning is introduced, to define this relationship

**is a {**

*Υ**N*×

*N*} diagonal matrix of the square of the response frequencies,

*ω*

^{2}

_{ri}, and

**Δ**= (

**−**

*Λ***)/**

*Υ**ε*represents the frequency detuning values which are assumed to be small, i.e. of order

*ε*

^{1}. The purpose of this detuning parameter is to represent the amplitude dependent nature of the response frequencies,

*ω*

_{ri}, of the system.

Now making the substitution of equations (2.13), (2.15a) and (2.16) into equation (2.8) and eliminating the similar terms using equation (2.14) with the substitution of equation (2.15b) gives

*ε*

^{2}are included. Although in principle higher order

*ε*may also be derived, typically, it is assumed that these terms are vanishingly small. Balancing the terms associated with different orders

*ε*in equation (2.18) gives a set of homological equations, written as

*a*

*b*

*c*

*a*

*b*

**n**

_{(2)}is dependent on

**h**

_{(1)}which in turn is determined by

**n**

_{(1)}.

To satisfy equation (2.19a), it can be obtained that

**u**, i.e.

*u*

_{i}, responds at a single response frequency

*ω*

_{ri}, and its trial solution may be assumed to be

where *U*_{i}, *ϕ*_{i} and *ω*_{ri} are the displacement amplitude, phase lag and response frequency of *u*_{i}, respectively. Using the left-hand side of equation (2.22), a vector **u***_{(j)} (of length ${\mathcal{L}}_{j}$) containing all the polynomial combinations of *u*_{ip} and *u*_{im} present in **n**_{(j)} may be defined, and then **n**_{(j)}, **h**_{(j)} and **n**_{u(j)} are expressed as,

*n*

_{(j)}], [

*h*

_{(j)}] and [

*n*

_{u(j)}] are $\{N\times {\mathcal{L}}_{j}\}$ matrices of corresponding coefficients. Substituting equation (2.23) into the

*j*th homological equation gives

**u***

_{(j)}are functions of

*u*

_{ip}and

*u*

_{im}, which may be written as

where $\ell =1,\dots ,{\mathcal{L}}_{j}$, and *s*_{pj,ℓ,i}, *s*_{mj,ℓ,i}, *m*_{pj,ℓ} and *m*_{mj,ℓ} are exponents of *u*_{ip}, *u*_{im}, *r*_{p} and *r*_{m} in the ℓth element of **u***_{(j)}, respectively. Using the right-hand side of equation (2.22) and equation (2.25), the second time derivative of *u**_{(j)ℓ} can be derived to be

*ω**

_{(j)ℓ}is the response frequency of

*u**

_{(j)ℓ}, i.e.

**u***

_{(j)}may be written as

*ω*

^{2}

_{(j)ℓ}. So, finally equation (2.28) is substituted into equation (2.24) to obtain

*Hadamard*product operator and

*β*_{(j)}is a $\{N\times {\mathcal{L}}_{j}\}$ coefficients matrix, whose {

*i*, ℓ}

^{th}element is

**n**

_{u(j)}is expected to be as simple as possible, all terms in it are set to zero, except where there is a zero in

*β*_{(j)}, for which case [

*n*

_{u(j)}]

_{i, ℓ}= [

*n*

_{(j)}]

_{i, ℓ}in order for equation (2.29) to be satisfied. The criteria for the determination of the resonant terms and the computation of the coefficients of the non-resonant terms is given as

*a*

and

*b*

Once the coefficient matrices [*n*_{u(j)}] are determined to a specific level of accuracy, i.e. *ε*^{j}, the vector of resonant nonlinear terms can be formed using equation (2.15b). Then equation (2.14) can be rearranged, with the substitution of the trial solution equation (2.22), to give

*χ*

_{n}and ${\stackrel{~}{\chi}}_{n}$ are time-independent complex conjugates. By inspection of equation (2.32), a solution can be found by setting

*χ*

_{i}= 0 (or equivalently ${\stackrel{~}{\chi}}_{i}=0$) which then finally leads to the expressions describing the relationship between the fundamental response amplitudes and frequencies.

### 3. An order-reduced variant of the nonlinear transform

As described in §1, previous authors have considered reduction of system size via a transform—for example [37,38]. For the direct normal form method, an essential step within the nonlinear transform is to construct the vector containing all relevant combinations of polynomial nonlinear terms, **u***_{(j)}, and its corresponding coefficient matrices, [*n*_{(j)}]. The elements in **u***_{(j)} are defined by the structure of the nonlinear terms in the original problem, e.g. the number of states originally in **N**_{x}, but then transformed into **N**_{q}, **N**_{v} and **n**_{(j)}. In fact, the number of polynomial terms in **u***_{(j)} will be exponentially increasing with the degrees-of-freedom of the considered systems. For example, considering nonlinear systems of quadratic and cubic nonlinearities, ${\mathcal{L}}_{1}={\mathrm{C}}_{2N+2}^{1}+({\mathrm{C}}_{2N+2}^{1}+{\mathrm{C}}_{2N+2}^{1})+({\mathrm{C}}_{2N+2}^{1}+{\mathrm{C}}_{2N+2}^{2}{\mathrm{C}}_{2}^{1}+{\mathrm{C}}_{2N+2}^{3})$ for an *N*-degree-of-freedom system. Correspondingly, the growth of the system size would significantly complicate the determination of the nonlinear coefficient matrices, [*n*_{(j)}], the resonance criteria coefficient matrices, *β*_{(j)} and the non-resonant coefficient matrices, [*h*_{(j)}]. This issue will be compounded when a higher-level of accuracy, e.g. *ε*^{2}, is required to be considered. For example, to find **n**_{(2)}, large matrix manipulations may need to be performed if ${\mathcal{L}}_{j}$ is large, which is typically the case via equation (2.20b).

One approach is to use symbolic manipulation methods to solve such problems, but here we consider another natural solution which is to reduce the size of **u***_{(j)} without significantly compromising the accuracy. This is based on the observation that for many applications, the majority of the **u***_{(j)} are non-responding and their associated coefficients are zero. For example, for nonlinear systems under a harmonic excitation, only a small number of eigenmodes of nearly commensurable natural frequencies to the external forcing frequency would respond significantly. Based on this, a sparse resonant state vector **u**_{sp} is introduced whose entries are set to be zero except for *R* specifically selected terms where *R* < *N*. These selected non-zero-response states will be either the externally resonant ones for the external forcing situation or the internally resonant ones for the free vibration.

When making the substitution of **u**_{sp}, equation (2.20) becomes

*a*

and

*b*

*R*× 1} vector containing the response of the

*R*dominant states only. Now due to the reduction of the number of terms in $\hat{\mathbf{u}}$ relative to that of the original state vector,

**u**, the basis vector,

**u***

_{(j)}, and coefficients matrices, [

*n*

_{(j)}], [

*n*

_{u,(j)}], [

*h*

_{(j)}] will be simplified significantly.

Furthermore, as a number of the states are now set to zero in advance, the final normal form of, equation (2.14), will become

*R*resonant equations and, then the number of time-independent equations related to response amplitudes and frequencies also shrinks to

*R*. Here ${\hat{\mathbf{P}}}_{u}$ is the {

*R*× 2} matrix of amplitude of the effective modal forces. Therefore, this variant of the nonlinear transform also performs a model order reduction.

### 4. Example 1: a 2 d.f. discrete mass-spring system

The first example system considered is shown in figure 1. This system has been previously studied in [9,24,42], and for the purpose of making a comparison we follow the format introduced by these previous authors. The system consists of a mass supported by a vertical and a horizontal spring. The equation of motion of this example system, with *L* = 1, is given by

*f*

_{i}cos(

*Ωt*), for

*i*= 1, 2 have been added after the conservative equations (i.e. with

*ζ*

_{i}=

*f*

_{i}= 0) are derived from figure 1, [9]. For the purpose of applying the direct normal form methods, the terms in equation (4.2) are assumed to be the order

*ε*

^{1}terms. It should be noted that increasing the

*ζ*

_{i}damping values can lead to increased complexity in the dynamic behaviour, as shown Touzé & Amabili [10].

As the conservative form of equation (4.1) is naturally linearly decoupled, it can be described in the matrix form of equation (2.4) by setting $\mathbf{q}={[{q}_{1},\hspace{0.17em}{q}_{2}]}^{\u22ba}={[{x}_{1},\hspace{0.17em}{x}_{2}]}^{\u22ba}$, where

*a*

*b*

Here, a non-internal-resonant case is considered, i.e. *ω*_{r1}≠*nω*_{r2} and *ω*_{r2}≠*nω*_{r1}, where *n* = 1, 2, …. We assume that the system is only forced in the horizontal direction, i.e. *f*_{2} = 0, thus only the horizontal response is significant. Therefore, the vector of reduced state vector is $\hat{\mathbf{u}}=[{u}_{1}]$. Substituting $\mathbf{q}={[{q}_{1},\hspace{0.17em}{q}_{2}]}^{\u22ba}={[{u}_{1p}+{u}_{1m},\hspace{0.17em}0]}^{\u22ba}$ into equation (4.3b), we obtain the vector of nonlinear terms of *ε*^{1} significance as

Expressing ${\mathbf{n}}_{(1)}(\hat{\mathbf{u}})$ in the matrix form, i.e. **n**_{(1)} = [*n*_{(1)}]**u***_{(1)}, gives

**u***

_{(1)}has been reduced significantly. Based on equation (4.5),

*β*_{(1)}and [

*h*

_{(1)}] are found, using equations (2.30) and (2), to be

*ω*

_{r2}=

*ω*

_{2}is used as

*u*

_{2}= 0 is assumed. Note that the coefficients of resonant terms are identified by having boxes around them in [

*n*

_{(1)}], equation (4.5). Using equation (2.23), the vectors of resonant and non-resonant terms of

*ε*

^{1}significance are obtained as

*a*

*b*

*ε*

^{2}must be considered for capturing the nonlinear behaviour, as will be explained below.

Therefore, the nonlinear terms of *ε*^{2} significance for this example are now determined. Considering equation (3.1b), **Δ**, ${\dot{\mathbf{h}}}_{(1)}(\hat{\mathbf{u}})$, $(\mathrm{\partial}{\mathbf{N}}_{v}/\mathrm{\partial}\mathbf{v})|{\hspace{0.17em}}_{{\mathbf{u}}_{sp},{\dot{\mathbf{u}}}_{sp}}$, and $(\mathrm{\partial}{\mathbf{N}}_{v}/\mathrm{\partial}\dot{\mathbf{v}})|{\hspace{0.17em}}_{{\mathbf{u}}_{sp},{\dot{\mathbf{u}}}_{sp}}$ are required to be derived first, which in this example gives

*a*

*b*

*c*

*d*

*u*

_{ip}

*u*

_{im}will go to zero when the assumed solution equation (2.22) is substituted, and so to maintain the simplest form of the equations, we eliminate them at this stage. Now it can be determined that,

*a*

*b*

It is noteworthy that including up to cubic nonlinear terms should provide an accurate solution for this example, thus the nonlinear terms vector **n**_{(2)} is truncated at $\mathcal{O}({\hat{\mathbf{u}}}^{4})$. This assumption will be verified by the simulation results shown below. By comparing equations (4.7b) and (4.9a), it can be determined that **u***_{(2)} = **u***_{(1)}. Now, projecting ${\mathbf{n}}_{(2)}=\mathbf{\Delta}{\mathbf{h}}_{(1)}+(\mathrm{\partial}{\mathbf{N}}_{v}/\mathrm{\partial}\mathbf{v})|{\hspace{0.17em}}_{\mathbf{u},\dot{\mathbf{u}}}{\mathbf{h}}_{(1)}+(\mathrm{\partial}{\mathbf{N}}_{v}/\mathrm{\partial}\dot{\mathbf{v}})|{\hspace{0.17em}}_{\mathbf{u},\dot{\mathbf{u}}}{\dot{\mathbf{h}}}_{(1)}$ in terms of **u***_{(2)}, the coefficient matrix [*n*_{(2)}] is found to be,

*a*

*b*

*n*

_{(2)}]

_{1}and [

*n*

_{(2)}]

_{2}are the first and second rows of [

*n*

_{(2)}], respectively. Similarly, the coefficients of non-resonant terms of

*ε*

^{2}significance are computed, leading to

*h*

_{(2)}]

_{1}and [

*h*

_{(2)}]

_{2}are the first and second rows of [

*h*

_{(2)}], respectively. Then, the vector representing the

*ε*

^{2}contribution is found to be

*ε*

^{2}part is,

**N**

_{u}=

**n**

_{u(1)}+

**n**

_{u(2)}, the reduced resonant equation of motion of the first resonance is

*u*

_{1p}= ((

*U*

_{1}/2)e

^{−iϕ1})e

^{iωr1t}and

*u*

_{1m}= ((

*U*

_{1}/2)e

^{iϕ1})e

^{−iωr1t}into equation (4.14) and arranging the result as equation (2.32) gives

*χ*

_{1}= 0 and eliminating the phase-related term gives

*ζ*

_{i}= 0 and

*f*

_{1}= 0, gives the backbone curve expression of the system first mode to be

*et al*. [9]. Once the fundamental response amplitudes and frequencies are obtained by solving equation (4.16) or equation (4.17), the physical displacement response may be computed using the corresponding reverse transform described in §2, so that

The reason for needing the *ε*^{2} terms in the direct normal form, is based on the following. First, the quadratic terms of the type found in this example generate cubic terms in the nonlinear normal mode decomposition, as pointed out by Touzé & Amabili [10]. However, in the direct normal form method, these terms are captured in the *ε*^{2} expansion not the *ε*^{1}, and most previous applications of the method were restricted to the *ε*^{1} version. So for example, when the comparison was carried out by Breunung & Haller [24], the *ε*^{1}-order approximation was used, and as a result, the softening nonlinear effects were not captured appropriately. As we show in the results below, (and for the full system in [43]) this is corrected by the addition of the *ε*^{2} terms. It should be noted that the direct normal form method relies on the nonlinear terms being small in the sense that they should be significantly smaller than the *ω*^{2}_{ni} values. However, in this example, the nonlinear coefficients are of the same order as the *ω*^{2}_{ni} values, and yet despite this, it is not the determining factor for the discrepancy observed by Breunung & Haller [24]—see [43] for further details.

Now specific values are chosen to compute the analytical forced response curves defined by, equation (4.16), and backbone curves, equation (4.17), of the example mass-spring system. In order to verify the analytical results, both the backbone curves and the resonant response curves are computed using the continuation Matlab toolbox, COCO [44,45].

The simulation uses the same undamped parameters as [9,24], while the forcing and damping parameters are slightly different—as will be explained below. Specifically, two sets of damping and forcing parameters are employed, i.e. *ζ*_{1} = *ζ*_{2} = 0.001, *f*_{1} = [0.001, 0.002], *f*_{2} = 0 and *ζ*_{1} = *ζ*_{2} = 0.01, *f*_{1} = [0.01, 0.02], *f*_{2} = 0 and the results are shown in figure 2. The plots are presented in the projection of the response amplitude of the physical coordinates, |*x*_{i}|, against the external forcing frequency, *Ω*. In figure 2, panels (*a*,*c*) show the |*x*_{1}| response, and panels (*b*,*d*) show the |*x*_{2}|. In each panel, there are two forcing cases shown, and the small damping (*ζ*_{1} = *ζ*_{2} = 0.001) case is shown in panels (*a*,*b*) with the higher damping case (*ζ*_{1} = *ζ*_{2} = 0.01) shown in panels (*c*,*d*). It can be seen that the order *ε*^{2} (reduced-order) backbone curves correctly predict the softening dynamics of the example system which is consistent with the findings in [9,10] and [24]. Furthermore, the backbone curves and forced-damped response curves computed using COCO are also shown in figure 2. It can be seen that the analytical and numerical curves are in very close agreement, particularly for the smaller forcing cases. As would be expected, the level of agreement is less good for the higher forcing curves, and in general, agreement is better when amplitudes of excitation are smaller. Specifically, it can be seen that as amplitude increases, there is the increasing difference between the analytically computed curves and the numerically computed curves. This is due to the truncation of the analytical approximation at order *ε*^{2}.

It is noteworthy that as the horizontal mode is assumed to be the dominant one here, and as a result, the vertical response, |*x*_{2}|, is only made up of contributions from the non-resonant terms (see equation (4.18)). This implies that the reduced-order normal form technique is also able to correctly estimate the response amplitudes of the non-resonant terms of the forced response of the system. The forcing parameters were chosen to be similar to those used by Touzé *et al*. [9], except here we are primarily interested in how accurate the undamped backbone curves are. Therefore, we have reduced *ζ*_{2} to the same (small) level as *ζ*_{1}, whereas in [10] the authors discussed the coupling effect of when *ζ*_{2} was at a range of different levels.

It is important to remember this type of model order reduction will be limited by the fact that a specific part of the system has been neglected. One instance, in this example, is that the backbone curve expression in equation (4.17) is an implicit function of *ω*_{r1}, but with both modes included, it would be an implicit function of both *ω*_{r1} and *ω*_{r2}. Given that we would expect *ω*_{r2} to vary with amplitude, using the assumption that *ω*_{r2} = *ω*_{2} (as in equation (4.17)) is an obvious limitation of the reduced-order approach.

### 5. Example 2: a continuous plate system

A pinned-pinned rectangular thin plate of large-amplitude vibration is now considered, as shown schematically in figure 3. Unlike the previous example, in this case, we focus on the fact that the system has an internal resonance. As a result, we have selected the example to have only cubic nonlinear terms and hence it only requires the *ε*^{1} order normal form transformation to find the backbone curves. In fact, the case of 1 : 1 internal resonance with both quadratic and cubic nonlinearity, has previously been investigated in [10,37,38].

This example has already been studied in detail in [46], and here it is used to consider the application of the reduced-order-direct normal form technique for a system with internal resonance. Following the approach in [46], the nonlinear partial differential equations governing the dynamic behaviour of the plate of large-amplitude responses can be projected onto an orthogonal set of linear eigenmodes, from which a set of ordinary differential equations can be obtained of the form

*i*= 1, 2, …,

*N*. The detailed model derivation process is given, for example, in [41]. Due to the geometry and boundary conditions of the plate, only cubic nonlinearity terms appear in equation (5.1). In this paper, we consider the same plate investigated in [46]. The geometrical and material properties of the plate and the values of the natural frequencies and nonlinear terms coefficients are given in appendix A.

For this example plate, only its backbone curves are investigated, and also, it is assumed that the system is excited at a frequency in the vicinity of its second (and third) natural frequency. Considering just the first four eigenfrequencies (*N* = 4), there may exist a 1 : 1 resonant interaction between the second and third mode of the plate system because *ω*_{n2} ≈ *ω*_{n3}. Therefore, it may not be possible to discern which of these two modes will be the dominant one, in order to construct the sparse vector of resonant states. Here, we may set *u*_{2}≠0 and *u*_{3}≠0 to represent the situation that an initial motion is initiated in either the second or third mode.

Following the approach in [46], and using $\hat{\mathbf{u}}={[{u}_{2},{u}_{3}]}^{\u22ba}$, which reduces the system from a four-mode to a two-mode approximation, the vector of nonlinear terms is found to be

**n**

_{(1)}associated with the second and third modes. Note that due to the specific value of the nonlinear coefficients, see table 2, ${{\mathbf{n}}_{(1)}(\hat{\mathbf{u}})}_{[1]}={{\mathbf{n}}_{(1)}(\hat{\mathbf{u}})}_{[4]}=0$ thus they are omitted. From equation (5.2), the vector of polynomials terms can be defined and its associated

**matrix is computed as**

*β*where *ω*_{r} = *ω*_{r2} = *ω*_{r3} is used as the 1 : 1 internal resonance is assumed. The coefficient matrix of **n**_{(1)} is found in which the terms corresponding to the resonant ones are identified with boxes and then the coefficients of non-resonant terms are computed, as

*a*

*b*

*U*

_{3}= 0 and

*U*

_{2}= 0 results in the expressions of two single-mode (i.e. non-internally resonant) backbone curves given by

*a*

*b*

^{i(|ϕ3−ϕ2)|}= 1, the phase relationships

*ϕ*

_{3}−

*ϕ*

_{2}= 0 or

*π*, are assumed, which means there are two double-mode backbone curves with identical response frequencies and amplitudes but of different modal phase differences denoted by the superscript of

*D*

^{±}

_{23}. See also [22] for a comparison of these results.

Figure 4 shows the analytically approximated backbone curve results of the example plate using the coefficient values in table 2. Again, the results from numerical continuation serve as the benchmark solution. The plate is assumed to be forced in the third (linear) mode-shape only. The frequency is set to be close to the second/third eigenfrequency, and the forcing amplitude is specifically chosen to make the resonant amplitude close to the thickness of the plate, i.e. *f*_{3} = 0.02, which will ensure that the large-amplitude response assumption is satisfied.

Figure 4 shows that the double-mode backbone curves, *D*^{±}_{23}, bifurcate from the single-mode backbone curve, *S*_{3}. Note that as the results are projected in the modal coordinates, the two double-mode backbone curves overlap (as opposed to being in the *u*_{2}, *u*_{3} coordinate space). As we are forcing the system only in mode 3, the *q*_{3} response is driven by the external forcing, and has a hardening response curve (figure 4*b*) as would be expected for a system with cubic nonlinearity. The *q*_{2} response (figure 4*a*) is only triggered after the *D*^{±}_{23} curves bifurcate from the *S*_{3} curve at an amplitude of approximately *q*_{3} = 2.9 × 10^{−4}. In contrast with the previous example, here the amplitudes are 10^{−3} smaller, and as a result, the analytically computed backbone curve matches very accurately with the numerical solution, for both the single- and double-mode responses.

### 6. Conclusion

In this paper, we introduced a technique with the ability to perform a simultaneous direct normal form transform and model order reduction for nonlinear dynamic problems. This method is derived via its application to a general *N*-degree-of-freedom system of coupled oscillators with polynomial type nonlinearities. The main advantages of this new method are that (i) it is directly applicable to systems of second-order ordinary differential equations which are the natural description form of mechanical structures and (ii) it significantly simplifies the derivation process of the normal form transformation when applied to larger-scale systems. However, it should also be noted that model order reduction methods of this type will, by definition, neglect part of the system behaviour (often called the residual), and therefore care should be exercised to ensure that all the required dynamic behaviour is appropriately captured by the reduced-order model. In our case, this was done by verifying the model using numerical continuation simulations of the forced damped system.

Following on from the general derivations, we demonstrated the performance of the proposed technique using two examples. The first example was a 2 d.f. system that has both quadratic and cubic nonlinearities, but no resonant interactions. Due to the fact that the direct normal form method captures cubic terms generated by the presence of quadratic nonlinearities at order *ε*^{2}, the application of the normal form transformation has to be performed at least order *ε*^{2} in order to capture the key dynamic behaviour. However, as we demonstrated, this can be done while simultaneously reducing the system from 2 to 1 d.f. As a result, all the intermediate matrices constructed during the nonlinear transform and the final resulting equation of motion were significantly simplified without sacrificing the accuracy of the solution. The approximated backbone curves results were then shown to be able to accurately predict the softening behaviour and the resonant response of the example system when harmonically forced, by comparison with numerical continuation simulations of the forced damped system.

Finally, the application of the new method to the example of a thin, rectangular, simply supported plate was investigated. Due to the specific geometry and boundary conditions of the plate, there exists a 1 : 1 internal resonance between the second and third modes of the system. In previous work, a four-mode reduced-order model was used to explain this behaviour. However, by using the reduced-order-direct normal form, only the two interacting modes were required to capture the key dynamic behaviour. Accordingly, the dynamics of the plate system were reduced to a system of two resonantly interacting oscillators. The associated numerical continuation results show a high degree of consistency with the analytically obtained backbone curve solutions for both single-mode and resonant double-mode responses.

### Data accessibility

The datasets supporting this article have been uploaded as part of the electronic supplementary material.

### Authors' contributions

Both authors contributed equally to this work.

### Competing interests

We declare we have no competing interests.

### Funding

The authors acknowledge the support of the Engineering and Physical Science Research Council. This work was started while D.J.W. was supported by EPSRC grant no. EP/K003836/2, and finished under the grant no. EP/R006768/1. X.L. was supported by a Department of Mechanical Engineering studentship and also EPSRC grant no. EP/J016942/1.

## Acknowledgments

The authors also thank Ayman Nasir for assistance in proof reading the manuscript.

## Appendix A. The parameters of the example plate

The geometrical and material properties of the example plate considered in §5 are shown in table 1. Note that SI units are used throughout this example. The natural frequencies and coefficients of nonlinear terms related to the first four mode of the example plate are shown table 2.

length | width | thickness | density | Young's modulus | Poisson's ratio |
---|---|---|---|---|---|

0.38 m | 0.48 m | 0.0005 m | 2700 kg m^{−3} |
70 GPa | 0.31 |

ω_{n1} = 58.9 |
α^{[1]}_{111} = 5.45 × 10^{9} |
ω_{n2} = 143.9 |
α^{[2]}_{112} = 2.36 × 10^{10} |

α^{[1]}_{122} = 2.36 × 10^{10} |
α^{[2]}_{222} = 3.14 × 10^{10} |
||

α^{[1]}_{133} = 2.27 × 10^{10} |
α^{[2]}_{233} = 6.51 × 10^{10} |
||

ζ_{1} = 0.01 |
α^{[1]}_{144} = 2.44 × 10^{10} |
ζ_{2} = 0.01 |
α^{[2]}_{244} = 1.24 × 10^{11} |

α^{[1]}_{234} = 7.43 × 10^{10} |
α^{[2]}_{134} = 7.43 × 10^{10} |
||

ω_{n3} = 150.8 |
α^{[3]}_{113} = 2.27 × 10^{10} |
ω_{n4} = 235.8 |
α^{[4]}_{114} = 2.44 × 10^{10} |

α^{[3]}_{223} = 6.51 × 10^{10} |
α^{[4]}_{224} = 1.24 × 10^{11} |
||

α^{[3]}_{333} = 3.14 × 10^{10} |
α^{[4]}_{334} = 1.32 × 10^{11} |
||

ζ_{3} = 0.01 |
α^{[3]}_{344} = 1.32 × 10^{12} |
ζ_{4} = 0.01 |
α^{[4]}_{444} = 5.58 × 10^{10} |

α^{[3]}_{124} = 7.43 × 10^{10} |
α^{[4]}_{123} = 7.43 × 10^{10} |

### References

- 1.
Poincaré H (translated by D. Goroff). 1993 New methods of celestial mechanics (English translation of:*Les méthodes nouvelles do la méchanique celeste*, 1892). College Park, MD: American Institute of Physics, USA. Google Scholar - 2.
Guckenheimer J, Holmes P . 1983**Nonlinear oscillations, dynamical systems, and bifurcations of vector fields**. New York, NY: Springer. Crossref, Google Scholar - 3.
Arnold Vi . 1988**Geometrical methods in the theory of ordinary differential equations**. Berlin, Germany: Springer. Google Scholar - 4.
- 5.
Chow SN, Li C, Wang D . 1994**Normal forms and bifurcation of planar vector fields**. Cambridge, UK: Cambridge University Press. Crossref, Google Scholar - 6.
Kahn PB, Zarmi Y . 2000 Nonlinear dynamics: a tutorial on the method of normal forms.**Am. J. Phys.**, 907–919. (doi:10.1119/1.1285895) Crossref, Web of Science, Google Scholar**68** - 7.
Murdock J . 2002**Normal forms and unfoldings for local dynamical systems**. Berlin, Germany: Springer. Google Scholar - 8.
Jezequel L, Lamarque CH . 1991 Analysis of nonlinear dynamic systems by the normal form theory.**J. Sound Vib.**, 429–459. (doi:10.1016/0022-460X(91)90446-Q) Crossref, Web of Science, Google Scholar**149** - 9.
Touzé C, Thomas O, Chaigne A . 2004 Hardening/softening behaviour in non-linear oscillations of structural systems using non-linear normal modes.**J. Sound Vib.**, 77–101. (doi:10.1016/j.jsv.2003.04.005) Crossref, Web of Science, Google Scholar**273** - 10.
Touzé C, Amabili M . 2006 Nonlinear normal modes for damped geometrically nonlinear systems: application to reduced-order modelling of harmonically forced structures.**J. Sound Vib.**, 958–981. (doi:10.1016/j.jsv.2006.06.032) Crossref, Web of Science, Google Scholar**298** - 11.
- 12.
Neild SA, Wagg DJ . 2011 Applying the method of normal forms to second-order nonlinear vibration problems.**Proc. R. Soc. A**, 1141–1163. (doi:10.1098/rspa.2010.0270) Link, Google Scholar**467** - 13.
Cammarano A, Hill TL, Neild SA, Wagg DJ . 2014 Bifurcations of backbone curves for systems of coupled nonlinear two mass oscillator.**Nonlinear Dyn.**, 311–320. (doi:10.1007/s11071-014-1295-3) Crossref, Web of Science, Google Scholar**77** - 14.
Hill TL, Cammarano A, Neild SA, Wagg DJ . 2015 Interpreting the forced responses of a two-degree-of-freedom nonlinear oscillator using backbone curves.**J. Sound Vib.**, 276–288. (doi:10.1016/j.jsv.2015.03.030) Crossref, Web of Science, Google Scholar**349** - 15.
Neild SA, Champneys AR, Wagg DJ, Hill TL, Cammarano A . 2015 The use of normal forms for analysing nonlinear mechanical vibrations.**Phil. Trans. R. Soc. A**,**373**20140404 . (doi:10.1098/rsta.2014.0404) Link, Web of Science, Google Scholar - 16.
Liu X, Cammarano A, Wagg DJ, Neild SA, Barthorpe RJ . 2016*N*− 1 modal interactions of a three-degree-of-freedom system with cubic elastic nonlinearities.**Nonlinear Dyn.**, 497–511. (doi:10.1007/s11071-015-2343-3) Crossref, Web of Science, Google Scholar**83** - 17.
Hill TL, Cammarano A, Neild SA, Barton DAW . 2017 Identifying the significance of nonlinear normal modes.**Proc. R. Soc. A**,**473**20160789 . (doi:10.1098/rspa.2016.0789) Link, Google Scholar - 18.
Touzé C, Thomas O, Chaigne A . 2002 Asymmetric non-linear forced vibrations of free-edge circular plates. Part 1: Theory.**J. Sound Vib.**, 649–676. (doi:10.1006/jsvi.2002.5143) Crossref, Web of Science, Google Scholar**258** - 19.
Thomas O, Touzé C, Chaigne A . 2003 Asymmetric non-linear forced vibrations of free-edge circular plates. Part II: experiments.**J. Sound Vib.**, 1075–1101. (doi:10.1016/S0022-460X(02)01564-X) Crossref, Web of Science, Google Scholar**265** - 20.
Vakakis AF, Manevitch LI, Mikhlin YV, Pilipchuk VN, Zevin AA . 1996**Normal modes and localization in nonlinear systems**. New York, NY: Wiley. Crossref, Google Scholar - 21.
Kerschen G, Worden K, Vakakis AF, Golinval JC . 2006 Past, present and future of nonlinear system identification in structural dynamics.**Mech. Syst. Signal Process.**, 505–592. (doi:10.1016/j.ymssp.2005.04.008) Crossref, Web of Science, Google Scholar**20** - 22.
Manevitch AI, Manevitch LI . 2003 Free oscillations in conservative and dissipative symmetric cubic two-degree-of-freedom systems with closed natural frequencies.**Meccanica**, 335–348. (doi:10.1023/A:1023362112580) Crossref, Web of Science, Google Scholar**38** - 23.
Hill TL, Cammarano A, Neild SA, Wagg DJ . 2015 Out-of-unison resonance in weakly nonlinear coupled oscillators.**Proc. R. Soc. A**,**471**20140659 . (doi:10.1098/rspa.2014.0659) Link, Google Scholar - 24.
Breunung T, Haller G . 2018 Explicit backbone curves from spectral submanifolds of forced-damped nonlinear mechanical systems.**Proc. R. Soc. A**,**474**20180083 . (doi:10.1098/rspa.2018.0083) Link, Google Scholar - 25.
Lamarque CH, Touzé C, Thomas O . 2012 An upper bound for validity limits of asymptotic analytical approaches based on normal form theory.**Nonlinear Dyn.**, 1931–1949. (doi:10.1007/s11071-012-0584-y) Crossref, Web of Science, Google Scholar**70** - 26.
Temam R . 1988**Infinite-dimensional dynamical systems in mechanics and physics**. New York, NY: Springer. Crossref, Google Scholar - 27.
Dowell EH, Tang D . 2003**Dynamics of very high dimensional systems**. Hackensack, NJ: World Scientific. Crossref, Google Scholar - 28.
Manneville P . 1990**Dissipative structures and weak turbulence**. New York, NY: Academic Press. Google Scholar - 29.
Lacarbonara W, Nayfeh AH, Kreider W . 1998 Experimental validation of reduction methods for nonlinear vibrations of distributed-parameter systems: analysis of a buckled beam.**Nonlinear Dyn.**, 95–117. (doi:10.1023/A:1008389810246) Crossref, Web of Science, Google Scholar**17** - 30.
Steindl A, Troger H . 2001 Methods for dimension reduction and their application in nonlinear dynamics.**Int. J. Solids Struct.**, 2131–2147. (doi:10.1016/S0020-7683(00)00157-8) Crossref, Web of Science, Google Scholar**38** - 31.
Apiwattanalunggarn P, Shaw SW, Pierre C, Jiang D . 2003 Finite-element-based nonlinear modal reduction of a rotating beam with large-amplitude motion.**J. Vib. Control**, 235–263. (doi:10.1177/1077546303009001749) Crossref, Web of Science, Google Scholar**9** - 32.
Schilders WH, Van der Vorst HA, Rommes J . 2008**Model order reduction: theory, research aspects and applications**,**vol. 13**. Berlin, Germany: Springer. Crossref, Google Scholar - 33.
Rutzmoser J, Rixen D, Tiso P . 2014 Model order reduction using an adaptive basis for geometrically nonlinear structural dynamics. In*Conf. on noise and vibration engineering, Leuven, Belgium, 15–17 September*, pp. 20–25. Heverlee, Belgium: KU Leuven. Google Scholar - 34.
Rand RH . 1971 Nonlinear normal modes in two-degree-of-freedom systems.**J. Appl. Mech.**, 561–561. (doi:10.1115/1.3408826) Crossref, Web of Science, Google Scholar**38** - 35.
Jiang D, Pierre C, Shaw SW . 2005 Nonlinear normal modes for vibratory systems under harmonic excitation.**J. Sound Vib.**, 791–812. (doi:10.1016/j.jsv.2005.01.009) Crossref, Web of Science, Google Scholar**288** - 36.
Kuether RJ, Allen MS . 2014 A numerical approach to directly compute nonlinear normal modes of geometrically nonlinear finite element models.**Mech. Syst. Signal Process.**, 1–5. (doi:10.1016/j.ymssp.2013.12.010) Crossref, Web of Science, Google Scholar**46** - 37.
Amabili M, Touzé C . 2007 Reduced-order models for nonlinear vibrations of fluid-filled circular cylindrical shells: comparison of POD and asymptotic nonlinear normal modes methods.**J. Fluids Struct.**, 885–903. (doi:10.1016/j.jfluidstructs.2006.12.004) Crossref, Web of Science, Google Scholar**23** - 38.
Touzé C, Amabili M, Thomas O . 2008 Reduced-order models for large-amplitude vibrations of shells including in-plane inertia.**Comput. Methods Appl. Mech. Eng.**, 2030–2045. (doi:10.1016/j.cma.2008.01.002) Crossref, Web of Science, Google Scholar**197** - 39.
Elliott AJ, Cammarano A, Neild S, Hill T, Wagg DJ . 2018 Comparing the direct normal form and multiple scales methods through frequency detuning.**Nonlinear Dyn.**, 2919–2935. (doi:10.1007/s11071-018-4534-1) Crossref, PubMed, Web of Science, Google Scholar**94** - 40.
Xin ZF, Neild SA, Wagg DJ, Zuo ZX . 2013 Resonant response functions for nonlinear oscillators with polynomial type nonlinearities.**J. Sound Vib.**, 1777–1788. (doi:10.1016/j.jsv.2012.09.022) Crossref, Web of Science, Google Scholar**332** - 41.
Wagg DJ, Neild SA . 2015**Nonlinear vibration with control - for flexible and adaptive structures**, 2nd edn. Berlin, Germany: Springer. Google Scholar - 42.
Bellizzi S, Bouc R . 2005 A new formulation for the existence and calculation of nonlinear normal modes.**J. Sound Vib.**, 545–569. (doi:10.1016/j.jsv.2004.11.014) Crossref, Web of Science, Google Scholar**287** - 43.
Liu X, Wagg DJ . 2019*ε*^{2}-order normal form analysis for a two-degree-of-freedom nonlinear coupled oscillator. In*Proc. of the First Int. Nonlinear Dynamics Conf. Rome, Italy, 17–19 February*. Google Scholar - 44.
Dankowicz H, Schilder F . 2013**Recipes for continuation**,**vol. 11**. Philadelphia, PA: SIAM. Crossref, Google Scholar - 45.
Schilder F, Dankowicz H . 2017 Continuation Core and Toolboxes (COCO). https://sourceforge.net/projects/cocotools/. Google Scholar - 46.
Liu X, Wagg DJ, Neild SA . 2017 An explanation for why natural frequencies shifting in structures with membrane stresses using backbone curve models.**Nonlinear Dyn.**, 9–19. Google Scholar**1**