## Abstract

Vibration problems are naturally formulated with second-order equations of motion. When the vibration problem is nonlinear in nature, using normal form analysis currently requires that the second-order equations of motion be put into first-order form. In this paper, we demonstrate that normal form analysis can be carried out on the second-order equations of motion. In addition, for forced, damped, nonlinear vibration problems, we show that the invariance properties of the first- and second-order transforms differ. As a result, using the second-order approach leads to a simplified formulation for forced, damped, nonlinear vibration problems.

## 1. Introduction

This paper concerns a class of discrete nonlinear vibration problems for which the equations of motion can be written as
1.1
where **x** is the *N*×1 displacement vector, *M*, *C* and *K* are the *N*×*N* real symmetric mass, damping and stiffness matrices, respectively, *Γ*_{x} is an *N*×1 vector of nonlinear terms—which are assumed to be small—*N* is the number of degrees of freedom in the system and an overdot represents differentiation with respect to time. The forcing is assumed to be sinusoidal and is represented by *P*_{x}**r**, where *P*_{x} is an *N*×2 forcing amplitude matrix and **r** = {*r*_{p} *r*_{m}}^{T} is a 2×1 forcing vector with *r*_{p} = e^{iΩt} and *r*_{m} = e^{−iΩt}, where *Ω* is the forcing frequency, and the subscripts *p* and *m* indicate the sign of the exponential term, plus and minus, respectively.

The motivation for this paper comes from the problem of how to transform equation (1.1) into a simpler form. In particular, the linearized, unforced, version of equation (1.1), with *Γ*_{x} = 0 and *P*_{x} = 0, can be transformed (with constraints) into a series of uncoupled equations relating to the *normal modes of vibration* of the system. Each of the resulting modes of vibration relates to the dominant physical configuration of the system at the corresponding natural frequency. This link with the physical problem allows the linear modal model to be used in combination with experimental techniques, particularly for what has now become known as *modal testing* (Ewins 2000).

Historically, for the nonlinear problem, there has been considerable work carried out on trying to identify *nonlinear* normal modes of vibration—a recent summary of this work is given by Kerschen *et al.* (2006)—see also Boivin *et al.* (1995), Nayfeh *et al.* (1999) and Touze & Amabili (2006). An alternative, but related approach is to search for a simplifying transformation of equation (1.1) using the *method of normal forms*. Comprehensive discussions of normal form theory can be found, for example, in Arnold (1988), Murdock (2002) and Sanders *et al.* (2007), and a survey of recent developments is given by Stolovitch (2009). Normal forms have been used to treat a range of dynamics problems, see for example Hsu (1983), Fredriksson & Nordmark (2000), Pelinovsky & Yang (2002), Leung & Zhang (2003) and Wang & Bajaj (2009). Mathematical aspects of the normal form technique have been studied by many authors, including Zhuravlev (1983), Walcher (1993), Bruno & Walcher (1994), Mayer *et al.* (2004), Stolovitch (2009), Meyer *et al.* (2009) and references therein.

In vibration problems, the relationship between system resonances and external forcing is particularly important. The method of normal forms has already been applied to forced nonlinear vibration problems, with the system equations expressed in the first-order form (Jezequel & Lamarque 1991; Nayfeh 2000). In this paper, we apply the method of normal forms to forced nonlinear vibration problems with the system equations expressed in the second-order form. We demonstrate that the invariance properties of the first- and second-order transforms differ. The invariance properties of the second-order transform leads to simpler transformed equations, which are therefore easier to solve. Two other modifications are described: (i) a matrix formulation is used to simplify the normal form process and (ii) a new formula is derived for computing the terms in the Lie bracket.

These latter two modifications have already been applied to first-order systems (Wagg & Neild 2009) and are summarized in §2, which describes the first-order normal form approach. In §3, the method of normal forms for equations in the second-order form is presented. Section 4 compares the normal form techniques (for equations in the first- and second-order form) for a single-degree-of-freedom system and then in §5 the method of normal forms, using equations in the second-order form, is applied to a two degree-of-freedom system. Conclusions are drawn in §6.

## 2. Method of normal forms for first-order systems

### (a) General formulation

The first step is to rewrite the second-order equation of motion, equation (1.1), in first-order form, using (**y** has size 2*N*×1), to give
2.1
where
2.2

Now a linear modal transformation **y** = *Φ***q** can be applied, where *Φ* is the 2*N*×2*N* eigenvector matrix relating to *A*. The corresponding eigenvalues of *A* are placed in the leading diagonal of a matrix forming the diagonal matrix *Λ* where, *Λ* = *Φ*^{−1}*AΦ*. Note that the leading diagonal values of *Λ* will be termed . Applying this linear modal transformation results in an equation of motion of the form
2.3
where the nonlinear and direct forcing terms are given by
2.4

The aim of the normal form method is to apply transformations resulting in an equation of the form , where the nonlinear and forcing terms (*Γ*_{u}(**u**,**r**) and *P*_{u}**r**, respectively) have a simpler form than the original equations. In addition, we seek the property that with the selection of a simple trial solution for **u**, *u*_{n} = *U*_{n}e^{iωrnt} for *u*_{n}, , leads to total cancellation of the time-dependent exponential terms. This results in a time-independent, spatial, relationship that is usually referred to as the vibration modes of the system.

### (b) Dealing with forcing terms

To simplify the forcing terms, we apply a transform of the form
2.5
where [*e*] has size 2*N*×2. Substituting this transformation into equation (2.3) gives
2.6
where *W* is a 2×2 diagonal matrix with the first and second leading diagonal terms being i*Ω* and −i*Ω*, respectively. By grouping the **r** terms, equations (2.6) can be rewritten as
2.7
where
2.8
which represents the relationship between the forcing matrices and the transformation matrix [*e*]. Equation (2.8) can be simplified by writing , where the element in the *n*th row () and *k*th column (*k* = 1,2) of , can be related to the corresponding element in [*e*] (*e*_{n,k}) using . It can be seen that when *Ω* is close to the *n*th eigenvalue, resulting in −*λ*_{n} + *i*^{2k−1}*Ω* ≈ 0, *e*_{n,k} becomes large, and these are called near-resonance terms.

To deal with the resonant (or near-resonant) terms, all those in [*e*] are set to zero, with the corresponding terms in *P*_{q} set to equal those in *P*_{v}, see for example Jezequel & Lamarque (1991). All the non-resonance terms in [*e*] are set such that for these terms with the corresponding terms in *P*_{v} set to zero. As a result, only resonant terms remain in the forcing term *P*_{v}**r** and, as will be seen in the example discussed in §4 of the paper, setting the response frequencies *ω*_{ri} to match the relevant resonant forcing frequency leads to the desired cancellation of the time-varying exponential terms.

Now the nonlinear term in the linear normal form representation of the equation of motion, equation (2.7), is partitioned into a series of functions with reduced levels of significance
2.9
Here *ϵ* is a small parameter and may be thought of as a book-keeping aid—it allows the tracking of the significance of each term.

The fourth step is to apply a further nonlinear transform, from **v** to **u**, given by
2.10
This is often referred to as the *near-identity transform* (e.g. Nayfeh 1993). Applying this, as-yet unspecified, transform results in the dynamic equation
2.11
Ideally, *Γ*_{u}(**u**,**r**) = 0 such that the dynamic equation is linear. However, owing to the fact that the transformation is a near-identity transform, this is not usually possible, and the aim is to simplify equation (2.11) as much as possible without invalidating the assumption that the transformation is near-identity, such that **h**(**u**,**r**) is of order *ϵ*^{1}.

Eliminating **v** from equation (2.9) using equation (2.10), and then replacing using equation (2.11), gives
2.12
Now, equating the zeroth- and first-order powers of *ϵ* gives
2.13
and
2.14
where a Taylor series expansion has been applied to **f**_{1}(**u**,**r**) such that . Conditions for the higher order *ϵ* terms can also be derived; however, these are not considered here.

In considering the *ϵ*^{1} equation, it is assumed that the response frequency for each state, (which may be positive or negative) for states , respectively, is close to the natural frequency of the state. A 2*N*×2*N* diagonal response matrix *Υ* is defined in which the *n*th diagonal element is i*ω*_{rn} for , note that . This allows the eigenvalue matrix *Λ* to be written in terms of the response frequency matrix *Υ* and a detuning matrix Δ; *Λ* = *Υ* + Δ, where the *n*th diagonal element of Δ is *λ*_{n}−i*ω*_{rn}.

The system is assumed to be weakly nonlinear, so that for small damping values *λ*_{n} ≈ i*ω*_{nn}, where *ω*_{nn} is the undamped natural frequency of the *n*th mode, and thus Δ is small—order *ϵ*^{1}, which on substitution in equation (2.14) leads to an *ϵ*^{2} term. This allows equation (2.14), the *ϵ*^{1} equation, to be written as
2.15
Now a vector **u***, which contains all the combinations of *u*_{n} (), *r*_{p} and *r*_{m} terms that are present in **f**_{1}(**u**,**r**) and is of length *L*, is defined. This allows the **f**_{1}(**u**,**r**), **g**_{1}(**u**,**r**) and **h**_{1}(**u**,**r**) terms to be expressed in matrix form
2.16
where matrices [*f*], [*a*] and [*b*] are of size 2*N*×*L*. Substituting these expressions into equation (2.15) gives
2.17

### (c) A formula for the Lie bracket terms

Consider the ℓth element in **u*** to have the form
2.18
where the *s*_{ℓn}, *m*_{ℓp} and *m*_{ℓm} constants indicate the power of the *u*_{n}, *r*_{p} and *r*_{m} terms, respectively, in *u**_{ℓ}. The ℓth element of the derivative vector in equation (2.17) may be written as
2.19
2.20
where the relationship has been used to write d*u*_{n}/d*t* = i*ω*_{rn}*u*_{n} for 1 ≤ *n* ≤ 2*N*. Since the derivative of with respect to time is linearly related to , from equation (2.20), the derivative of **u*** may be written as
2.21
where is a diagonal matrix with, using equation (2.20), the ℓth diagonal element being given by
2.22
Using this information, equation (2.17) may now be rewritten as
2.23
Considering non-zero **u***(**u**,**r**) solution to this equation gives
2.24
where is the Lie bracket. Using equation (2.22), the element in the *n*th row and ℓth column of matrix , may be written as
2.25
where *b*_{n,ℓ} is the element in the *n*th row and ℓth column of matrix [*b*].

It is desirable for **g**(**u**,**r**) to be as simplified as possible; therefore, [*a*] should contain as many zeros as possible. The restriction on this is that the transformation must be near-identity, and so the terms in **h**(**u**,**r**) and therefore [*b*] must be small. In the ideal case, [*a*] = [0] is selected as an initial assumption. Then equations (2.24) and (2.25) can be used to find the *b*_{n,ℓ} elements using *b*_{n,ℓ} = *f*_{n,ℓ}/*β*_{n,ℓ}. However, this is not possible if a *β*_{n,ℓ} value is zero or approximately zero for the system parameters under consideration. In both these cases, *b*_{n,ℓ} is set to zero and, to satisfy equation (2.24), *a*_{n,ℓ} = *f*_{n,ℓ} is selected.

The key feature of the dynamic equation
2.26
where *ϵ*^{2} terms and higher have been assumed to be insignificant, is that only resonant terms exist within **g**_{1} and *P*_{u} and therefore the equation can be solved exactly using a trial solution where each term in **u** takes the form *u*_{n} = *U*_{n}e^{iωrnt}. This will be seen when an example is considered in §4.

## 3. Normal forms applied to second-order nonlinear vibration problems

In the following analysis, the same notation has been adopted as for the first-order normal form to highlight the similar steps used in the two methods; however, it must be noted that the content and size of the various matrices and functions are different.

### (a) General formulation

For the second-order problem, standard linear modal theory can be used to first transform equation (1.1) into a form in which the undamped linear terms are decoupled. This is done using the undamped linear modes, which are the modes based on the equation . This may be rewritten in eigen form as , where **p**_{n} is the *n*th modeshape (and eigenvector of *M*^{−1}*K*) and is the corresponding square of the natural frequency (and eigenvalue). Now the *N*×*N* matrix of modeshapes *Φ*, in which the *n*th column is **p**_{n}, and the diagonal eigenvalue matrix *Λ*, in which the *n*th diagonal element is , are defined. Equation (1.1) may now be written as
3.1
in which **x** = *Φ***q**. Rearranging, and noting that by definition *M*^{−1}*KΦ* = *ΦΛ*, gives
3.2
where
3.3
The damping terms, which are assumed to be small, are included within the nonlinear function *Γ*_{q}.

Now the transformation designed to simplify the forcing terms, in a similar way to that used to generate the Jordan form for the first-order method, is applied. Substituting the transform (2.5), where [*e*] now has size *N*×2, into equation (3.2) gives
3.4
where, as with the first-order method, *W* is a 2×2 diagonal matrix with the first and second diagonal values being i*Ω* and −i*Ω*, respectively. This may be written in the form
3.5
where the relationship between the nonlinear terms is
3.6
and the relationship between the forcing and the transformation matrix is
3.7
where the *n*th row () and *k*th column (*k* = 1,2) of may be written in terms of the corresponding element in [*e*] using (recalling that ). As with the first-order method, when the forcing is close to a natural frequency, then the corresponding terms in [*e*] become large. These terms, which are near-resonant terms, are set to zero in [*e*] and the corresponding terms in *P*_{v} are set to equal those in *P*_{q} to satisfy equation (3.7). All the other terms in [*e*] are set such that for these terms and the corresponding terms in *P*_{v} are set to zero.

### (b) Near-identity transform for second-order case

Now the near-identity transform can be applied using
3.8
3.9and
3.10
where equation (3.8) is the current dynamic equation, equation (3.9) is the transformation and equation (3.10) is the resulting dynamic equation. Combining these equations to eliminate **v** gives
3.11
Applying a Taylor series expansion to **f**_{1} and equating zeroth- and first-order powers of *ϵ* gives
3.12
and
3.13
Equations for higher order *ϵ* may also be derived; however, this is beyond the scope of the current discussion.

To satisfy the *ϵ*^{1} equation, the form of the response of the states needs to be considered. As with the first-order normal form approach, the transform removes non-resonant nonlinear terms, resulting in a response for each state at a single response frequency, . However, since the differential equation in **u** is second order, the trial solutions for the states must consist of both a positive and a negative complex exponential term. The state vector **u** is therefore split into components **u** = **u**_{p} + **u**_{m}, allowing the trial solutions for the *n*th state to be
3.14
for 1 ≤ *n* ≤ *N*, where *ω*_{rn} are positive. The amplitudes of the *u*_{nm} states, *U*_{n}e^{iθn}/2, are the complex conjugates of the *u*_{np} states to ensure that **u** is real. This is required to ensure that the *ϵ*^{0} representation of **q** and therefore **x** are real for all time (noting that for the second-order representation, the modeshape matrix, *Φ*, is real, in contrast to the first-order representation, and so **q** must be real to ensure real **x**). For convenience, the trial solutions have been written such that . The diagonal, *N*×*N*, response matrix *Υ* is now defined in which the *n*th diagonal element is i*ω*_{rn}. This allows the time derivatives of **u** to be written as and .

Now a frequency detuning matrix, Δ, is introduced to simplify the analysis close to resonance. It is assumed that the response frequencies are close to the linear undamped natural frequencies of the modes, *ω*_{rn} ≈ *ω*_{nn}. The *n*th diagonal elements in the eigenvalue matrix *Λ* and the square of the response matrix *ΥΥ* are and , respectively, which allow the detuning matrix to be written as *Λ* = −*ΥΥ* + Δ, where the *n*th diagonal element of Δ is . Since *ω*_{rn} ≈ *ω*_{nn} for all *n*, Δ is small—order *ϵ*^{1} and so the *ϵ*^{1} order equation, (3.13), may be written as
3.15

As with the first-order method, a vector **u*** (of length *L*) is defined; however, here it contains all the combinations of *u*_{np}, *u*_{nm} (1 ≤ *n* ≤ *N*), *r*_{p} and *r*_{m} terms that are present in **f**_{1}(**u**,**r**). Using **u** = **u**_{p} + **u**_{m} and allows the following matrix expressions to be defined:
3.16
where [*f*], [*a*] and [*b*] are of size *N*×*L*. Substituting equation (3.16) into equation (3.13) gives
3.17

### (c) Formula for the Lie bracket terms in the second-order case

Considering the ℓth element in **u***, and expressing it in terms of its states, gives
3.18
Using and , the time derivative of the ℓth element of vector **u*** may be written as
3.19

Since the derivative of with respect to time is linearly related to , from equation (3.19), the second derivative of **u*** may be written as
3.20
where is a diagonal matrix of size *N*×*N* in which, using equation (3.19), the ℓth diagonal element is given by
3.21
Using this information, equation (3.17) may now be rewritten as
3.22
Considering non-zero **u***(**u**,**r**), solution to this equation gives
3.23
Using equation (3.21), the element in the *n*th row and ℓth column of matrix , may be written as
3.24
where *b*_{n,ℓ} is the element in the *n*th row and ℓth column of matrix [*b*].

As with the first-order formulation, [*a*] and [*b*] can now be selected by considering the size of the *β*_{n,ℓ} terms. It is desirable for [*a*] to contain as many zeros as possible so that the dynamic equation in **u** is as simple as possible. Therefore, most of the terms in [*a*] are set to zero and equation (3.23) is satisfied by setting *b*_{n,ℓ} = *f*_{n,ℓ}/*β*_{n,ℓ} for each element. The exception to this is where *β*_{n,ℓ} is zero or close to zero. In these cases, the corresponding [*b*] values are set to zero and [*a*] = [*f*] is selected.

## 4. A single-degree-of-freedom Duffing oscillator

As a first example, a single-degree-of-freedom Duffing oscillator is considered. The equation of motion is given by
4.1
where the system has light damping, the nonlinearity is small and the forcing is close to resonance, *Ω* ≈ *ω*_{n}. The forcing may be rewritten in the form shown in equation (1.1), where and *p* = *F*e^{iϕ}/2 (where is the complex conjugate of *p*).

### (a) First-order normal forms

The first-order normal form for equation (4.1) has already been derived by Jezequel & Lamarque (1991). Using either their approach, or the modified approach set out in §2, the near-identity transformation equation can be used to find the response *x*. The component of this response at the forcing frequency, *x*_{Ω}, is given by
4.2
This is based on the relationships calculated from the equations for the dynamics
4.3
and
4.4
which can be solved numerically to give the response amplitude *U* and relative phase *ϕ* for given forcing frequencies *Ω*. Note that for all these equations, the trial solutions *u*_{1} = (*U*/2)e^{iωt} and *u*_{2} = (*U*/2)e^{−iωt} have been used. We note also that
4.5
where are the eigenvalues of the first-order linearized system.

### (b) Second-order normal forms

As a one-degree-of-freedom system, the Duffing oscillator represented in equation (4.1) is already in a modal form. The transformation **x** = *Φ***q** is therefore a unity transformation (*Φ* = 1 with the corresponding eigenvalue ), and equation (4.1) may be written in the form of equation (3.2), , where
4.6
Now the transform **q** = **v** + [*e*]**r**, equation (2.5), is applied. Using (*n* = 1, *k* = 1,2), it can be seen that, with the forcing close to resonance, both terms in [*e*] become large compared with indicating resonance terms. Therefore, [*e*] = [00] is selected and *P*_{v} = *P*_{q} is used to ensure that equation (3.7) is satisfied. Since [*e*] = [00], equation (3.6) gives .

Provided the damping and nonlinearity are small, using equation (3.9) allows to be defined, giving 4.7

where to indicate small terms the notation and have been introduced. The vector **u*** is defined such that the terms contained in **u*** allow to be written as , equation (3.16). So [*f*] and **u*** may be written as
4.8
Note that in deriving these expressions, the relationship is used. Now may be calculated using equation (3.24) to give
4.9
Note that if the detuning approximation from equations (3.13) to (3.15) had not been made, the terms would be , indicating the potential for 1 : 3 resonance. To satisfy , equation (3.23), while ensuring that the near-identity transform is small, results in [*a*] and [*b*] matrices (defined in equation (3.16))
4.10
Using equations (3.9), (3.10) and (3.16) allows the second-order differential equation in **u** () and the transform (**v** = **u** + *ϵ*[*b*]**u***) to be written as
4.11
and
4.12
respectively, noting that *P*_{u} = *P*_{v} (equation (2.13)) and . Adopting the trial solution *u*_{1p} = (*U*/2)e^{iωt} = (*U*/2)*r*_{p} and *u*_{1m} = (*U*/2)e^{−iωt} = (*U*/2)*r*_{m} where *U* is real, these may be written as
4.13
and
4.14
where *x* = **q** = **v** = *v*_{1}. Finally, equation (4.13) can be harmonically balanced *exactly*—as the higher frequency terms have been removed in the transform from *v* to *u* to give the time-independent equations
4.15
and
4.16
The phase between the forcing and the response can be eliminated by squaring and adding equations (4.15) and (4.16) to give a quadratic in *Ω*^{2}
4.17
This equation can be used to find valid solutions for forcing frequency given a range of response amplitudes *U* and a forcing amplitude *F*, allowing the calculation of the frequency–response curve for *U*.

### (c) Comparison of the two solutions

Figure 1*a* shows the amplitude of response at frequency *Ω* calculated using a fourth-order Runge–Kutta time-stepping simulation and the first- and second-order normal forms. The forcing frequency of the simulation was slowly increased to follow the upper branch of the solution and then slowly decreased from the maximum frequency to follow the lower branch. For the first-order normal form solution, the frequency at which the peak response occurs is approximately 0.8 per cent different from that calculated using the time-stepping. However, the second-order normal form solution gives an even closer approximation to the time-stepping simulation. Figure 1*b* shows that the second-order normal form solution also captures the response at three times the forcing frequency with a high degree of accuracy.

The difference in the first- and second-order methods comes from a difference in their invariance properties. We explain the mechanism for this with the following remarks:

### Remark 4.1

In first-order form, the eigenvalue decomposition results in the natural frequencies *ω*_{n} being split into two eigenvalue components i*ω*_{n} and −i*ω*_{n} () if damping is included). In contrast, the second-order form has eigenvalues of for each *ω*_{n}.

### Remark 4.2

The result of this is that the trial solutions for the first-order form applied to a one-degree-of-freedom system are *u*_{n} = (*U*/2)e^{iωrnt} for *n* = 1,2 and *ω*_{r1} = −*ω*_{r2} = *Ω*, whereas for the second-order form it is *u*_{n} = (*U*/2)(e^{iωrnt} + e^{−iωrnt}) for *n* = 1 and *ω*_{r1} = *Ω*.

### Remark 4.3

To understand the difference in the properties of the first- and second-order formulations, consider an unforced, undamped system with *Γ*_{x} consisting of just a small linear detuning term
4.18
where *α* is order *ϵ*^{1}. In this case, *N* = 1, *L* = 2 and
4.19
So, for the second-order approach, the near-identity transform is invariant. This is because for the near-identity transform **v** = **u** + *ϵ*[*b*]**u*** to be invariant, *b*_{i} = 0, ∀*i*. As a result, the response frequency is exact; *ω*_{r1} = *ω*_{n}(1 + *γ*) where *α* = (2*γ* + *γ*^{2}).

### Remark 4.4

The reason for the difference in these invariance properties is that for the second-order form, the response of the *n*th mode at the *n*th response frequency is completely captured by the *n*th dynamic equation in **u** (the *n*th row of the matrix equation (3.10)), none of the response is contained in the transform. Specifically, when transforming from **v** to **u**, the terms that remain in the *n*th dynamic equation are those that have a time-dependent component of the form e^{iωrnt} or e^{−iωrnt} once the trial solutions have been made. These terms can be exactly harmonically balanced owing to the form of the trial solution (as stated in remark 4.2). This can be seen by considering the equation defining resonant terms (i.e. terms that remain in the dynamic equation), which is from equation (3.24). In particular, note the ± on the right-hand side.

### Remark 4.5

In contrast, even for this example, which contains no sub- or super-harmonics, the first-order near-identity transform is not invariant. For the first-order form, during the **v** to **u** transform, the terms that remain in the *n*th dynamic equation are those that have a time-dependent component of the form e^{iωrnt} once the trial solutions have been made. As a result, the terms with an e^{−iωrnt} time-dependent component, once the trial solutions have been made, become part of the transform. In line with this, the equation defining the resonant terms is from equation (2.25). In contrast to the expression in remark 4.4, note the lack of ± on the right-hand side.

### Remark 4.6

The first-order transform for the linear detuning example results in
4.20
where *ω*_{r1} = −*ω*_{r2} is the natural frequency including small linear detuning. From the dynamic equation, the response frequency is predicted to be *ω*_{r1} = *ω*_{n}(1 + *γ* + *γ*^{2}/2), which is accurate to order *ϵ*^{1} but not exact. Taking a detuning of *γ* = 0.075, which is similar to that observed in the Duffing example, results in a frequency error of 0.26 per cent, which is of the same order as the error in figure 1.

This indicates where the two transform techniques lose equivalence, leading to a difference in accuracy when simulating resonant behaviour. We note that the second-order method has further advantages:

### Remark 4.7

A key operation in the transformation is to define . This is done using equations (2.24), (2.25) for the first-order and equations (3.23), (3.24) for the second-order. The dimension of is 2*N*×*L* in the first-order case but only *N*×*L* in the second-order case.

### Remark 4.8

The frequency–response curve for *U* in the second-order case can be found from the analytical expression (4.17). However, in the case of the first-order normal form, the relative phase cannot be easily eliminated since it is contained within the *e*_{1,2} and *e*_{2,1} terms as well as the terms on the right-hand side of equations (4.3) and (4.4). The frequency–response curve for *U*, for a given *F*, must therefore be solved numerically using these equations.

## 5. Two-degree-of-freedom oscillator

Consider the following two-degree-of-freedom oscillator, written in modal form
5.1
and
5.2
The forcing is close to the first-mode natural frequency of the linearized system, *Ω* ≈ *ω*_{n}, and for the linearized system, the second mode has twice the natural frequency of the first mode. It is assumed that the nonlinear and damping terms are small.

This can be written as the matrix modal equation used for the second-order normal forms, given by equation (3.2), by defining
5.3
with *r*_{p} = e^{iωt} and *r*_{m} = e^{−iωt}.

Now the forcing transform given in equation (2.5), **q** = **v** + [*e*]**r**, is applied to give a dynamic equation in the form of equation (3.5), . The matrix [*e*] and the forcing function in equation (3.5), *P*_{v}**r**, are selected using equation (3.7) and the relationship . Noting that *e*_{1,1} and *e*_{1,2} would be near-resonant terms if the corresponding terms in *P*_{v} were set to zero since the forcing is close to the first natural frequency, [*e*] and *P*_{v} are set to
5.4
where . To finish the forcing transformation, the nonlinear term in the resulting dynamic equation, (3.5), must be derived using equation (3.6) giving
5.5

The near-identity transform, equation (3.9), is now derived. Firstly, since the damping and nonlinearity are small, using equation (3.9) gives where
5.6
where to indicate small terms the book-keeping notation , and have been used. Using the substitution **u** = **u**_{p} + **u**_{m} allows the nonlinear term to be written in the form given in equation (3.16), where [*f*] and **u*** are
5.7

By considering each element in **u*** taking the form given in equation (3.18), the 2×14 matrix can be calculated using equation (3.24), giving
5.8

Now matrices [*a*] and [*b*], which represent the order *ϵ*^{1} terms in the dynamic equation in terms of *u* and the transform from *v* to *u*, respectively (as defined in equation (2.16)), can be derived. Recall that the expression , equation (2.24), must be satisfied with the desire that as many elements in [*a*] as possible are zero, while ensuring that the terms in [*b*] are small, non-resonant terms. In deriving these matrices, the response frequencies must be selected. Firstly, since the forcing is close to the first natural frequency, the first response frequency may be written as *ω*_{r1} = *Ω*. Secondly, since the undamped linearized natural frequency for the second mode is twice that of the first mode, and recalling that the response frequencies are assumed to be close to the natural frequencies, the second response frequency may be written as *ω*_{r2} = 2*ω*_{r1} = 2*Ω*. Noting that [*b*] is related to via equation (5.8), it can be seen that the terms [1,1],[1,2],[1,5],[1,6],[1,8],[1,11],[2,3],[2,4],[2,9],[2,14] are resonant and so for these terms *a*_{n,l} must be set to *f*_{n,l}, giving
5.9

With reference to equations (3.10) and (3.16), the dynamic equations in terms of *u* may now be written as
5.10
and
5.11
where the relationships with the *n*th diagonal element in matrix *Υ* being i*ω*_{rn}, *ω*_{r1} = *Ω*, *ω*_{r2} = 2*Ω*, , and have been used.

Substituting the trial solutions, given in equation (3.14), into these equations and then balancing the harmonic terms, and for equation (5.10) and and for equation (5.11), give exact relationships
5.12
and
5.13
for equation (5.10) and
5.14
and
5.15
for equation (5.11). These equations can be solved numerically, or recognizing that the response in *u*_{1} will be significantly larger than that in *u*_{2} (since the frequency of *u*_{1} is *ω*_{r1} = *Ω*, whereas it is *ω*_{r2} = 2*Ω* for *u*_{2}) an approximate solution can be found in two steps. Firstly, a crude approximation to the response of mode 1 is found by ignoring the nonlinear terms in equations (5.12) and (5.13), and eliminating the phase lag *θ*_{1} to give the amplitude of response
5.16
For the second mode, where the nonlinear terms can be viewed as a forcing, the approximate amplitude of response can be found using
5.17
along with equation (5.16). In addition, for this case, it can be shown that 2*θ*_{1}−*θ*_{2} = *π*. Using these approximate relationships, a more accurate expression for the *U*_{1} amplitude may be found by reconsidering equations (5.12) and (5.13), including the nonlinear terms. In this case, the expression for *U*_{1} is
5.18
where *U*_{2} is calculated from equation (5.17).

The response of the two modes at frequencies other than *ω*_{n} can be calculated by considering the transformation from **u** to **q**. Using equations (2.5), (3.9) and (3.16), this may be written as
5.19
Using this equation, it can be seen that the response of mode *n* at the response frequency, *ω*_{rn}, is represented entirely by *u*_{n}, as noted in remark 4.4, so that
5.20
where the first subscript of *q* indicates the mode number and the second subscript indicates that the expression is just the portion of the total modal response containing the terms at this frequency. Whereas the response of mode 2 at, for example, the forcing frequency comes from [*b*] and [*e*]
5.21
the response of mode 1 at twice the forcing frequency comes from [*b*]
5.22

Figures 2 and 3 show the predicted response amplitudes of modes 1 and 2, respectively, using the normal form technique compared with the results of time-stepping simulations. It can be seen that the agreement is very good.

### Remark 5.1

In addition to demonstrating the accuracy of the normal form technique, this example shows that although the transformed dynamic equation in **u** contains just the responses at the response frequencies, *ω*_{rn} for the *n*th mode, the response of the second mode is predicted accurately. Compare this to the averaging technique (see Tondl *et al.* (2000) and Verhulst (1996) with its implementation on a multi-mode system demonstrated in Gonzalez-Buelga *et al.* (2008)), for example: averaging assumes a response of each mode just at the response frequencies of the form
and
5.23
where *ω*_{r1} = *Ω* and *ω*_{r2} = 2*Ω* and results in a prediction of zero response of the second mode in the steady state. The reason for this discrepancy is that the mechanism for generating a response in the second mode at twice the forcing frequency is a frequency mixing of the first-mode response at the forcing frequency with the second-mode response at the forcing frequency—which, using averaging, is assumed to be zero (the response is taken to be entirely at twice the forcing frequency). This mechanism is captured in the normal form technique owing to the forcing transformation **q** = **v** + [*e*]**r**, equation (2.5), which was applied prior to the normal form transformation. Using this transformation, the first mode was unaffected, *q*_{1} = *v*_{1}; however, the second-mode response was altered, *q*_{2} = *v*_{2} + *E*(*r*_{p} + *r*_{m}), see equation (5.4). The term *E*(*r*_{p} + *r*_{m}), where , represents the linear undamped response of the second mode at the forcing frequency. Inclusion of this response at the forcing frequency resulted in the forcing term −*βE*(*u*_{1p}*r*_{p} + *u*_{1m}*r*_{m}) present in the equation for the dynamics of the second mode in **u**, see equation (5.11), and hence results in the response of the second mode at twice the forcing frequency.

## 6. Conclusions

When a vibration problem is nonlinear in nature, equation (1.1), using normal form analysis has, up until now, required that the second-order equations of motion be put into first-order form. In this paper, we have demonstrated how a normal form analysis can be carried out for vibration problems of the type expressed as equation (1.1). We have shown that the invariance properties of the first- and second-order transforms differ, because in the first-order approach terms are split. As a result, using the second-order approach leads to a simpler form and improves accuracy when simulating resonant behaviour.

- Received May 27, 2010.
- Accepted September 22, 2010.

- This journal is © 2011 The Royal Society