A Transformation Method for Solving Conservative Nonlinear Two-Degree-of-Freedom Systems

Recently, we have introduced in [1, 2] a procedure that transforms the original equation of motion of single degree of freedom systems, having nonlinear restoring forces, into nonlinear ordinary differential equations of the Duffing type. Furthermore, this procedure has been used in conjunction with the Chebyshev polynomials to derive the equivalent representation form of nonlinear system with dissipative and sinusoidal driving forces [2] in which the numerical comparisons between the equivalent and the original equations of motion demonstrate the effectiveness of this method in predicting the corresponding system dynamics response. Motivated by these results, we herein expand the application of this technique to obtain the equivalent representation form of two-degree-of-freedom homogeneous, undamped, and ordinary differential equations of the form


Introduction
Recently, we have introduced in [1,2] a procedure that transforms the original equation of motion of single degree of freedom systems, having nonlinear restoring forces, into nonlinear ordinary differential equations of the Duffing type.Furthermore, this procedure has been used in conjunction with the Chebyshev polynomials to derive the equivalent representation form of nonlinear system with dissipative and sinusoidal driving forces [2] in which the numerical comparisons between the equivalent and the original equations of motion demonstrate the effectiveness of this method in predicting the corresponding system dynamics response.Motivated by these results, we herein expand the application of this technique to obtain the equivalent representation form of two-degree-of-freedom homogeneous, undamped, and ordinary differential equations of the form that models the dynamics response of systems that arise in physics and engineering applications [3][4][5][6], where  1 and  2 are the masses of the system,  11 ,  12 ,  21 , and  22 are stiffness parameters,  1 ( 1 ,  2 ) and  2 ( 1 ,  2 ) are nonlinear restoring forces, and  1 (0) =  10 , ẋ 1 (0) = ẋ 10 ,  2 (0) =  20 , and ẋ 2 (0) = ẋ 20 are the system initial conditions.
The method is based on replacing first the system restoring forces by polynomial expressions and, then, we use a transformation technique to replace the resulting equations by two uncoupled nonlinear differential equations of the Duffing type whose solutions are based on Jacobi elliptic functions.It is shown that the appropriate derivation of one of these solutions can improve the numerical estimates of the equivalent equations of motion even at larger and feasible nonlinear parameter values.

A Nonlinear Transformation Approach
The equivalent representation form of (1) is obtained by replacing it by two uncoupled, nonlinear ordinary differential equations of the Duffing type by following an approach similar to the one described in [1], This procedure yields 2

Mathematical Problems in Engineering
Here we assume that the coefficients   are determined by the following expressions: where and the values of   are fitted to satisfy (4) and ( 5) [7].Notice, in this case, that the system (2) can be solved in closed form by using Jacobi elliptic functions.This yields [8]  1 =  10  ( 1 , However, an alternative approach to find the equivalent representation form of (1) can be derived if we use the form (2) 1 or (2) 2 and then substitute its exact solution into the remaining equation in (1).In this case, the influence of the nonlinear effects of  1 or  2 into the system dynamics response will be taken into account by the resulting coupled equation of motion whose approximate solution could be based on the usage of Jacobi elliptic functions.
Other equivalent representation forms of (1) are possible since these depend on the polynomial order use to write the equivalent form of the system restoring forces [9][10][11].For instance, let us suppose that the nonlinear conservative forces of (1) are of the cubic type and rewrite (1) in its canonical representation form [12]: with initial conditions given as  1 (0) =  10 , u 1 (0) = u 10 ,  2 (0) =  20 , and u 2 (0) = u 20 .Here the dots denote the time derivative,  1 and  2 are known as the system normal coordinates,  is a nonlinear parameter, and  1 ,  2 , and  1 through  8 are corresponding system parameters that are defined in accordance with the physics of the system.Now, let us replace (7) by their corresponding equivalent equations of motion of the Duffing type which are valid for the complete range of oscillation amplitudes [1], by using the relations where ] 1 , ] 2 , ] 11 , and ] 22 are fitted constants that satisfy (8), and the parameters   are determined by solving the equations This yields Then, (7) can be rewritten in equivalent form as whose exact solutions are similar to those given by (6).We shall next explore the applicability of the nonlinear transformation method in obtaining the equivalent representation form of (1) and, then, we will compare the numerical integration solutions of five dynamics systems, having nonlinear restoring forces with rational or irrational terms, with respect to their equivalent representation forms [13][14][15][16][17][18][19][20].First, let us consider the case for which the restoring forces are of the cubic type.
Table 2: RMSE values computed by considering the canonical representation form of ( 13) with ] 1 = −9, ] 2 = −0.215,] 11 = −2.325,and  1, that the computed values of the RMSE become smaller for increasing values of the nonlinear parameter  which shows that our solution procedure can be used to obtain solutions that follow the numerical integration solutions of the original equations of motion even at larger values of .For illustrative purposes, we have plotted in Figure 1 the amplitude-time response curves obtained from (13), and those provided by substituting (14) 1 into (13) 2 .We have selected the value of  = 20 which provides a highly nonlinear dynamic system response.Notice from Figure 1 that both solutions are almost the same.In fact, the RMSE values do not exceed 0.1179, for  1 , and 0.0233, for  2 .In Figure 1, the gray solid lines describe the numerical integration solutions of ( 13), while the blue dashed lines represent the numerical solutions obtained from our derived equivalent nonlinear equations of motion.Table 2 shows the RMSE values obtained by using the canonical representation form of (13).One can notice from Table 2 that the smallest RMSE values were computed at decreasing values of the nonlinear parameter  when the equivalent form (11) 1 was substituted into (7) 2 (the fifth and the sixth columns).Also, one can notice that when  = 0, (7) describes the oscillatory behavior of the resulting linear dynamics system; however, this is not the case when we use the equivalent representation form of (13) given by ( 14) because the terms  2 and  4 are different from zero.This situation increases the RMSE errors as listed in Table 1.In an attempt to further quantify the accuracy of our proposed procedure, we next introduce the global error estimation based on a technique similar to the ones developed in [14][15][16] in which the system global error, GE, can be determined from the following expression: where  1 and  2 represent the values obtained by numerically integrating the equivalent equations of motion and  1 and  2 are the numerical values computed by using the original system equations of motion.Figure 2 displays the global error curve plotted versus time by considering the numerical integration solutions of both, the resulting equation obtained by substituting (14) 1 into (13) 2 and the original equations of motion (13).One can notice from Figure 2, that the global error tends to grow linearly at increasing time values, as predicted by Calvo and Hairer in [14].
To further assess the accuracy of our proposed nonlinearization method, we will next examine a hyperelastic shear suspension system.

Example 2: Hyperelastic Shear Suspension System.
Here the equations of motion of a linear absorber attached to a rigid body that is supported symmetrically by incompressible, homogeneous, and isotropic hyperelastic shear blocks are given as with initial conditions given as (0) =  10 , σ (0) = 0, (0) =  10 , and ż (0) = 0.Here  and  denote, respectively, the simple shear deformation of the load and the motion of the linear absorber system,  describes the nonlinear material response effects,  2 is a tuning system parameter, and ,  2 , and  1 are defined as The details of the derivation of the equations of motion given by ( 19) are provided in [2].Next, we set  1 = / 10 and  2 = / 10 and rewrite (12) in its normalized form with  1 (0) = 1, ẏ 1 (0) = 0,  2 (0) = 1, and ẏ 2 (0) = 0. Notice that the canonical representation form of (21) can be obtained by using the transformation (16); this yields exactly (7), where the system parameters are given by the following relations: Since (21) 2 is only coupled with  1 through the linear term − 2 2  10 /( 10 ) 1 , we slightly modify our procedure and assume that the equivalent representation form of ( 21) is given as .
where the parameters   , computed from (4), are given as  21) is considered, their approximate solutions, that can be obtained from (7), are given by (6).Tables 3 and 4 show the RMSE values computed from the numerical integration solutions of (21) and from the derived equivalent representation forms.Here, we have used as initial conditions the values of  10 (0) = 1, σ (0) = 0,  10 (0) = −1, and ż (0) = 0, with  = 1/5 and  2 = 1.One can notice from Table 3 that the lowest RMSE values are tabulated in the fifth and the sixth columns at values of  > 0.25.In fact, as the value of  tends to increase, the RMSE values tend to decrease.Figure 3 shows the amplitude-time response curves computed by substituting (23) 1 into (21) 2 .In this case, we have considered that the nonlinear parameter  has the value of 50.Notice that the numerical integration solutions of (21) are almost indistinguishable from those computed by using the equivalent transformation forms.Here the gray solid lines represent the numerical integration solutions of (21), while the blue dashed lines described the numerical integration solutions obtained from the equivalent equations of motion.Therefore, we can conclude that our equivalent transformation approach provides accurate solutions to (21) if the form (23) 1 is used to derive the approximate solution of  2 .
Also notice from Table 4 that the smallest RMSE values, computed when the canonical form (7) of ( 19) is used to derive its equivalent representation form, are on 0 <  ≤ 0.25, that is, when the equivalent form (11) 1 is substituted into (7) 2 .Therefore, if one wants to have solutions that describe well the numerical integration solutions of (21), we need to use their canonical representation form at values of  on the interval 0 <  ≤ 0.25.For larger nonlinear system parameter values, good predictions of the original equation of motion are Here, the gray solid lines describe the numerical integration solutions of (21), while the blue dashed lines represent the numerical solutions obtained from the derived equivalent nonlinear equations of motion.In this case, the computed RMSE values were 0.1006, for  1 , and 0.0102, for  2 .

Example 3:
A System with Three Nonlinear Springs.We now derive the equivalent representation form of two-mass system with three nonlinear springs introduced in [17] and whose differential equation of motion is given by with initial conditions  1 (0) =  10 , ẋ 1 (0) = 0,  2 (0) =  20 , and ẋ 2 (0) = 0.If we introduce the transformations  1 = ( 1 +  0 )/2 and  2 = ( 0 −  1 )/2, thus (25) can be written as with initial conditions given as  0 (0) =  00 , ẏ 0 (0) = ẏ 00 ,  1 (0) =  10 , and ẏ 1 (0) = ẏ 10 .Notice that (26) has the same form of (7) with system parameters given as Therefore and in accordance with our transformation approach, (26) are transformed into (11) in which To assess the accuracy of the equivalent equations of motion when compared to the original ones (25), let us consider the parameter values of  1 =  2 =  3 =  4 = 1,  = 1, and  = −2.9 and initial conditions  0 (0) = 0.5, ẏ 0 (0) = 0,  1 (0) = 1.5, and ẏ 1 (0) = 0. Figure 4 shows the amplitude-time response curves obtained by numerical integration of the corresponding equations of motion by using the Runge-Kutta method.One can notice from Figure 4 that the amplitudetime response curves obtained from (11) provide periodic solutions that cannot capture the high frequency components exhibited by the numerical amplitude-time curves plotted from (25).However, if we substitute the exact solution of (11) 1 into (7) 2 and plot its corresponding numerical integration solution, the high frequency components are captured as illustrated in Figure 5. Therefore, we can conclude that if during the dynamic analysis processes the motion of the coordinate  1 exhibits harmonic behavior at lower angular frequency values, then the dynamics response behavior of this coordinate solution can be used to obtain the high frequency components exhibited by the remaining coordinate.In fact and for the nonlinear parameter value of  = −2.9, the RMSE values, for the time interval shown in Figure 5, do not exceed 0.3159 and 0.1792 for  0 and  1 , respectively.
We next increase the system nonlinearity at the value of  = −3.1 and plot the corresponding amplitude versus time response curves.Notice from Figure 6 that when  ≥ 4, the numerical integration solution of the equations of motion (25) obtained by using the Runge-Kutta numerical method fails in spite of using the several solver options provided by  Figure 4: Amplitude-time response curves computed from the numerical integration solutions of (26) and those provided by (11).The parameter values used to obtain these plots were  = 1,  1 =  2 =  3 =  4 = 1, and  = −2.9, with  0 (0) = 0.5, ẏ 0 (0) = 0,  1 (0) = 1, ẏ 1 (0) = 0, ] 11 = 1.1935, and ] 2 = 1.1935.Here, the gray solid lines describe the numerical integration solutions of (26), while the blue dashed lines represent the numerical solutions obtained from the derived equivalent nonlinear equations of motion.In this case, the computed RMSE values were 0.3159, for  0 , and 0.0358, for  1 .Figure 5: Amplitude-time response curves computed from the numerical integration solutions of (26) and those provided by substituting (11) 1 into (7) 2 .The parameter values used to obtain these plots were  = 1,  1 =  2 =  3 =  4 = 1, and  = −2.9, with  0 (0) = 0.5, ẏ 0 (0) = 0,  1 (0) = 1, ẏ 1 (0) = 0, ] 11 = 1.1935, and ] 2 = 1.1935.Here, the gray solid lines describe the numerical integration solutions of (26), while the blue dashed lines represent the numerical solutions obtained from the derived equivalent nonlinear equations of motion.In this case, the computed RMSE values were 0.3159, for  0 , and 0.1792, for  1 .Mathematica 9.0 or the MATLAB V.2012a computer packages.In an attempt to capture the dynamics system response with this nonlinear value, we have numerically solved (25) by using the Enhanced Multistage homotopy perturbation method (EMHPM) introduced in [18].This technique also fails at values of  ≥ 4, as illustrated in Figure 6 by the black and the purple dots.However, the predicted amplitude versus time response curves computed by substituting (11) 2 into (7) 1 exhibits periodic behavior for the time interval shown in Figure 6.We can conclude that this transformation approach can help to obtain the numerical integration solutions of nonlinear dynamics systems in which the numerical integrations of the original equations of motion could fail because of the convergence of the corresponding numerical method.

Example 4: Finite Extensibility Nonlinear Elastic
Absorber.We now focus our attention on finding the approximate solution of a finite extensibility nonlinear elastic (FENE) absorber attached to a nonlinear primary system [19], whose equations of motions are of the form and the nonlinear restoring forces are given by Here, the gray solid lines and the black and the purple dots describe the numerical integration solutions of (26) by using, respectively, the Runge-Kutta and the EMHPM methods, while the blue dashed lines representing the numerical solutions obtained from the derived equivalent nonlinear equations of motion.

Example 5:
A Two-Degree-of-Freedom System with Irrational Restoring Force.As a final example, let us consider the system shown in Figure 9 having two masses,  1 and  2 which are connected by four elastic springs of stiffness  1 and  2 having an undeformed length,  0 .Furthermore, the element with mass  1 is moving on a smooth horizontal bearing surface whose motion is restricted by an elastic linear spring of stiffness  3 .By using Newton second law, it is easy to show that the equations of motion are where the irrational restoring forces are given by Here The system (41) can be cast in different form by setting  1 = ( 1 −  2 )/2 and  2 = ( 2 +  1 )/2, and after some algebraic manipulations, we get that  with initial conditions  1 (0) =  10 , ẇ 1 (0) = 0,  2 (0) =  20 , and ẇ 2 (0) = 0.If we set  1 =  1 / 10 and  2 =  2 / 20 , then (43) becomes where with  1 (0) = 1, ż 1 (0) = 0, ż 2 (0) = 1, and  2 (0) = 0.By using the Chebyshev polynomials of the first kind, the cubic-like representation form of the irrational restoring force becomes where Here (− 2 ) and (− 2 ) are the complete elliptic integral of the first and second kind for the modulus , respectively.Thus, (44) can be rewritten as Application of the proposed nonlinear transformation approach to (48) yields where the   are given as Figure 10 illustrates the numerical predictions obtained by substituting (49) 1 into (48) 2 , as well as those computed from the numerical integration solution of (48).In Example 5, we have used the following system parameter values:

Conclusions
We have introduced a nonlinear transformation approach to determine the equivalent representation form of two degreeof-freedom oscillatory systems with nonlinear restoring forces.It has been shown that this nonlinear transformation approach decoupled the system equation into two Duffing equations.It has been found that the numerical integration of the equivalent equations of motion, of the five dynamic systems examined here, follows well the computed numerical values obtained from the equations of motion.Furthermore, it has been shown that by substituting the decouple equation related to the  1 ( 1 ) coordinate system into the original equation of motion described by  2 ( 2 ), it provides the lower RMSE values.However, when the canonical representation form of the original equations of motion is used to decouple the corresponding equations, the smallest RMSE values were found at small values of the nonlinear parameter .We also have found in Example 3 that our equivalent representation form can capture high frequency system components if the exact solution of the cubic-like representation form of the coordinate  1 is used to find the dynamics response of the  0 coordinate.Besides, the corresponding equivalent equations of motion can be numerically integrated at the value of  = −3.1 while the numerical integration solution of the original equations of motion (25), by using the Runge-Kutta, fails in spite of using the several solver options provided by Mathematica 9.0 or the MATLAB V.2012a computer packages.We next applied the EMHPM to find the numerical solution of (25) but it was found that this technique provides numerical results that diverge when  ≥ 4.However, the numerical integration of our equivalent representation form provides solutions that describe the system dynamics periodic behavior.In the case of Examples 4 and 5, we believe that the RMSE and the global error magnitudes are mainly due to two factors: (a) the Chebyshev polynomial representation of the system restoring forces and (b) the equivalent representation form of the original equations of motion that are based on cubic-like Duffing's equations.These two factors affect the accuracy achieved during the solution process; however, the quantitative and qualitative system response are captured with good degree of accuracy.
Therefore, it can be concluded, in accordance with some of the systems studied here, that if one wants to have the smallest RMSE values when using the proposed nonlinear transformation approach, the canonical representation form has to be used for weak nonlinear system; however, for larger values of , the accurate results are obtained when the original equations of motion were decoupled.Based on the numerical simulation results, it is clear that our method is a promising technique that could provide equivalent representation forms of dynamic systems with two or more degrees of freedom that can follow the numerical integration solutions of the original equations of motion, at larger nonlinear parameter values.

Figure 9 :
Figure 9: Schematic of a two-degree-of-freedom system with irrational restoring force.

Table 5 :
(8)E values computed by using the equivalent representation form of (31) with  1 = −1.991, 2 = −1.115,11=−0.7,and22=−4.6. 2 follow the amplitude-time response curves computed from (48) with RMSE values of 0.2095 for  1 and 0.1318 for  2 , respectively.The system global error as illustrated in Figure11tends to linearly grow as time increases.Of course, other values of the system initial conditions can be used to study the quantitative and qualitative system dynamics response; however, when the system parameter changes, it could be necessary to adjust the values of   and   in order to satisfy(8).