A combined FD-HB approximation method for steady-state vibrations in large dynamical systems with localised nonlinearities

The approximation of steady-state vibrations within non-linear dynamical systems is well-established in academics and is becoming increasingly important in industry. However, the complexity and the number of degrees of freedom of application-oriented industrial models demand efficient approximation methods for steady-state solutions. One possible approach to that problem are hybrid approximation schemes, which combine advantages of standard methods from the literature. The common ground of these methods is their description of the steady-state dynamics of a system solely based on the degrees of freedom affected directly by non-linearity—the so-called non-linear degrees of freedom. This contribution proposes a new hybrid method for approximating periodic solutions of systems with localised non-linearities. The motion of the non-linear degrees of freedom is approximated using the Finite Difference method, whilst the motion of the linear degrees of freedom is treated with the Harmonic Balance method. An application to a chain of oscillators showing stick-slip oscillations is used to demonstrate the performance of the proposed hybrid framework. A comparison with both pure Finite Difference and Harmonic Balance method reveals a noticeable increase in efficiency for larger systems, whilst keeping an excellent approximation quality for the strongly non-linear solution parts.


Introduction and literature review
The literature provides a variety of methods for the direct numerical approximation of steady-state solutions of nonlinear systems. These methods typically form an algebraic equation system (AES) of the form F(X) = 0 that can be solved e.g. by Newton-like schemes. To name some representatives, the Harmonic Balance (HB), Finite Difference (FD) and Shooting (SH) method [5,15,18] are capable of approximating periodic motions and can also be extended to quasi-periodic solutions [1,6,16]. These approaches based on an AES are suitable for both stable and unstable solutions, Hartmut Hetzler hetzler@uni-kassel.de 1 Engineering Dynamics, University of Kassel, Mönchebergstrasse 7, 34125 Kassel, Germany but can differ significantly in their numerical convergence speed depending on the problem. The analysis of steady-state motions is also of particular interest for practitioners and engineers. With this use and application case in mind, the present publication focusses on a class of mechanical structures, where non-linearities are limited to smaller (local) areas, e.g. contacts or joints causing frictional forces. Here, this leads to a few degrees of freedom (DoF) whose equations of motion involve strong non-linear terms-the so-called non-linear degrees of freedom. These potentially exhibit locally strong non-linear behaviour. In such cases, it is necessary to obtain a detailed picture of the steady-state dynamics and it is consequently essential to resolve the DoF with non-linear behaviour more finely. There are two ways to locally improve the accuracy of a corresponding method: On the one hand, this can be handled within the computational method itself, as it was realized in [9] for the HB method. Here, a harmonic selection technique is utilized to enhance accuracy by increasing the number of higher harmonics of the Fourier series for DoFs showing a strong non-linear dynamic behaviour. For an FD scheme, it is Bold indicates the proposed method (FD/HB-method) in the presented framework possible to increase the discrete resolution in time for certain degrees of freedom, as proposed in [13] for time-domain methods.
On the other hand, approximation methods typically have different strengths making them more suitable for certain problems. As an example, the HB is very efficient if only a few harmonics are needed to approximate a solution. This is typically the case for weakly non-linear problems [19, p. 131]. The FD may be better suited for strongly non-linear problems [1]. The SH method on the other hand becomes problematic when dealing with highly repulsive solutions. 1 For the case of systems with strong local non-linearities discussed here, the combination of two different computation methods into a hybrid approximation framework is a promising approach. Strengths can be combined and weaknesses mitigated thereby increasing computational efficiency. The basic idea is a partitioning of the DoF of the system into a non-linear and a linear node set [6,17,21], where different approximation methods will be chosen for each set. The nodes coupling the linear and non-linear set need to be addressed by both methods. The resulting AES is then solved by Newton-like methods. This approach of partitioning the system into a linear and a non-linear part, and combining individual methods to hybrid approximations schemes is common in literature [4,6,21,22,24]. Surveying different approaches reveals that vibrations of the linear subsystem are almost always approximated using HB (cf. Table 1). For the linear part, the HB-method yields the frequency response matrix for a defined frequency range and thus allows for an efficiently determination of the steady-state response to an excitation. In particular, this frequency response matrix can be explicitly stated in analytical form. In this context, the excitation of the linear sub-systems may either stem from external forces or from the adjacent non-linear subsystem, which is coupled in the sense of a Master-Slave reduction. In a way, this reduction 2 is similar to dynamic condensation, which is often employed in the context of model order reduction [14,24]. Although the linear and non-linear equations of motion are coupled, there are no general restrictions on the selection of the approximation method for the non-linear subsystem. Scheme 1 [21,22], which is based on combining a SH method for the non-linear node set with a HB method for the linear set. The main advantage of this hybrid method is the ability of selecting a specific numerical time integration (NTI) scheme. This is beneficial when dealing with non-smooth dynamical systems involving e.g. non-smooth friction or impacts. However, the Jacobian matrix is evaluated numerically, which increases the computational effort for systems with a high numbers of non-linear nodes [22]. Scheme 2 [4,6,24] is based on combining two HB schemes for the linear as well as for the non-linear subsystem. The resulting AES depends only on the Fourier coefficients of the non-linear nodes which are chosen as masternodes while the linear nodes are treated as dependent slave nodes. Thus, no costly transformations from time to frequency domain of the adjacent non-linear nodes are needed. Nevertheless, the resulting AES for the Fourier coefficients of the non-linear subsystem may show poor conditioning regarding numerical convergence for strong non-linear behaviour. One reason can be seen in the global ansatz functions of the HB scheme, resulting in extensive coupling of the algebraic equations and thus a fully populated Jacobian matrix. Here, changing of one single entry of the solution vector causes various shifts in the residual equations which may lead to poor convergence. Scheme 3 In this paper, a new hybrid approximation scheme for steady-state vibrations is presented which combines the HB and the FD method (cf. also [12]).
Here, the FD method is used to solve the non-linear node set, yielding a band structured Jacobian matrix that can be specified with adequate knowledge of the non-linear forces analytically. In contrast to the SH method-that follows the flow of the ODE within the NTI scheme, with periodicity explicitly enforced by the associated AESthe periodicity and thus the corresponding boundary value problem is fulfilled a priori by the FD approximation. Hence, FD is less sensitive to initial conditions of a Newton-like method when approximating highly repulsive solutions [19, p. 129]. Furthermore, both FD and HB can handle differential-algebraic equations. 3 without much effort [15, p. 27]. In conjunction with the analytical evaluation of the linear nodes' motion, these properties are promising prerequisites for an efficient and fast converging computational method.
This contribution is divided into three parts. The first part deals with the identification and separation of the linear and non-linear node sets. In addition, the local non-linearities considered within this contribution are defined. The second part is dedicated to the proposed Finite Difference/Harmonic Balance (FD/HB) method itself with a more detailed view on how the steady-state dynamics of a mechanical system 4 can be condensed to the non-linear node set. Therefore, a general treatment of the linear domain is suggested with a transformation in frequency domain by employing the Fourier transform. This section proceeds with the application to periodic solutions and finally the derivation of the algebraic equation system of the proposed method. The third part gives an application to a self-excited system showing periodic stick-slip vibrations. A reference solution is defined and followed by demonstrating the performance of the proposed method compared to both FD and HB approximation applied to the whole system. Finally, the results are discussed and a conclusion and outlook are given.

Mechanical systems with localised non-linearities
The basic idea of the proposed framework is to increase efficiency and accuracy by applying different numerical approximation schemes to the non-linear and linear parts of the equations of motion (EOM) of a mechanical system. In a first step, the DoF subjected to only linear terms are grouped into a set (L) denoted as the linear nodes. The remaining DoF are collected in the set (N) and referred to as non-linear nodes. 5 Now the vector u(t) containing all DoF of the system can be split into two coordinate vectors u N (t) ∈ (N) and u L (t) ∈ (L) in which all N N non-linear and N L linear DoF are collected separately. This subdivision is illustrated in Fig. 1 and can efficiently be introduced by means of a permutation of coordinates in u(t) 6 . Inserting this permutation into the EOM and partitioning the matrices results in Here, all linear forces are expressed by corresponding matrices while f N (u N ,u N ) collects all remaining non-linear forces. Please note that there are no non-linear couplings between the subsystems by definition: all coupling terms are linear and are completely described by the off-diagonal block matrices indexed with NL and LN. Before embedding Eq. (1) into the proposed framework, some further comments on the notation of the localised nonlinear forces shall be given. In general, there is no generic mathematical structure for the non-linear forces f N (u N ,u N ). However, since this contribution is focussed on localised nonlinearities, a closed analytical form for this class of nonlinear forces can be found, which was adopted from Krack and Gross [15, p. 141]. For example, if the relative motion of an i-th and j-th DoF causes a cubic restoring force f cubic (u N ) between each other, this force has to be considered only within the i-th and j-th equation of the momentum balance. Since no other DoF is (directly) affected by this cubic restor- 6 This means a reassembling of the degrees of freedom with u N u L = Bu =ũ. The Matrix B denotes a boolean permutation matrix containing exactly one unit entry in each row and each column. The remaining entries are zero and since ũ = u the relation B −1 = B holds. Note that applying this permutation only the sequence of the DoF within u(t) is changed that has no impact on the system dynamics.
ing force, the force vector may be written where γ is some parameter and a ∈ R N N ×1 only contains zeros except a i = 1 and a j = −1. Here, the actual values of u N are projected onto the vector a, so that only the relative deflection of the i-th and j-th entry is used to calculate the local cubic restoring force f (u) = γ u 3 . The evaluated force law is then multiplied by a to get the non-linear force vector that can be added within the global momentum balance Eq. (1). Hence, local non-linear effects influence only the DoF located in the same (small) area. As a popular example, friction in joints or contacts act locally since in general the friction force is not directly affected by state variables outside the contact zone. Since e.g. friction is (mostly modelled as) velocity dependent, the example shown in Eq. (2) has to be extended. A more general notation of N nl local non-linear forces summing up to f N is given by where a k , b k , c k ∈ R N N ×1 are vectors containing integers.
Here, a k selects the affected rows of the momentum balance in Eq. (1) and the b k and c k the affected DoF. The local non-linear behaviour is described by a scalar function This mathematical structure of f N (u N ,u N ) is not mandatory for the hybrid approximation method discussed below. However, it offers significant advantages since it allows for an analytical derivation of the Jacobian matrix and, thus, an efficient numerical solution of the non-linear AES.

Hybrid approximation method
The proposed hybrid approximation method intends to combine the advantages of HB and FD. To this end, the HB is applied to the linear subsystem in order to transform the problem to the frequency domain: this yields an algebraic relation between the amplitudes of the excitation and the stationary frequency response of linear part. Using this relation allows for expressing the linear DoFs (slave coordinates) as functions of the non-linear master DoFs and thus reduce the dimension of the problem. In a second step, the FD method is used for approximating the periodic motion of the entire system described by the non-linear master DoFs u N (t). Since both approaches apply to different domains, spectral-temporal transformations will be necessary.

System dynamics expressed by the non-linear node set
Introducing the abbreviations f (N) C and f (L) C into equation (1) yields of both the non-linear and linear set u N (t) and u L (t) and of the external forces f (L) EX (t). By inserting into Eq. (5b), a linear dependence from U N ( jω) to U L ( jω) is given by being the dynamical stiffness matrix that describes the influence from the j-th on the i-th domain. The continuous frequency axis of the frequency domain is denoted byω. Multiplication with the inverse of G LL ( jω) gives the Fourier transform of the linear subset u L (t) that is directly dependent on the FT of the non-linear subset u N (t) and the forcing f (L) EX (t), see Eq. (6). Re-transforming into time domain by the inverse Fourier transformation (iFT) gives the deflection u L (t) of the linear subsystem. Now that the interaction between the non-linear and linear subsystems has been identified, the feedback f (N) C of the linear to the non-linear subset is sought. Therefore, the deflection, velocity and acceleration of the linear subset u L (t) are given by with Inserting them into the linear feedback forces f (N) With Eq. (8), the coupling force f (N) C is given by showing a superposition of the feedback of the linear struc- holds. Please note, that only those non-linear DoF u (C) N (t) ∈ ∂ (N) which actually couple with the linear domain are transformed into frequency domain [cf. Eq. (5b), where M LN , P LN , C LN are typically sparse]. This is illustrated in Fig. 2.
Hence, the equations of motion cf. Eq. (1) can be condensed to the reduced order motion equations (ROME) This procedure corresponds to a Master-Slave reduction approach, where the linear node set corresponds to the slave coordinates that are only related to the non-linear master nodes and a given external excitation. Note that so far neither the type of the solution, e.g. periodic or quasi-periodic motion, nor the approximation method for u N (t) and thus u(t) are restricted. Consequently, these considerations apply in general, provided that the solution to be approximated is steady-state.

Application to periodic oscillations
The approach from the previous section may be applied to general stationary solutions with discrete frequency spectra (equilibria, quasi-/periodic solutions), whereas this publication focusses on periodic motions. Such a steady-state, periodic solution u(t) of Eq. (1) fulfils a boundary value problem (BVP) in time. The boundary condition is given by where T denotes the period of the oscillation with ω = 2π T being the base frequency. In order to apply the approach from the previous section to periodic solutions, the Fourier transforms must be specified for this BVP. Assuming u N (t), u L (t) and f (L) EX (t) are periodic functions in time, their Fourier transforms can be expressed by where δ(·) denotes the Dirac delta function and ωÛ N ( jω), ωÛ L ( jω) and ωF where the coefficients of δ( jω − jkω) must be zero to fulfil the equation. For k ∈ Z, this leads tô which means that steady-state dynamics of the linear DoF u L (t) are determined in frequency domain. In accordance to Eq. (12), the coupling forces can be calculated in time domain. Here, the convolution properties of the Dirac delta function are considered and the individual weights of the Fourier transforms are given bŷ with k ∈ Z. Finally, the steady-state solution for the linear subsystem related to t and u (C) N (t) is given by with the Fourier transform given in Eq. (17). Henceforth, the relations for the linear node set given in Eqs. (18) and (20) are understood as complex Fourier series with complex Fourier coefficientsÛ L,k = 1 Truncating the Fourier series gives the approximation with the harmonic truncation order H . This corresponds to the classical HB approach. As was shown above, the solution of the linear subsystem is derived analytically and hence, an approximation of the solution of the linear subsystem is given by Eq. (21a). Note that the approximation consists of (2 H + 1) global ansatz functions in time domain and that the maximal accessible truncation order H depends on the time respective frequency resolution of the non-linear DoF and the external forces. A detailed discussion of the selection of H is given in the following subsection.

Approximation of the non-linear subsystem
This section proceeds with a brief discussion of the FD method that is chosen for the approximation of u N (t) in time domain. In order to do so, the time derivatives of u N (t) occurring in the ROME, cf. Eq. (13), are approximated by the corresponding finite difference quotients on discrete time values on a prescribed grid t i ∈ T . General forms of these difference quotients are given by cf. [5,10] Within this contribution it is assumed that T = {t 0 +i t} i∈Z is a set of equidistant grid points and t is the time step. Here, α (1) k m , α (2) k m are the weights of the finite difference with index 8 k m ∈ M ⊂ Z. Further on, the approximations for the time derivations in Eq. (22) are substituted into the BVP solution framework given in Eq. (14). The periodicity of u N (t) is set with u N t N FD = u N (t 1 ). For the approximation method, only one period T of the steady-state solution is considered.
where 1 ≤ i ≤ N FD , θ = ω t is the step size and k m ∈ M. One advantage of this non-dimensionalization is that the angle grid T θ is independent of the period T and hence, is a fixed grid bounded by 0 and 2π . The approximations oḟ u N (t i ) andü N (t i ) are now substituted into Eq. (13). Due to the approximation of the derivatives and neglect of higher , evaluating the modified ROME at all values on the grid T θ gives a residual function

Residual equations of the Finite Difference/Harmonic Balance method
The  EX ( jkω) and N FD θ = 2π , the two integrals reduce tô with N FD being a prescribed number of steps. In accordance to the Nyquist Shannon sampling theorem, the number N FD has to be sufficiently high to ensure a correct depiction of the frequency domain [20, p. 124]. The sampling frequency f S = N FD T must at least be higher than two times the maximum occurring frequency in the solution, so that no aliasing effects will occur. Here, spectral leakage does not occur, since the evaluation domain is a multiple of the period duration of the lowest frequency.
Next   The global residual equation is build up by utilizing the Kronecker product 9 ⊗, by which the sums of both the DFT and iDFT can be written as a global matrix multiplication  For a global notation of the residual function, the Kronecker product provides again a notation for both the dynamical stiffness matrices G i j ( jkω) and the difference quotients. The global dynamical stiffness matrices read 9 The definition of the Kronecker product relevant in the present context can be found in Appendix A.
Here, the Kronecker product allocates the equations for the coefficients for the corresponding kth Fourier order in a global manner. Introducing the matrices α (1) , α (2) containing the difference weight factors, a global notation is analogously found for the difference quotients where I N ∈ R N N ×N N and X corresponds to the globally stored approximated velocities and X are the approximated accelerations, see Appendix C. Finally, the closed form algebraic equation system for the suggested FD/HB method reads where ω is the base frequency of the periodic solution and witĥ . . .
being the external forces acting on both the linear and nonlinear structure. The non-linear forces are given by Hence, with Eq. (33) an algebraic equation system was derived that is non-linear within the solution vector X, since F N X is non-linear in general.
Finding X 0 , being a root of the given global residual function R X , is a common problem in numerical calculation and is achieved, for example, by utilizing the Newton-Raphson algorithm. This scheme needs the Jacobian matrix J (n) = ∂ R ∂ X n that is ideally sought in an analytical form to avoid costly approximations of J (n) . For the suggested FD/HB method the Jacobian matrix can directly be derived from Eq. (33) with The crucial point is the evaluation of the derivatives of the non-linear forces . If no information of these derivatives is available, it is calculated numerically by finite difference approximations, which become costly for larger dimensions of X. To circumvent the numerical approximation, a detailed prescription of the non-linearities is reasonable. Provided that the local formulation for the nonlinear forces in Eq. (3) can be used, the derivation of the non-linear forces w.r.t. X (n) is Here, for the actual solution vector X (n) into the global equations and D θ = ω θ (α (1) ⊗ I N ) is the operator for the first order derivation of u N (t), cf. Eq. (32a).
Concluding this section, a hybrid approximation method was presented, combining the benefits of FD and HB. The scheme is valid for systems showing steady-state vibrations. Under this assumption, it was shown that the dynamics are purely representable by the non-linear node set in terms of a Master-Slave reduction, cf. Eq. (13). By applying this approach to periodic oscillations, it was shown that the deflection of the linear DoF is analytically derived with EX . When analysing self-excited vibrations, the angular frequency ω is in general unknown. Hence, an additional equation has to be added that is denoted as phase condition. Some common variants are listed in [23, p. 304].

Numerical investigation
This section focusses on the application of the proposed hybrid method. Here, an important aspect is the inclusion of the analytical solution of the linear DoF. Thus, the system dynamics can be described solely in terms of the non-linear DoF, which can significantly reduce the system size. Hence, this section aims to determine the scalability of accuracy in conjunction with relative computational effort for an increasing number N L of linear DoF. This is done by investigating a chain of oscillators with a single non-linearity allowing a straight forward increase of the linear DoF. On this Basis, FD/HB is compared with both the classical HB 10 and pure FD method. As reference solution the Shooting method is selected and a nominal error metric is defined.

Stick-slip vibrations of a chain of oscillators
As a numerical example, self-excited (periodic) vibrations of a chain of oscillators with local non-linearity (cf. Sect. 2) are studied, see Fig. 3. Although this is a strongly abstracted example, it exhibits essential properties that also occur in large FE systems.
Here, the corresponding system matrices are band structured, since the state variables are often not coupled over 10 In the present context, for numerically performing the HB method the AFT (alternating frequency-time) scheme of Cameron and Griffin is used, cf. [3,15]. In contrast to the implementation suggested in [15], the discrete Fourier transformations in the AFT scheme are performed by the matrices W (DFT) and W (iDFT) and not by means of the FFT (Fast Fourier Transformation) in order to achieve comparability of results. Here, N fft = 2 10 time samples are considered. larger spatial distances. Depending on the FE-shape functions, only the neighbouring nodes affect each other. Within the current minimal example, the non-linearity is located at the right end of the chain of oscillators being a frictional contact between the corresponding DoF and a belt moving with velocity v B . The friction force f μ is modelled by an exponentially decaying Stribeck friction characteristic that causes self-excited oscillations [11]. Since the discussed FD approach is applicable for smooth differential equations, 11 the sticking domain of the friction law is softened by a regularization where μ = μ ST − μ SL and k and are the regularization parameters, μ SL is the sliding friction coefficient for v rel → ∞ and μ ST N is the maximal sticking force corresponding to Coulombs law. In accordance to the local force representation in Eq. (3), the corresponding vectors are a = c = −1 and relative velocity is v rel = v B + c u N . The parameters of the friction curve are given in Table 2.
For a variable number of linear DoF, the dimensionless mass and stiffness parameter m and k of each oscillator are selected in accordance to the discretization of an elastic rod from which follows that to ensure comparable dynamics when adding linear DoF. The given dependence of mass and stiffness properties can be interpreted as a finer discretization of the linear domain of the underlying continuous system. In the following, viscous modal damping is assumed with a damping ratio δ = N N +N L 100 and global damping matrix P = δ k C. Since the resulting equations of motion [cf. Eq. (1)] are autonomous, the base frequency of the limit cycle is a variable and thus, an additional unknown in the AES. Consequently, the integral phase condition from Doedel et al. [7] is used to close the equation system.
In summary, due to the (regularised) decaying Stribeck friction characteristic, the rest position of the system becomes unstable for low belt velocities v B . 12 However, the amplitudes of the oscillations occurring are limited by the non-linearity of the system. This results in an asymptotically stable limit cycle in the sense of Ljapunow, which is approximated in the following for different numbers of linear DoF.

Quantitative comparison in accuracy and computation time
To highlight the abilities of the proposed hybrid framework, a quantitative comparison between the FD/HB method, pure FD and pure HB method is carried out under the following assumptions: • the comparison is exemplarily shown for N FD = 60 time samples • for FD/HB, the Fourier series approximation in frequency domain is done for four truncation orders H ∈ {3, 5, 10, H max }, where H max is the maximum number of harmonics in accordance to the Nyquist-Shannon theorem • for FD and FD/HB, the same difference quotients were used (see below) • for the HB method, the maximum number of harmonics is H max .
For all three methods, the error to reference ε M and computation time are compared. A sufficiently accurate approximation of the limit cycles of the system shown in Fig. 3, computed via the Shooting method, is used as a reference for the comparison. 13 For the FD based schemes, the time derivatives of first and second order are approximated using the difference quotientṡ Since the approximation of the first derivativeu i corresponds to an extended formulation for a second order upwind 12 However, the viscous damping must be sufficiently low so that the eigenvalues of the system linearised around the rest position have a positive real part. 13 For the SH method as reference, the trapezoid rule in its implementation as ode23t in Matlab©is chosen for the numerical time integration that has no numerical damping. For finding a root of the corresponding SH residual function, a trust-region dogleg algorithm within the fsolve implementation is utilized. For the technical realization of the SH method, the implementation proposed in [16, p. 301] was used. scheme, the second derivativeü i is discretized by a second order central difference scheme. 14 First, the occurring stick-slip limit cycle is analysed for FD/HB, HB and FD at v B = 0.3.
For the mechanical benchmark system consisting of one non-linear and one linear DoF, the results are shown in Figs. 4 14 The individual treatment of the first and second derivative prevents oscillations within the values at the grid points due to high local Péclet numbers. This may occur for problems with strong convective character, if the first derivative is approximated using central differences at a coarse time grid [25]. This comparison shows that all three methods meet the reference with no large deviations. In addition, no significant differences can be seen for different approximation orders H for the linear part of the FD/HB method.
This statement does not hold for an increasing number of DoF. This is exemplary shown for the the limit cycles of a chain of oscillators with 19 linear DoF, see Figs. 6 and 7.
Here, the motion of the non-linear DoF is shown in fig. 6 and the motion of the first linear DoF u L,1 being most distant from the non-linear subset in Fig. 7.
As indicated by these plots, the truncation order of the Fourier series for the FD/HB has a significant impact on the approximation results. Especially for the linear DoF, the velocity is much smaller than the reference predicts, if an insufficient number of harmonics is considered, see H ∈ {3, 5}. In addition, oscillations for the non-linear DoF occur within the velocity nearu N (t i ) ≈ −v B , cf. Fig. 6. Here, the frequency resolution of the linear node set and, thus, the coupling forces is too low. That causes the additional oscillations within the non-linear DoF.
To further examine the capabilities of the proposed method, two aspects are examined, comparing FD/HB to classical FD and HB: the first aspect, is a quantitative comparison via a nominal error to reference. Here, the error in the non-linear state z N (t) = (u N (t),u N (t)) is compared over one period T . The corresponding error norm is defined by  15 More precisely, u (SH) N is a numerically via NTI calculated solution to an initial value problem in the time interval from 0 to T with the SH solution as initial condition. The NTI is performed with the trapezoid rule in its implementation as ode23t in Matlab©. 16 Initial conditions are estimated using the tangent ∂ X ext 0 at the current solution point at v 0 B = 0.3. The corresponding initial guess for the Newton-Raphson iteration at v 1 B = 0.31 is given with where s denotes the arc-length of the implicit curve defined by the residual function R ext (v B , X ext ) = 0. For a more detailed discussion on prediction and continuation of implicitly defined curves, the interested reader is referred to Allgower and Georg [8] and Marx and Vogt [16]. As expected, the comparison shows that the FD/HB becomes more accurate for higher truncation orders H within the coupling forces, which on the other hand leads to an increase in computation time. However, the FD/HB is computationally more efficient than pure FD for all H for systems with more than 10 linear DoF.
Interestingly, for the maximum harmonic order H max , the computation time becomes lower than the time needed for pure FD (at about N L = 10), whilst the relative approximation errors ε M of FD/HB and FD are almost identical, see Fig. 8. One possible explanation is that the information content of the motion of the linear DoF is the same for both methods. The only difference is that this information And since the approximation for the linear DoF is solved analytically within the hybrid approach, the AES is significantly smaller for a larger number of linear DoF than for pure FD and therefore more efficient. This effect also scales with the number of time grid points N FD . For higher values of N L , the FD/HB is more efficient than the pure FD is at lower numbers of N FD and vice versa. As expected, the HB provides sufficiently accurate results. For all investigated numbers of linear DoF the nominal error is on a par with that of FD and FD/HB with H max , see Fig. 8. However, the HB method requires more computing time compared to the other methods. This may be caused by the AFT algorithm itself, where two discrete Fourier transformations 17 are carried out.
A closer look at the error diagram in Fig. 8 shows that the approximation error for truncation orders H < H max seems to saturate, when increasing the number of linear DoF. In order to be able to explain this observation, the Fourier coefficients of the non-linear deflection u N (t) of the SH solution are shown in Fig. 10. The absolute value of the kth Fourier coefficientÛ N,k of the non-linear DoF is plotted against the corresponding harmonic number k and the number of linear DoF N L . As expected, it can be seen that the absolute values of the Fourier coefficients decrease with the harmonic order, but increase for higher numbers of N L within the considered range. Therefore, the nominal error ε M will increase for low H since essential frequency content of the non-linear DoF is not considered in the calculation of the lin- 17 The implementation of the DFT and iDFT via the Fast Fourier Transformation promises a more efficient calculation [15], but was not used for reasons of comparability between the methods. ear DoF and the coupling forces f (N) C (u L ,u L ,ü L ). Following fig. 10, the next higher Fourier orders of H ∈ {3, 5, 10} form a plateau for increasing N L that may cause the saturation of the nominal approximation error in Fig. 8.
As a summary for this section, one can state that the proposed FD/HB scheme can be more efficient than FD and HB, if a certain number of linear DoF is surpassed. Compared to both methods, the error of the proposed method is in the same order of magnitude if a sufficient number of harmonics is used. The decisive factor is that this method requires less computing time as the number of linear degrees of freedom increases. For the HB, a high sensitivity regarding initial conditions of the Newton iteration was observed here -especially for small deviations regarding the autonomous base frequency. Furthermore, the results show that the harmonic order H has a high impact on the accuracy of the approximation and has to be selected in accordance to the Nyquist-Shannon theorem to gain maximum accuracy. When taking the maximum number H max of harmonics, there is a characteristic point depending on N FD and N L , where it is more efficient to use the proposed framework than pure FD.

Conclusion
This contribution presents a hybrid approximation method, which combines the Finite Difference (FD) and Harmonic Balance (HB) method, for investigating steady-state solutions of non-linear dynamical systems. The method is derived for, but not limited to mechanical systems. The basic idea starts by separating the equations of motion into two sets containing non-linear and solely linear equations. The nonlinearities are assumed to act only locally and the two sets are bi-directionally connected, by the so-called coupling forces. The goal is now to describe the steady-state dynamics of the entire dynamical system just in terms of the degrees of freedom (DoF) that occur in the non-linearities -the so-called non-linear DoF. This is possible since the linear equations of motion can be solved analytically by applying a Fourier transform. Consequently, a closed form expression for the linear DoF solely depending on the -yet unknownnon-linear DoF and the corresponding excitation forces can be provided. For the special case of periodic motions, the Finally, the corresponding algebraic equation system for the FD/HB method is given and the Jacobian matrix is analytically derived. Compared to the HB method, the Jacobian matrix is band structured for FD/HB method, which improves convergence. In addition, this method can deal with highly repulsive steady-state solutions in contrast to hybrid shooting (SH) methods.
In order to study the performance of the proposed method, periodic stick-slip limit cycles of a chain of oscillators with varying numbers of DoF are approximated. The results are compared to solutions calculated by the FD method and the HB method each applied to both the non-linear and linear DoF. For all methods, the computational time is compared and the limit cycles are validated against a SH reference solution. The results of the comparison show that for systems with a small number of non-linear DoF, the hybrid approximation method can improve computational efficiency for larger systems. The level of accuracy is directly depending on the number of considered harmonics H of the closed form expression for the linear DoF and the number of FD time grid points N FD . The maximum accuracy possible is limited to the truncation order H max according to the Nyquist-Shannon theorem, which depends on the number N FD . As a major advantage of the hybrid framework, the accuracy of that approximation with H max is on an equal level with the approximation gained by pure FD and pure HB. However, it takes significantly less time for systems with a higher number of linear DoF due to the efficient approximation of the linear DoF via HB within the proposed framework.
For future research, a harmonic selection technique will be integrated into the hybrid framework. Here, an individual truncation order for the Fourier series within every linear DoF can be selected. On the basis of the general approach via Fourier transformation, the method will be extended to quasi-periodic oscillations.
Funding Open Access funding enabled and organized by Projekt DEAL. No funding was received to assist with the preparation of this manuscript. The authors have no competing interests to declare that are relevant to the content of this article.
Availability of data, code and materials Not applicable for this contribution.

Declarations
Consent to participate Not applicable for this contribution.

Consent for publication Not applicable for this contribution.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.