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

Abstract

A nonlinear transformation approach based on a cubication method is developed to obtain the equivalent representation form of conservative two-degree-of-freedom nonlinear oscillators. It is shown that this procedure leads to equivalent nonlinear equations that describe well the numerical integration solutions of the original equations of motion.

1. 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 𝑚₁00𝑚₂𝑥̈₁𝑥̈₂+𝑘₁₁𝑘₁₂𝑘₂₁𝑘₂₂𝑥₁𝑥₂+𝑓₁𝑥₁,𝑥₂𝑓₂𝑥₁,𝑥₂=0 (1) that models the dynamics response of systems that arise in physics and engineering applications [3–6], where 𝑚₁ and 𝑚₂ are the masses of the system, 𝑘₁₁, 𝑘₁₂, 𝑘₂₁, and 𝑘₂₂ are stiffness parameters, 𝑓₁(𝑥₁,𝑥₂) and 𝑓₂(𝑥₁,𝑥₂) are nonlinear restoring forces, and 𝑥₁(0)=𝑥₁₀, 𝑥̇₁(0)=𝑥̇₁₀, 𝑥₂(0)=𝑥₂₀, and 𝑥̇₂(0)=𝑥̇₂₀ 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.

2. 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 𝑥̈₁+𝑎₁𝑥₁+𝑎₂𝑥³₁=0;𝑥̈₂+𝑎₃𝑥₂+𝑎₄𝑥³₂=0. (2) Here we assume that the coefficients 𝑎_𝑖 are determined by the following expressions: 𝜕𝐹₁𝑎₁,𝑎₂𝜕𝑎₁=0,𝜕𝐹₁𝑎₁,𝑎₂𝜕𝑎₂=0,𝜕𝐹₂𝑎₃,𝑎₄𝜕𝑎₃=0,𝜕𝐹₂𝑎₃,𝑎₄𝜕𝑎₄=0, (3) where 𝐹₁𝑎₁,𝑎₂=^𝑢2₀^𝑢1₀𝑘₁₁𝑥₁+𝑘₁₂𝑥₂+𝑓₁𝑥₁,𝑥₂𝑚₁−𝑎₁𝑥₁−𝑎₂𝑥³₁²𝑑𝑥₁𝑑𝑥₂⟶min,𝐹₂𝑎₃,𝑎₄=^𝑢4₀^𝑢3₀𝑘₂₁𝑥₁+𝑘₂₂𝑥₂+𝑓₂𝑥₁,𝑥₂𝑚₂−𝑎₃𝑥₂−𝑎₄𝑥³₂²𝑑𝑥₂𝑑𝑥₁⟶min (4)(5) 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] 𝑥₁=𝑥₁₀𝑐𝑛𝜔₁𝑡,𝜅²₁,𝜔₁=𝑎₁+𝑎₂𝑥²₁₀,𝜅²₁=𝑎₂𝑥²₁₀2𝜔²₁,𝑥₂=𝑥₂₀𝑐𝑛𝜔₂𝑡,𝜅²₂,𝜔₂=𝑎₃+𝑎₄𝑥²₂₀,𝜅²₂=𝑎₃𝑥²₂₀2𝜔²₂. (6)

However, an alternative approach to find the equivalent representation form of (1) can be derived if we use the form (2)₁ or (2)₂ and then substitute its exact solution into the remaining equation in (1). In this case, the influence of the nonlinear effects of 𝑥₁ or 𝑥₂ 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–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]: 𝑢̈₁𝑢̈₂+𝜔²_𝑛100𝜔²_𝑛2𝑢₁𝑢₂𝜑+𝜀¹𝑢³₁+𝜑₂𝑢²₁𝑢₂+𝜑₃𝑢₁𝑢²₂+𝜑₄𝑢³₂𝜑₅𝑢³₁+𝜑₆𝑢²₁𝑢₂+𝜑₇𝑢₁𝑢²₂+𝜑₈𝑢³₂=0, (7) with initial conditions given as 𝑢₁(0)=𝑢₁₀, 𝑢̇₁(0)=𝑢̇₁₀, 𝑢₂(0)=𝑢₂₀, and 𝑢̇₂(0)=𝑢̇₂₀. Here the dots denote the time derivative, 𝑢₁ and 𝑢₂ are known as the system normal coordinates, 𝜀 is a nonlinear parameter, and 𝜔_𝑛1, 𝜔_𝑛2, and 𝜑₁ through 𝜑₈ 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 𝐹₁𝑏₁,𝑏₂=^𝜈2₀^𝜈1₀𝜔²_𝑛1𝑢₁+𝜀𝜑₁𝑢³₁+𝜑₂𝑢²₁𝑢₂+𝜑₃𝑢₁𝑢²₂+𝜑₄𝑢³₂−𝑏₁𝑢₁−𝑏₂𝑢³₁²𝑑𝑢₁𝑑𝑢₂𝐹⟶min,₂𝑏₃,𝑏₄=^𝜈11₀^𝜈22₀𝜔²_𝑛2𝑢₂+𝜀𝜑₅𝑢³₁+𝜑₆𝑢²₁𝑢₂+𝜑₇𝑢₁𝑢²₂+𝜑₈𝑢³₂−𝑏₃𝑢₂−𝑏₄𝑢³₂²𝑑𝑢₂𝑑𝑢₁⟶min, (8) where 𝜈₁, 𝜈₂, 𝜈₁₁, and 𝜈₂₂ are fitted constants that satisfy (8), and the parameters 𝑏_𝑖 are determined by solving the equations 𝜕𝐹₁𝑏₁,𝑏₂𝜕𝑏₁=0,𝜕𝐹₁𝑏₁,𝑏₂𝜕𝑏₂=0,𝜕𝐹₂𝑏₃,𝑏₄𝜕𝑏₃=0,𝜕𝐹₂𝑏₃,𝑏₄𝜕𝑏₄=0. (9)

This yields 𝑏₁=−−30𝜑₂𝜈₁𝜈₂−64𝜑₃𝜈₁𝜈²₂−135𝜑₄𝜈³₂−192𝜈₁𝜔²_𝑛1192𝜈₁,𝑏₂=−−192𝜑₁𝜈³₁−70𝜑₂𝜈²₁𝜈₂+105𝜑₄𝜈³₂192𝜈³₁,𝑏₃=−−135𝜑₅𝜈³₁₁−64𝜑₆𝜈²₁₁𝜈₂₂−30𝜑₇𝜈₁₁𝜈²₂₂−192𝜈₂₂𝜔²_𝑛2192𝜈₂₂,𝑏₄=−−105𝜑₅𝜈³₁₁−70𝜑₇𝜈₁₁𝜈²₂₂−192𝜑₈𝜈³₂₂192𝜈³₂₂. (10) Then, (7) can be rewritten in equivalent form as 𝑢̈₁+𝑏₁𝑢₁+𝑏₂𝑢³₁=0,𝑢̈₂+𝑏₃𝑢₂+𝑏₄𝑢³₂=0, (11) 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–20]. First, let us consider the case for which the restoring forces are of the cubic type.

2.1. Example 1: A Cubic Nonlinear System

Let us consider the conservative system built with a nonlinear spring of the cubic type. The corresponding equations of motion are given as [13] 𝑚₁00𝑚₂𝑥̈₁𝑥̈₂+𝑘₁+𝑘₂−𝑘₂−𝑘₂2𝑘₂𝑥₁𝑥₂+𝜀𝑥³₁0=0, (12) with initial conditions 𝑥₁(0)=𝑥₁₀, 𝑥̇₁(0)=0, 𝑥₂(0)=𝑥₂₀, and 𝑥̇₂(0)=0. If we introduce the transformations 𝑦₁=𝑥₁/𝑥₁₀ and 𝑦₂=𝑥₂/𝑥₂₀ so that 𝑦₁(0)=1, 𝑦̇₁(0)=0, 𝑦₂(0)=1, 𝑦̇₂(0)=0, thus (12) can be written in normalized form as 10𝑦̈01𝑦̈¹₂+𝑘₁+𝑘₂𝑚₁−𝑘₂𝑥₂₀𝑚₁𝑥₁₀−𝑘₂𝑥₁₀𝑚₂𝑥₂₀2𝑘₂𝑚₂𝑦₁𝑦₂+𝜀𝑥²₁₀𝑦³₁𝑚₁0=0. (13) If we use our proposed nonlinear transformation approach then the equivalent representation form of (13) is given as 𝑦̈₁+𝑎₁𝑦₁+𝑎₂𝑦³₁=0;𝑦̈₂+𝑎₃𝑦₂+𝑎₄𝑦³₂=0, (14) where the 𝑎_𝑖, determined from (4) and (5), are given as 𝑎₁=−−32𝑘₁𝑢₁𝑥₁₀−32𝑘₂𝑢₁𝑥₁₀+45𝑘₂𝑢₂𝑥₂₀32𝑚₁𝑢₁𝑥₁₀,𝑎₂=−−32𝜀𝑢³₁𝑥³₁₀−35𝑘₂𝑢₂𝑥₂₀32𝑚₁𝑢³₁𝑥₁₀,𝑎₃=−45𝑘₂𝑢₁₁𝑥₁₀−32𝑘₂𝑢₂₂𝑥₂₀32𝑚₂𝑢₂₂𝑥₂₀,𝑎₄=35𝑘₂𝑢₁₁𝑥₁₀32𝑚₂𝑢³₂₂𝑥₂₀, (15) and 𝑢₁, 𝑢₂, 𝑢₁₁, and 𝑢₂₂ are the corresponding fitting parameters described earlier.

To apply our nonlinear transformation approach to the canonical representation form of (13), we first use the following linear coordinate transformation [12]: 𝑥₁𝑥₂=𝑅₁𝑅₂𝑅₁𝑓₁𝑅₂𝑓₂𝑢₁𝑢₂, (16) which yields the canonical form (7). Here, the system parameters of (7) are described by the relations 𝜔²_𝑛1=12𝑚₁𝑘₁+3𝑘₂−𝑘²₁−2𝑘₁𝑘₂+5𝑘²₂,𝜔²_𝑛2=12𝑚₁𝑘₁+3𝑘₂+𝑘²₁−2𝑘₁𝑘₂+5𝑘²₂,𝜑₁≡−𝜀𝑅⁴₁𝑥²₁₀,𝜑₂≡−3𝜀𝑅³₁𝑅₂𝑥₁₀𝑥₂₀,𝜑₃≡−3𝜀𝑅²₁𝑅²₂𝑥²₂₀,𝜑₄≡−𝜀𝑅₁𝑅³₂𝑥³₂₀𝑥₁₀,𝜑₅≡−𝜀𝑅³₁𝑅₂𝑥³₁₀𝑥₂₀,𝜑₆≡−3𝜀𝑅²₁𝑅²₂𝑥²₁₀,𝜑₇≡−3𝜀𝑅₁𝑅³₂𝑥₁₀𝑥₂₀,𝜑₈≡−𝜀𝑅⁴₂𝑥²₂₀,𝑓₁≡𝑘₁+𝑘₂−𝑚₁𝜔²_𝑛1𝑘₂,𝑓₂≡𝑘₁+𝑘₂−𝑚₁𝜔²_𝑛2𝑘₂,𝑅²₁≡1𝑚₁1+𝑓²₁,𝑅²₂≡1𝑚₁1+𝑓²₂. (17) Thus, the approximate solutions of (13), by considering their canonical representation form, are given by expressions similar to those of (6).

Tables 1 and 2 show the root-mean-square error (RMSE) values computed by comparing the numerical integration solutions of the original equations of motion (13), with respect to numerical results computed from their derived equivalent representation forms by considering the initial conditions values of 𝑥₁(0)=2, 𝑥̇₁(0)=0, 𝑥₂(0)=−1, and 𝑥̇₂(0)=0, the parameter values of 𝑚₁=𝑚₂=1, 𝑘₁=3, and 𝑘₂=1, and different values of 𝜀. In the second and third columns of Table 1 are tabulated the RMSE values computed by using the equivalent representation forms of (13) given by (14), while the fourth and the fifth columns provide the RMSE values obtained by substituting (14)₁ into (13)₂. The sixth and the seventh columns of Table 1 contain the RMSE values computed by substituting (14)₂ into (13)₁. As we can see from Table 1, the lowest RMSE values are tabulated in the fifth and the sixth columns which implies that the most accurate solutions of (13) are obtained when the equivalent representation form (14)₁ is used to derive the approximate solution of 𝑥₂. Also notice from Table 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)₁ into (13)₂. 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 𝑥₁, and 0.0233, for 𝑥₂. 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 1 

RMSE values computed from the equivalent representation forms of (13) with 𝑢₁=1.55, 𝑢₂=1.5, 𝑢₃=1.25, and 𝑢₄=1.6. Here 𝑏₁₁=(𝑘₁+𝑘₂)/𝑚₁, 𝑏₁₂=𝑏₂₁=−𝑘₂𝑥₂₀/(𝑚₁𝑥₁₀), 𝑏₂₂=2𝑘₂/𝑚₂, and 𝑑₁=𝜀𝑥²₁₀/𝑚₁.

Table 2 

RMSE values computed by considering the canonical representation form of (13) with 𝜈₁=−9, 𝜈₂=−0.215, 𝜈₁₁=−2.325, and 𝜈₂₂=−8.5.

(a)

(b)

Figure 1 

Amplitude-time response curves computed from the numerical integration solutions of (13) and those provided by substituting (14)1 into (13)2. The parameter values used to obtain these plots were 𝑚1=𝑚2=1, 𝑘1=3, 𝑘2=1, and 𝜀=20, with 𝑥1(0)=2, 𝑥̇1(0)=0, 𝑥2(0)=−1, and 𝑥̇2(0)=0, and 𝑢1=1.55, 𝑢2=1.5, 𝑢3=1.25, and 𝑢4=1.6. Here, the gray solid lines describe the numerical integration solutions of (13), while the blue dashed lines represent the numerical solutions obtained from the 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)₁ was substituted into (7)₂ (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 𝑎₂ and 𝑎₄ 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–16] in which the system global error, GE, can be determined from the following expression: GE=𝑦₁(𝑡)−𝑝₁(𝑡)+𝑦₂(𝑡)−𝑝₂(𝑡), (18) where 𝑦₁ and 𝑦₂ represent the values obtained by numerically integrating the equivalent equations of motion and 𝑝₁ and 𝑝₂ 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)₁ into (13)₂ 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].

Figure 2 

Global error behavior.

To further assess the accuracy of our proposed nonlinearization method, we will next examine a hyperelastic shear suspension system.

2.2. 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 𝜎̈𝑧̈+100𝛽1−𝛼²₂−𝛼²₂𝛼²₂𝜎𝑧−𝜀𝑏₁𝜎³0=0, (19) with initial conditions given as 𝜎(0)=𝜎₁₀, 𝜎̇(0)=0, 𝑧(0)=𝑧₁₀, 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, 𝑟₂ is a tuning system parameter, and 𝛽, 𝛼₂, and 𝑏₁ are defined as 𝑚𝛽=𝑀,𝑏₁=𝛼²₂−1,𝛼²₂=𝛽𝑟²₂𝛽𝑟²₂.+1 (20) The details of the derivation of the equations of motion given by (19) are provided in [2]. Next, we set 𝑦₁=𝜎/𝜎₁₀ and 𝑦₂=𝑧/𝑧₁₀ and rewrite (12) in its normalized form 10𝑦̈01𝑦̈¹₂+𝛼1−²²𝑧₁₀𝜎₁₀−𝛼²₂𝜎₁₀𝛽𝑧₁₀𝛼²₂𝛽𝑦₁𝑦₂−𝜀𝜎²₁₀𝑦³₁0=0 (21) with 𝑦₁(0)=1, 𝑦̇₁(0)=0, 𝑦₂(0)=1, and 𝑦̇₂(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: 𝜔²_𝑛1=12𝛼²₂𝛽+1−𝛼²₂𝛽+1²+4𝛼²₂𝑏₁𝛽,𝜔²_𝑛2=12𝛼²₂𝛽+1+𝛼²₂𝛽+1²+4𝛼²₂𝑏₁𝛽,𝜑₁≡−𝑏₁𝜀𝑅⁴₁𝜎²₁₀,𝜑₂≡−3𝑏₁𝜀𝑅³₁𝑅₂𝜎₁₀𝑧₂₀,𝜑₃≡−3𝑏₁𝜀𝑅²₁𝑅²₂𝑧²₂₀,𝜑₄𝑏≡−¹𝜀𝑅₁𝑅³₂𝑧³₂₀𝜎₁₀,𝜑₅𝑏≡−¹𝜀𝑅³₁𝑅₂𝜎³₁₀𝑧₂₀,𝜑₆≡−3𝑏₁𝜀𝑅²₁𝑅²₂𝜎²₁₀,𝜑₇≡−3𝑏₁𝜀𝑅₁𝑅³₂𝜎₁₀𝑧₂₀,𝜑₈≡−𝑏₁𝜀𝑅⁴₂𝑧²₂₀,𝑓₁≡1−𝜔²_𝑛1𝛼²₂,𝑓₂≡1−𝜔²_𝑛2𝛼²₂,𝑅²₁≡11+𝑓²₁𝛽,𝑅²₂≡11+𝑓²₂𝛽. (22) Since (21)₂ is only coupled with 𝑦₁ through the linear term −𝛼²₂𝜎₁₀/(𝛽𝑧₁₀)𝑦₁, we slightly modify our procedure and assume that the equivalent representation form of (21) is given as 𝑦̈₁+𝑎₁𝑦₁+𝑎₂𝑦³₁=0;𝑦̈₂+𝑎₃𝑦₂=0, (23) where the parameters 𝑎_𝑖, computed from (4), are given as 𝑎₁=−−32𝑢₁𝜎₁₀+45𝛼²₂𝑢₂𝑧₁₀32𝑢₁𝜎₁₀,𝑎₂=−−32𝑏₁𝜀𝑢³₁𝜎³₁₀−35𝛼²₂𝑢₂𝑧₁₀32𝑢³₁𝜎₁₀,𝑎₃=−𝛼²₂𝑢₁₁𝜎₁₀+𝛼²₂𝑢₂₂𝑧₁₀𝛽𝑢₂₂𝑧₁₀, (24) and 𝑢₁, 𝑢₂, 𝑢₁₁, and 𝑢₂₂ are the corresponding fitting parameters. The assumption that (21)₂ can be transformed into a linear equation of the form (23)₂ is based on the fact that in (21)₂, there is not 𝑦₂(𝑡) terms of the cubic type.

Next, if the canonical representation form of (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 𝜎₁₀(0)=1, 𝜎̇(0)=0, 𝑧₁₀(0)=−1, and 𝑧(̇0)=0, with 𝛽=1/5 and 𝑟₂=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)₁ into (21)₂. 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)₁ is used to derive the approximate solution of 𝑦₂.

Table 3 

RMSE values computed from the equivalent representation forms of (21) with 𝑢₁=−0.4, 𝑢₂=0.1, 𝑢₃=−0.01, and 𝑢₄=−0.315. Here 𝑏₁₁=1, 𝑏₁₂=𝑏₂₁=−𝛼²₂𝜎₁₀/(𝛽𝑧₁₀), 𝑏₂₂=𝛼²₂/𝛽, and 𝑑₁=−𝜀𝜎²₁₀.

Table 4 

RMSE values computed by considering the canonical representation form of (21) with 𝜈₁=0.85, 𝜈₂=0.45, 𝜈₁₁=−0.01, and 𝜈₂₂=0.1.

(a)

(b)

Figure 3 

Amplitude-time response curves computed from the numerical integration solutions of (21) and those provided by substituting (23)1 into (21)2. The parameter values used to obtain these plots were 𝛽=1/5 and 𝑟2=1, 𝜀=50, with 𝜎10(0)=1, 𝜎̇(0)=0, 𝑧10(0)=−1, 𝑧(̇0)=0, and 𝑢1=−0.4, 𝑢2=0.1, 𝑢3=−0.01, and 𝑢4=−0.315. 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.

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)₁ is substituted into (7)₂. 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 obtained when the equivalent form (23)₁ is substituted into (21)₂.

2.3. 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 𝑚0𝑥̈0𝑚𝑥̈¹₂+𝑘₁+𝑘₂−𝑘₂−𝑘₂𝑘₁+𝑘₂𝑥₁𝑥₂+𝑘₃𝑥₁−𝑥₂³+𝜀𝑘₄𝑥³₁𝑘₃𝑥₂−𝑥₁³+𝜀𝑘₄𝑥³₂=0, (25) with initial conditions 𝑥₁(0)=𝑥₁₀, 𝑥̇₁(0)=0, 𝑥₂(0)=𝑥₂₀, and 𝑥̇₂(0)=0. If we introduce the transformations 𝑥₁=(𝑦₁+𝑦₀)/2 and 𝑥₂=(𝑦₀−𝑦₁)/2, thus (25) can be written as 𝑦̈₀𝑦̈₁+𝜔²_𝑛100𝜔²_𝑛2𝑦₀𝑦₁+𝜀𝑘₄(4𝑚)𝑦³₀+𝜀𝑘₄(4𝑚)𝑦₀𝑦²₁𝜀𝑘₄(4𝑚)𝑦²₀𝑦₁+𝛽+𝜀𝑘₄()𝑦4𝑚³¹=0, (26) with initial conditions given as 𝑦₀(0)=𝑦₀₀, 𝑦̇₀(0)=𝑦̇₀₀, 𝑦₁(0)=𝑦₁₀, and 𝑦̇₁(0)=𝑦̇₁₀. Notice that (26) has the same form of (7) with system parameters given as 𝜔²_𝑛1=𝑘₁𝑚,𝜔²_𝑛2=𝑘₁+2𝑘₂𝑚,𝜑₁=𝜑₃=𝜀𝑘₄(),𝜑4𝑚₂=𝜑₄=0,𝜑₆=𝜀𝑘₄(),𝜑4𝑚₈=𝛽+𝜀𝑘₄4𝑚,𝜑₅=𝜑₇=0,𝛽=2𝑘₃𝑚. (27) Therefore and in accordance with our transformation approach, (26) are transformed into (11) in which 𝑏₁=𝜑₃𝜈²₂+3𝜔²_𝑛13;𝑏₂=𝜑₁;𝑏₃=𝜑₆𝜈²₁₁+3𝜔²_𝑛23;𝑏₄=𝜑₈. (28) 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, 𝑚=1, and 𝜀=−2.9 and initial conditions 𝑦₀(0)=0.5, 𝑦̇₀(0)=0, 𝑦₁(0)=1.5, and 𝑦̇₁(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 amplitude-time 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)₁ into (7)₂ 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 𝑦₁ 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 𝑦₀ and 𝑦₁, respectively.

(a)

(b)

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.

(a)

(b)

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.

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 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)₂ into (7)₁ 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.

(a)

(b)

Figure 6 

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 𝜀=−3.1, 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 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.

2.4. 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 𝑚₁00𝑚₂𝑥̈₁𝑥̈₂+𝑘₁0𝑥00𝑥¹₂+𝐹₁𝑥₁,𝑥₂𝐹₂𝑥₁,𝑥₂=0, (29) and the nonlinear restoring forces are given by 𝐹₁𝑥₁,𝑥₂=𝜀𝑥³₁−𝑘_𝑓𝑥₂−𝑥₁1−𝑥₂−𝑥₁/𝑥₀²,𝐹₂𝑥₁,𝑥₂=𝑘_𝑓𝑥₂−𝑥₁1−𝑥₂−𝑥₁/𝑥₀², (30) with initial conditions 𝑥₁(0)=𝑥₁₀, 𝑥̇₁(0)=0, 𝑥₂(0)=𝑥₂₀, and 𝑥̇₂(0)=0. Here, 𝑚₁ is the mass of the primary system, 𝑚₂ represents the mass of the FENO system, 𝑘₁ and 𝜀 are, respectively, the linear and nonlinear stiffness parameters of the main system, and 𝑘_𝑓 and 𝑥₀ represent the coupling stiffness and the maximum possible extension of the FENO system, respectively. By introducing the transformations 𝑦₁=𝑥₁/𝑥₁₀ and 𝑦₂=𝑥₂/𝑥₂₀ into (29), we get that 𝑦̈₁𝑦̈₂+𝜔²_𝑛10𝑦00𝑦¹₂+𝜑₁𝑥³₁−𝜑₂𝐴𝑦₂−𝑦₁𝑥₁₀/𝑥₂₀1−𝐴²𝑦₂−𝑦₁𝑥₁₀/𝑥₂₀²𝜑₃𝐴𝑦₂−𝑦₁𝑥₁₀/𝑥₂₀1−𝐴²𝑦₂−𝑦₁𝑥₁₀/𝑥₂₀²=0, (31) with initial conditions given as 𝑦₁(0)=1, 𝑦̇₁(0)=0, 𝑦₂(0)=1, and 𝑦̇₂(0)=0, and 𝐴=𝑥₂₀/𝑥₀, 𝜑₁=𝜀𝑥²₁₀/𝑚₁, 𝜑₂=𝑘_𝑓𝑥₀/(𝑚₁𝑥₁₀𝑥₀), and 𝜑₃=𝑘_𝑓𝑥₀/(𝑚₂𝑥₂₀). We next use the Chebyshev polynomials of the first kind [9, 11] to derive the cubic-like representation form of the nonlinear rational restoring force of (31). This procedure yields 𝐹𝑦=𝐴𝑦₂−𝑦₁𝑥₁₀/𝑥₂₀1−𝐴²𝑦₂−𝑦₁𝑥₁₀/𝑥₂₀²=𝐴𝑦1−𝐴²𝑦²=𝑏₁−3𝑏₃𝑦+4𝑏₃𝑦³, (32) where 𝑦=(𝑦₂−𝑦₁𝑥₁₀/𝑥₂₀), and by definition 𝑏_2𝑛+1(𝐴)=2𝜋¹₋₁1−𝑦^2−1/2𝐹𝑦𝑇_2𝑛+1𝑦𝑑𝑦, (33) and 𝑇_2𝑛−1 are the Chebyshev polynomials of the first kind defined as 𝑇_2𝑛+1(𝑥)=cos(2𝑛+1)cos⁻¹(𝑥,[])𝑥∈−1,1,𝑛=0,1,2,…. (34) By using (33), we obtain that 𝑏₁=2𝐴11−𝐴²,𝑏−1₃=2𝐴³41−𝐴²+𝐴²31−1−𝐴²−4. (35) Thus, the FENO nonlinear restoring force 𝐹(𝑦) can be written in equivalent form as 𝐹𝑦≡𝛼𝑦+𝛽𝑦³, (36) in which 1𝛼=𝐴211−𝐴²−6−1𝐴²41−𝐴²+𝐴²31−1−𝐴²,8−4𝛽=𝐴³41−𝐴²+𝐴²31−1−𝐴²−4. (37) Based on (37), it is clear that the equivalent representation form (36) of the restoring force (32) holds if and only if 𝑥₂₀<𝑥₀. Thus, (31) can be written as 𝑦̈₁𝑦̈₂+𝜔²_𝑛10𝑦00𝑦¹₂+𝜑₁𝑥³₁−𝜑₂𝛼𝑦₂−𝑦₁𝑥₁₀𝑥₂₀+𝛽𝑦₂−𝑦₁𝑥₁₀𝑥₂₀³𝜑₃𝛼𝑦₂−𝑦₁𝑥₁₀𝑥₂₀+𝛽𝑦₂−𝑦₁𝑥₁₀𝑥₂₀³=0. (38) We next apply the aforementioned nonlinear transformation approach to (38), to obtain that 𝑦̈₁+𝑎₁𝑦₁+𝑎₂𝑦³₁=0;𝑦̈₂+𝑎₃𝑦₂+𝑎₄𝑦³₂=0, (39) where the 𝑎_𝑖 are given as 𝑎₁=𝜔²_𝑛1+𝜑₂2𝛼𝑥₂₀32𝑥₁₀𝑢₁−45𝑥₂₀𝑢₂+𝛽𝜎₂64𝑥₁₀𝑥₂₀𝑢₁𝑢₂−30𝑥²₁₀𝑢²₁−45𝑥²₂₀𝑢²₂164𝑥²₂₀𝑢₁,𝑎₂=𝜑₁+𝜑₂×70𝛼𝑥³₂₀𝑢₂+𝛽64𝑥³₁₀𝑢³₁−70𝑥²₁₀𝑥₂₀𝑢²₁𝑢₂+35𝑥³₂₀𝑢³₂×164𝑥³₂₀𝑢³₁,𝑎₃=𝜑₃64𝑥³₂₀𝑢₂₂×𝛼−90𝑥₁₀𝑥₂₀𝑢₁₁+64𝑥³₂₀𝑢₂₂+𝛽𝑥₁₀𝑢₁₁×−45𝑥²₁₀𝑢²₁₁+64𝑥₁₀𝑥₂₀𝑢₁₁𝑢₂₂−30𝑥²₂₀𝑢²₂₂,𝑎₄=𝜑₃×6464𝛽+352𝛼𝑥₁₀𝑥₂₀𝑢₁₁+𝛽𝑥³₁₀𝑢³₁₁𝑥³₂₀𝑢³₂₂−70𝛽𝑥₁₀𝑢₁₁𝑥₂₀𝑢₂₂. (40) To assess the degree of accuracy attaining by (39), let us consider the parameter values of 𝑚₁=1, 𝑚₂=1/2, 𝑘₁=1, 𝑘_𝑓=1, 𝑥₀=1, and 𝜀=1, with 𝑥₀(0)=1/2, 𝑥̇₀(0)=0, 𝑥₁(0)=0.99, and 𝑥̇₁(0)=0. Figure 7 shows the amplitude-time response curves of (31) compared to those obtained by substituting (39)₁ into (38)₂. Here, the computed parameter values are 𝛼=−15.4771, 𝛽=37.0369, 𝑎₁=0.3305, 𝑎₂=0.5512, 𝑎₃=−107.549, and 𝑎₄=96.7734, with 𝑢₁=−1.991, 𝑢₂=−1.115, 𝑢₁₁=−0.7, and 𝑢₂₂=−4.6. In this case, the computed RMSE values for 𝑦₁ and 𝑦₂ were, respectively, 0.1965 and 0.2443. Table 5 shows the RMSE computed at different values of 𝑥₂₀. Notice that when 𝑥₂₀ tends to 𝑥₀, the RMSE values tend in general, to decrease, while the global error exhibits a linear increasing behavior, as illustrated in Figure 8.

Table 5 

RMSE values computed by using the equivalent representation form of (31) with 𝑢₁=−1.991, 𝑢₂=−1.115, 𝑢₁₁=−0.7, and 𝑢₂₂=−4.6.

(a)

(b)

Figure 7 

Amplitude-time response curves computed from the numerical integration solutions of (31) and those provided by substituting (39)1 into (38)2. The parameter values used to obtain these plots were 𝑚1=1, 𝑚2=1/2, 𝑘1=1, 𝑘𝑓=1, 𝑥0=1, and 𝜀=1, with 𝑥0(0)=1/2, 𝑥̇0(0)=0, 𝑥1(0)=0.99, 𝑥̇1(0)=0, 𝑢1=−1.991, 𝑢2=−1.115, 𝑢11=−0.7, and 𝑢22=−4.6. Here, the gray solid lines describe the numerical integration solutions of (31), while the blue dashed lines represent the numerical solutions obtained from the derived equivalent nonlinear equations of motion.

Figure 8 

Global error behavior plotted versus time.

2.5. 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, 𝑚₁ and 𝑚₂ which are connected by four elastic springs of stiffness 𝑘₁ and 𝑘₂ having an undeformed length, 𝑙₀. Furthermore, the element with mass 𝑚₁ is moving on a smooth horizontal bearing surface whose motion is restricted by an elastic linear spring of stiffness 𝑘₃. By using Newton second law, it is easy to show that the equations of motion are 10𝑥̈01𝑥̈¹₂+𝑘₃𝑚₁0𝑥00𝑥¹₂+−2𝑘₁𝑚₁0+𝐵02𝑘₁𝑚₂+𝐵₁𝑥₂−𝑥₁𝑥₂−𝑥₁+𝐹₁𝑥₁,𝑥₂𝐹₂𝑥₁,𝑥₂=0, (41) where the irrational restoring forces are given by 𝐹₁𝑥₁,𝑥₂=𝐶𝑥₂−𝑥₁1+𝐴²𝑥₂−𝑥₁²,𝐹₂𝑥₁,𝑥₂=𝐶₁𝑥₂−𝑥₁1+𝐴²𝑥₂−𝑥₁². (42) Here 𝐵=−2𝑘₂/𝑚₁, 𝐵₁=2𝑘₂/𝑚₂, 𝐶=2𝑘₂/𝑚₁, and 𝐶₁=−2𝑘₂/𝑚₂. The system (41) can be cast in different form by setting 𝑥₁=(𝑤₁−𝑤₂)/2 and 𝑥₂=(𝑤₂+𝑤₁)/2, and after some algebraic manipulations, we get that 10𝑤̈01𝑤̈¹₂+𝑘₃2𝑚₁𝑘₁₁𝑘₂₂−𝑘₃2𝑚₁𝑤₁𝑤₂+𝑤₂𝐶+𝐶₁21+𝐴²𝑤²₂−𝑤₂𝐶₁−𝐶21+𝐴²𝑤²₂=0, (43) with initial conditions 𝑤₁(0)=𝑤₁₀, 𝑤̇₁(0)=0, 𝑤₂(0)=𝑤₂₀, and 𝑤̇₂(0)=0. If we set 𝑧₁=𝑤₁/𝑤₁₀ and 𝑧₂=𝑤₂/𝑤₂₀, then (43) becomes 10𝑧̈01𝑧̈¹₂+𝑘₃2𝑚₁𝑘₁₁𝑤₂₀𝑤₁₀𝑘₂₂−𝑘₃𝑤₁₀2𝑚₁𝑤₂₀𝑧₁𝑧₂+𝑧₁𝐶+𝐶₁𝑤₂₀𝑤₁₀1+𝐴²𝑤²₂𝑤²₂₀−𝑧₁𝐶₁−𝐶𝑤₂₀𝑤₁₀1+𝐴²𝑤²₂𝑤²₂₀=0, (44) where 𝑘₁₁=12−𝑘₃/2+2𝑘₁𝑚₁+2𝑘₁𝑚₂+𝐵+𝐵₁,𝑘₂₂=12𝑘₃/2+2𝑘₁𝑚₁+2𝑘₁𝑚₂+𝐵₁−𝐵, (45) with 𝑧₁(0)=1, 𝑧̇₁(0)=0, 𝑧̇₂(0)=1, and 𝑧₂(0)=0. By using the Chebyshev polynomials of the first kind, the cubic-like representation form of the irrational restoring force becomes 𝑧₁1+𝐴²₁𝑤²₂𝑤²₂₀≡𝛼𝑧₁+𝛽𝑧³₁, (46) where 4𝛼=𝜋𝐴²𝐸−𝐴²−𝐾−𝐴²−1𝐴²×8+5𝐴²𝐾−𝐴²−8+𝐴²𝐸−𝐴²,𝛽=163𝜋𝐴⁴×8+5𝐴²𝐾−𝐴²−8+𝐴²𝐸−𝐴². (47) Here 𝐾(−𝐴²) and 𝐸(−𝐴²) are the complete elliptic integral of the first and second kind for the modulus 𝐴, respectively. Thus, (44) can be rewritten as 10𝑧̈01𝑧̈¹₂+𝑘₃2𝑚₁𝑘₁₁𝑤₂₀𝑤₁₀𝑘₂₂−𝑘₃𝑤₁₀2𝑚₁𝑤₂₀𝑧₁𝑧₂+𝐶+𝐶₁𝑤₂₀𝛼𝑧₁+𝛽𝑧³₁𝑤₁₀−𝐶₁−𝐶𝛼𝑧₁+𝛽𝑧³₁𝑤₂₀/𝑤₁₀=0. (48) Application of the proposed nonlinear transformation approach to (48) yields 𝑧̈₁+𝑎₁𝑧₁+𝑎₂𝑧³₁=0;𝑧̈₂+𝑎₃𝑧₂+𝑎₄𝑧³₂=0, (49) where the 𝑎_𝑖 are given as 𝑎₁=164𝑚₁𝑤₁₀𝑢₁32𝑘₃𝑤₁₀𝑢₁+90𝛼𝑚₁𝑤₁₀𝑢₂+90𝑘₁₁𝑚₁𝑤₂₀𝑢₂+45𝛽𝑚₁𝑤₁₀𝑢³₂,𝑎₂=−3564𝑚₁𝑤₁₀𝑢³₁2𝛼𝑤₁₀𝑢₂+2𝑘₁₁𝑤₂₀𝑢₂+𝛽𝑤₁₀𝑢³₂,𝑎₃=−45𝑘₃𝑤₁₀𝑢₁₁−64𝛼𝑚₁𝑤₂₀𝑢₂₂−64𝑘₂₂𝑚₁𝑤₂₀𝑢₂₂64𝑚₁𝑤₂₀𝑢³₂₂,𝑎₄=−−35𝑘₃𝑤₁₀𝑢₁₁−64𝛽𝑚₁𝑤₂₀𝑢³₂₂64𝑚₁𝑤₂₀𝑢³₂₂. (50)

Figure 9 

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

Figure 10 illustrates the numerical predictions obtained by substituting (49)₁ into (48)₂, as well as those computed from the numerical integration solution of (48). In Example 5, we have used the following system parameter values: 𝑚₁=1, 𝑚₂=1/2, 𝑘₁=1, 𝑘₂=1, 𝑘₃=1, and 𝐴=1, with initial conditions given as 𝑤₁(0)=4, 𝑤̇₁(0)=0, 𝑤₂(0)=2, and 𝑤̇₂(0)=0. The predicted parameter values are 𝛼=4.7082, 𝛽=−2.1374, 𝑎₁=−0.5442, 𝑎₂=1.7564, 𝑎₃=8.5590, and 𝑎₄=0.1652, with 𝑢₁=−0.68, 𝑢₂=0.5, 𝑢₁₁=−0.3, and 𝑢₂₂=0.55. Notice from Figure 10 that our numerical integration predictions computed from (49)₁ and (48)₂ follow the amplitude-time response curves computed from (48) with RMSE values of 0.2095 for 𝑧₁ and 0.1318 for 𝑧₂, respectively. The system global error as illustrated in Figure 11 tends 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).

(a)

(b)

Figure 10 

Amplitude-time response curves computed from the numerical integration solutions of (48) and those provided by substituting (49)1 into (48)2. The parameter values used to obtain these plots were 𝑚1=1, 𝑚2=1/2, 𝑘1=1, 𝑘2=1, 𝑘3=1, 𝐴=1, 𝑤1(0)=4, 𝑤̇1(0)=0, 𝑤2(0)=2, 𝑤̇2(0)=0, 𝛼=4.7082, 𝛽=−2.1374, 𝑎1=−0.5442, 𝑎2=1.7564, 𝑎3=8.5590, and 𝑎4=0.1652, with 𝑢1=−0.68, 𝑢2=0.5, 𝑢11=−0.3, and 𝑢22=0.55. Here, the gray solid lines describe the numerical integration solutions of (31), while the blue dashed lines represent the numerical solutions obtained from the derived equivalent nonlinear equations of motion.

Figure 11 

3. Conclusions

We have introduced a nonlinear transformation approach to determine the equivalent representation form of two degree-of-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 original equations of motion. Furthermore, it has been shown that by substituting the decouple equation related to the 𝑥₁(𝑦₁) coordinate system into the original equation of motion described by 𝑥₂(𝑦₂), 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 𝑦₁ is used to find the dynamics response of the 𝑦₀ 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was funded by Tecnológico de Monterrey-Campus Monterrey, through the Research Chair in Nanomaterials for Medical Devices and the Research Chair in Intelligent Machines. Additional support was provided from the European Commission Project, IREBID (FP7-PEOPLE-2009-IRSES-247476), and from Consejo Nacional de Ciencia y Tecnología (Conacyt), México.