Comparing the direct normal form and multiple scales methods through frequency detuning

Approximate analytical methods, such as the multiple scales (MS) and direct normal form (DNF) techniques, have been used extensively for investigating nonlinear mechanical structures, due to their ability to offer insight into the system dynamics. A comparison of their accuracy has not previously been undertaken, so is addressed in this paper. This is achieved by computing the backbone curves of two systems: the single-degree-of-freedom Duffing oscillator and a non-symmetric, two-degree-of-freedom oscillator. The DNF method includes an inherent detuning, which can be physically interpreted as a series expansion about the natural frequencies of the underlying linear system and has previously been shown to increase its accuracy. In contrast, there is no such inbuilt detuning for MS, although one may be, and usually is, included. This paper investigates the use of the DNF detuning as the chosen detuning in the MS method as a way of equating the two techniques, demonstrating that the two can be made to give identical results up to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon ^2$$\end{document}ε2 order. For the examples considered here, the resulting predictions are more accurate than those provided by the standard MS technique. Wolfram Mathematica scripts implementing these methods have been provided to be used in conjunction with this paper to illustrate their practicality. Electronic supplementary material The online version of this article (10.1007/s11071-018-4534-1) contains supplementary material, which is available to authorized users.


Introduction
In recent years, there has been substantial interest in the study of backbone curves, due to their utility in studying lightly damped nonlinear vibrations in multi-degree-of-freedom (MDOF) mechanical structures. The motivation for this paper comes from observations made by the authors when comparing backbone curves found using the multiple scales (MS) method (see, for instance, [1]) and those found using the normal form method, defined in [2].
The normal form method in [2] was developed as a technique that can be applied directly to systems of weakly coupled second-order nonlinear differential equations. This concept is not entirely uncommon, having previously been proposed in [3], but it is the matrix formulation proposed in [2] that is considered particu-larly beneficial to the current work. We will call this the "direct" normal form (DNF) method 1 in order to differentiate it from the "classical" method described, for example, by Jezequel and Lamarque [4], Arnold [5], Murdock [6], Kahn and Zarmi [7] and Nayfeh [8]; the latter is not investigated here, as similar comparisons have previously been made, for example, in [2].
In recent years, the DNF method (and other normal form methods similar to this) has been used extensively to capture the responses of nonlinear systems. This includes, but is not limited to, describing modal interactions and bifurcations in backbone curves [9][10][11][12], recognising out-of-unison resonance in a taut cable [13], reduced-order modelling [14], nonlinear system identification [15,16], investigating aeroelastic systems under fluid flow [17,18], exploring applicability conditions for nonlinear superposition [19], and quantifying the significance of nonlinear normal modes [20]. In contrast with the recent development of the DNF method, the MS method is well established in the literature, with thorough discussions regarding its development readily available, for example, in [21][22][23][24][25][26].
Perturbation methods require the repeated application of a number of steps, building up an increasingly accurate solution by addressing smaller terms in each repetition. In the practical application of these methods, the steps can require significant computational effort and produce increasingly complex expressions, which can, arguably, hide the mathematical insight gained from employing such a technique. In this paper, we consider the "accuracy" of these methods by assessing the result after one or two repetitions of their respective steps. It is generally recognised that these techniques converge to the correct solution with many repetitions, so can ultimately be considered as precise as each other.
A contributing factor in the accuracy of the DNF method, as shown in [27], is the frequency detuning which arises in its formulation. In physical terms, this can be interpreted as a series expansion around the natural frequency of the underlying linear system. This is not naturally present in the MS technique; however, several examples of alternative detunings, applied to MS technique, can be found in the literature [28][29][30][31][32][33]. The attempt that most closely resembles the detuning of the DNF method is found in [32], although this proposed detuning is only employed in a small number of papers, such as [34,35]. In [32], an ε-expansion is applied, not only to time, as is standard, but also to frequency. The paper presents the updated frequencyamplitude relationships and suggests that they appear more accurate, although it was not possible for this to be verified with numerical data. The motivation for expanding the frequency is solely to remove the secular terms in the response, and so the technique lacks the physical motivation that is present in the DNF method, as described in detail in the current paper.
Further attempts to detune the MS method have been proposed, though a number of these focus on the forced case in which it is common practice to perturb the forcing frequency [28,29]. A more thorough investigation is given in [30], and a comparison of the MS method and the generalised method of averaging can be found in [31]. Additionally, the detuning applied in [32] has also been applied to the Lindstedt-Poincaré method of strained parameters and the generalised method of averaging, with these detuned methods producing identical truncated results [33].
In this paper, a comparison on the DNF and MS techniques is provided, with emphasis placed on the detuning used. Specifically, in Sect. 2, the two techniques are briefly outlined and compared using the Duffing oscillator as an example system, a system which is adopted in [27,[31][32][33]. The two techniques are equated by introducing a detuning step, which is physically interpreted as a perturbation about the response frequency rather than the linear frequency, into the MS technique in Sect. 3. The detuning approach employed in the DNF method will be applied in the MS method, and it will be shown that doing so allows the two methods to be equated. By considering a more general detuning, it is shown that using MS both the fundamental and the harmonic response predictions are affected by the detuning. This is in contrast to the DNF technique, in which only the harmonic response changes. In Sect. 4, the techniques are compared for a two-mode system, where it is shown that the techniques give the same results if the MS method is modified to include the detuning. Conclusions are drawn in Sect. 5.

Approximate methods
This section introduces the DNF and MS techniques, giving an overview of how they are applied to a single-degree-of-freedom (SDOF) oscillator of the following form Here, x denotes the displacement, ω n represents the linear natural frequency, and n x (x) is a nonlinear term. For both techniques, the nonlinear term is assumed to be small. Here, this is indicated by ε, which may be thought of as a bookkeeping parameter that allows the relative size of terms to be tracked [8]. As such, ε is taken to have a value of unity, such that it does not alter the equations. The application of the techniques is described as a series of steps, with the Duffing oscillator (n x (x) = αx 3 ) being used as an example.

Direct normal form
The normal form approach is typically used to find periodic solutions to the equation of motion of a system. The objective of this approach is to apply a transform to the equation of motion to give a resonant form, in terms of transformed coordinate u, that can be solved exactly by using the following form for the solution, which assumes that system will respond as a single harmonic where u p and u m are used to denote the positive and negative parts of the exponents, respectively, and A c and ω r represent the initial amplitude and nonlinear natural frequency, respectively. Time is denoted by t and φ 0 denotes the phase of the response. Once u has been found, the harmonics of the response can be recovered using the transform equation. The DNF approach is applied to equations of motion that are expressed in the linear modal coordinates, q, where q = x for SDOF systems. This means that x could be used instead of q in the following equations. However, q has been kept to allow easier comparison with the MDOF case discussed in Sect. 4. The transform may be summarised as Here, n q (q) represents the small nonlinear terms of the untransformed equation and, as q = x for a one degree-of-freedom system, n q (q) = n x (q). In addition, the detuning ω 2 n = ω 2 r + εδ, which will later be utilised in the MS method, is applied. The harmonics are now captured by the product of h 1 , a 1 × vector of coefficients, and u * 1 , an × 1 vector consisting of all the combinations of u p and u m that arise in n q (u p + u m ); these harmonics are also assumed to be small. The method of finding the harmonics, h 1 u * 1 , and the transformed nonlinear terms n u1 u * 1 require three steps. These are now discussed, along with their application to the Duffing oscillator. An explanation of how both the steps and the detuning are derived is given in Appendix A, together with an indication of how they may be modified for MDOF systems.
By eliminating q from the original differential equation in Eq. (3) using the transform and then simplifying using the transformed equation of motion, the ε i balance equation is given by the homological equation: The excitation of these equations, which defines the vectors u * i , is given as where D{n q (u)} represents the Jacobian of n q (u) and arises from the Taylor expansion of n(q) = n(u + εh 1 u * 1 + ε 2 h 2 u * 2 ). These equations are solved using the following steps (which can be followed in Online Resource 1), first for the ε 1 equation, as illustrated below, and then for the ε 2 terms by making the necessary modification to the first step.
Step 1 N F The substitution q = u = u p +u m is made in the nonlinear term to give n q (q) = n q (u p + u m ) = n e1 u * 1 . Here, n e1 contains coefficient values and u * 1 is defined above.
For the Duffing oscillator, n q (u p + u m ) = α(u p + u m ) 3 , giving Step 2 N F Using Eq. (2), the variables u p and u m in u * 1 are written as a series of complex exponentials in time. The resulting vector is double differentiated with respect to time. The second derivative with respect to time can be expressed as a Hadamard product (•); d 2 u * 1 /dt 2 = −dd • u * 1 . Further details on this are given in Appendix A.
For the Duffing example, using Eqs. (2) and (7), u * 1 may be written as and so Step 3 N F Now, h 1 and n u1 may be found using where 1 1, is a 1 × row vector with every element being one. This expression is derived in Eq. (A.7) in Appendix A. For each nonzero element in the bracketed term, the corresponding value in n u1 is set to zero and the value in h 1 is selected to satisfy the equation.
For the zero elements in the bracketed term, the corresponding terms in n u1 are set to match those in n e1 and the h 1 terms are set to zero. The result of this is a series of coefficients representing resonant terms in n u1 and harmonic (i.e. non-resonant) terms in h 1 .
For the Duffing oscillator, Eq. (10) becomes This allows us to find the required vectors Note that the zeros on the left-hand side of Eq. (11) correspond to the resonant terms in Eq. (7) being set to zero, a feature that will also be observed in the MS method. Furthermore, it is important to note that there is some freedom of choice between the h 1 and n u1 coefficients in Eq. (11). However, one of the advantages of this method is that the non-resonant terms, and only the non-resonant terms, in u * are removed from the transformed equation of motion.
The near-identity transform to order ε 1 may now be written as From Eq. (3), along with Eq. (12), the transformed equation of motion may be written as To get the frequency-amplitude relationship for the backbone curve, we substitute the base solutions for u p and u m into Eq. (14) and then exactly balance either the e i(ω r t−φ 0 ) or e −i(ω r t−φ 0 ) terms (there are no nonresonant terms as these have been removed) to give This solution can be refined by repeating these steps, addressing the terms with increasing powers of ε in turn. While each repetition leads to a more refined solution, they becoming increasingly onerous to perform algebraically. Thus, it is desirable to approach the true solution in the smallest possible number of iterations. This basis will be used to compare the DNF and MS methods in later in the paper.
To illustrate this refinement, if the ε 2 terms are included in the near-identity transform by repeating the steps a second time, the following, more precise, solution can be obtained: As a result, the frequency-amplitude relationship will now be given by

Multiple scales
The method of multiple scales is an established technique that is discussed at length in the literature (for example, see [22,23,26,31,32] and references therein), and here we provide a brief summary of this technique to form a basis on which modifications can be discussed later.
Following this review of the method, in Sect. 2.3, solutions found using the frequency detuning proposed in [32] will be presented; a more thorough investigation is given in Sect. 3, in which a comparison will be made between this detuning and that used in the DNF method.
The approach builds on the standard perturbation method in which the response is split into a series of terms with reducing significance x = X 0 + ε X 1 + ε 2 X 2 + · · · . In MS, each of these time-dependent components are treated as functions of multiple timescales.
If these timescales are used, this response is assumed to be of the form Here, the prescribed timescales are fast time over which oscillations occur, τ = ωt, a slower time over which the amplitudes evolve, T = εt, and a timescale which is slower still, given by T s = ε 2 t. This definition of τ , which incorporates frequency, is more typically associated with the Lindstedt-Poincaré method, but is applied here to allow a simpler comparison with the DNF method. These times, τ , T , and T s , are treated as independent variables, such that derivatives with respect to t can be expressed Note that fast time-frequency, ω is typically set to the linear natural frequency ω n , such that τ = ωt = ω n t. It is this selection of fast time that is now considered, and which gives the result listed in Table 1. Substituting Eq. (19) into a general representation of an undamped, unforced nonlinear oscillator and using ω = ω n , gives where substituting Eq. (18) into the right-hand equation in Eq. (20), removing the terms of order ε 3 and higher, and balancing for ε lead to To find the solution for the components of x, firstly the ε 0 order balance in Eq. (21) is solved to give where A(T, T s ) and φ(T, T s ) are slow time-varying amplitude and phase functions, respectively, which are defined by the initial conditions of the system. This allows the ε 1 equation of Eq. (21) to be written as From Eq. (23), A(T, T s ), φ(T, T s ), and X 1 may be calculated using the following steps (as demonstrated in Online Resource 2), written at ε 1 order. These steps may then be reapplied to the ε 2 balance in Eq. (21), to find the ε 2 solution.
Step 1 M S The resonant terms, i.e. those that respond at τ = ω n t in Eq. (23), are removed and equated, writing where Res{n x (X 0 )} represents the resonant terms in n x (X 0 ).

This equation is then solved to find A(T, T s ) and φ(T, T s ).
For the Duffing oscillator example, we can write Balancing the sin(τ + φ) and cos(τ + φ) terms gives respectively. These can be solved to give where φ c is an integration constant representing a phase offset at t = 0. Hence, using Eq. (22), we can write where we have recalled that τ = ω n t and T = εt such Step 2 M S The remaining terms in Eq. (23), where NRes{n x (X 0 )} represents the non-resonant terms in n x (X 0 ) are now considered. Here the righthand side may be viewed as an "excitation" of a linear dynamic system in X 1 which can be solved to generate harmonic responses terms in x. For the Duffing oscillator example, we have where Eq. (27) has been used. Solving this linear differential equation gives Hence, the order ε 1 solution, x = X 0 + ε X 1 , is given by Here, we have written A c (T s ) = A c and φ c (T s ) = φ c as the timescale T s is not used to find the ε 1 frequencyamplitude relationship.
As with the DNF methods, these steps can be repeated for higher-order ε terms. Applying these to the ε 2 terms, the refined solution is given by Again, A c and φ c are now a constants, though these would be functions of higher-order timescales if a higher ε-order solution was being sought.

Duffing oscillator backbone curves
It is now possible to compare the expressions for the frequency-amplitude relationship derived using two repetitions of the steps in each method. These are given by At this point, it is apparent that the DNF method detunes around the square of the response frequency, whereas the MS method directly detunes ω r . The corresponding higher harmonic response amplitudes are given in Table 1. Figure 1 presents the fundamental and third harmonic backbone curves for the Duffing oscillator at ε 1 -and ε 2 -order, along with the results found using the numerical continuation software Auto 07p [36]. In this figure, the influence that the type of detuning has on the results can be clearly seen. Both at ε 1 -and at ε 2 -order the DNF curve remains close to that of the numerical solution, whereas the backbone derived from the MS method diverges from this at higher amplitudes. Considering panels (a) and (c), it is evident that the ε 2order solution remains close to the numerical curve for a greater range of amplitudes, though this introduces a hardening-to-softening behaviour at higher values of A 1 . Further, the harmonic components are poorly captured by the MS method in both panels (b) and (d).
The results displayed in Fig. 1 demonstrate the differences that can occur when a detuning is applied to the square of ω r , as opposed to directly to the linear term, and provide motivation for the application of the DNF detuning in the MS method, as described in Sect. 3.2. In particular, in contrast to the explicit form for the  Fig. 1 Comparison of first-order accurate (ε 1 ) response curves found using approximate methods and numerical continuation for the undamped Duffing oscillator in terms of a the fundamental amplitude, b the third harmonic and c other harmonics, using ω n = 1 and α = 0.5 MS relationship, the DNF method gives an implicit equation in ω 2 r . This can be easily rearranged to give a quadratic equation in ω 2 r which is easily solved and square rooted to give an explicit equation for ω r . This process becomes more complicated at higher orders of ε, at which point it is possible that either a Taylor expansion or numerical continuation would need to be used. That being said, the accuracy of the curves in Fig. 1 suggests that it is unlikely that these higher orders would be necessary to obtain a strong approximation of the true solution.

Equating the techniques
In this section, we compare the derivations of the DNF and MS approaches. To do this, we first consider frequency detuning. The importance of this step for the DNF method was assessed in [27], in which the Duffing oscillator was used to demonstrate that it is this detun-ing which increases the accuracy of the technique in comparison with the classical normal form method. In light of the fact that perturbation methods repeat a specific set of steps to find a solution, as demonstrated in Sect. 2, we consider whether the same approach may be used in the MS method to improve the agreement with the DNF method at the same number of repetitions.
It should be noted that it is possible to introduce the intrinsic time-dependent amplitudes of the MS method to the DNF technique to allow transient behaviour to be captured. This is not investigated further here, as this paper focuses on the unforced, undamped behaviour of systems.

Detuning the MS method
In the derivation of the DNF technique, a frequency detuning is applied, in which the square of the natural frequency is assumed to be detuned from the square of the response frequency such that the substitution ω 2 n = ω 2 r +εδ can be made, where δ is introduced as a detuning parameter. This is discussed in Appendix A where, for multiple degrees of freedom, the equation is written = + εΔ. This allows the linear natural frequency to be replaced with the response frequency, ω r , and a detuning term, δ, in the ε 1 relationship, Eq. (A.3), and results in coefficients in (dd − ω 2 r 1 1, ) expression that are exactly zero, see Step3 N F . This detuning has been discussed in [37], where it was shown that the detuning does not affect the frequency-amplitude relationship, but does improve the prediction of the third harmonic. The physical interpretation of this is associated with how the underlying linear system is defined-normally we consider the Duffing oscillator to have a linear stiffness term ω 2 n x (and hence a natural frequency of ω n ), but the same result can be achieved by treating the linear stiffness term as ω 2 r x and modifying the nonlinear term to compensate for this, giving αx 3 + δx, where δ is a detuning parameter. Adopting the second approach can result in a smaller nonlinear term which more closely meets the key assumption that the non-linearity is of order ε 1 . Note that this interpretation of the detuning does not specifically rely on the assumption that δ is small, provided the new nonlinear term, αx 3 +δx, remains small.

Detuned multiple scales
Now let us consider how frequency detuning, which we will view as a means to express the equation of motion in terms of ω r , may be used in a MS approach, resulting in the detuned multiple scales approach (dMS). Firstly, when selecting the timescales we set the fast time as ω = ω r and hence τ = ω r t. The result of this is that Eq. (20) is modified tö where, now, τ = ω r t, whereas previously, in Eq. (21), τ = ω n t. Now, we apply a frequency tuning to remove the ω n terms. This tuning can take a number of forms, but let us select the same detuning as used in the DNF approach, with its link to modifying the linear and nonlinear stiffness terms, and use ω 2 n = ω 2 r + εδ. Substituting this and Eq. (18) into Eq. (33) and balancing for ε i give which may be compared with Eq. (21) for the standard MS approach. The solution to the ε 0 equation is the same as before, namely Eq. (22), although note that now τ = ω r t, whereas, previously, τ = ω n t. The ε 1 equation may be solved using the steps outlined previously. Step1 M S involves balancing the resonant term using the modified equation and for the Duffing oscillator this results in These can be solved to give Note here, that the frequency shift is captured using δ, as δ is defined as the detuning parameter, and so φ(T ) is set to a constant. Using this and recalling that ω 2 n = ω 2 r + εδ result in a frequency response equation given by This is identical to the expression obtained by the DNF, as shown in Eq. (15) and Table 1.
Step 2 M S is applied to find the harmonics captured by X 1 . With the resonant terms removed, the ε 1 balance may be expressed as where X 0 = A c cos(τ + φ c ) and τ = ω r t has been used. Solving this differential equation gives Hence the order ε 1 solution, x = X 0 + ε X 1 , is given by This is identical to the response predicted using the DNF approach, see Table 1. As previously mentioned, a similar detuning of the MS technique is considered in [32], which introduces an ε perturbation, ω 2 = ω 2 0 + εω 1 + ε 2 ω 2 + · · · , to resolve the issues of secular terms in the response. 2 Once truncated to order ε 1 , this expansion can be seen to be the same as that in the DNF method, though without the physical interpretation of a series expansion about the underlying natural frequency. It should be noted that, in [32], the first term is given as a square simply because it is convenient.
Note that the steps for the dMS method are illustrated in Online Resource 3.

Comparison of detuned multiple scales and direct normal form
It has been shown that the predicted response using the DNF method can be matched by the dMS method. Now we compare these two techniques in more detail for the case where the amplitude of response is assumed to be fixed, i.e. A(T ) = A c . As with all the discussions up to this point, we will consider the ε 1 accuracy case for a SDOF system. The form of the response for the methods may be written as where the ε 0 term on the right-hand side of both equations represents the resonant response and the ε 1 term the harmonic response. The resonant response takes the form = A c cos(ω r t + φ c ). 2 Interestingly, it can be seen that ω 0 = ω n , even though the ω i terms are arbitrary in [32].
Here we have set φ(T ) = φ c as discussed in the previous section. In both cases, the expression for the response frequency is derived by considering the resonant terms in the ε 1 equation.

For MS, for the case where A(T ) = A c and φ(T ) = φ c , Eq. (35) can be simplified to give
and hence, applying this in the dMS method and using ω 2 n = ω 2 r + εδ, we get In the case of the DNF approach, the transformed dynamic equation isü + ω 2 n u + εn u1 u * 1 = 0 where n u1 u * 1 is determined using Step3 N F . This step solves (dd − ω 2 r 1 1, ) • h 1 = n q1 − n u1 by considering the elements term by term. For elements where (dd −ω 2 r 1 1, ) = 0, the resonant terms, the equation is satisfied by setting the corresponding elements in n u1 equal to those in n q1 . This is equivalent to stating that n u1 u * 1 = Res{n q u * 1 } = Res{n q (q)}. By substituting this, along with the solution for u into the transformed equation of motion and noting that, for a SDOF system, n q = n x , we can write to obtain to the same expression as dMS, Eq. (43). Now, considering the harmonic contribution, from Eq. (39), we have Recalling for the dMS technique that fast time is defined as τ = ω r t, the double derivative of X 1 with respect to τ may be written as In the case of the DNF method, the harmonic terms are found in Step 3 N F where (dd − ω 2 r 1 1, ) • h 1 = n e1 − n u1 is considered. For the non-resonant, or harmonic, elements this equation is satisfied by setting the lefthand side of the equation equal to the values in n e1 on an element-by-element basis. From the derivation in Appendix A, it can be seen that this solution arises from Eq. (A.3) and may be expressed as Recalling that, for a SDOF system, Γ = ω 2 r and that h 1 is a coefficient matrix, we can rewrite Eq. (48) as DNF: It can be seen that this is the same as the harmonic expression for the direct MS approach, Eq. (46), by following the relationship in Eq. (42) and setting u = X 0 and h 1 u * 1 = X 1 . From this, we can conclude that, at an accuracy level of ε 1 , the prediction of periodic oscillations using the DNF and MS methods can be made the same. This requires the MS technique to use τ = ω r t, as in the DNF method, for fast time and to remove ω n from the equations of motion using the frequency tuning ω 2 n = ω 2 r + εδ. As discussed in [37], this is justified based on the idea that the system can be linearised about a stiffness ω 2 r x rather than ω 2 n x to potentially reduce the size of the nonlinear term. This may be substituted into Eq. (41) to give the full solution to order ε 1 .

Alternative frequency tunings
So far in this section, we have shown that the MS and DNF techniques are equivalent, to order ε 1 , under the special conditions that the fast time is set to τ = ω r t and the stiffness term, ω 2 n x, in the equation of motion is rewritten as (ω 2 r + εδ)x, where δ can still be viewed as a detuning parameter. However, this frequency tuning approach raises the question about the predicted response when a different detuning parameter is selected.
For the case of the DNF method, this has been addressed in [37] for both single-and multi-degree-offreedom systems. Consider the arbitrary frequency tuning ω 2 n = ω 2 d + δ d , where ω d is the detuned frequency with ω d = ω r for the standard technique described in 2.1. In the MDOF notation used in Appendix A, the equivalent expression is = d + Δ d . In [37], it was shown that the prediction of the response at the fundamental frequency is independent of the chosen detuning at order ε 1 . The reason for this is that the only change to the ε 1 balance, Eq. (A.3), is that d , a diagonal matrix of ω 2 ri terms, is replaced by a diagonal matrix of ω 2 di terms. The result is that, in Step3 N F , h 1 and n u1 are now found using Eq. (10). When satisfying Eq. (10) using the arbitrary frequency tuning, we apply the rule defined in Step3 N F to entries that are approximately zero in the brackets, rather than looking for values which are exactly zero. The corresponding terms in n u1 are still set to those in n q1 . As these terms are the same as those for the case where ω d = ω r , the resulting nonlinear function in u, n u1 , also remains the same. Hence, the ε 1 order equation of motion in u, and the subsequent response at the resonant frequency, is independent of the selection of ω d . However, the harmonic response prediction is affected, as each term in this is dependent on the nonnear-zero value of the bracketed term in Eq. (10). The result is that, for the Duffing oscillator, the vector for h 1 becomes which may be compared to Eq. (12) for the case where the standard frequency tuning, ω d = ω r , is used. Thus, in contrast to the fundamental frequency response, which is not a function of the detuning parameter, the prediction for the third harmonic amplitude is dependent on the choice of detuning frequency and is given by .
(51) Figure 2 shows the DNF prediction of the response of the Duffing oscillator in terms of the first and third harmonics for a range of frequency tuning frequencies, ω d = ω r + (ω n − ω r )γ , from γ = 0, corresponding to the standard detuning used in DNF (i.e. ω d = ω r ) to γ = 1, where no detuning is used (i.e. ω d = ω n ). Figure 2a shows that the prediction of the response at the resonant frequency is robust to the choice of detuning parameter; however, the third harmonic response is affected by its choice and is better captured using the standard DNF detuning (γ = 0) than with no detuning (γ = 1). Now, considering MS, we have already seen that the selection of the fast time-frequency and the subsequent frequency tuning equation (for the case where this frequency did not match ω n ) does affect both the resonant response and that of the harmonics at order ε 1 accuracy. In general, if we write the fast time as τ = ω d t, where for dMS ω d = ω r , the ε balance equations [see Eq. (21) for MS and Eq. (34) for dMS] become where the ε i ω 2 n X i terms have been rewritten as ε i ω 2 d X i + ε i+1 δ d X i prior to the balancing using the ω 2 n = ω 2 d + δ d frequency tuning expression. In addition, the Taylor series expansion ω n = ω d +εδ d /(2ω d ) has been used to remove ω n in the slow dynamics term 2ω n X † ‡ 0 . Following this approach results in X 0 = A(T ) cos(τ +φ(t)) = A(T ) cos(ω d t +φ(t)) and, using the ε 1 equation, we find that, for the Duffing oscillator, which may be compared to Eqs. (25) and (36) for the MS and dMS techniques, respectively. Solving the differential equations in φ(T ) and A(T ) and substituting the solutions into the X 0 expression give where A c and φ c are the values of A(T ) and φ(T ) at t = 0. Recalling that δ d is defined in ω 2 n = ω 2 d + δ d , this gives the response frequency ω r = ω 2 n /(2ω d ) Writing ω d = ω r + (ω n − ω r )γ results in the response frequency equation and, from solving the ε 1 expression, the resulting harmonic response amplitude is Figure 3 demonstrates that varying γ from 0 to 1 transforms the response from the DNF/dMS to the standard MS response. For the MS technique, δ d = 0 and, hence, the frequency shift away form ω n is captured by φ(T ). However, for the dMS technique, ω d = ω r and so φ(T ) = φ c represents the fact that the X 0 response is at response frequency ω r . These represent two special cases, for a general frequency tuning with fast time τ = ω d t; the thin, green curves in Fig. 3 represent a continuum between these two cases. Note that the accuracy of the DNF method is only reached when the detuning from that method is used. Interestingly, the fundamental response is independent of the detuning for the DNF method, whereas this is not the case for MS.

Example: non-symmetric, two-mass oscillator
A 2DOF system is considered in this section, allowing the two methods to be compared using a more complex system, as well as examining the robustness of the frequency tuning methods.
The system under consideration is the same as that in [38]. It consists of a two-mass oscillator with a symmetric underlying system of linear springs; two cubic nonlinear springs are added in parallel with the corresponding linear springs, one grounding the first mass and one connecting the two masses. Therefore, the force-deflection equation for the grounding of the first mass is F = k 1 (Δx) + κ 1 (Δx) 3 and the relationship is similar for the springs connecting the masses, given by F = k 2 (Δx) + κ 2 (Δx) 3 . Here, k i and κ i are the spring constants of the linear and nonlinear springs, respectively. As in [38], both techniques are applied directly to the modal equations of motion to ease the comparison of solutions. These are given bÿ where is a diagonal matrix of the squared natural frequencies of the underlying linear system, ω 2 ni , and with β = 16 κ 2 κ 1 . The application of the methods is largely the same as for the Duffing oscillator considered in previous sections, so only a brief overview is given below. For brevity, solutions will only be considered to order ε 1 . In addition, we provide scripts, as supplementary material, that allow the equations to be derived symbolically using Wolfram Mathematica.

Multiple scales
Standard perturbations are again implemented, giving the two modal coordinates as The notation τ i = ω ni t has been introduced to ensure that the fast time for each mode corresponds to the appropriate natural frequency. It should be noted that, for this model, we consider the case τ 1 ≈ τ 2 . Implementing this perturbation, as well as the corresponding adaptation of the derivative given in Eq. (19), results in zeroth-and first-order perturbation equations that take the same form as in Eq. (21), and hence, the former can be solved to give where A i, j denotes the amplitude of the jth harmonic in the ith mode. These solutions can be applied to the first-order equation [equivalent to the SDOF equation in Eq. (21)] to give the ε 1 equations The nonlinear terms, n qi (Q 10 , Q 20 ), are lengthy and contain fundamental and harmonic terms. Therefore, they are not shown here.  Collecting the resonant terms allows the amplitudes and phases to be calculated, as in the SDOF case, so that 2ω n1 A 1,1 (T ) ‡ sin(τ 1 + φ 1 (T )) These equations can now be solved to give: where β = 16κ 2 κ 1 . Note that the expressions for phase enforce the condition that neither fundamental amplitude can be equal to zero. Therefore, recalling that τ 1 = ω n1 t results in This leads to the following compatibility condition This expression can now be used to find A c1 in terms of A c2 , or vice versa. However, the explicit solution is non-trivial and is not shown here. The resulting backbone curves from Eq. (64) are given in Fig. 4 and discussed in Sect. 4.3. Due to the involved process required to find the harmonics, analytical solutions for these are not given, but have been derived using Wolfram Mathematica and solved numerically to allow comparison between the techniques; this is discussed in Sect. 4.3.

Direct normal form
This technique also closely mirrors its SDOF counterpart, so only a brief description of the key differences is given. The resonant equations of motion are once again found in terms of u, with Steps1 N F -3 N F are followed as previously described and are not shown here due to the large size of the matrices involved. Full workings are shown in [38]. The frequency-amplitude relationships that arise are given by These results are comparable with Eq. (64) and the resulting backbone curves are, again, displayed in Fig. 4. Similarly, the equations for the harmonics are algebraically complex and are solved numerically.

Detuned multiple scales
The key difference when applying the dMS in two degrees of freedom (2DOF) is that separate frequency tunings are to be applied to each mode Again, the resonant equations are used to find the amplitude, phase, and now detuning parameter. For the 2DOF case under consideration, these are given by Thus, substituting these values into Eq. (68) gives the frequency-amplitude equations as Comparing these with Eq. (67) demonstrates that the results from the dMS method, once again, match those from DNF. It should be made clear that, in Figs. 4 and 5, the curve for the dMS method has not been printed as it is coincident with the DNF curve.
As with the SDOF case, the final forms of q i are identical, though this is not shown here for reasons of brevity.

Comparison of the techniques
The fundamental backbone curves for the first and second modal responses are given in Fig. 4. Four backbone curves are shown for each technique. Panels (a) and (b) correspond to the first backbone curve of the system, that is, the curve which initiates at the first natural frequency of the underlying linear system, ω n1 ; panels (c) and (d) represent the second backbone curve. These results are comparable to those for the Duffing oscillator in Fig. 1, with the MS curve underestimating the numerical continuation results and the DNF/dMS results again remain closer to the numerical continuation results. The difference between the methods grows significantly with increasing amplitude. In particular, the MS results diverge noticeably from the numerical and DNF/dMS counterparts at higher amplitudes. As verified in [38], this is the result of the loss of influence of the higher-order terms during the linearisation of the system.
Interestingly, the third harmonic components of the backbone curves in Fig. 5 are qualitatively different from the equivalent curve for the Duffing oscillator. While the amplitudes of the third harmonics from the MS method in the SDOF case were greater than those from numerical continuation, Fig. 5 shows that the opposite is true for the 2DOF responses. This inconsistency suggests that the MS method is less robust to changes in the system compared to the DNF and dMS methods, which remains consistent across the two cases, although higher-order cases have not been considered in this study.

Conclusions and discussion
This paper presents a comparison between the multiple scales and direct normal form techniques and investigates whether the two methods can produce equivalent results. In particular, the detuning used in the DNF method was applied in the MS method to investigate  whether a similar level of accuracy could be achieved. The frequency detuning, which can be physically interpreted as a way of reducing the amplitude of the nonlinear term based on adapting the effective linear stiffness, is inherent in the DNF method and has been shown to improve the prediction of the harmonic response content. In applying this detuning in the MS method, it was shown that the two methods could be equated, giving identical solutions up to ε 2 order.
The DNF is advantageous insofar as a natural detuning approach is intrinsic in its formulation, whereas this is not the case for the MS technique. It is, therefore, the decision of the user as to whether a detuning is utilised to increase the accuracy of the method. Furthermore, it has been demonstrated that the fundamental response prediction is robust to changes in detuning in the DNF method. Since this is not the case for the MS technique, we observe that there is room for further optimisation of the detuning to be applied, which could further increase the accuracy of the method.
To aid the understanding of these methods, as well as the differences in their implementation, Wolfram Mathematica files for the 2DOF case have been provided as open access data files. These closely follow the steps defined in Sect. 2 and are designed to be used in conjunction with this paper to give a practical understanding of each procedure. Here, Λ is a diagonal matrix with the i th diagonal element being the square of the i th linear natural frequency, ω 2 ni . This reduces to Eq. (3) for the case where a single mode is considered, with q = q, u = u, n = n and Λ = ω 2 n1 = ω 2 n to indicate that the terms are now scalar quantities or functions.
To find vector u * 1 , we make use of the fact that in the transformed equation of motion the harmonic terms have been removed such that the response in the i th coordinates u i may be expressed as A ci 2 e i(ω ri t−φ 0i ) + A ci 2 e −i(ω ri t−φ 0i ) .

(A.2)
Using this the nonlinear function n q expressed in terms of u may be written as n q (u) = n q (u p + u m ) = n e1 u * 1 . Here u * is an × 1 vector consisting of all combinations of u pi and u mi generated by n q (u p + u m ) and n e1 is a n × matrix (for an n degree-of-freedom system) that contains coefficient values. The subscript e in n e1 indicates that this term can be thought of as an excitation term in the following discussion. Now, considering Eq. (A.1), the transform can be substituted into the equation of motion in q. The resulting expression can first be simplified using the resonant equation of motion and then balanced in terms of ε to give where Γ is a diagonal matrix with the i th diagonal element being the square of the i th response frequency, ω 2 ri . Here, a Taylor series expansion has been used on the nonlinear function: n(q) = n(u) + O(ε). In addition, a form of frequency tuning, similar to that used in [37], is employed: we write Λ = Γ + εΔ when considering the h 1 u * 1 terms. Now, observing that each element in u * 1 is made up of u p and u m elements which themselves are complex exponentials in time, if follows that each term in u * 1 may be written as complex exponentials in time. Therefore, when differentiating u * twice with time, each element maps onto a scaled version of itself. So, we may representü * 1 as d 2 u * 1 d t 2 = −dd • u *

(A.4)
where dd is a vector of length × 1 and • is the Hadamard product (element-wise matrix multiplication). Using this, we can write the first term in Eq. It can be seen that this reduces to the equation given in Step3 N F of the normal form description for the SDOF case (n = 1). Using Step3 N F , the equation can be solved to find h 1 and n u1 and hence identify the transform and transformed equation of motion, respectively.