Lyapunov Equivalent Representation Form of Forced, Damped, Nonlinear, Two Degree-of-Freedom Systems

Featured Application: The qualitative and quantitative dynamic behavior of two degree-of-freedom nonlinear systems can be studied by using their corresponding decoupled one degree of freedom Dufﬁng type equivalent representation forms in the sense of Lyapunov, with the advantage of capturing amplitude-dependent nonlinear mode shapes. Abstract: The aim of this paper focuses on ﬁnding equivalent representation forms of forced, damped, two degree-of-freedom, nonlinear systems in the sense of Lyapunov by using a nonlinear transformation approach that provides decoupled, forced, damped, nonlinear equations of the Dufﬁng type, under the assumption that the driving frequency and the external forces are equal in both systems. The values of Lyapunov characteristic exponents (LCEs), Lyapunov largest exponents (LLE), and time-amplitude and frequency-amplitude curves computed from numerical integration solutions, indicate that the decoupled Dufﬁng-type equations are equivalent, in the sense of Lyapunov, to the original dynamic system, since both set of motion equations tend to have the same qualitative and quantitative behaviors.


Introduction
Minorsky [1], Caughey [2], Iwan [3,4], Sinha and Srinivasan [5], Agrwal and Denman [6], Yuste and Sánchez [7,8], Cai [9], Farzaneh and Tootoonchi [10], and many others have reported approaches that replace linear and nonlinear dynamic systems by equivalent ones with known solutions that are closed to the original system, which produce the same oscillations that appear in the original equations of motion. . There, the authors focused on using linearization, weighted mean-square, and least squares techniques to determine the equivalent expressions that provide solutions with the same quantitative and qualitative original system dynamics response behavior.
On the other hand, a nonlinear transformation approach that is based on a cubication method, was recently introduced to obtain the equivalent representation form of conservative two degree-offreedom nonlinear oscillators [11]. There, the authors developed an approach to replace a two degree-of-freedom homogeneous, undamped system by another equivalent system with known solutions that were closed to the original one. In that approach, they first replaced the system restoring forces by polynomial expressions and then used a transformation technique to replace the resulting x 2 + k 11 k 12 k 21 k 22 into equivalent decoupled, nonlinear equations of the Duffing type [16,17] in the Lyapunov sense, with similar quantitative and qualitative dynamic behaviors. In Equation (1), m 1 and m 2 are the masses of the system, c ij and k ij are the damping and stiffness system coefficients, respectively, f 1 (x 1 , x 2 ) and f 2 (x 1 , x 2 ) are even, nonlinear restoring forces, and Q * 1 (t) and Q * 2 (t) are driving periodic forces. The system initial conditions are assumed to be given by x 1 (0) = x 10 , x 2 (0) = x 20 , .
x 20 . To achieve such equivalence in the sense of Lyapunov, first, the system (1) is written in its normal canonical form in an attempt to predict the system's dynamic behavior. Then, from the normal form representation of the system (1) it is assumed that the modal system's generalized coordinates can be equivalently expressed as a power series expansion [18][19][20][21] to decouple each normal mode equation into a forced, damped, nonlinear differential equation of the Duffing type, whose approximate solution has been widely discussed in the literature [22][23][24].
It is conjectured that the equivalent form representation of the decoupled equations will have the same Lyapunov characteristic exponents value (LCEs) as those for the normal canonical form of the original equations of motion. In other words, it can be conjectured that the resulting decoupled equations are equivalent in the sense of Lyapunov if the solutions of the decoupled expressions are solutions of the normal canonical form of the original equations of motion, and vice versa [25,26]. If the decoupled Duffing equations of a two-degree-of-freedom dynamical system are equivalent in the sense of Lyapunov to the canonical form of the original equations of motion, then the Lyapunov characteristic exponents are the same, i.e., λ i = λ ie , with i = 1, 2, 3, 4, where λ i , and λ ie are the canonical form, and the equivalent Lyapunov characteristic exponents, respectively. It is believed that the decoupled form, in the sense of Lyapunov, of the normal canonical representation of Equation (1) can be used as a valuable tool for understanding the significance of nonlinear normal modes of several nonlinear systems. In fact, the generalization of conservative nonlinear normal modes and its computation has led to the use of these in a diverse range of applications in structural dynamics [27]. Furthermore, nonlinear normal modes have been used to identify localized modes in micro-electromechanical devices [28] for the development of a new kind of low frequency acoustic absorber [29] and for the analysis of nonlinear vibrations for double-walled carbon nanotubes [30]. Also, the nonlinear normal modes have helped in the understanding of the dynamic response behavior of full-scale aircrafts and satellites [31,32] as a tool for sub-structuring [33], for the construction of reduced-order models [34], and for damage detection in engineering structures [35], among others uses.
Therefore, the aim of this paper is to verify the previously mentioned conjecture by introducing a nonlinear approach to replace the canonical normal form of a dynamical two-degree-of-freedom system by two equivalent decoupled expressions of the Duffing type. A few examples are presented to outline the basic ideas behind the determination of the equivalent representation form in the sense of Lyapunov of nonlinear, forced, damped, two degree-of-freedom systems. The steps involved in the proposed transformation approach to decouple nonlinear equations are introduced next.

Transformation Technique
Before the decoupling procedure to find the equivalent representation form of Equation (1) in the sense of Lyapunov is introduced, Equation (1) is first written into its canonical, normal mode by applying the following linear transformation [36] x 1 which yields: where u 1 and u 2 are the normal coordinates of the linear, undamped, free vibration system, κ is the modal matrix which consists of the characteristic vectors that represent the natural modes of the linear system obtained from (1), and these are given by κ = [R 1 {κ} 1 , R 2 {κ} 2 ], where R i are scale factors that can be determined from √ 1/M ii , in which M ii = κ T i Mκ i is the generalized mass and {κ} 1 and {κ} 2 denote column eigenvectors. If the scale factors are chosen so that M ii = κ T i Mκ i = 1, i = 1, 2, then the normal modes are mass-orthonormal, and thus, ω 2 ni = κ T i Mκ i , i = 1, 2. Furthermore, ϕ 1 through ϕ 8 , and P 1 (t) and P 2 (t) are parameters that are defined in accordance with the physics of the system. Here, the initial conditions are assumed to be given as u 1 (0) = u 10  u 20 . It is further assumed that the linear transformation (3) preserves the Lyapunov exponents of the original system (1); this implies that Lyapunov exponents of the original system (1) are equal to those computed from (4) and hold for any solution of the original system and the corresponding solution of the transformed one, as stated by Barabanov [37].
On the other hand, one can notice that the nonlinear Equation (4) does not have known exact solutions and is as difficult to solve as the original equations of motion (1). In order to avoid these difficulties, it is now assumed that the modal system (4) can be expressed, for the system generalized coordinates (displacements or rotations), in a third-order power series expansion [18][19][20][21]. Thus, Equation (4) becomes .. ..
Expressions (5) and (6) provide approximate decoupled representation forms of the original equations of motion whose accuracy depends on the unknown coefficients: a i , b i , and µ i . Here, the external forces and the driving frequencies are considered to be the same in both systems [16,38]. Once these coefficients are determined, the approximate solutions of (5) and (6) can be derived by using perturbation or numerical techniques [39], and their LCEs can be numerically computed by using the procedure discussed in [40,41].
The following remarks are important to set the number of terms of the power series expansion of Equations (5) and (6) that provides an invariant polynomial expression which describes the shapes of the invariant manifolds that are related to the system's nonlinear normal modes [21,38,42]: (a) Odd terms arise in the right-hand term (RHT) of the invariant polynomial expressions of Equations (5) and (6) if the restoring forces of the dynamic system are described by an invariant odd polynomial expression.
(b) If the restoring forces of the original system have mixed-parity nonlinearities, then even and odd terms arise in the RHT invariant polynomial expressions of Equations (5) and (6). (c) Velocity-dependent terms arise mainly in the RHT of the invariant polynomial expressions of Equations (5) and (6) in order to capture not only the effective trend of the system's nonlinearities responsible for displaying amplitude-dependent nonlinear mode shapes, but also to take into account decay rate effects.

Determination of a i , b i , and µ i
To determine the coefficients a i , b i , and µ i , a minimization procedure analogous to the one followed by Caughley in [2] is assumed. Here, it is assumed that the square of the difference between the acceleration of the originally coupled system and that of the simplified cubic one, given by Expressions (5) and (6), could be minimized for each system by considering that the weighted mean square errors, U 1 and U 2 , can be determined by using the following expressions [17]: Furthermore, the minimization of U 1 and U 2 with respect to the coefficients, a i , b i and µ i , can be achieved from [9] ∂U 1 ∂a i = 0 with i = 1, 2, 3, Equations (9) and (10) yield, after performing some algebraic computer calculations with the help of Mathematica 11.1.1.0 symbolic computer package, the following expressions: Here, η i , η ii , V i , and V ii describe integration constants whose values must minimize Expressions (7) and (8). The general procedure to compute the integration constants η i , η ii , V i , and V ii is discussed in the following section.

Computation of η i , η ii , V i , and V ii Values
To determine the above integration constants, Equations (11) and (12) are substituted into Equations (7), and (13)- (16) are substituted into Equation (8). This step yields, after performing the corresponding integrations, the following expressions: Notice that the resulting polynomial expressions, (17) and (18), provide convex, continuous functions in the relative interior of their domains which depend on the values of η i , η ii , V i , and V ii [43,44]. These functions, S i , have a unique minimum value within their bounded domain intervals, as discussed in [45,46] and references cited therein. Thus, the values of η i , η ii , V i , and V ii that provide the minimum values of U i could be computed by following a classical minimization algorithm in which the partial derivatives of Equations (17) and (18) must be determined, respectively, with respect to the unknowns: η i , η ii , V i , and V ii . This step provides the following equations:    which are set equal to zero to find the critical points at which Equations (7) and (8) are minimized. Notice from Equation (9) through (16) that the denominator terms depend, respectively, on η 1 , V 1 , η 22 , and V 22 . Therefore, these coefficients cannot be zero. Thus, ∂S 1 /∂η 1 = 0, ∂S 1 /∂η 2 = 0, ∂S 1 /∂V 1 = 0, ∂S 1 /∂V 2 = 0, ∂S 2 /∂η 11 = 0, ∂S 2 /∂η 22 = 0, ∂S 2 /∂V 11 = 0, ∂S 2 /∂V 22 = 0 are identically satisfied, if, and only if, η 2 = 0, V 2 = 0, η 11 = 0, and V 11 = 0. Then, the expressions for a i and b i that minimize the weighted mean square error of using Expressions (7) and (8) to determine the coefficients, a i , b i , and µ i that describes the equivalent representation form of the original equation of motion, into the approximate forms (5) and (6), simplify to Thus, the above approach and the usage of power series expansion up to the third-order to replace the modal restoring forces of Equations (5) and (6) provide the following uncoupled Duffing-type equivalent equations: .. ..
It is important to point out that such equivalence in this two degree-of-freedom (DOF), non-integrable system (4) is only equivalent, in an approximate form, to two single degree-of-freedom integrable equations. Furthermore, in the conservative case, almost periodic solutions can be observed in two-DOF nonlinear system of coupled equations, but they are absent in the uncoupled, single-DOF, nonlinear Equation (29). However, the linear combination of their approximate solutions, in accordance with Equation (3), could exhibit, as in the case of the original system of differential Equation (1), almost periodic solutions, as will be demonstrated in Example 2. In accordance with Shaw and Pierre [20], the advantage of having decoupled modal equations is related to the possibility of obtaining information about the system modal amplitudes and thus, the numerical integration solution of these equivalent expressions, which could provide valuable information about the qualitative and quantitative behaviors of the system dynamics.
On the other hand, if the restoring forces, f 1 (x 1 , x 2 ) and f 2 (x 1 , x 2 ), of the dynamic system (1) contain even, nonlinear effects, then, their equivalent representation forms in the sense of Lyapunov could also be expanded to power series in which decay terms must be considered. Section 5 focuses on studying the effects of these new terms in the decoupled form of a nonlinear, dynamic system.

Numerical Validation
The applicability of the proposed approach of replacing, in the sense of Lyapunov, a two-degreeof-freedom system with equivalent expressions of the Duffing type will be examined by considering three nonlinear dynamic systems.

Example 1: Dynamic System with Cubic Nonlinearities
To study the accuracy attained by applying the proposed approach to decouple nonlinear differential equations and to validate the proposed conjecture, first, the frequency-amplitude response curves of the original unforced, undamped system (4) will be compared to the backbone curves obtained from the Duffing equations described by Equation (29). To accomplish this goal, the nonlinear dynamic systems discussed in [20,21] are considered, in which the corresponding equations of motion are described by  .. where f 1 is the magnitude of the driving force, and ω f represents the corresponding driving frequency.
The system initial conditions are assumed to be given as x 1 (0) = x 10 , .
(35) Figure 1 shows the simulation performed to verify the accuracy of the proposed methodology that decouples the nonlinear normal mode's differential equations. The system parameter values in dimensionless units are m = 1, k = 1, ε = 0.5 with c = 0 and f 1 = 0; these are similar to those used in [20]. To plot the frequency-amplitude curve of the first mode, initial conditions of u 1 , u 2 , . u 1 , . u 2 = (1, 0, 0, 0) were considered, and for the second mode, u 1 , u 2 , . u 1 , . u 2 = (0, 1, 0, 0) was used. It is clear from the results exhibited in Figure 1 that the frequency-amplitude curves of the original system compare well with those of the uncoupled equations. Therefore, it is concluded that the replacement of the original equation of motion for two equivalent equations of the Duffing-type describes the qualitative and quantitative dynamic system behaviors well, as shown in Figure 2, in which the frequency-amplitude backbone curves for several values of k and ε have been plotted.
Next, the damped nonlinear case was considered with parameter values of m = 1, c = 0.3, k = 1 and ε = 0.5 to compare the numerical integration of Equation (30) with respect to those given by expressions (29). The results of the simulations in physical coordinates are depicted in Figure 3 in which the amplitude-time curves have been plotted by using the system's initial conditions of x 2 = (0, 0, 2, 0). It is evident from Figure 3, that the decoupled Equation (29) is in agreement with the numerical integration solutions of the original equations of motion (30). In fact, the LCEs of the original system have mean values of λ 1 = −0.0756, λ 2 = −0.0802, λ 3 = −0.3698, and λ 4 = −0.3743 bits/second, while the LCEs computed from the equivalent models are λ 1e = −0.0749, λ 2e = −0.0750, λ 3e = −0.3726, and λ 4e = −0.3773 bits/second. In spite of having some discrepancies in the computed numerical values, the order of magnitude of the LCEs of the equivalent expressions tend to be similar to those computed from the original equations of motion with root mean square error (RMSE) values of 0.1117 and 0.0788 for the first and second mode, respectively.

of 22
To further assess the accuracy of the proposed nonlinear decoupling procedure, the frequencyamplitude steady-state response curves of the forced, damped nonlinear system (30) were plotted with parameter values of m = 1, k = 1, c = 0.3, ε = 0.5 and f 1 = 0.25.  Equations (33) to (35). The parameter values used to obtain these plots for different values of and were 1 = 2 = 1 , with = 0 and 1 = 0 . Here, the black solid lines describe the solutions obtained from Equations (34) and (35), while the colored, dashed lines represent the backbone curves computed from Equation (33).  Equations (33) to (35). The parameter values used to obtain these plots for different values of and were 1 = 2 = 1 , with = 0 and 1 = 0 . Here, the black solid lines describe the solutions obtained from Equations (34) and (35), while the colored, dashed lines represent the backbone curves computed from Equation (33).  (29) and (30). The parameter values used to obtain these plots were 1 = 2 = 1 , = 1 , = 0.5 with = 0.3 , with initial conditions given by ( 1 , 2 ,̇1,̇2) = (0, 0, 2, 0) . Here, the black lines describe the numerical integration solutions of Equation (30), while the blue, dashed lines represent the numerical solutions obtained from the nonlinear equations of motion (29).
x 2 = (0, 0, 2, 0). Here, the black lines describe the numerical integration solutions of Equation (30), while the blue, dashed lines represent the numerical solutions obtained from the nonlinear equations of motion (29).
It can be seen from Figure 4 that the numerical integration solutions of the decoupled Equation (29) provide curves that slightly differ close to the resonant region to those obtained from the numerical integration solutions of the original equations of motion, but its qualitative and quantitative system dynamic predictions are, in general, good. The corresponding Largest Lyapunov Characteristic Exponent (LLE) curves plotted versus the driving frequency, ω f , are exhibited in Figure 5. It is evident from Figure 5 that the LLE curve computed from the equivalent uncoupled expressions (29) is closed to that estimated from Equation (30). In this case, the RMSE does not exceed the value of 0.0056. This confirms that the proposed conjecture is numerically true.

Example 2: A Forced System with Cubic Nonlinearities
Here, the nonlinear dynamic system introduced in [21] was examined, which has the form with initial conditions given by 1 (0) = 10 , ̇1(0) = 0, 2 (0) = 20 , and ̇2(0) = 0. To find the canonical representation form of Equation (36), the following transformation approach is applied [36]: which yields the canonical form (4). In this case, Oscillation amplitude, Oscillation amplitude, As a second example to investigate the validity of the proposed conjecture, a forced, damped nonlinear dynamic system with cubic nonlinearities was next examined.

Example 2: A Forced System with Cubic Nonlinearities
Here, the nonlinear dynamic system introduced in [21] was examined, which has the form with initial conditions given by x 1 (0) = x 10 , .
x 2 (0) = 0. To find the canonical representation form of Equation (36), the following transformation approach is applied [36]: which yields the canonical form (4). In this case, where and ..
where a i , b i , and µ i , are defined by Equations (27) and (28).
To assess the precision of the proposed approach, the frequency amplitude response curves were next computed using the following system parameter values: m 1 = 1 kg, m 2 = 1.5 kg, k 1 = 2 N/m, k 2 = 3.5 N/m, k 3 = 5 N/m, c 1 = 0.066 N-s/m, c 2 = 0.099 N-s/m, ε 1 = 1 N/m 3 , ε 2 = 1 N/m 3 and Q 1 = 0.2 N. The computation is performed starting with ω f = 0, with the initial conditions given by x 2 = (0, 0, 0, 0) and the driving frequency, ω f , increased gradually at small incremental driving frequency step values of ∆ω f = 0.05. The steady-state vibration amplitude of the previous solution is used as the initial condition to obtain the corresponding numerical integration solution at the next ω f value. Figure 6 shows the frequency-amplitude response diagrams from the numerically computed form, Equation (36), and from the corresponding equivalent representation forms (42). Notice that the numerical integration of (42) slightly differs from the numerical integration solutions of Equation (36) near to the second resonance region; however, the nonlinear modal equations capture the dynamics of the original system well. The LLE average values for Expressions (36) and (42) are λ LLE = −0.0229 bits/second, and λ LLEe = −0.0226 bits/second, respectively. These computed LLE average values are almost the same, which confirms the validity of the proposed conjecture.
Next the undamped and unforced case was considered with the above system parameter values and x 1 , .
x 2 = (0, −0.1, 0, 0). Then, the transformation relationships (42) were used to plot the time-amplitude response curves. Notice from Figure 7 that the linear combination solution (37) describes the qualitative and quantitative system dynamics well. In fact, almost periodic solutions can be observed when the modal displacements are transformed back into the original coordinate systems by using the relationships (37). Therefore, it is concluded that the equivalent representation form of Equation (36) validates the conjecture, since the numerical integration solutions follow the dynamic behavior, observed during the numerical integration of the original equations of motion, well.
As a final example, the nonlinear absorber system introduced by Ji in [50] was examined.

Example 3: A Nonlinear Absorber System
The following equations of motion,

Example 3: A Nonlinear Absorber System
The following equations of motion,

Example 3: A Nonlinear Absorber System
The following equations of motion, ..
describe the dynamics of a nonlinear absorber attached to a single degree-of-freedom nonlinear oscillator [50]. Here, x 1 , m 1 , k 1 , k 2 and c 1 represent the displacement, mass, linear and nonlinear stiffnesses, and damping coefficient of the primary system, respectively, while x 2 , m 2 , k 3 , k 4 and c 2 are system parameters related to the nonlinear secondary absorber system, and f 0 and ω f are the external excitation force and driving frequency, respectively. It is easy to show that system (43) has the normal canonical modal representation form (4) with To obtain the frequency-amplitude response curves of the nonlinear absorber system, the parameter values of m 1 = 10 kg, m 2 = 0.6 kg, c 1 = 0.1 Ns/m, c 2 = 0.08 Ns/m, k 1 = 44 N/m, k 2 = 8 N/m 3 , k 3 = 2 N/m, k 4 = −0.15 N/m 3 , and f 0 = 0.37 N were considered, which are similar to those used in [50]. Figure 8 illustrates the frequency-amplitude response curves obtained from the numerical integrations of Equations (29) and (43), in which the equivalent damping coefficients, µ 1 and µ 2 , are bigger than zero. Figure 8 shows that the decoupled solutions describe the dynamic behavior of the original equations of motion well. The LLE numerically-calculated curves are illustrated in Figure 9. As before, these computed LLE curves agree well along the driving frequency range, since the estimated RMSE value is lower than 0.0017.
Next, the following system parameter values are used to plot the frequency-amplitude curves shown in Figure 10: m 1 = 10 kg, m 2 = 0.8 kg, c 1 = 0.1 Ns/m, c 2 = 0.08 Ns/m, k 1 = 44 N/m, k 3 = 2 N/m, k 4 = −0.65 N/m 3 , f 0 = 0.37 N, and k 2 = 0 and 8 N/m 3 . As can be seen from Figure 10, the frequency-amplitude curves obtained using the corresponding decoupled expressions follow closely, in spite of having an absorber device that can have linear or nonlinear stiffness effects; these curves are obtained by numerically integrating the original expressions (43). In this case, the corresponding LLE-computed curves for the original and equivalent equations in the sense of Lyapunov are shown in Figure 11.
To further assess the accuracy of the proposed approach, the values of k 3 = 9.016 N/m, and 20.828 N/m with m 2 = 0.6 kg were used to excite internal resonances of the types, ω n2 = 2ω n1 and ω n2 = 3ω n1 , which are second-order and third-order internal resonance relationships, respectively. Figure 12 illustrates the frequency-amplitude curves obtained from the numerical integration solutions of Equations (29) and (43).
It is evident that the equivalent equations in the sense of Lyapunov capture both types of internal resonances without any additional assumptions or simplifications in the method. The estimated LLE curves for the original and equivalent equations are shown in Figure 13. From Figure 13, it is very interesting to see the agreement in the computed LLE curves which is also confirmed by the computed RMSE values and by the bifurcation diagrams shown in Figures 14 and 15. As a consequence, one can conclude that the proposed conjecture is true, i.e., the decoupled Duffing Equations (29) are equivalent, in the sense of Lyapunov, to the normal canonical form of the original equations of motion (43) if the system's nonlinearities are small and the oscillation amplitudes are moderate. In other words, although our proposed approach is straightforward and easy to apply, this uncoupling process could be only justifiable for small deviations in the system's linear behavior because, for large amplitudes, coupling and nonlinear terms are important. However, in spite of these drawbacks, the proposed decoupling technique can capture internal resonances, something that is cumbersome when applying other techniques [38,51]. Of course, further investigation into the applicability of this approach to the prediction of internal resonances in other dynamic systems and its potential advantages or disadvantages are beyond the scope of this paper, and these issues will be addressed in a forthcoming article. curves for the original and equivalent equations are shown in Figure 13. From Figure 13, it is very interesting to see the agreement in the computed LLE curves which is also confirmed by the computed RMSE values and by the bifurcation diagrams shown in Figures 14 and 15. As a consequence, one can conclude that the proposed conjecture is true, i.e., the decoupled Duffing Equations (29) are equivalent, in the sense of Lyapunov, to the normal canonical form of the original equations of motion (43) if the system's nonlinearities are small and the oscillation amplitudes are moderate. In other words, although our proposed approach is straightforward and easy to apply, this uncoupling process could be only justifiable for small deviations in the system's linear behavior because, for large amplitudes, coupling and nonlinear terms are important. However, in spite of these drawbacks, the proposed decoupling technique can capture internal resonances, something that is cumbersome when applying other techniques [38,51]. Of course, further investigation into the applicability of this approach to the prediction of internal resonances in other dynamic systems and its potential advantages or disadvantages are beyond the scope of this paper, and these issues will be addressed in a forthcoming article.    Figure 13. From Figure 13, it is very interesting to see the agreement in the computed LLE curves which is also confirmed by the computed RMSE values and by the bifurcation diagrams shown in Figures 14 and 15. As a consequence, one can conclude that the proposed conjecture is true, i.e., the decoupled Duffing Equations (29) are equivalent, in the sense of Lyapunov, to the normal canonical form of the original equations of motion (43) if the system's nonlinearities are small and the oscillation amplitudes are moderate. In other words, although our proposed approach is straightforward and easy to apply, this uncoupling process could be only justifiable for small deviations in the system's linear behavior because, for large amplitudes, coupling and nonlinear terms are important. However, in spite of these drawbacks, the proposed decoupling technique can capture internal resonances, something that is cumbersome when applying other techniques [38,51]. Of course, further investigation into the applicability of this approach to the prediction of internal resonances in other dynamic systems and its potential advantages or disadvantages are beyond the scope of this paper, and these issues will be addressed in a forthcoming article.       Amplitude, x 1 [m] Amplitude, x 1 [m] Amplitude, x 2 [m] Amplitude,     n2 = 2  n1 , k 3 = 9.016 N/m  n2 = 2  n1 , k 3 = 9.016 N/m

Fifth-Order Power Series Expansion
It is important to bear in mind that the predictions obtained in the above examples consider a truncated power series to replace the nonlinear restoring forces, = ( 1 , 2 ), by uncoupled cubic polynomial expressions, which, if inaccurate, would cause the conjecture to be false. In an attempt to

Fifth-Order Power Series Expansion
It is important to bear in mind that the predictions obtained in the above examples consider a truncated power series to replace the nonlinear restoring forces, = ( 1 , 2 ), by uncoupled cubic polynomial expressions, which, if inaccurate, would cause the conjecture to be false. In an attempt to

Fifth-Order Power Series Expansion
It is important to bear in mind that the predictions obtained in the above examples consider a truncated power series to replace the nonlinear restoring forces, f i = (x 1 , x 2 ), by uncoupled cubic polynomial expressions, which, if inaccurate, would cause the conjecture to be false. In an attempt to avoid this situation, the modal system (4) was next replaced by a truncated fifth-order power series expansion that contained nonlinear damping terms [13,38,[51][52][53]. As usual, driving forces were assumed to remain constant during the transformation approach. The uncoupled dynamic equations of motion are now written as .. ..
The Expressions (53) and (54) provide equivalent decoupled representation forms, in the sense of Lyapunov, of the original equations of motion whose accuracy depends on the unknown coefficients: a i , b i , and µ i . By following the procedure described in Section 3, the coefficients, a i , b i , and µ i , were found to be The values of η i , η ii , V i , and V ii were obtained by minimizing the following expressions: For comparison purposes, only the dynamic system examined in Example 1 was considered, since the oscillators studied in Examples 2 and 3 had small errors when the equivalent expressions in the sense of Lyapunov were compared to the original equations of motion. First, we focused on studying the unforced dynamic system (30) and plotted the amplitude-time curves by using the equivalent cubic, and quintic expressions: (29), (53) Figure 16 shows the time-amplitude curves for the two modes of the system. Notice that as a consequence of the fifth-order and nonlinear decay terms of the quintic approach, the equivalent representation forms (53) and (54) provide an improvement in the RMSE values. Next, the driving force magnitude of f 1 = 0.25 was considered and then, the frequency-amplitude response, as well as the LLE curves were plotted. As can be seen in Figure 17, the quintic equivalent representation form in the sense of Lyapunov provided a better approximation to the exact numerical values. An improvement was achieved in the frequency-amplitude response curves obtained from the quintic approach since the estimated curves match the numerical results well; however, one can notice from Figure 18, that when the quintic approach is used to compute the LLE values, the RMSE value is slightly higher than that computed by using the cubic model, which is mainly due to the numerical procedure used to compute these values rather than to the decoupling proposed approach.  0.0109 , 1 = 2.9995 , 2 = 0.1228 , 3 = 0.0008 , 1 = 0.0807 , 2 = −0.0206 , 3 = 0.0066 , 4 = 0.375, 5 = 0, 6 = 0, with fitting parameter values of 1 = 1, 2 = 0.3, 1 = 1, 2 = −0.15, 11 = −0.01, 22 = 1, 11 = 0, 22 = 1, and weighted mean square error values of U1 = −2.16 × 10 −5 , U2 = 0. Here, the black, solid lines describe the numerical integration solutions of Equation (30), while the blue, dashed lines represent the numerical solutions obtained from the nonlinear equations of motion, (53) and (54).

Conclusions
A transformation approach has been proposed to find the equivalent representation form in the sense of Lyapunov of forced, damped, nonlinear, two degree-of-freedom dynamic systems. This proposed approach of finding equivalent expressions has the main advantage of providing simple algebraic relations whose coefficient values are determined by minimizing the error of replacing the original restoring forces by equivalent ones. By studying the dynamic responses of three nonlinear dynamic systems, the validity of the proposed conjecture has been numerically examined by comparing the LCEs of the original equations with those obtained from their equivalent uncoupled Based on these results, one can conclude that the improvement of the numerical prediction obtained from the equivalent representation forms will depend on the physical system, its nonlinearities, and its initial condition values.

Conclusions
A transformation approach has been proposed to find the equivalent representation form in the sense of Lyapunov of forced, damped, nonlinear, two degree-of-freedom dynamic systems. This proposed approach of finding equivalent expressions has the main advantage of providing simple algebraic relations whose coefficient values are determined by minimizing the error of replacing the original restoring forces by equivalent ones. By studying the dynamic responses of three nonlinear dynamic systems, the validity of the proposed conjecture has been numerically examined by comparing the LCEs of the original equations with those obtained from their equivalent uncoupled expressions. From the numerical predictions obtained, it is concluded that the LCEs and the LLE of the equivalent expression are similar to those obtained from the original equations of motion. In fact, numerical results have shown that the proposed procedure predicts the quantitative behavior of the nonlinear systems examined here well; however, there are some discrepancies which could be removed if additional terms in the power series expansion of the dynamical system's restoring forces are considered. Furthermore, when the system's restoring forces were replaced by a fifth-order truncated power series expansion in which decay rate terms were considered, numerical predictions computed from the nonlinear dynamic system of Example 1, indicated that the RMSE values were lower than those of the cubic truncated power series expansion, since, in this quintic approach, decay rate terms were considered. In this case, both the qualitative and the quantitative dynamic system behaviors in physical coordinates were predicted well in spite of having uncoupled modal equivalent oscillators, as shown in Figures 16  and 17.
Finally, and contrary to other proposed approaches, the one introduced here is capable of describing the dynamic behavior of nonlinear systems when internal resonance exists between the two modes, and thus, it is possible to uncouple, in modal coordinates, the resonant modes, as illustrated in Figure 12. Therefore, it can be concluded, in accordance with the numerical evidence of the dynamic systems examined here, that it is possible to have two decoupled Duffing-type equivalent representation forms of the original equations of motion in the sense of Lyapunov that can predict the dynamics of the original system, and accuracy can be improved if additional terms in the truncated power series are considered, not only to describe elastic forces, but also to take into account decay rate effects.
Of course, the system's restoring forces could be replaced by other forms rather than Duffing-type equivalent representation expressions.
However, the applicability of the proposed approach to other nonlinear dynamic systems by using truncated power series expansions or alternative forms to replace the system restoring forces could be considered only for small deviations of the system's linear behavior, because for large amplitudes, coupling and nonlinear terms are important.