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

Authors: Alex Elías-Zúñiga; Daniel Olvera Trejo; Inés Ferrer Real; Oscar Martínez-Romero
Journal: Mathematical Problems in Engineering
 (2014)

---

## 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.


---

## Body

## 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
(1){m100m2}[x¨1x¨2]+{k11k12k21k22}[x1x2]+{f1(x1,x2)f2(x1,x2)}=0
that models the dynamics response of systems that arise in physics and engineering applications [3–6], where m1 and m2 are the masses of the system, k11, k12, k21, and k22 are stiffness parameters, f1(x1,x2) and f2(x1,x2) are nonlinear restoring forces, and x1(0)=x10, x˙1(0)=x˙10, x2(0)=x20, and x˙2(0)=x˙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.


## 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
(2)x¨1+a1x1+a2x13=0;x¨2+a3x2+a4x23=0.
Here we assume that the coefficients ai are determined by the following expressions:
(3)∂F1(a1,a2)∂a1=0,∂F1(a1,a2)∂a2=0,∂F2(a3,a4)∂a3=0,∂F2(a3,a4)∂a4=0,
where
(4)F1(a1,a2)=∫0u2∫0u1(k11x1+k12x2+f1(x1,x2)m1hhhhhhhhhh-a1x1-a2x13k11x1+k12x2+f1(x1,x2)m1)2dx1dx2⟶min⁡,(5)F2(a3,a4)=∫0u4∫0u3(k21x1+k22x2+f2(x1,x2)m2hhhhhhhhhh-a3x2-a4x23k21x1+k22x2+f2(x1,x2)m2)2dx2dx1⟶min⁡
and the values of ui 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]
(6)x1=x10cn(ω1t,κ12),ω1=a1+a2x102,κ12=a2x1022ω12,x2=x20cn(ω2t,κ22),ω2=a3+a4x202,κ22=a3x2022ω22.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 x1 or x2 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]:
(7){u¨1u¨2}+{ωn1200ωn22}[u1u2]+ɛ{φ1u13+φ2u12u2+φ3u1u22+φ4u23φ5u13+φ6u12u2+φ7u1u22+φ8u23}=0,
with initial conditions given as u1(0)=u10, u˙1(0)=u˙10, u2(0)=u20, and u˙2(0)=u˙20. Here the dots denote the time derivative, u1 and u2 are known as the system normal coordinates, ɛ is a nonlinear parameter, and ωn1, ωn2, 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
(8)F1(b1,b2)=∫0ν2∫0ν1(ωn12u1+ɛ(φ1u13+φ2u12u2+φ3u1u22+φ4u23)hhhhhhhhhmm-b1u1-b2u13)2du1du2⟶min⁡,F2(b3,b4)=∫0ν11∫0ν22(ωn22u2+ɛ(φ5u13+φ6u12u2+φ7u1u22+φ8u23)hhhhhhhhhhmm-b3u2-b4u23)2du2du1⟶min⁡,
where ν1, ν2, ν11, and ν22 are fitted constants that satisfy (8), and the parameters bi are determined by solving the equations
(9)∂F1(b1,b2)∂b1=0,∂F1(b1,b2)∂b2=0,∂F2(b3,b4)∂b3=0,∂F2(b3,b4)∂b4=0.This yields(10)b1=--30φ2ν1ν2-64φ3ν1ν22-135φ4ν23-192ν1ωn12192ν1,b2=--192φ1ν13-70φ2ν12ν2+105φ4ν23192ν13,b3=--135φ5ν113-64φ6ν112ν22-30φ7ν11ν222-192ν22ωn22192ν22,b4=--105φ5ν113-70φ7ν11ν222-192φ8ν223192ν223.
Then, (7) can be rewritten in equivalent form as
(11)u¨1+b1u1+b2u13=0,u¨2+b3u2+b4u23=0,
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]
(12){m100m2}[x¨1x¨2]+{k1+k2-k2-k22k2}[x1x2]+{ɛx130}=0,
with initial conditions x1(0)=x10, x˙1(0)=0, x2(0)=x20, and x˙2(0)=0. If we introduce the transformations y1=x1/x10 and y2=x2/x20 so that y1(0)=1, y˙1(0)=0, y2(0)=1, y˙2(0)=0, thus (12) can be written in normalized form as
(13){1001}[y¨1y¨2]+{(k1+k2)m1-k2x20(m1x10)-k2x10(m2x20)2k2m2}[y1y2]+{ɛx102y13m10}=0.
If we use our proposed nonlinear transformation approach then the equivalent representation form of (13) is given as
(14)y¨1+a1y1+a2y13=0;y¨2+a3y2+a4y23=0,
where the ai, determined from (4) and (5), are given as
(15)a1=--32k1u1x10-32k2u1x10+45k2u2x2032m1u1x10,a2=--32ɛu13x103-35k2u2x2032m1u13x10,a3=-45k2u11x10-32k2u22x2032m2u22x20,a4=35k2u11x1032m2u223x20,
and u1, u2, u11, and u22 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)[x1x2]={R1R2R1f1R2f2}[u1u2],
which yields the canonical form (7). Here, the system parameters of (7) are described by the relations
(17)ωn12=12m1[(k1+3k2)-k12-2k1k2+5k22],ωn22=12m1[(k1+3k2)+k12-2k1k2+5k22],φ1≡-ɛR14x102,φ2≡-3ɛR13R2x10x20,φ3≡-3ɛR12R22x202,φ4≡-ɛR1R23x203x10,φ5≡-ɛR13R2x103x20,φ6≡-3ɛR12R22x102,φ7≡-3ɛR1R23x10x20,φ8≡-ɛR24x202,f1≡k1+k2-m1ωn12k2,f2≡k1+k2-m1ωn22k2,R12≡1m1(1+f12),R22≡1m1(1+f22).
Thus, the approximate solutions of (13), by considering their canonical representation form, are given by expressions similar to those of (6).Tables1 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 x1(0)=2, x˙1(0)=0, x2(0)=-1, and x˙2(0)=0, the parameter values of m1=m2=1, k1=3, and k2=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)1 into (13)2. The sixth and the seventh columns of Table 1 contain the RMSE values computed by substituting (14)2 into (13)1. 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)1 is used to derive the approximate solution of x2. 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)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 x1, and 0.0233, for x2. 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 u1=1.55, u2=1.5, u3=1.25, and u4=1.6. Here b11=(k1+k2)/m1, b12=b21=-k2x20/(m1x10), b22=2k2/m2, and d1=εx102/m1.
ε
y
¨
1
+
a
1
y
1
+
a
2
y
1
3
=
0
{
∑
}, y¨2+a3y2+a4y23=0.
y
¨
1
+
a
1
y
1
+
a
2
y
1
3
=
0
{
∑
}, y¨2+b21y1+b22y2=0.
y
¨
2
+
a
3
y
2
+
a
4
y
2
3
=
0
{
∑
}, y¨1+b11y1+b12y2+d1y13=0.
RMSEy1(x)
RMSEy2(x)
RMSEy1(x)
RMSEy2(x)
RMSEy1(x)
RMSEy2(x)
0.01
0.2992
0.9138
0.2992
0.3062
1.1924
0.9138
0.25
0.3248
0.8863
0.3248
0.4004
0.1797
0.8863
0.75
0.1511
0.8222
0.1511
0.2873
0.2167
0.8222
1.50
0.0398
0.8373
0.0398
0.1989
0.1902
0.8373
5
0.1488
0.9068
0.1488
0.0818
0.0930
0.9068
20
0.1179
0.9490
0.1179
0.0233
0.0082
0.9490
100
0.0697
0.9635
0.0697
0.0047
0.0011
0.9635Table 2
RMSE values computed by considering the canonical representation form of (13) with ν1=-9, ν2=-0.215, ν11=-2.325, and ν22=-8.5.
ε
u
¨
1
+
b
1
u
1
+
b
2
u
1
3
=
0
{
∑
},u¨2+b3u2+b4u23=0.
u
¨
1
+
b
1
u
1
+
b
2
u
1
3
=
0
{
∑
},u¨2+ωn22u2+ε{φ5u13+φ6u12u2+φ7u1u22+φ8u23}=0.
u
¨
2
+
b
3
u
2
+
b
4
u
2
3
=
0
{
∑
}, u¨1+ωn12u1+ε{φ1u13+φ2u12u2+φ3u1u22+φ4u23}=0.
RMSEu1(x)
RMSEu2(x)
RMSEu1(x)
RMSEu2(x)
RMSEu1(x)
RMSEu2(x)
0.01
0.06010
0.0277
0.0044
0.0109
0.0606
0.0248
0.25
1.5351
0.6069
0.3611
0.2661
1.5824
0.4956
0.75
2.0910
0.8049
1.3988
0.6512
2.2351
0.7348
1.50
2.1222
0.8349
2.1151
0.8220
2.5577
1.1983
5
1.9552
0.8733
1.8608
0.8739
3.1331
8.8874
20
1.9116
0.9044
1.9889
0.9090
3.1549
8.7513
100
1.9495
0.9136
1.9704
0.9191
3.4887
8.1949Amplitude-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 m1=m2=1, k1=3, k2=1, and ɛ=20, with x1(0)=2, x˙1(0)=0, x2(0)=-1, and x˙2(0)=0, and u1=1.55, u2=1.5, u3=1.25, and u4=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.
(a)
(b)Table2 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 a2 and a4 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:
(18)GE=(y1(t)-p1(t))+(y2(t)-p2(t)),
where y1 and y2 represent the values obtained by numerically integrating the equivalent equations of motion and p1 and p2 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].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(19){100β}[σ¨z¨]+{1-α22-α22α22}[σz]-{ɛb1σ30}=0,
with initial conditions given as σ(0)=σ10, σ˙(0)=0, z(0)=z10, and z˙(0)=0. Here σ and z denote, respectively, the simple shear deformation of the load and the motion of the linear absorber system, ɛ describes the nonlinear material response effects, r2 is a tuning system parameter, and β, α2, and b1 are defined as
(20)β=mM,b1=α22-1,α22=βr22βr22+1.
The details of the derivation of the equations of motion given by (19) are provided in [2]. Next, we set y1=σ/σ10 and y2=z/z10 and rewrite (12) in its normalized form
(21){1001}[y¨1y¨2]+{1-α22z10σ10-α22σ10(βz10)α22β}[y1y2]-{ɛσ102y130}=0
with y1(0)=1, y˙1(0)=0, y2(0)=1, and y˙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:
(22)ωn12=12[(α22β+1)-(α22β+1)2+4α22b1β],ωn22=12[(α22β+1)+(α22β+1)2+4α22b1β],φ1≡-b1ɛR14σ102,φ2≡-3b1ɛR13R2σ10z20,φ3≡-3b1ɛR12R22z202,φ4≡-b1ɛR1R23z203σ10,φ5≡-b1ɛR13R2σ103z20,φ6≡-3b1ɛR12R22σ102,φ7≡-3b1ɛR1R23σ10z20,φ8≡-b1ɛR24z202,f1≡1-ωn12α22,f2≡1-ωn22α22,R12≡1(1+f12β),R22≡1(1+f22β).
Since (21)2 is only coupled with y1 through the linear term -α22σ10/(βz10)y1, we slightly modify our procedure and assume that the equivalent representation form of (21) is given as
(23)y¨1+a1y1+a2y13=0;y¨2+a3y2=0,
where the parameters ai, computed from (4), are given as
(24)a1=--32u1σ10+45α22u2z1032u1σ10,a2=--32b1ɛu13σ103-35α22u2z1032u13σ10,a3=-α22u11σ10+α22u22z10βu22z10,
and u1, u2, u11, and u22 are the corresponding fitting parameters. The assumption that (21)2 can be transformed into a linear equation of the form (23)2 is based on the fact that in (21)2, there is not y2(t) 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 σ10(0)=1, σ˙(0)=0, z10(0)=-1, and z˙(0)=0, with β=1/5 and r2=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 y2.Table 3
RMSE values computed from the equivalent representation forms of  (21) with u1=-0.4, u2=0.1, u3=-0.01, and u4=-0.315. Here b11=1, b12=b21=-α22σ10/(βz10), b22=α22/β, and d1=-εσ102.
ε
y
¨
1
+
a
1
y
1
+
a
2
y
1
3
=
0
{
∑
}, y¨2+a3y2=0.
y
¨
1
+
a
1
y
1
+
a
2
y
1
3
=
0
{
∑
}, y¨2+b21y1+b22y2=0.
y
¨
2
+
a
3
y
2
=
0
{
∑
}, y¨1+b11y1+b12y2+d1y13=0.
RMSEy1(x)
RMSEy2(x)
RMSEy1(x)
RMSEy2(x)
RMSEy1(x)
RMSEy2(x)
0.01
0.9448
1.3409
0.9448
2.3366
1.7724
1.3409
0.25
0.5880
1.2351
0.5880
1.1768
0.7684
1.2351
0.75
0.1609
1.0394
0.1609
0.1774
0.3309
1.0394
1.5
0.2402
0.7685
0.2402
0.2079
0.3157
0.7685
5
0.3637
0.4461
0.3637
0.1376
0.0910
0.4461
20
0.2489
0.3267
0.2489
0.0414
0.0125
0.3267
100
0.1226
0.2995
0.1226
0.0083
0.0011
0.2995Table 4
RMSE values computed by considering the canonical representation form of  (21) with ν1=0.85, ν2=0.45, ν11=-0.01, and ν22=0.1.
ε
u
¨
1
+
b
1
u
1
+
b
2
u
1
3
=
0
{
∑
}, u¨2+b3u2+b4u23=0.
u
¨
1
+
b
1
u
1
+
b
2
u
1
3
=
0
{
∑
}, u¨2+ωn22u2+ε{φ5u13+φ6u12u2+φ7u1u22+φ8u23}=0.
u
¨
2
+
b
3
u
2
+
b
4
u
2
3
=
0
{
∑
}, u¨1+ωn12u1+ε{φ1u13+φ2u12u2+φ3u1u22+φ4u23}=0.
RMSEu1(x)
RMSEu2(x)
RMSEu1(x)
RMSEu2(x)
RMSEu1(x)
RMSEu2(x)
0.01
0.0077
0.0168
0.0066
0.0174
0.0028
0.0050
0.25
0.1953
0.3769
0.1325
0.3596
0.1516
0.2393
0.75
0.8844
1.2761
0.7365
1.1609
0.9850
0.9016
1
1.0156
1.3704
0.9953
1.4230
1.4544
3.8230
1.5
0.8727
1.1634
0.9922
1.3683
1.5028
4.3001
5
0.9301
1.2149
0.9678
1.4201
0.8801
4.0894
20
0.9010
1.2577
0.9729
1.5285
1.2238
3.8974
100
0.9282
1.270
1.0162
1.5621
1.2479
3.8438Amplitude-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 r2=1, ɛ=50, with σ10(0)=1, σ˙(0)=0, z10(0)=-1, z˙(0)=0, and u1=-0.4, u2=0.1, u3=-0.01, and u4=-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 y1, and 0.0102, for y2.
(a)
(b)Also notice from Table4 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 obtained when the equivalent form (23)1 is substituted into (21)2.


### 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
(25){m00m}[x¨1x¨2]+{k1+k2-k2-k2k1+k2}[x1x2]+{k3(x1-x2)3+ɛk4x13k3(x2-x1)3+ɛk4x23}=0,
with initial conditions x1(0)=x10, x˙1(0)=0, x2(0)=x20, and x˙2(0)=0. If we introduce the transformations x1=(y1+y0)/2 and x2=(y0-y1)/2, thus (25) can be written as
(26){y¨0y¨1}+{ωn1200ωn22}[y0y1]+{ɛk4(4m)y03+ɛk4(4m)y0y12ɛk4(4m)y02y1+(β+ɛk4(4m))y13}=0,
with initial conditions given as y0(0)=y00, y˙0(0)=y˙00, y1(0)=y10, and y˙1(0)=y˙10. Notice that (26) has the same form of (7) with system parameters given as
(27)ωn12=k1m,ωn22=k1+2k2m,φ1=φ3=ɛk4(4m),φ2=φ4=0,φ6=ɛk4(4m),φ8=β+ɛk44m,φ5=φ7=0,β=2k3m.
Therefore and in accordance with our transformation approach, (26) are transformed into (11) in which
(28)b1=φ3ν22+3ωn123;b2=φ1;b3=φ6ν112+3ωn223;b4=φ8.
To assess the accuracy of the equivalent equations of motion when compared to the original ones (25), let us consider the parameter values of k1=k2=k3=k4=1, m=1, and ɛ=-2.9 and initial conditions y0(0)=0.5, y˙0(0)=0, y1(0)=1.5, and y˙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 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)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 y1 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 y0 and y1, respectively.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 m=1, k1=k2=k3=k4=1, and ɛ=-2.9, with y0(0)=0.5, y˙0(0)=0, y1(0)=1, y˙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 y0, and 0.0358, for y1.
(a)
(b)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 m=1, k1=k2=k3=k4=1, and ɛ=-2.9, with y0(0)=0.5, y˙0(0)=0, y1(0)=1, y˙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 y0, and 0.1792, for y1.
(a)
(b)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 t≥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 t≥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.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 m=1, k1=k2=k3=k4=1, and ɛ=-3.1, with y0(0)=0.5, y˙0(0)=0, y1(0)=1, y˙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.
(a)
(b)


### 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
(29){m100m2}[x¨1x¨2]+{k1000}[x1x2]+{F1(x1,x2)F2(x1,x2)}=0,
and the nonlinear restoring forces are given by
(30)F1(x1,x2)=ɛx13-kf(x2-x1)1-((x2-x1)/x0)2,F2(x1,x2)=kf(x2-x1)1-((x2-x1)/x0)2,
with initial conditions x1(0)=x10, x˙1(0)=0, x2(0)=x20, and x˙2(0)=0. Here, m1 is the mass of the primary system, m2 represents the mass of the FENO system, k1 and ɛ are, respectively, the linear and nonlinear stiffness parameters of the main system, and kf and x0 represent the coupling stiffness and the maximum possible extension of the FENO system, respectively. By introducing the transformations y1=x1/x10 and y2=x2/x20 into (29), we get that
(31){y¨1y¨2}+{ωn12000}[y1y2]+{φ1x13-φ2A(y2-y1x10/x20)1-A2(y2-y1x10/x20)2φ3A(y2-y1x10/x20)1-A2(y2-y1x10/x20)2}=0,
with initial conditions given as y1(0)=1, y˙1(0)=0, y2(0)=1, and y˙2(0)=0, and A=x20/x0, φ1=ɛx102/m1, φ2=kfx0/(m1x10x0), and φ3=kfx0/(m2x20). 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
(32)F(y)=A(y2-y1x10/x20)1-A2(y2-y1x10/x20)2=Ay1-A2y2=(b1-3b3)y+4b3y3,
where y=(y2-y1x10/x20), and by definition
(33)b2n+1(A)=2π∫-11(1-y2)-1/2F(y)T2n+1(y)dy,
and T2n-1 are the Chebyshev polynomials of the first kind defined as
(34)T2n+1(x)=cos⁡[(2n+1)cos⁡-1(x)],hx∈[-1,1],n=0,1,2,….
By using (33), we obtain that
(35)b1=2A(11-A2-1),b3=2A3(41-A2+A2(1-31-A2)-4).
Thus, the FENO nonlinear restoring force F(y) can be written in equivalent form as
(36)F(y)≡αy+βy3,
in which
(37)α=1A{2(11-A2-1)hhhhh-6A2(41-A2+A2(1-31-A2)-4)},β=8A3{41-A2+A2(1-31-A2)-4}.
Based on (37), it is clear that the equivalent representation form (36) of the restoring force (32) holds if and only if x20<x0. Thus, (31) can be written as
(38){y¨1y¨2}+{ωn12000}[y1y2]+{φ1x13-φ2{α(y2-y1x10x20)+β(y2-y1x10x20)3}φ3{α(y2-y1x10x20)+β(y2-y1x10x20)3}}=0.
We next apply the aforementioned nonlinear transformation approach to (38), to obtain that
(39)y¨1+a1y1+a2y13=0;y¨2+a3y2+a4y23=0,
where the ai are given as
(40)a1=ωn12+φ2(2αx20(32x10u1-45x20u2)+βσ2{64x10x20u1u2-30x102u12-45x202u22})164x202u1,a2=φ1+φ2×(70αx203u2hhh+β(64x103u13-70x102x20u12u2+35x203u23))×164x203u13,a3=φ364x203u22×(×(-45x102u112+64x10x20u11u22-30x202u222))α(-90x10x20u11+64x203u22)+βx10u11hhhh×(-45x102u112+64x10x20u11u22-30x202u222)),a4=φ364×(64β+35(2αx10x20u11+βx103u113)x203u223-70βx10u11x20u2235(2αx10x20u11+βx103u113)x203u223).
To assess the degree of accuracy attaining by (39), let us consider the parameter values of m1=1, m2=1/2, k1=1, kf=1, x0=1, and ɛ=1, with x0(0)=1/2, x˙0(0)=0, x1(0)=0.99, and x˙1(0)=0. Figure 7 shows the amplitude-time response curves of (31) compared to those obtained by substituting (39)1 into (38)2. Here, the computed parameter values are α=-15.4771, β=37.0369, a1=0.3305, a2=0.5512, a3=-107.549, and a4=96.7734, with u1=-1.991, u2=-1.115, u11=-0.7, and u22=-4.6. In this case, the computed RMSE values for y1 and y2 were, respectively, 0.1965 and 0.2443. Table 5 shows the RMSE computed at different values of x20. Notice that when x20 tends to x0, 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 u1=-1.991, u2=-1.115, u11=-0.7, and u22=-4.6.
x
20
RMSEy1(x)
RMSEy2(x)
0.75
0.9981
1.1200
0.80
1.0245
1.0467
0.85
1.0330
0.9560
0.90
0.5851
0.5029
0.95
1.0122
0.8019
0.99
0.1965
0.2443Amplitude-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 m1=1, m2=1/2, k1=1, kf=1, x0=1, and ɛ=1, with x0(0)=1/2, x˙0(0)=0, x1(0)=0.99, x˙1(0)=0, u1=-1.991, u2=-1.115, u11=-0.7, and u22=-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.
(a)
(b)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 Figure9 having two masses, m1 and m2 which are connected by four elastic springs of stiffness k1 and k2 having an undeformed length, l0. Furthermore, the element with mass m1 is moving on a smooth horizontal bearing surface whose motion is restricted by an elastic linear spring of stiffness k3. By using Newton second law, it is easy to show that the equations of motion are
(41){1001}[x¨1x¨2]+{k3m1000}[x1x2]+{-2k1m1+B002k1m2+B1}[(x2-x1)(x2-x1)]+{F1(x1,x2)F2(x1,x2)}=0,
where the irrational restoring forces are given by
(42)F1(x1,x2)=C(x2-x1)1+A2(x2-x1)2,F2(x1,x2)=C1(x2-x1)1+A2(x2-x1)2.
Here B=-2k2/m1, B1=2k2/m2, C=2k2/m1, and C1=-2k2/m2. The system (41) can be cast in different form by setting x1=(w1-w2)/2 and x2=(w2+w1)/2, and after some algebraic manipulations, we get that
(43){1001}[w¨1w¨2]+{k3(2m1)k11k22-k3(2m1)}[w1w2]+{w2(C+C1)21+A2w22-w2(C1-C)21+A2w22}=0,
with initial conditions w1(0)=w10, w˙1(0)=0, w2(0)=w20, and w˙2(0)=0. If we set z1=w1/w10 and z2=w2/w20, then (43) becomes
(44){1001}[z¨1z¨2]+{k3(2m1)k11w20w10k22-k3w10(2m1w20)}[z1z2]+{z1(C+C1)w20w101+A2w22w202-z1(C1-C)w20w101+A2w22w202}=0,
where
(45)k11=12(-k3/2+2k1m1+2k1m2+B+B1),k22=12(k3/2+2k1m1+2k1m2+B1-B),
with z1(0)=1, z˙1(0)=0, z˙2(0)=1, and z2(0)=0. By using the Chebyshev polynomials of the first kind, the cubic-like representation form of the irrational restoring force becomes
(46)z11+A12w22w202≡αz1+βz13,
where
(47)α=4πA2{E(-A2)-K(-A2)-1A2hhhhhhh×((8+5A2)K(-A2)-(8+A2)E(-A2))1A2},β=163πA4×{(8+5A2)K(-A2)-(8+A2)E(-A2)}.
Here K(-A2) and E(-A2) are the complete elliptic integral of the first and second kind for the modulus A, respectively. Thus, (44) can be rewritten as
(48){1001}[z¨1z¨2]+{k3(2m1)k11w20w10k22-k3w10(2m1w20)}[z1z2]+{(C+C1)w20{αz1+βz13}w10-(C1-C){αz1+βz13}w20/w10}=0.
Application of the proposed nonlinear transformation approach to (48) yields
(49)z¨1+a1z1+a2z13=0;z¨2+a3z2+a4z23=0,
where the ai are given as
(50)a1=164m1w10u1(u2332k3w10u1+90αm1w10u2hhhhhhhhhhhh+90k11m1w20u2+45βm1w10u23),a2=-3564m1w10u13(2αw10u2+2k11w20u2+βw10u23),a3=-45k3w10u11-64αm1w20u22-64k22m1w20u2264m1w20u223,a4=--35k3w10u11-64βm1w20u22364m1w20u223.Figure 9
Schematic of a two-degree-of-freedom system with irrational restoring force.Figure10 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: m1=1, m2=1/2, k1=1, k2=1, k3=1, and A=1, with initial conditions given as w1(0)=4, w˙1(0)=0, w2(0)=2, and w˙2(0)=0. The predicted parameter values are α=4.7082, β=-2.1374, a1=-0.5442, a2=1.7564, a3=8.5590, and a4=0.1652, with u1=-0.68, u2=0.5, u11=-0.3, and u22=0.55. Notice from Figure 10 that our numerical integration predictions computed from (49)1 and (48)2 follow the amplitude-time response curves computed from (48) with RMSE values of 0.2095 for z1 and 0.1318 for z2, 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 ui and uii in order to satisfy (8).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 m1=1, m2=1/2, k1=1, k2=1, k3=1, A=1, w1(0)=4, w˙1(0)=0, w2(0)=2, w˙2(0)=0, α=4.7082, β=-2.1374, a1=-0.5442, a2=1.7564, a3=8.5590, and a4=0.1652, with u1=-0.68, u2=0.5, u11=-0.3, and u22=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.
(a)
(b)Figure 11


## 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]
(12){m100m2}[x¨1x¨2]+{k1+k2-k2-k22k2}[x1x2]+{ɛx130}=0,
with initial conditions x1(0)=x10, x˙1(0)=0, x2(0)=x20, and x˙2(0)=0. If we introduce the transformations y1=x1/x10 and y2=x2/x20 so that y1(0)=1, y˙1(0)=0, y2(0)=1, y˙2(0)=0, thus (12) can be written in normalized form as
(13){1001}[y¨1y¨2]+{(k1+k2)m1-k2x20(m1x10)-k2x10(m2x20)2k2m2}[y1y2]+{ɛx102y13m10}=0.
If we use our proposed nonlinear transformation approach then the equivalent representation form of (13) is given as
(14)y¨1+a1y1+a2y13=0;y¨2+a3y2+a4y23=0,
where the ai, determined from (4) and (5), are given as
(15)a1=--32k1u1x10-32k2u1x10+45k2u2x2032m1u1x10,a2=--32ɛu13x103-35k2u2x2032m1u13x10,a3=-45k2u11x10-32k2u22x2032m2u22x20,a4=35k2u11x1032m2u223x20,
and u1, u2, u11, and u22 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)[x1x2]={R1R2R1f1R2f2}[u1u2],
which yields the canonical form (7). Here, the system parameters of (7) are described by the relations
(17)ωn12=12m1[(k1+3k2)-k12-2k1k2+5k22],ωn22=12m1[(k1+3k2)+k12-2k1k2+5k22],φ1≡-ɛR14x102,φ2≡-3ɛR13R2x10x20,φ3≡-3ɛR12R22x202,φ4≡-ɛR1R23x203x10,φ5≡-ɛR13R2x103x20,φ6≡-3ɛR12R22x102,φ7≡-3ɛR1R23x10x20,φ8≡-ɛR24x202,f1≡k1+k2-m1ωn12k2,f2≡k1+k2-m1ωn22k2,R12≡1m1(1+f12),R22≡1m1(1+f22).
Thus, the approximate solutions of (13), by considering their canonical representation form, are given by expressions similar to those of (6).Tables1 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 x1(0)=2, x˙1(0)=0, x2(0)=-1, and x˙2(0)=0, the parameter values of m1=m2=1, k1=3, and k2=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)1 into (13)2. The sixth and the seventh columns of Table 1 contain the RMSE values computed by substituting (14)2 into (13)1. 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)1 is used to derive the approximate solution of x2. 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)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 x1, and 0.0233, for x2. 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 u1=1.55, u2=1.5, u3=1.25, and u4=1.6. Here b11=(k1+k2)/m1, b12=b21=-k2x20/(m1x10), b22=2k2/m2, and d1=εx102/m1.
ε
y
¨
1
+
a
1
y
1
+
a
2
y
1
3
=
0
{
∑
}, y¨2+a3y2+a4y23=0.
y
¨
1
+
a
1
y
1
+
a
2
y
1
3
=
0
{
∑
}, y¨2+b21y1+b22y2=0.
y
¨
2
+
a
3
y
2
+
a
4
y
2
3
=
0
{
∑
}, y¨1+b11y1+b12y2+d1y13=0.
RMSEy1(x)
RMSEy2(x)
RMSEy1(x)
RMSEy2(x)
RMSEy1(x)
RMSEy2(x)
0.01
0.2992
0.9138
0.2992
0.3062
1.1924
0.9138
0.25
0.3248
0.8863
0.3248
0.4004
0.1797
0.8863
0.75
0.1511
0.8222
0.1511
0.2873
0.2167
0.8222
1.50
0.0398
0.8373
0.0398
0.1989
0.1902
0.8373
5
0.1488
0.9068
0.1488
0.0818
0.0930
0.9068
20
0.1179
0.9490
0.1179
0.0233
0.0082
0.9490
100
0.0697
0.9635
0.0697
0.0047
0.0011
0.9635Table 2
RMSE values computed by considering the canonical representation form of (13) with ν1=-9, ν2=-0.215, ν11=-2.325, and ν22=-8.5.
ε
u
¨
1
+
b
1
u
1
+
b
2
u
1
3
=
0
{
∑
},u¨2+b3u2+b4u23=0.
u
¨
1
+
b
1
u
1
+
b
2
u
1
3
=
0
{
∑
},u¨2+ωn22u2+ε{φ5u13+φ6u12u2+φ7u1u22+φ8u23}=0.
u
¨
2
+
b
3
u
2
+
b
4
u
2
3
=
0
{
∑
}, u¨1+ωn12u1+ε{φ1u13+φ2u12u2+φ3u1u22+φ4u23}=0.
RMSEu1(x)
RMSEu2(x)
RMSEu1(x)
RMSEu2(x)
RMSEu1(x)
RMSEu2(x)
0.01
0.06010
0.0277
0.0044
0.0109
0.0606
0.0248
0.25
1.5351
0.6069
0.3611
0.2661
1.5824
0.4956
0.75
2.0910
0.8049
1.3988
0.6512
2.2351
0.7348
1.50
2.1222
0.8349
2.1151
0.8220
2.5577
1.1983
5
1.9552
0.8733
1.8608
0.8739
3.1331
8.8874
20
1.9116
0.9044
1.9889
0.9090
3.1549
8.7513
100
1.9495
0.9136
1.9704
0.9191
3.4887
8.1949Amplitude-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 m1=m2=1, k1=3, k2=1, and ɛ=20, with x1(0)=2, x˙1(0)=0, x2(0)=-1, and x˙2(0)=0, and u1=1.55, u2=1.5, u3=1.25, and u4=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.
(a)
(b)Table2 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 a2 and a4 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:
(18)GE=(y1(t)-p1(t))+(y2(t)-p2(t)),
where y1 and y2 represent the values obtained by numerically integrating the equivalent equations of motion and p1 and p2 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].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(19){100β}[σ¨z¨]+{1-α22-α22α22}[σz]-{ɛb1σ30}=0,
with initial conditions given as σ(0)=σ10, σ˙(0)=0, z(0)=z10, and z˙(0)=0. Here σ and z denote, respectively, the simple shear deformation of the load and the motion of the linear absorber system, ɛ describes the nonlinear material response effects, r2 is a tuning system parameter, and β, α2, and b1 are defined as
(20)β=mM,b1=α22-1,α22=βr22βr22+1.
The details of the derivation of the equations of motion given by (19) are provided in [2]. Next, we set y1=σ/σ10 and y2=z/z10 and rewrite (12) in its normalized form
(21){1001}[y¨1y¨2]+{1-α22z10σ10-α22σ10(βz10)α22β}[y1y2]-{ɛσ102y130}=0
with y1(0)=1, y˙1(0)=0, y2(0)=1, and y˙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:
(22)ωn12=12[(α22β+1)-(α22β+1)2+4α22b1β],ωn22=12[(α22β+1)+(α22β+1)2+4α22b1β],φ1≡-b1ɛR14σ102,φ2≡-3b1ɛR13R2σ10z20,φ3≡-3b1ɛR12R22z202,φ4≡-b1ɛR1R23z203σ10,φ5≡-b1ɛR13R2σ103z20,φ6≡-3b1ɛR12R22σ102,φ7≡-3b1ɛR1R23σ10z20,φ8≡-b1ɛR24z202,f1≡1-ωn12α22,f2≡1-ωn22α22,R12≡1(1+f12β),R22≡1(1+f22β).
Since (21)2 is only coupled with y1 through the linear term -α22σ10/(βz10)y1, we slightly modify our procedure and assume that the equivalent representation form of (21) is given as
(23)y¨1+a1y1+a2y13=0;y¨2+a3y2=0,
where the parameters ai, computed from (4), are given as
(24)a1=--32u1σ10+45α22u2z1032u1σ10,a2=--32b1ɛu13σ103-35α22u2z1032u13σ10,a3=-α22u11σ10+α22u22z10βu22z10,
and u1, u2, u11, and u22 are the corresponding fitting parameters. The assumption that (21)2 can be transformed into a linear equation of the form (23)2 is based on the fact that in (21)2, there is not y2(t) 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 σ10(0)=1, σ˙(0)=0, z10(0)=-1, and z˙(0)=0, with β=1/5 and r2=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 y2.Table 3
RMSE values computed from the equivalent representation forms of  (21) with u1=-0.4, u2=0.1, u3=-0.01, and u4=-0.315. Here b11=1, b12=b21=-α22σ10/(βz10), b22=α22/β, and d1=-εσ102.
ε
y
¨
1
+
a
1
y
1
+
a
2
y
1
3
=
0
{
∑
}, y¨2+a3y2=0.
y
¨
1
+
a
1
y
1
+
a
2
y
1
3
=
0
{
∑
}, y¨2+b21y1+b22y2=0.
y
¨
2
+
a
3
y
2
=
0
{
∑
}, y¨1+b11y1+b12y2+d1y13=0.
RMSEy1(x)
RMSEy2(x)
RMSEy1(x)
RMSEy2(x)
RMSEy1(x)
RMSEy2(x)
0.01
0.9448
1.3409
0.9448
2.3366
1.7724
1.3409
0.25
0.5880
1.2351
0.5880
1.1768
0.7684
1.2351
0.75
0.1609
1.0394
0.1609
0.1774
0.3309
1.0394
1.5
0.2402
0.7685
0.2402
0.2079
0.3157
0.7685
5
0.3637
0.4461
0.3637
0.1376
0.0910
0.4461
20
0.2489
0.3267
0.2489
0.0414
0.0125
0.3267
100
0.1226
0.2995
0.1226
0.0083
0.0011
0.2995Table 4
RMSE values computed by considering the canonical representation form of  (21) with ν1=0.85, ν2=0.45, ν11=-0.01, and ν22=0.1.
ε
u
¨
1
+
b
1
u
1
+
b
2
u
1
3
=
0
{
∑
}, u¨2+b3u2+b4u23=0.
u
¨
1
+
b
1
u
1
+
b
2
u
1
3
=
0
{
∑
}, u¨2+ωn22u2+ε{φ5u13+φ6u12u2+φ7u1u22+φ8u23}=0.
u
¨
2
+
b
3
u
2
+
b
4
u
2
3
=
0
{
∑
}, u¨1+ωn12u1+ε{φ1u13+φ2u12u2+φ3u1u22+φ4u23}=0.
RMSEu1(x)
RMSEu2(x)
RMSEu1(x)
RMSEu2(x)
RMSEu1(x)
RMSEu2(x)
0.01
0.0077
0.0168
0.0066
0.0174
0.0028
0.0050
0.25
0.1953
0.3769
0.1325
0.3596
0.1516
0.2393
0.75
0.8844
1.2761
0.7365
1.1609
0.9850
0.9016
1
1.0156
1.3704
0.9953
1.4230
1.4544
3.8230
1.5
0.8727
1.1634
0.9922
1.3683
1.5028
4.3001
5
0.9301
1.2149
0.9678
1.4201
0.8801
4.0894
20
0.9010
1.2577
0.9729
1.5285
1.2238
3.8974
100
0.9282
1.270
1.0162
1.5621
1.2479
3.8438Amplitude-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 r2=1, ɛ=50, with σ10(0)=1, σ˙(0)=0, z10(0)=-1, z˙(0)=0, and u1=-0.4, u2=0.1, u3=-0.01, and u4=-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 y1, and 0.0102, for y2.
(a)
(b)Also notice from Table4 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 obtained when the equivalent form (23)1 is substituted into (21)2.


## 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
(25){m00m}[x¨1x¨2]+{k1+k2-k2-k2k1+k2}[x1x2]+{k3(x1-x2)3+ɛk4x13k3(x2-x1)3+ɛk4x23}=0,
with initial conditions x1(0)=x10, x˙1(0)=0, x2(0)=x20, and x˙2(0)=0. If we introduce the transformations x1=(y1+y0)/2 and x2=(y0-y1)/2, thus (25) can be written as
(26){y¨0y¨1}+{ωn1200ωn22}[y0y1]+{ɛk4(4m)y03+ɛk4(4m)y0y12ɛk4(4m)y02y1+(β+ɛk4(4m))y13}=0,
with initial conditions given as y0(0)=y00, y˙0(0)=y˙00, y1(0)=y10, and y˙1(0)=y˙10. Notice that (26) has the same form of (7) with system parameters given as
(27)ωn12=k1m,ωn22=k1+2k2m,φ1=φ3=ɛk4(4m),φ2=φ4=0,φ6=ɛk4(4m),φ8=β+ɛk44m,φ5=φ7=0,β=2k3m.
Therefore and in accordance with our transformation approach, (26) are transformed into (11) in which
(28)b1=φ3ν22+3ωn123;b2=φ1;b3=φ6ν112+3ωn223;b4=φ8.
To assess the accuracy of the equivalent equations of motion when compared to the original ones (25), let us consider the parameter values of k1=k2=k3=k4=1, m=1, and ɛ=-2.9 and initial conditions y0(0)=0.5, y˙0(0)=0, y1(0)=1.5, and y˙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 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)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 y1 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 y0 and y1, respectively.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 m=1, k1=k2=k3=k4=1, and ɛ=-2.9, with y0(0)=0.5, y˙0(0)=0, y1(0)=1, y˙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 y0, and 0.0358, for y1.
(a)
(b)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 m=1, k1=k2=k3=k4=1, and ɛ=-2.9, with y0(0)=0.5, y˙0(0)=0, y1(0)=1, y˙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 y0, and 0.1792, for y1.
(a)
(b)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 t≥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 t≥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.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 m=1, k1=k2=k3=k4=1, and ɛ=-3.1, with y0(0)=0.5, y˙0(0)=0, y1(0)=1, y˙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.
(a)
(b)


## 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
(29){m100m2}[x¨1x¨2]+{k1000}[x1x2]+{F1(x1,x2)F2(x1,x2)}=0,
and the nonlinear restoring forces are given by
(30)F1(x1,x2)=ɛx13-kf(x2-x1)1-((x2-x1)/x0)2,F2(x1,x2)=kf(x2-x1)1-((x2-x1)/x0)2,
with initial conditions x1(0)=x10, x˙1(0)=0, x2(0)=x20, and x˙2(0)=0. Here, m1 is the mass of the primary system, m2 represents the mass of the FENO system, k1 and ɛ are, respectively, the linear and nonlinear stiffness parameters of the main system, and kf and x0 represent the coupling stiffness and the maximum possible extension of the FENO system, respectively. By introducing the transformations y1=x1/x10 and y2=x2/x20 into (29), we get that
(31){y¨1y¨2}+{ωn12000}[y1y2]+{φ1x13-φ2A(y2-y1x10/x20)1-A2(y2-y1x10/x20)2φ3A(y2-y1x10/x20)1-A2(y2-y1x10/x20)2}=0,
with initial conditions given as y1(0)=1, y˙1(0)=0, y2(0)=1, and y˙2(0)=0, and A=x20/x0, φ1=ɛx102/m1, φ2=kfx0/(m1x10x0), and φ3=kfx0/(m2x20). 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
(32)F(y)=A(y2-y1x10/x20)1-A2(y2-y1x10/x20)2=Ay1-A2y2=(b1-3b3)y+4b3y3,
where y=(y2-y1x10/x20), and by definition
(33)b2n+1(A)=2π∫-11(1-y2)-1/2F(y)T2n+1(y)dy,
and T2n-1 are the Chebyshev polynomials of the first kind defined as
(34)T2n+1(x)=cos⁡[(2n+1)cos⁡-1(x)],hx∈[-1,1],n=0,1,2,….
By using (33), we obtain that
(35)b1=2A(11-A2-1),b3=2A3(41-A2+A2(1-31-A2)-4).
Thus, the FENO nonlinear restoring force F(y) can be written in equivalent form as
(36)F(y)≡αy+βy3,
in which
(37)α=1A{2(11-A2-1)hhhhh-6A2(41-A2+A2(1-31-A2)-4)},β=8A3{41-A2+A2(1-31-A2)-4}.
Based on (37), it is clear that the equivalent representation form (36) of the restoring force (32) holds if and only if x20<x0. Thus, (31) can be written as
(38){y¨1y¨2}+{ωn12000}[y1y2]+{φ1x13-φ2{α(y2-y1x10x20)+β(y2-y1x10x20)3}φ3{α(y2-y1x10x20)+β(y2-y1x10x20)3}}=0.
We next apply the aforementioned nonlinear transformation approach to (38), to obtain that
(39)y¨1+a1y1+a2y13=0;y¨2+a3y2+a4y23=0,
where the ai are given as
(40)a1=ωn12+φ2(2αx20(32x10u1-45x20u2)+βσ2{64x10x20u1u2-30x102u12-45x202u22})164x202u1,a2=φ1+φ2×(70αx203u2hhh+β(64x103u13-70x102x20u12u2+35x203u23))×164x203u13,a3=φ364x203u22×(×(-45x102u112+64x10x20u11u22-30x202u222))α(-90x10x20u11+64x203u22)+βx10u11hhhh×(-45x102u112+64x10x20u11u22-30x202u222)),a4=φ364×(64β+35(2αx10x20u11+βx103u113)x203u223-70βx10u11x20u2235(2αx10x20u11+βx103u113)x203u223).
To assess the degree of accuracy attaining by (39), let us consider the parameter values of m1=1, m2=1/2, k1=1, kf=1, x0=1, and ɛ=1, with x0(0)=1/2, x˙0(0)=0, x1(0)=0.99, and x˙1(0)=0. Figure 7 shows the amplitude-time response curves of (31) compared to those obtained by substituting (39)1 into (38)2. Here, the computed parameter values are α=-15.4771, β=37.0369, a1=0.3305, a2=0.5512, a3=-107.549, and a4=96.7734, with u1=-1.991, u2=-1.115, u11=-0.7, and u22=-4.6. In this case, the computed RMSE values for y1 and y2 were, respectively, 0.1965 and 0.2443. Table 5 shows the RMSE computed at different values of x20. Notice that when x20 tends to x0, 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 u1=-1.991, u2=-1.115, u11=-0.7, and u22=-4.6.
x
20
RMSEy1(x)
RMSEy2(x)
0.75
0.9981
1.1200
0.80
1.0245
1.0467
0.85
1.0330
0.9560
0.90
0.5851
0.5029
0.95
1.0122
0.8019
0.99
0.1965
0.2443Amplitude-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 m1=1, m2=1/2, k1=1, kf=1, x0=1, and ɛ=1, with x0(0)=1/2, x˙0(0)=0, x1(0)=0.99, x˙1(0)=0, u1=-1.991, u2=-1.115, u11=-0.7, and u22=-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.
(a)
(b)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 Figure9 having two masses, m1 and m2 which are connected by four elastic springs of stiffness k1 and k2 having an undeformed length, l0. Furthermore, the element with mass m1 is moving on a smooth horizontal bearing surface whose motion is restricted by an elastic linear spring of stiffness k3. By using Newton second law, it is easy to show that the equations of motion are
(41){1001}[x¨1x¨2]+{k3m1000}[x1x2]+{-2k1m1+B002k1m2+B1}[(x2-x1)(x2-x1)]+{F1(x1,x2)F2(x1,x2)}=0,
where the irrational restoring forces are given by
(42)F1(x1,x2)=C(x2-x1)1+A2(x2-x1)2,F2(x1,x2)=C1(x2-x1)1+A2(x2-x1)2.
Here B=-2k2/m1, B1=2k2/m2, C=2k2/m1, and C1=-2k2/m2. The system (41) can be cast in different form by setting x1=(w1-w2)/2 and x2=(w2+w1)/2, and after some algebraic manipulations, we get that
(43){1001}[w¨1w¨2]+{k3(2m1)k11k22-k3(2m1)}[w1w2]+{w2(C+C1)21+A2w22-w2(C1-C)21+A2w22}=0,
with initial conditions w1(0)=w10, w˙1(0)=0, w2(0)=w20, and w˙2(0)=0. If we set z1=w1/w10 and z2=w2/w20, then (43) becomes
(44){1001}[z¨1z¨2]+{k3(2m1)k11w20w10k22-k3w10(2m1w20)}[z1z2]+{z1(C+C1)w20w101+A2w22w202-z1(C1-C)w20w101+A2w22w202}=0,
where
(45)k11=12(-k3/2+2k1m1+2k1m2+B+B1),k22=12(k3/2+2k1m1+2k1m2+B1-B),
with z1(0)=1, z˙1(0)=0, z˙2(0)=1, and z2(0)=0. By using the Chebyshev polynomials of the first kind, the cubic-like representation form of the irrational restoring force becomes
(46)z11+A12w22w202≡αz1+βz13,
where
(47)α=4πA2{E(-A2)-K(-A2)-1A2hhhhhhh×((8+5A2)K(-A2)-(8+A2)E(-A2))1A2},β=163πA4×{(8+5A2)K(-A2)-(8+A2)E(-A2)}.
Here K(-A2) and E(-A2) are the complete elliptic integral of the first and second kind for the modulus A, respectively. Thus, (44) can be rewritten as
(48){1001}[z¨1z¨2]+{k3(2m1)k11w20w10k22-k3w10(2m1w20)}[z1z2]+{(C+C1)w20{αz1+βz13}w10-(C1-C){αz1+βz13}w20/w10}=0.
Application of the proposed nonlinear transformation approach to (48) yields
(49)z¨1+a1z1+a2z13=0;z¨2+a3z2+a4z23=0,
where the ai are given as
(50)a1=164m1w10u1(u2332k3w10u1+90αm1w10u2hhhhhhhhhhhh+90k11m1w20u2+45βm1w10u23),a2=-3564m1w10u13(2αw10u2+2k11w20u2+βw10u23),a3=-45k3w10u11-64αm1w20u22-64k22m1w20u2264m1w20u223,a4=--35k3w10u11-64βm1w20u22364m1w20u223.Figure 9
Schematic of a two-degree-of-freedom system with irrational restoring force.Figure10 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: m1=1, m2=1/2, k1=1, k2=1, k3=1, and A=1, with initial conditions given as w1(0)=4, w˙1(0)=0, w2(0)=2, and w˙2(0)=0. The predicted parameter values are α=4.7082, β=-2.1374, a1=-0.5442, a2=1.7564, a3=8.5590, and a4=0.1652, with u1=-0.68, u2=0.5, u11=-0.3, and u22=0.55. Notice from Figure 10 that our numerical integration predictions computed from (49)1 and (48)2 follow the amplitude-time response curves computed from (48) with RMSE values of 0.2095 for z1 and 0.1318 for z2, 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 ui and uii in order to satisfy (8).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 m1=1, m2=1/2, k1=1, k2=1, k3=1, A=1, w1(0)=4, w˙1(0)=0, w2(0)=2, w˙2(0)=0, α=4.7082, β=-2.1374, a1=-0.5442, a2=1.7564, a3=8.5590, and a4=0.1652, with u1=-0.68, u2=0.5, u11=-0.3, and u22=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.
(a)
(b)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 thex1(y1) coordinate system into the original equation of motion described by x2(y2), 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 coordinatey1 is used to find the dynamics response of the y0 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 t≥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.



---

*Source: 237234-2014-05-07.xml*