Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences
You have accessResearch articles

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

X. Liu

X. Liu

Nuclear Advanced Manufacturing Research Centre, University of Sheffield, Rotherham S60 5WG, UK

Google Scholar

Find this author on PubMed

D. J. Wagg

D. J. Wagg

Department of Mechanical Engineering, University of Sheffield, Sheffield S1 3JD, UK

[email protected]

Google Scholar

Find this author on PubMed


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 [27]. 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 [1317] 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 [2633], 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,3436]. 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 2N 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

where x is the {N × 1} vector of physical displacements, i.e. x=[x1,,xN], M and K are the {N × N} mass and stiffness matrices respectively, Nx 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 Px and a {2 × 1} forcing frequency related vector r, written as
where Ω 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 Nx are polynomial functions of x, 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. Mx¨+Kx=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. q=[q1,,qN], and Φ is the {N × N} matrix of linear modeshapes found from the eigenvalue problem ΦΛ = M−1KΦ, and Λ is a {N × N} diagonal matrix containing the square of the natural frequencies of the vibration modes, i.e. ω2ni 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 Φ gives

Equation (2.3) is then multiplied by (ΦMΦ)1 and, rearranged using M−1KΦ = ΦΛ, leading to
where Nq(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

where 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
where 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 Nv is an {N × 1} vector of nonlinear and damping terms in v and Pv 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

where Pv,i and Pq,i denote the ith row of Pv and Pq, respectively. Now making use of equation (2.11) and equation (2.10), the terms in the force transform matrix, e, can be determined using,
where ei is the ith 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 Nv 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

where u and H(u,u˙,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 H(u,u˙,r) is denoted to be of order ε1. The resulting normal form after the near-identity transform has been applied may be expressed as
which includes only the resonant terms, such that the ith state, ui 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 H(u,u˙,r) vector in equation (2.13) and the Nu(u,u˙,r) vector in equation (2.14) are approximated by a series of functions of reducing levels of significance, i.e.

where h(i)(u,u˙,r) and nu(i)(u,u˙,r) are {N × 1} vectors in the same form as H(u,u˙,r) and Nu(u,u˙,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 [v,v˙]=[u,u˙], we obtain
Furthermore, in recognition of the fact that the ωri values are close to, but not equal to ωni in general, a frequency detuning is introduced, to define this relationship
where Υ is a {N × N} diagonal matrix of the square of the response frequencies, ω2ri, 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

Note that only the terms of up to order ε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
Note that the sequential calculations, typical of the normal form approach is represented in equation (2.20b) because 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

Then, following the approach of [4], we consider that the response of each of the individual states in u, i.e. ui, responds at a single response frequency ωri, and its trial solution may be assumed to be

where Ui, ϕi and ωri are the displacement amplitude, phase lag and response frequency of ui, respectively. Using the left-hand side of equation (2.22), a vector u*(j) (of length Lj) containing all the polynomial combinations of uip and uim present in n(j) may be defined, and then n(j), h(j) and nu(j) are expressed as,

where [n(j)], [h(j)] and [nu(j)] are {N×Lj} matrices of corresponding coefficients. Substituting equation (2.23) into the jth homological equation gives
Note that the individual terms of u*(j) are functions of uip and uim, which may be written as

where =1,,Lj, and spj,ℓ,i, smj,ℓ,i, mpj,ℓ and mmj,ℓ are exponents of uip, uim, rp and rm 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

where ω*(j)ℓ is the response frequency of u*(j)ℓ, i.e.
Therefore, the second time derivative of the vector u*(j) may be written as
where Υ˘ is a {Lj×Lj} diagonal matrix of ω2(j)ℓ. So, finally equation (2.28) is substituted into equation (2.24) to obtain
where ° denotes the Hadamard product operator and β(j) is a {N×Lj} coefficients matrix, whose {i, ℓ}th element is
As nu(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 [nu(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


More details on this process including the different types of resonance obtained are given in [41].

Once the coefficient matrices [nu(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

where χn and χ~n are time-independent complex conjugates. By inspection of equation (2.32), a solution can be found by setting χi = 0 (or equivalently χ~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 Nx, but then transformed into Nq, Nv 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, L1=C2N+21+(C2N+21+C2N+21)+(C2N+21+C2N+22C21+C2N+23) 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 Lj 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 usp 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 usp, equation (2.20) becomes



where u^ is a {R × 1} vector containing the response of the R dominant states only. Now due to the reduction of the number of terms in u^ relative to that of the original state vector, u, the basis vector, u*(j), and coefficients matrices, [n(j)], [nu,(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

which only has R resonant equations and, then the number of time-independent equations related to response amplitudes and frequencies also shrinks to R. Here 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

where the coefficients of the nonlinear terms are
and the damping terms 2ζiωix˙i and forcing terms ficos(Ωt), for i = 1, 2 have been added after the conservative equations (i.e. with ζi = fi = 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].
Figure 1.

Figure 1. The example 2 d.f. system. For previous discussions of this example see [9,24,42].

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 q=[q1,q2]=[x1,x2], where


Here, a non-internal-resonant case is considered, i.e. ωr1r2 and ωr2r1, where n = 1, 2, …. We assume that the system is only forced in the horizontal direction, i.e. f2 = 0, thus only the horizontal response is significant. Therefore, the vector of reduced state vector is u^=[u1]. Substituting q=[q1,q2]=[u1p+u1m,0] into equation (4.3b), we obtain the vector of nonlinear terms of ε1 significance as


Expressing n(1)(u^) in the matrix form, i.e. n(1) = [n(1)]u*(1), gives

Here because of introducing u^, the size of 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
where ωr2 = ω2 is used as u2 = 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
As the example system has quadratic nonlinearity, the level of accuracy ε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), Δ, h˙(1)(u^), (Nv/v)|usp,u˙sp, and (Nv/v˙)|usp,u˙sp are required to be derived first, which in this example gives

Note that in equation (4.8b) derivatives of terms containing uipuim 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,

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 O(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 n(2)=Δh(1)+(Nv/v)|u,u˙h(1)+(Nv/v˙)|u,u˙h˙(1) in terms of u*(2), the coefficient matrix [n(2)] is found to be,

where resonant terms have boxes around them, and
where [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
where [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
and the vector of non-resonant terms for the order ε2 part is,
Now using Nu = nu(1) + nu(2), the reduced resonant equation of motion of the first resonance is
Substituting u1p = ((U1/2)e−iϕ1)eiωr1t and u1m = ((U1/2)eiϕ1)e−iωr1t into equation (4.14) and arranging the result as equation (2.32) gives
Setting χ1 = 0 and eliminating the phase-related term gives
which describes the forced response curve. Furthermore, substituting the forcing and damping related terms to zero, i.e. ζi = 0 and f1 = 0, gives the backbone curve expression of the system first mode to be
Note that expressions analogous to this have been previously obtained, for example by Touzé 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 ω2ni values. However, in this example, the nonlinear coefficients are of the same order as the ω2ni 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, f1 = [0.001, 0.002], f2 = 0 and ζ1 = ζ2 = 0.01, f1 = [0.01, 0.02], f2 = 0 and the results are shown in figure 2. The plots are presented in the projection of the response amplitude of the physical coordinates, |xi|, against the external forcing frequency, Ω. In figure 2, panels (a,c) show the |x1| response, and panels (b,d) show the |x2|. 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.

Figure 2.

Figure 2. The backbone curves and resonance curves of the 2 d.f. example system described in equation (4.1) for the case where its horizontal mode is dominant. The solid-black and solid-blue lines denote the analytically approximated backbone curves and forced response curves, respectively. The dashed-green and dashed-red lines represent the numerically computed backbone curves and forced response curves (from COCO), respectively. The black dots indicate branch points, where the backbone curves start. Parameters: ω1 = 2, ω2 = 4.5, (a,b): ζ1 = ζ2 = 0.001, f1 = [0.001, 0.002] and f2 = 0; (c,d): ζ1 = ζ2 = 0.01, f1 = [0.01, 0.02] and f2 = 0. The stability of the solution curves is not indicated on this figure. (Online version in colour.)

It is noteworthy that as the horizontal mode is assumed to be the dominant one here, and as a result, the vertical response, |x2|, 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].

Figure 3.

Figure 3. The plate system example. For further details, see the discussion by Liu et al. [46]. (Online version in colour.)

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

for 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 u2≠0 and u3≠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 u^=[u2,u3], which reduces the system from a four-mode to a two-mode approximation, the vector of nonlinear terms is found to be

where the subscript [2;3] indicates the nonlinear terms in n(1) associated with the second and third modes. Note that due to the specific value of the nonlinear coefficients, see table 2, n(1)(u^)[1]=n(1)(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

Using equations (5.3) and (5.4), the vectors resonant and non-resonant nonlinear terms are determined as
Finally, the time-invariant equations governing the response frequencies and amplitudes of the second and third modes are found to be,
Successively setting U3 = 0 and U2 = 0 results in the expressions of two single-mode (i.e. non-internally resonant) backbone curves given by
Combining the bracketed contents in equation (5.7) and rearranging gives the expression of the double-mode (i.e. internally resonant) backbone curves, written as
where, for balancing the complex terms, ei(|ϕ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. f3 = 0.02, which will ensure that the large-amplitude response assumption is satisfied.

Figure 4.

Figure 4. The backbone curves and resonance curves of the example plate system. The solid-red and solid-green lines denote the analytically approximated backbone curves S3 and D±23, respectively. The dashed-black and solid-blue lines are the numerically computed backbone curves and forced responses (using COCO), respectively. The red and green dots denote the branch points. Note that the information on the stability of the solution curves is not indicated on this figure. Also note that the unit of the response amplitude, |qi|, is metre and the plate thickness is 5 × 10−4m. The units of ωr is radians/second. (Online version in colour.)

Figure 4 shows that the double-mode backbone curves, D±23, bifurcate from the single-mode backbone curve, S3. Note that as the results are projected in the modal coordinates, the two double-mode backbone curves overlap (as opposed to being in the u2, u3 coordinate space). As we are forcing the system only in mode 3, the q3 response is driven by the external forcing, and has a hardening response curve (figure 4b) as would be expected for a system with cubic nonlinearity. The q2 response (figure 4a) is only triggered after the D±23 curves bifurcate from the S3 curve at an amplitude of approximately q3 = 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.


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.


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.

Table 1. The properties of the example plate.

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

Table 2. The values of coefficients for the lowest four modes of the example plate.

ωn1 = 58.9 α[1]111 = 5.45 × 109 ωn2 = 143.9 α[2]112 = 2.36 × 1010
α[1]122 = 2.36 × 1010 α[2]222 = 3.14 × 1010
α[1]133 = 2.27 × 1010 α[2]233 = 6.51 × 1010
ζ1 = 0.01 α[1]144 = 2.44 × 1010 ζ2 = 0.01 α[2]244 = 1.24 × 1011
α[1]234 = 7.43 × 1010 α[2]134 = 7.43 × 1010
ωn3 = 150.8 α[3]113 = 2.27 × 1010 ωn4 = 235.8 α[4]114 = 2.44 × 1010
α[3]223 = 6.51 × 1010 α[4]224 = 1.24 × 1011
α[3]333 = 3.14 × 1010 α[4]334 = 1.32 × 1011
ζ3 = 0.01 α[3]344 = 1.32 × 1012 ζ4 = 0.01 α[4]444 = 5.58 × 1010
α[3]124 = 7.43 × 1010 α[4]123 = 7.43 × 1010


Electronic supplementary material is available online at

Published by the Royal Society. All rights reserved.