# Accurate Solutions of Conservative Nonlinear Oscillators by the Enhanced Cubication Method

Authors: Alex Elías-Zúñiga; Oscar Martínez-Romero
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

The enhanced cubication method is applied to develop
approximate solutions for the most common nonlinear oscillators found
in the literature. It is shown that this procedure leads to amplitude-time
response curves and angular frequency values with maximum relative
errors lower than those found by previously developed approximate solutions.


---

## Body

## 1. Introduction

In this paper, we use the enhanced-cubication method to develop the approximate solution of some conservative nonlinear oscillators of the form(1)x¨+f(x)=0;x(0)=A,x˙(0)=0,
where f(x) is the system restoring force which could have rational or irrational elastic terms, and A is the oscillation amplitude. We recall that the enhanced-cubication method uses the Chebyshev polynomial expansion to replace the system restoring force by an equivalent one. Then, the nonlinearization method is used to find and equivalent Duffing equation for (1) that is valid for the complete range of oscillation amplitudes [1]. Elías-Zúñiga and coworkers used this technique to derive the approximate solution of the Duffing-harmonic oscillator of the form
(2)x¨+f(x)=0;f(x)=ɛx3(B+Cx2),
where B, C, and ɛ are constant system parameters [1]. They found by using this method that the approximate angular frequency values agree well with the exact ones for the whole range of oscillation amplitude values with a maximum relative error value that is less than 0.055%.Here in this paper, we examine the application of the enhanced-cubication method to obtain approximate solutions of the general pendulum equation, the generalized power-form elastic term oscillator, the general restoring force term oscillators, and of the singular oscillator that arises in plasma physics. To assess the accuracy of the enhanced-cubication method, we compare our derived approximate angular frequency values with those obtained by numerical integration. It is shown that in most cases, our derived solutions have maximum relative errors that are lower than those obtained by using other solution techniques.We will next start with the application of our enhanced cubication approach by deriving the approximate solution of the general pendulum equation.


## 2. A General Pendulum Equation

We first consider the general pendulum equation of the form [2]
(3)d2ydt2+f(y)=0,f(y)=-by+asiny,
with initial conditions y(0)=A and y˙(0)=0, and system parameters a and b. Next, we use the following transformation x=y/A and rewrite (3) as
(4)d2xdt2+f(x)=0,f(x)=-bx+aAsin(Ax);x(0)=1,x˙(0)=0.
Then, we write f(x) in equivalent form by using the Chebyshev polynomials of the first kind [3–6]:
(5)f(x)=∑n=0N‍b2n+1(x10)T2n+1(x),
where
(6)b2n+1=2π∫-1+1‍11-x2f(x)T2n+1(x)dx,(7)T1(x)=x;T3(x)=4x3-3x;T5(x)=16x5-20x3+5x.
Thus, the equivalent restoring force f(x) can be written as
(8)f(x)≡b1(q)T1(y)+b3(q)T3(y)+b5(q)T5(y)≈α(A)x+β(A)x3+γ(A)x5,
where the coefficients α(A), β(A), and γ(A) depend on f(x), the amplitude of oscillation, A, and the Chebyshev coefficient terms. Therefore, the nonlinear differential equation (4) can be replaced by an equivalent cubic-quintic Duffing equation of the form
(9)d2xdt2+α(A)x+β(A)x3+γ(A)x5≈0.
Since the cubication procedure requires a cubic polynomial representation of the restoring force, we shall next transform the restoring force term of (9) into an equivalent cubic polynomial
(10)f(x)=α(A)x+β(A)x3+γ(A)x5≡δ(A)x+ϵ(A)x3,
where δ(A) and ϵ(A) can be found by using the equivalent nonlinearization method [2]. In this method, we replace the terms α(A)x+β(A)x3+γ(A)x5 by the cubic polynomial δ(A)x+ϵ(A)x3 that satisfies
(11)F(δ,ϵ)=∫0σ‍(αx+βx3+γx5-δx-ϵx3)2dx⟶min,∂F∂δ=0,∂F∂ϵ=0.
By following this procedure, the coefficients δ(A) and ϵ(A) can be determined by using the expressions
(12)δ(A)=121(21α-5γσ),ϵ(A)=19(9β+10γσ2),
where the value of σ is fitted to satisfy (11). Thus, the equivalent representation form of (4) is given by
(13)d2xdt2+δ(A)x+ϵ(A)x3≈0
whose exact angular frequency-amplitude relationship is given by
(14)ΩC=πωC2K(kC2),
where
(15)ωC=δ(A)+ϵ(A),(16)kC=ϵ(A)2(δ(A)+ϵ(A)).
Here, the coefficients α(A), β(A), and γ(A) are obtained by using (6)–(8). It is easy to show that these coefficients are given as
(17)α(A)=6aA4(A(A2-80)J1(A)-16(A2-20)J2(A))-b,β(A)=32aA4(-A(A2-60)J1(A)+2(7A2-120)J2(A)),γ(A)=32aA4(A(A2-48)J1(A)-12(A2-16)J2(A)),
where J1(A) and J2(A) are the first and second order Bessel functions of the first kind, and K(-A2) and E(-A2) are the complete elliptic integral of the first and second kinds for the modulus A, respectively. Figure 1 illustrates the relative errors achieved by comparing the approximate angular frequency values computed from (14) with the exact ones, obtained by using Radhakrishnan et al. procedure introduced in [7]. As we may see from Figure 1 and for 0∘≤A≤113∘, the relative error values are not bigger than 0.2%.Figure 1
Prediction of the error attained by computing the general pendulum angular frequency oscillator values by the enhanced cubication method. Notice that the maximum relative error attained between the approximate and the exact angular frequency values does not exceed of 0.2% for0∘≤A≤113∘.To further assess the applicability of the enhanced cubication approach, we next derive the solution of an oscillator with fractional restoring form.


## 3. A Generalized Power-Form Elastic Term Oscillator

In this case, the restoring force term in (3) has the form
(18)f(y)=hsgn(y)|y|α.
The restoring force given by (18) is similar in form to the one proposed by Pilipchuk [9], and due to its oddness with respect to time, the exponent α can continuously take any nonnegative real value, such as odd, even, rational, or irrational that is, 0≤α<∞. Notice that h is a constant parameter whose value depends on the physical system under consideration. We next use the transformation y=x/A and write the corresponding equation of motion in the form
(19)d2xdt2+f(x)=0,f(x)=c1sgn(x)|x|α;c1=hA(α-1);x(0)=1,x˙(0)=0.
Then, we use Chebyshev polynomial expansion and write the restoring force as
(20)f(x)≈α(A)x+β(A)x3+γ(A)x5,
where
(21)α(A)=3c1(m-5)(m-3)Γ[m/2+1]2πΓ[(7+m)/2],β(A)=-8c1(m-5)(m-1)Γ[m/2+1]πΓ[(7+m)/2],γ(A)=8c1(m-3)(m-1)Γ[m/2+1]πΓ[(7+m)/2].
Here, Γ[w] is the Euler gamma function. Table 1 shows the corresponding values of the fitting parameter σ that satisfies (11), for selected values of α exponent. In all cases, we found that the maximum relative errors attained by comparing the approximate angular frequency values computed from (14) with the exact ones never exceed 1E-06 for the interval 0≤A≤100. The amount of error attained by our enhanced cubication approach is remarkably lower than those previously found in the literature. See, for instance, [10–12] and references cited therein. This confirms the usefulness of our enhanced cubication approach to obtain approximate solutions of nonlinear oscillators with generalized power-form elastic terms. For illustrative purposes, Figure 2 shows the corresponding amplitude-time curve by assuming that α=2003, h=1, A=1, and σ=0.99608975. There, the solid line represents the numerical integration solution of (19) obtained by using the fourth-order Runge-Kutta algorithm provided by the Mathematica symbolic package, while the red dashed lines represent the enhanced cubication approximate solution. As we can see from Figure 2, both solutions agree well. In fact, the relative error between both solutions is less than 1E-06.Table 1
Fitting parameterσ for selected values of the exponent α.
Exponent value
Fitting parameter value
α
σ
1/3
1.0481977
2/3
1.045645
4/3
1.041253
5/3
1.039345
7/5
1.040857
2003
0.99608975Figure 2
Amplitude-time response curve by assuming thatα=2003, h=1, A=1, and σ=0.99608975. Here, the solid lines represents the numerical integration solution of (19), while the dashed lines represent the enhanced cubication approximate solution.We next develop the approximate solution of a most general restoring force-type that involves terms that describing some of the most common conservative nonlinear oscillators found in the literature.


## 4. A General Restoring Force Term

Here, we use the enhanced cubication approach to develop the approximate solution of a general oscillator described by (1) with a restoring force given by
(22)f(y)=λ1y+λ2sgn(y)|y|α1(1+α3y2)α2,
where α1, α2, α3, λ1, and λ2 are system parameter values. If we introduce the transformation x=y/A, then the general oscillator differential equation becomes
(23)d2xdt2+f(x)=0,f(x)=λ1x+λ2A(α1-1)sgn(y)|y|α1(1+α3A2y2)α2,
with initial conditions x(0)=1 and x˙(0)=0. We next use Chebyshev polynomials to find the equivalent representation of f(x) as
(24)f(x)≈α(A)x+β(A)x3+γ(A)x5,
where
(25)α(A)=λ1+([(4+α1)2,α2,(3+α1)2,-A2α3]Γ[(3+α1)2](-1+(-1)α1)Aα1-5λ2Γ[1+α12]×([(4+α1)2,α2,(3+α1)2,-A2α3]Γ[(3+α1)2]-(80(2+α1)(3+α1)+16A2(2+α1)(30+7α1-14α2)α3+35A4(3+α1-2α2)(5+α1-2α2)α32)×2F1[(2+α1)/2,α2,(3+α1)/2,-A2α3]Γ[(3+α1)/2]+16(2+α1)(1+A2α3)×(5(3+α1)+A2(15+2α1-4α2)α3)×2F1[(4+α1)/2,α2,(3+α1)/2,-A2α3]Γ[(3+α1)/2][(4+α1)2,α2,(3+α1)2,-A2α3]Γ[(3+α1)2]))×(π(3+α1-2α2)(5+α1-2α2)α32)-1,β(A)=((πΓ[(7+α1)2])4(-1+(-1)α1)Aα1-1×λ2Γ[1+α12]×([(6+α1)2,α2,(7+α1)2,-A2α3](5+α1)([(4+α1)2,α2,(5+α1)2,-A2α3]7(3+α1)×2F1[(2+α1)2,α2,(3+α1)2,-A2α3]-26(2+α1)×2F1[(4+α1)2,α2,(5+α1)2,-A2α3])+20(2+α1)(4+α1)×2F1[(6+α1)2,α2,(7+α1)2,-A2α3]))-1,γ(A)=-1((πΓ[(7+α1)2])4(-1+(-1)α1)×Aα1-1λ2Γ[1+α12]×([(6+α1)2,α2,(7+α1)2,-A2α3]5(3+α1)×([(6+α1)2,α2,(7+α1)2,-A2α3](5+α1)×2F1[(2+α1)2,α2,(3+α1)2,-A2α3]+4(2+α1)×(-5(5+α1)2F1[(4+α1)2,α2,(5+α1)2,-A2α3(4+α1)2])+5(4+α1)2F1[(6+α1)2,α2,(7+α1)2,-A2α3(6+α1)2](6+α1)2)))-1.
Here,
(26)2F1[a,b,c,z]=∑‍n=0∞(a)n(b)n(c)nn!zn
is the Gauss hypergeometric function. Thus, the approximate solution to (23) is obtained by using (13), (14), (15), and (16).To assess the accuracy of our enhanced cubication solution, let us consider first the case of a nonlinear dynamical oscillator with an irrational restoring force for whichα1=1, α2=1/2, α3=1, λ1=1, and λ2=-1 [13]. Figure 3 shows that the maximum error attained by comparing our approximate angular frequency values with the exact ones is less than 0.05%. This error is bigger than the error value of 0.04% computed by Lai and Xiang in [13] by using the first two order approximations of the generalized Senator Bapat perturbation technique. Both solutions as well as the numerical integration of (23) are plotted in Figure 4 by considering that A=10. As we can see from Figure 4, all solutions have the same qualitative and quantitative behaviors. Here, the solid black line represents the numerical integration solution, the red dashed line our proposed solution, while the black dots describe the solution derived by Lai and Xiang [13]. For convenience, we have plotted in Figure 5 the frequency ratio versus the amplitude A. We may notice that both approximate angular frequency values are close to the exact one. However, we found that Senator Bapat solution fail for small values of A, that is, for A<0.0012. Of course, we can improve the percentage error value computed from our enhanced cubication method by considering additional terms in the Chebyshev polynomial expansion of the restoring force term. In the present paper, we shall not consider these additional Chebyshev polynomial expansion terms.Figure 3
Prediction of the error attained by computing the elastic wire angular frequency oscillator values by the enhanced cubication method. Notice that the maximum relative error between the approximate and the exact angular frequency values does not exceed of 0.05%.Figure 4
Amplitude-time response curves of a nonlinear oscillator with an irrational restoring force for whichα1=1, α2=1/2, α3=1, λ1=1, and λ2=-1 with a fitting parameter value of σ=1.03875. Here, the red dashed line, the black dots, and the green solid line represent the enhance cubication solution, Senator Bapat perturbation solution, and the numerical integration solution of (23), respectively.Figure 5
Comparison of the frequency ratios versus the amplitudeA. Here, 0.1≤A≤20.We next study the applicability of the enhanced cubication method by deriving the approximate solution of the finite extensibility nonlinear oscillator for whichα1=1, α2=1, α3=-1, λ1=0, and λ2=1 [14, 15]. As we can see from Figure 6, the maximum relative error between the approximate angular frequency value and the exact one is not bigger than -3.451% at A=0.9953 with σ=0.9985. This value is lower than the error value of -3.53% for A=0.9 reported in [14, 15]. Figure 7 shows the amplitude-time response curves with an initial oscillation amplitude value of A=0.9953. Here, we compare our cubication solution (red dashed line) with respect to the approximate Beléndez et al. solution (black dots) and with the numerical integration solution of (23) (black solid line). As one can see from Figure 7, our derived solution is closer to the numerical integration one. Furthermore, we have computed the approximate angular frequency values of our derived solution for several values of A and compared these with respect to the exact values obtained from numerical integrations and with those values computed from the approximate solution derived by Beléndez and coworkers in [15]. Table 2 shows that in general, our enhanced angular frequency value is closer to the exact numerical one even for larger values of A. Finally, we have plotted the angular frequency ratio versus the amplitude A. It is clear from Figure 8 that both approximate solutions follow well the exact numerical value, however, Beléndez et al. solution shows higher deviation on 0.9<A<0.9953. This confirms the applicability of our enhanced-cubication approach to obtain accurate solutions of strongly nonlinear oscillators.Table 2
Comparison of the exact and the approximate angular frequency values for a finite extensibility nonlinear oscillator.
A
ω
exact
ω
Enhanced
/
ω
exact
ω
Belendez
/
ω
exact
0.5
1.1113
0.9988
1.0000
0.75
1.3369
0.9922
1.0004
0.85
1.5412
0.9848
0.9991
0.95
2.0355
0.9723
0.9707
0.98
2.4600
0.9776
0.9098
0.985
2.5901
0.9836
0.8861
0.9953
3.0856
1.0345
0.7875Figure 6
Finite extensibility nonlinear oscillator for whichα1=1, α2=1, α3=-1, λ1=0, and λ2=1. The maximum relative error value of -3.451% between the approximate angular frequency value and the exact one occurs at the oscillation amplitude of A=0.9953 with a fitting parameter value of σ=0.9985.Figure 7
Amplitude-time response curves of the finite extensibility nonlinear oscillator for whichα1=1, α2=1, α3=-1, λ1=0, λ2=1, and A=0.9953 with a fitting parameter value of σ=0.9985. Here, the red dashed line, the black dots, and the black solid line represent the enhance cubication solution, Beléndez et al. approximate solution, and the numerical integration solution of (23), respectively.Figure 8
Comparison of the frequency ratios versus the amplitudeA.Of course, other equivalent restoring force representations can be studied by assuming different parameters values ofα1, α2, α3, λ1, and λ2 in (22).


## 5. A Singular Oscillator

As a final example, we now focus our attention on deriving the approximate solution of the singular oscillator that describes the motion of injected electrons in plasma physics. In this case, the corresponding equation of motion is given as(27)d2xdt2+f(x)=0,f(x)=(1A2x)m,x(0)=1,x˙(0)=0,
where m is a positive entire number that usually takes the value of one [8, 16]. By using (6)–(8), we can write the approximate Chebyshev equivalent representation form of the restoring force as
(28)α(A)=3(-(-1/A2)m+(1/A2)m)(3+m)(5+m)Γ[1-m/2]4πΓ[(7-m)/2],β(A)=4(-(-1/A2)m-(1/A2)m)(1+m)(5+m)Γ[1-m/2]πΓ[(7-m)/2],γ(A)=4(-(-1/A2)m+(1/A2)m)(1+m)(3+m)Γ[1-m/2]πΓ[(7-m)/2].It was found in [8, 16, 17] that the exact angular frequency value of (27) when m=1 is given as
(29)ωex=1.2533131A.To find the approximate angular frequency value by using the enhanced cubication approach, we now substitute (28) into (14) and (16), to get that
(30)ΩC=π(5/14)(224σ2-48σ4-189)3AK[k2],(31)k2=(756-560σ2)(945-1120σ2+240σ4),
where σ is a fitting parameter that satisfies (11). In this particular case, we note that σ can be found by equating (29) with (30) since this provides an equation that depends only on σ. By numerically solving the resulting equation, we found two real roots for the fitting parameter σ. These real roots have the values of σ1=1.061800921 and σ2=1.3853631. In this way, we have found two possible approximate solutions to (27) since ΩC is equal to (29), that is, ΩC=1.2533131/A. Figure 9 shows the amplitude-time plots obtained by using the values of σ1, σ2, and the solution of (13)
(32)x(t)=cn(ωCt,kC2)
with initial conditions x(0)=A=1, and x˙(0)=0. Notice from Figure 9(a) that when we use the root value of σ1, the estimated amplitude-time curve follows surprisingly well the exact numerical solution curve. In this case, the modulus of the elliptic function, k, has the value of k=3.13864i, where i=-1. Also, we have plotted in Figure 9(a) the approximate solution of (27) derived in [8] for which its predicted numerical solution has a dispersion value of d=0.044 when compared to the numerical integration one. As we can see from Figure 9(a), these solutions are almost the same. In Figure 9(a), the solid line represents the derived enhanced-cubication solution, the dashed line represents the numerical integration solution, while the black dots are the estimated prediction values obtained by using the approximate solution derived in [8]. We next use the parameter value of σ2=1.3853631 and plot the amplitude-time curves by using (32). In this case, although our enhanced-cubication solution provides the same angular frequency value ΩC, the amplitude-time response curve showed in Figure 9(b) differs from the curve obtained by numerical integration. Therefore, we must conclude that the value of σ2 is not a valid solution to (27). Furthermore, one can numerically prove that this solution does not satisfy (27) at all times t. Thus, we can conclude that only the fitting parameter value of σ1=1.061800921 provides a solution that predicts well the trajectory of the singular oscillator. Finally, we must mention that during the derivation of the equivalent restoring force representation of (27), Chebyshev polynomial expansion provides a poor fit to this force term. Due to this finding, we ran several computer simulations by considering different initial condition values. Then, we plotted the amplitude-time response curve of our derived solution and compared this with the exact numerical one. We found in all cases that for A>0.001, our enhanced solution provided the same qualitative and quantitative behaviors than those of the numerical integration solutions of (27).Amplitude-time response curves of a singular oscillator with initial conditionsx(0)=1,  x˙(0)=0 (a) σ1=1.061800921, and (b) σ2=1.3853631. Here, the solid line represent the derived enhanced-cubication solution, the dashed lines represents the numerical integration solution, while the black dots represent the approximate solution derived in [8] by using Jacobi elliptic functions.
(a)
(b)


## 6. Conclusions

In this paper, we have used the enhanced cubication procedure introduced by Elías-Zúñiga et al. in [1] to develop approximate solutions of some strongly nonlinear oscillators that described the dynamical behavior response of several physical and engineering systems. During the solution processes, we used the enhanced cubication approach to replace the original equation of motion by the homogeneous Duffing equation. This transformation is achieved by computing the Chebyshev polynomial representation of the restoring force terms. Then, we used the homogenous Duffing equation to find the angular frequency expression since it is well known that the Duffing equation has an exact solution that depends on Jacobi elliptic functions. Based on these expressions, we have assessed the accuracy of the angular frequency equations for the complete range of oscillations amplitudes. For instance, we have examined the general pendulum response, given by (4), at the interval oscillation values of 0∘≤A≤113∘. In this case, we have found that the maximum relative error does not exceed 0.2%. Then, we used our enhanced cubication method to derive the approximate solution of a generalized power-form elastic term oscillator and found, in all cases, that the maximum error attained is less than 1E-06 which is remarkably lower when compared to other approximate solutions [10–12]. However, when we examined the nonlinear oscillator described by (23) for which α1=1, α2=1/2, α3=1, λ1=1, and λ2=-1, the percentage error value found by the enhanced cubication method is bigger than that obtained by using the generalized Senator Bapat perturbation technique [13]. Of course, further improvement could be achieved in our enhanced perturbation method if we consider more terms on the Chebyshev polynomial expansion (6).We next used the enhanced cubication approach and derived the approximate solution of the finite extensibility nonlinear oscillator and found a maximum error of about-3.451% for 0<A≤0.99853 which represents an improvement of the error values estimated by using other perturbation techniques [14, 15]. Finally, we examined the dynamical response of a singular oscillator (27) and found that our enhanced cubication technique provides two approximate solutions by using the same initial conditions. However, the solution related to the fitting parameter value of σ2=1.3853631 provides an amplitude-time curve whose trajectory differs from the numerical one, and then it was discarded as a possible solution to (27).Based on the previous results, we can conclude that our proposed enhanced-cubication method provides good approximate solutions to strongly nonlinear oscillators even for irrational restoring forces.In a forthcoming paper, we shall introduce some modifications to our enhanced cubication method to account for linear or nonlinear damping effects of strongly nonlinear oscillators.



---

*Source: 842423-2013-04-03.xml*