Using frequency detuning to compare analytical approximations for forced responses

It is possible to use numerical techniques to provide solutions to nonlinear dynamical systems that can be considered exact up to numerical tolerances. However, often, this does not provide the user with sufficient information to fully understand the behaviour of these systems. To address this issue, it is common practice to find an approximate solution using an analytical method, which can be used to develop a more thorough appreciation of how the parameters of a system influence its response. This paper considers three such techniques—the harmonic balance, multiple scales, and direct normal form methods—in their ability to accurately capture the forced response of nonlinear structures. Using frequency detuning as a method of comparison, it is shown that it is possible for all three methods to give identical solutions, should particular conditions be used.


Introduction
Modelling the dynamical behaviour of nonlinear structures remains an active area of research across the engineering disciplines, and one which becomes increasingly relevant as the push for more lightweight and efficient structures continues. While the complexity and variation in such structures are constantly increasing, the development of mathematical methods for assessing the equations of motion of these systems continues to be an important branch of research. As has been addressed in the literature [1][2][3][4][5][6], there are two main ways in which the free and harmonically forced response of a nonlinear system can be approximated. Numerical continuation methods [2] can be used to follow the entire set of periodic responses over a region of the parameter space. Alternatively, it is possible to apply analytical methods to develop an approximate solution for the dynamical response. In general, while the former can readily produce the forced and free response curves of the system [7][8][9], the analytical methods are able to provide additional insight into the nonlinear behavioural characteristics displayed by the structure (see, for instance, [10]). In addition, the numerical continuation technique is limited by the need for an initial solution. In particular, structures which exhibit behaviour such as isolas lack this initial point, hence benefit from the use of analytical techniques in conjunction with numerical methods, as shown in [11]. Moreover, the computational cost of response contin-uation can greatly increase for larger systems, which limits this numerical strategy to low-order applications.
In this work, the accuracy of three analytical methods in capturing the behaviour of a nonlinear system is assessed. Namely, these are the harmonic balance (HB), the multiple scales (MS), and direct normal form (DNF) techniques. These methods were chosen due to their wide application in the literature and, therefore, the results provide the user with the make a more informed decision regarding the accuracy and applicability of the methods. It should be noted that these techniques do not represent the entire field of analytical methods for approximating nonlinear systems; for example, the more mathematically rigorous spectral manifolds approach may also be applied (see [12]) A more complete review of the associated literature is presented in [3]. However, as motivation for their comparison, we present some prominent studies including each of the considered techniques. The HB balance method has been widely used to identify nonlinearities in mechanical structures [13][14][15], as well as calculating their nonlinear normal modes [16]. Furthermore, in [17], it has been shown that the relative simplicity of the method permits the use of large number of harmonics, allowing the accurate prediction of the behaviour of complex structures. There are a number of ways in which the MS method can be applied, as has been discussed in [18]. It has been widely used in the prediction of nonlinear dynamic behaviour, such as bifurcations in the frequency-amplitude relationship [19][20][21][22] and internal resonances [23][24][25]. These resonances have also been investigated using the DNF method [26][27][28][29], which has been further applied to explore the significance of nonlinear normal modes in relation to forced cases [30] and identify nonlinearity in structures [31,32].
This paper expands the recent discussion by the authors in [3], which considered the calculation of the free response of two nonlinear mass-spring systems using the MS and DNF techniques. The authors applied the frequency detuning from the DNF method, as explored in [33,34], in the MS technique to show that it is possible for both methods to give identical solutions. In the current study, this discussion is expanded in a number of ways, though the desire to achieve accurate results with minimal analytical effort remains a priority. First, the widely applied harmonic balance (HB) method is included in the discussion, providing a more exhaustive assessment of analytical approximation methods. Furthermore, the three methods are now applied to forced, damped systems, allowing the validity and applicability of the discussion to be expanded to situations which more closely represent real-world engineering applications. Finally, the variable nature in which the MS method can be applied is addressed by using the derivative expansion version [35,36], as opposed to the two-timing variant [1,37]. This paper is structured as follows: Sect. 2 provides an overview of the analytical techniques considered; these methods are then applied to a Duffing oscillator in Sect. 3; frequency detuning is used to compare the methods in Sect. 4; graphical results are presented in Sect. 5; and conclusions are drawn in Sect. 6.

Overview of techniques
In this paper, the techniques are initially presented with reference to a nonlinear, dynamical system, in which it is assumed that the damping, forcing, and nonlinear terms produce relatively small contributions to the system behaviour. The equations of motion of such a system are defined by where M, C, and K define the linear N ×N mass, damping, and stiffness matrices, respectively, x is an N × 1 vector of the nonlinear terms, and the dot notation represents derivatives with respect to time, t. In addition, r T = [r p , r m ] = [e +jΩt , e −jΩt ] represents the periodic nature of the forcing applied and P x = P x 2 ,P x 2 , wherê P x is an N × 1 vector of scalar terms which define the magnitude of this forcing. The bookkeeping parameter, ε-used to denote the relatively small nature of the damping, nonlinear, and forcing terms-is bracketed to denote the fact that it is not necessarily used in the HB method.
The system of equations in Eq. (1) expresses the dynamics in terms of the physical coordinates of the system. However, it is also possible to express them in terms of the modal coordinates, by defining the coordinate transform x = q, where is the matrix of mode shapes, found by solving the eigenanalysis problem defined by = M −1 K . This transformation is applied in this paper as it allows the comparison to be followed more simply and relates directly to the Duffing oscillator example, as discussed in Sect. 3.
Applying the modal transformation and premultiplying the resultant equation by T gives Therefore, by pre-multiplying Eq. (2) by ( T M ) −1 , the modal equations of motion can be simplified as where q is an N ×1 vector of modal contributions, = ( −1 M −1 K ) is a diagonal matrix of squared frequency terms, ω 2 n,i , and P q = −1 M −1 P x is the modal projection of the forcing amplitudes. q (q,q, r) = is a function that contains the nonlinear terms, similarly to x in Eq. (1), but is here expanded to include damping terms.

Harmonic balance method
Widely used across the literature [13][14][15][16][17]38,39], the HB method begins with the assumption that the ith displacement, q i , can be expressed in the form where n denotes the number of harmonic terms to be included in the trial solution, A i and φ i define the ith amplitude and phase, respectively, and the notation• represents the complex conjugate. ω r denotes the response frequency; note that this is the true frequency at which the system will oscillate, and this may differ from the forcing frequency, Ω. Thus, this assumption states that the time-varying displacement can be expressed as a series of harmonic terms which contain frequency content at integer multiples of ω r,i . The displacement expression defined in Eq. (4) can be applied in Eq. (2) to give a system of polynomials in terms of e +jΩt , effectively collecting the terms that respond at each frequency. The terms of these polynomial equations are then functions of A i and φ i . Since these equations are time-dependent, the coefficients of the harmonic terms must be balanced, as suggested by the name of this method. In fact, this balancing of harmonic terms is an integral step of both the MS and DNF methods, as both require the assumption that the coefficients of like harmonic terms must be equal.

Multiple scales technique
As previously mentioned, the derivative expansion version of the MS is used in this paper [1,37]. This version begins by first perturbing the time and displacement as respectively. From the latter of these, it is possible to define the derivatives with respect to time, t, as d dt where D n denotes the derivative with respect to T n . Implementing these definitions in Eq. (3) leads to the equation This step is explored further in [40]. To resolve the complex q term, it is useful to apply a Taylor expansion, so that the balanced exponents of Eq. (7) can be written as where ∂ ∂q 0 q denotes the vector of partial derivatives of q with respect to each of the elements of q 0 .
The treatment of this ε-balance is more thoroughly detailed in "Appendix A", with only the key steps given in this section. The initial step is to solve the ε 0 -order equation to give where A i and φ i denote the ith displacement and phase components, as they were in the HB method. However, unlike the previous technique, both of these parameters are now dependent on the faster timescales, T 1 , T 2 , . . .. Once this solution has been established, it can be implemented in the ε 1 -order equation. In doing so, it must be noted that the homogeneous form of this expression is identical to that used to find q 0 . Therefore, those terms which respond at ω r,i -referred to as secular terms-must be set to zero. The reason for doing so is that the homogeneous forms of the ε 0 -and ε 1 -order expressions are the same, so these secular terms would lead to a divergent solution in the latter; further details are given in "Appendix A". The solution of this expression leads to the frequency-amplitude relationship.
Then, only the non-secular terms are used to find a solution for q 1 by considering the updated expression where NRes{•} denotes the non-resonant component of •.
This process can then be repeated to find higherorder displacement expressions.

Direct normal form technique
Having applied the modal transform to the equations of motion in Eq. (3), two further transforms are applied as part of the DNF method; these can be summarised as [1,33]: -Forcing transform -The terms in the response which respond at frequencies close to the forcing frequency are isolated.
-Nonlinear near-identity transform -This transform is used to remove those nonlinear terms which are non-resonant.
The forcing transform takes the form where v is the transformed state of q and [e] is an N ×2 matrix used to isolate the resonant forcing terms, an explicit definition of which is given subsequently.
where W is a 2 × 2 diagonal matrix with entries +jΩ and −jΩ. and r are defined as in the previous section. This step is used to monitor which terms are close to resonance, which, for the kth mode, will be taken to mean ω n,i ≈ Ω. Thus, the kth row of P v can be explicitly written as It is desirable to write Eq. (11) in the form where any non-resonant forcing terms have been removed. This can be achieved by imposing the definitions Wr, r) and Applying the definition of W, the ith row in the latter of these definitions can be expressed as Combining this expression and Eq. (12), it is possible to explicitly define the kth row of [e] as Now, the nonlinear near-identity transform is applied, allowing the response to be separated into its fundamental (u) and harmonic (h) components. This is achieved by implementing the following expression for v: Similarly to the nonlinear and damping terms, the harmonics have been perturbed, formalising the assumption that higher-frequency harmonics will have less influence on the behaviour of the system. Furthermore, it is assumed that the expression for u is sinusoidal and will be written in exponential form, with separate vectors for the positive and negative exponents; this is expressed as u = u p +u m , where the subscripts denote the sign of the exponent. Therefore, the ith element of u is now written as where A i , φ i , and ω r,i denote the amplitude, phase, and response frequency of u i , respectively. As in Eq. (13), the desired form for the equations of motion in this step, with the non-resonant nonlinear terms removed, is given bÿ Following the procedure defined by [1], it is possible to implement Eq. (17) in Eq. (13), and then apply Eq. (19) to give the simplified, ε 1 -order expression This equation applies the perturbation u (u,u, r) = u,1 (u,u, r) + ε u,2 (u,u, r) + · · · and notes thaẗ u = ϒu, where ϒ is a diagonal matrix with ith term ω 2 r,i . Furthermore, the following detuning expression has been employed: a step which is typically applied in the DNF method and is explored further in [33,34]. The first equality here demonstrates the fact that instead of simply detuning around the natural frequency, we detune its square, as this is the form in which it arises in the equations of motion. The second expression makes use of the fact that although is not necessarily equal to ϒ, the dynamics are considered within some close neighbourhood of the natural frequency, so their difference will be small (and hence of order ε).
Balancing the ε 0 terms in Eq. (20), it can be seen that P u = P v . In the ε 1 -order balance, the u,i terms are used to manage resonant terms, which respond at ω r,i . These expressions typically comprise polynomials in terms of u p = {u pi }, u m = {u mi }, and r. The procedure for deriving the fundamental and harmonic responses from these polynomials is more thoroughly outlined in "Appendix B", but is briefly outlined here.
In the DNF method, and in contrast with the other methods discussed in this paper, the matrix notation u * i (u p , u m , r) is introduced. This expression is simply an N i × 1 vector including all the possible polynomial terms which could arise in the aforementioned expansion of u,i . This vector can then be pre-multiplied by a matrix of time-independent coefficients for each entry in u * i ; a more detailed explanation of this is given in "Appendix B".
Then, the difference between the frequency of these entries and the response frequency is captured by the introduction of the matrix β i , with element {k, } defined by Once this has been done, it is possible to express the frequency-amplitude relationship as as derived in "Appendix B". As well as finding this relationship, it is possible to calculate the harmonic response using Eq. (17), so that the physical response of the general nonlinear system is

Application to a Duffing oscillator
This section focuses on the application of these analytical approximation methods to the single-degree-offreedom (SDOF) Duffing oscillator, with equations of motion Although relatively simple in nature, the Duffing oscillator has been widely used in the discussion and comparison of analytical methods. In this study, it provides an example application of the discussion above that is easily accessible to the reader; this could be expanded to include more complex structures in future works, though the conclusions drawn would be similar to those in this paper. Given there is only one DOF, the system is automatically written in modal coordinates, allowing a cosmetic change to givë q + 2(ε)ζ ω nq + ω 2 n q + (ε)αq 3 = (ε)P cos(Ωt). (26) As such, it is possible to directly follow the modal application of the analytical techniques, as defined in Sect. 2.
In this comparison, the solutions will be found up to ε 1 -order. The ε 2 -order solutions are qualitatively similar to those shown in this section, but are more algebraically complicated and, therefore, are not presented.

Harmonic balance method
In this example, across all three methods, it will be assumed that only a single harmonic term will be necessary in the trial solution and that the response frequency will be equal to the forcing frequency. As such, Eq. (4) can be written as where A and φ denote the amplitude and phase of the mode, respectively. By implementing this displacement expansion in Eq. (26), the fully expanded form for the equations of motion can be expressed as It can be noted that the terms in the expansion which respond at 3Ω are assumed to be negligible and are, therefore, removed from the equation. This is a direct consequence of the use of a single harmonic in the HB trial solution. The solution of the Duffing oscillator using the HB technique is readily found in the literature [41], so is not expanded here. Instead, only the key results are presented: Free response

Multiple scales technique
In this section, the perturbed displacements, timescales, and derivatives are given by respectively, when truncated at ε 1 -order. Therefore, for the Duffing oscillator, the ε-expansion in Eq. (8), up to the first order, is given by Separating the secular and non-resonant terms in the ε 1 -order expression in Eq. (34) leads to the two expressions The first of these equations can now be used to find the frequency-amplitude relationship, and the second to find the displacement expression. Note that, by perturbing the displacement in the MS method, the implicit assumption of negligible harmonics seen in the HB technique has been removed. Instead, this is explicitly imposed by the inclusion of the bookkeeping parameter, ε. For Eq. (35) to represent steady-state dynamics, it is important to guarantee that there are no changes in amplitude and phase with respect to T 1 . Doing so is relatively straightforward for A(T 1 ), but the phase term must be considered in conjunction with a detuning parameter. Since the forcing frequency will be considered in some neighbourhood of the natural frequency, it is common practice to define Ω = ω n + εσ , where σ is an arbitrary detuning parameter. The inclusion of this detuning means that Ω T 0 = (ω n +εσ )t = ω n T 0 +σ T 1 .
Here, it must be noted that the use of a detuning parameter is inherent in this application of the MS method; this is utilised in later sections. Now, it becomes necessary to define the T 1 -dependent linear transformation of the phase angle Thus, the conditions for steady-state behaviour are Solving these equations, it can be concluded A(T 1 ) = A is constant, and that D 1 φ(T 1 ) = σ . Equation (35) can now be rewritten in terms of its real and imaginary components to give the system As with the HB method, these expressions can be squared and summed to eliminate the phase term, φ. Recalling that the detuning parameter is given by σ = 1 ε (Ω − ω n ), the final expression for the forced response is given by and the free vibration is defined by In addition, the phase in the forced case can be expressed as It can be further noted that the bookkeeping parameter remains present in the final expression, demonstrating that the relative contributions of the terms are monitored through the entire process. Furthermore, solving Eq. (36) allows the ε 1 -order displacement solution to be given by where it can be seen that an approximation to the third harmonic is captured.

Direct normal form technique
In the SDOF example, the three aforementioned transforms can be summarised as Here, the forcing transform is unity, as there is only a single mode to consider and the forcing is assumed to be near-resonant for this mode. However, the nearidentity transform still allows the separate consideration of the fundamental and harmonic components of the response. As shown in Eq. (18), It has been demonstrated, in the previous section and "Appendix A", that the frequency-amplitude relationship of the DNF method is defined by Eq. (23). Considering that the linear natural frequency is known, it is only necessary to calculate the nonlinear coefficients, Γ + u , and vector u * to define the system dynamics. Both of these are found in the expansion of Γ q (q,q, r) = Γ v (v,v, r), as found in Eq. (11). This is given by By applying Eq. (45) and converting to matrix notation, it is possible to express Eq. (46) as Recalling the expression in Eq. (22) and the subsequent separation of resonant and non-resonant terms, the vectors β, [Γ u,1 ], and [h 1 ] for the Duffing oscillator are given by These vectors can now be directly applied in Eq. (23).
Treating the real and imaginary parts separately and reconciling allows the system response to be approximated. This can be summarised as follows (52)

Comparison through frequency detuning
As previously mentioned, the discussion of the methods thus far has been provided so that they can be more easily compared. The system response predicted by the three methods is summarised in Table 1. By removing the bookkeeping parameter, it can be seen that the expressions for the free and forced responses in the HB and DNF method, in terms of the amplitude of the response at fundamental frequency, A, are identical to one another.
It has been highlighted, in [4], that the free response of the MS is a linearization of that found in the HB and DNF methods, as can be found via the use of a Taylor expansion. Here, it can be seen that this remains true in the forced case. Thus, the MS solution can be considered as an approximation to the DNF and HB expressions, a point which was demonstrated to lead to divergent backbone curves in [3].
However, considering the displacement expressions in Table 1, it can be seen that, by applying the same number of iterations of each method, the DNF and MS methods provide insight into the harmonic behaviour. This will not be provided by the HB method unless the trial solution is expanded to include a higher-order term.
In this study, the DNF frequency detuning, as given in Eq. (21), is applied in the MS method, resulting in the two methods giving identical expressions for the free response. In the present work, as well as considering a more general implementation of the MS method, the forced responses are considered, as outlined above.

General detuning comparison
As has already been shown, the use of a detuning is common practice in the MS method. Now, in this section, we introduce the specific detuning used in the DNF method, presented in Eq. (21), to the MS method. Recall that this is given by Applying this detuning in the general ε-expansion given in Eq. (8) leads to the updated expressions Once more, the secular terms must be set to zero and the trial solution in Eq. (9) is again implemented, giving Comparing this with the expression for the standard MS method (given in "Appendix A"), it can seen that, by removing the assumption that the response frequency is equal to the linear natural frequency, new terms arise in the ε 1 -order equation. The inclusion of the term (ω 2 n,i − ω 2 r,i ) A i e +j(ω n,i T 0 +φ i ) +e −j(ω n,i T 0 +φ i ) can be thought of as a detuning term that accounts for the influence that those terms which are close to resonance have on the vibration of the system. Collecting coefficients of e ±jω r,i T 0 , Eq. (55) can be written as Similarly to the DNF method, the terms in Eq. (56) inside the square brackets are complex conjugates, so both must be equal to zero to remove the secular terms. As such, the frequency-amplitude equation can be written as Recall that the equivalent expression in the DNF method, Eq. (23), is written as Further consideration is also given to Γ + u,i , which denotes the resonant terms of u , including both the damping and nonlinear terms. In the detuned multiple scales (dMS) case, the resonant terms of q are collected by the term However, it must be noted that there is a remaining additional term in Eq. (57), namely the term jω r,i ( Although the inclusion of an extra term may appear to nullify the equivalence of the two techniques, it can actually be noted that this expression is a damping term that corresponds to the higher-order derivatives. Therefore, it simply accounts for the perturbation of t, which is not present in the DNF method and the combination of these terms is directly equivalent to Γ +

u,i
The solution for q 1 can be determined in a similar way to the standard MS method, as in Eq. (10), by solving Recall that the corresponding ε 1 -order equation for the DNF method is It can be further noted thatḧ 1 = − h 1 and that v − n u1 = βh 1 = −NRes{ q }, so that the expression in Eq. (59) can be rewritten as It is immediately clear that the solutions of Eqs. (58) and (60) must be identical.

Duffing oscillator detuning comparison
These new steps can now be explored for the Duffing oscillator, beginning with the updated ε-balance, which is expressed as Note that the δ notation (the SDOF equivalent to in Eq. (21)) will be used until the calculation of the final solution, so that the influence of the detuning is more readily tracked. As above, the steady-state dynamics are of interest, so the conditions in Eq. (38) are applied. Separating the real and imaginary parts now leads to the following steady-state equations Therefore, by squaring and then summing these expressions, then introducing the definition of δ from Eq. (53), the free and forced responses are now given by respectively. Comparing these expressions with those in Eq. (49), it can be seen that the dMS and DNF forced responses are now identical, as is the case with the free response. The harmonic part of the displacement can now be calculated using Eq. (58), giving an expression identical to that of the DNF method, given in Eq. (52): Table 1; the bookkeeping parameter is simply removed from the expressions, so the HB curve is identical to that given by the DNF/dMS methods. The 'Numerical' backbone curve in Fig. 1 has been found using numerical continuation [2]. In this figure, the following parameter values have been used: ω n = 1, α = 0.6, with ζ = 0.005, P = 0.0015 in the first case and ζ = 0.0015, P = 0.005 in the second. These two forcing cases allow the influence of this detuning to be considered at both high and low amplitudes. In the first case, it can be seen that the variation between the techniques is very small. In fact, when the system is forced at low levels, it could be argued that the choice of method has a negligible influence on the response prediction.

Figure 1 shows graphical representations of the responses in
In any analytical approximation method that includes a perturbation or series expansion, it is possible to increase the number of step iterations or terms included so that the predicted solution will better agree with the true response. However, as has been previously discussed, the methods are compared here with the use of only a small number of iterations, as to add any more would create analytical complexity, which would add little further insight than the numerical technique. When the ratio of forcing to damping is increased, the differences between the solutions at ε 1 -order become much more noticeable. Both methods diverge from the numerical results, but the key difference is the rate at which they do so. While the DNF/dMS methods remain close to the numerical solution across the considered range, the MS method begins to diverge at approximately Ω = 1.1 Hz, which represents a shift in frequency of roughly 10%. When this point is reached, the MS method begins to underestimate the displacement. It has been shown that, by choosing the detuning found in the DNF method, it is possible to remove the underprediction of the forced response; in Fig. 1, the updated MS curve is identical to that of the DNF method.

Considerations for users
As has shown in the previous sections, by choosing a suitable frequency detuning, it is possible to achieve identical results when applying any of the techniques discussed to the example system. As such, it could be concluded that the user may choose the technique with which they are most familiar and ensure that the accuracy of the solutions will not be compromised. While there may be some truth to this assertion, the discussion provided in this study alludes to further considerations that could be taken.
The HB method has been demonstrated to be a much simpler technique than the others considered; it is arguably for this reason that it has been more widely applied, particularly in those cases that have seen a more algorithmic approach taken for larger systems. As a result of this simplicity, the method lacks the bookkeeping parameter used by the MS and DNF techniques. As has been previously discussed, this perturbation strategy not only measures the relative contribution of each term, but actually informs the user as to the neighbourhood of to linear natural frequency in which they may have confidence in the accuracy of the results. Thus, the HB should be seen as an effective method for gaining a quick understanding of the nonlinear behaviour, but is likely to be most useful when the exact values of the solution are not important.
The results presented in this paper, along with those of [3], have demonstrated that the MS and DNF are able to provide identical results. Given that the implementation of the two techniques requires a similar number of steps, this leaves the user to make a decision based on the nuances of the particular system they are considering. A key benefit of the MS has always been its ability to model transient behaviour. Since this is not affected by the frequency tuning, this asset still holds true. Meanwhile, the DNF technique utilises a matrix formulation to effectively track the resonant and nonresonant terms, a strategy that may prove useful in more complex structures with multiple harmonics. In summary, the results of the present work provide the user with two accurate options, which will suit a wide array of potential structures.

Conclusions
This paper presents and compares three analytical approximation methods, both in terms of a general nonlinear dynamical system and the Duffing oscillator, which has been widely used in the development and study of such techniques. In particular, the harmonic balance, multiple scales, and direct normal forms methods have been considered. In an initial comparison, it has been shown that the HB and DNF give identical solutions for the Duffing oscillator, although the inclusion of a bookkeeping parameter in the latter provides the potential for a simpler management of weak terms in more complex structures. Additionally, the matrix formulation of the DNF method also reveals the harmonic components of the behaviour, whereas this requires the introduction of a more complex trial solution in the HB method. In previous works, it has been shown that the free vibration expression given by the MS is a linearization of the HB/DNF solution. In this work, it has been shown that this property also holds true for the forced response.
To expand this comparison, the previously developed detuned multiple scales method is extended to include forced responses. Furthermore, an alternative form of the MS method is applied, demonstrating the broader applicability of the comparison. By considering the widely applied derivative expansion version of the MS method, it has been shown that, for a general nonlinear system, the analytical solutions of the forced, damped equations are, once more, identical in the DNF and dMS methods. This is an important result, as it means that, regardless of the method chosen, the user has the potential to capture complex nonlinear behaviour in forced structures (such as hysteresis and internal resonance) to the same level of accuracy. Particular attention has been given to the influence that the forcing and damping terms have on the divergence of the MS and DNF/dMS responses. In models with lower forcing and/or higher damping, it has been seen that the difference between these solutions is effectively negligible. However, as the forcing grows, the non-detuned MS method begins to under-predict the displacement, demonstrating that there is a region in which the linearization is valid and that the prediction may become inaccurate away from this.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A
Once the ε-expansion in Eq. (8) has been established, the ε 0 -order can be solved to give where A i and φ i denote the amplitude and phase of the fundamental response of the ith mode, respectively. This ε 0 -order solution can now be applied to the ε 1order terms in Eq. (8) to give e +jΩt + e −jΩt + P q,i r p .
Now, the secular terms-i.e., those that respond at ω n,i -must be set to zero. The reason for this is that the homogeneous form of Eq. (A.2) is given by (D 2 0 + )q 1 = 0. As such, the trial solution must be equal to the trial solution for q 0 . Therefore, a failure to remove the secular terms would lead to a divergent solution. Removing the secular terms in Eq. (A.2) leads to the expression jω n,i (D 1 A 0,i + A 0,i D 1 φ i )e +j(ω n,i T 0 +φ i ) A 0,i 2 e +j(ω n,i T 0 +φ i ) + e −j(ω n,i T 0 +φ i ) , jω n,i A 0,i 2 e +j(ω n,i T 0 +φ i ) − e −j(ω n,i T 0 +φ i ) , e +jΩt + e −jΩt − P q,i r p = 0, (A. 3) where Res{•} denotes the resonant terms of •. This secular equation can now be solved to find expressions for the amplitude and phase terms. This is done by separating the real and imaginary parts of the equation, then setting the coefficients of e ±j(ω n,i T 0 +φ i ) in each of those to zero. The solutions for A 0,i and φ i can be applied in the non-resonant equation: where NRes{•} denotes the non-resonant terms of •. This equation can then be solved to find q 1 . Now, the solutions for q 0 and q 1 can be applied to the ε 2 -order expression to find q 2 . The process is the same and, therefore, not shown here, but details can be found in [3].

Appendix B
Recall that Eq. (20) is given by (εḧ 1 + ε 2ḧ 2 + · · · ) + (εϒh 1 + ε 2 ϒh 2 + · · · ) + ε v (u + εh,u + εḣ, r) − (ε u,1 + ε 2 u,2 + · · · ) + P u r − P v r = 0. (B.1) Here, we apply a Taylor expansion to v (u + εh,u + εḣ, r), and simplify using the notation Then, the vector u * i , which captures all the possible polynomial terms which could arise in the expansion of u , is used to introduce the following matrix formu- By consideration of the u i vector, it is possible to find a solution for u. First, the th element of u * i is written as  where P + ui and P − ui denote elements {i, 1} and {i, 2} of P u , respectively. The variables Γ + ui and Γ − ui arise in the decomposition Γ ui = Γ + ui e +jω r,i t + Γ − ui e −jω r,i t . (B.10) The terms in Eq. (B.9) which are contained in square brackets are complex conjugates of one another and, for this equality to hold, it is necessary for both of these to be equal to zero. Therefore, the frequency-amplitude relationship is given by which can be solved to find the fundamental response of the system.