Accurate Solutions of Conservative Nonlinear Oscillators by the Enhanced Cubication Method

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.

1. Introduction

In this paper, we use the enhanced-cubication method to develop the approximate solution of some conservative nonlinear oscillators of the form 𝑥̈+𝑓(𝑥)=0;𝑥(0)=𝐴,𝑥̇(0)=0, (1) where 𝑓(𝑥) is the system restoring force which could have rational or irrational elastic terms, and 𝐴 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 𝑥̈+𝑓(𝑥)=0;𝑓(𝑥)=𝜀𝑥³𝐵+𝐶𝑥², (2) where 𝐵, 𝐶, 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] 𝑑²𝑦𝑑𝑡²+𝑓𝑦=0,𝑓𝑦=−𝑏𝑦+𝑎sin𝑦, (3) with initial conditions 𝑦(0)=𝐴 and 𝑦(̇0)=0, and system parameters 𝑎 and 𝑏. Next, we use the following transformation 𝑥=𝑦/𝐴 and rewrite (3) as 𝑑²𝑥𝑑𝑡²(𝑥)𝑓(𝑥)𝑎+𝑓=0,=−𝑏𝑥+𝐴sin(𝐴𝑥);𝑥(0)=1,𝑥̇(0)=0. (4) Then, we write 𝑓(𝑥) in equivalent form by using the Chebyshev polynomials of the first kind [3–6]: 𝑓(𝑥)=^𝑁_𝑛=0𝑏_2𝑛+1𝑥₁₀𝑇_2𝑛+1(𝑥), (5) where 𝑏_2𝑛+1=2𝜋⁺¹₋₁1√1−𝑥²𝑓(𝑥)𝑇_2𝑛+1(𝑥)𝑇𝑑𝑥,₁(𝑥)𝑇=𝑥;₃(𝑥)=4𝑥³𝑇−3𝑥;₅(𝑥)=16𝑥⁵−20𝑥³+5𝑥. (6)(7) Thus, the equivalent restoring force 𝑓(𝑥) can be written as 𝑓(𝑥)≡𝑏₁𝑞𝑇₁𝑦+𝑏₃𝑞𝑇₃𝑦+𝑏₅𝑞𝑇₅𝑦≈𝛼(𝐴)𝑥+𝛽(𝐴)𝑥³+𝛾(𝐴)𝑥⁵, (8) where the coefficients 𝛼(𝐴), 𝛽(𝐴), and 𝛾(𝐴) depend on 𝑓(𝑥), the amplitude of oscillation, 𝐴, and the Chebyshev coefficient terms. Therefore, the nonlinear differential equation (4) can be replaced by an equivalent cubic-quintic Duffing equation of the form 𝑑²𝑥𝑑𝑡²+𝛼(𝐴)𝑥+𝛽(𝐴)𝑥³+𝛾(𝐴)𝑥⁵≈0. (9) 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) where 𝛿(𝐴) and 𝜖(𝐴) can be found by using the equivalent nonlinearization method [2]. In this method, we replace the terms 𝛼(𝐴)𝑥+𝛽(𝐴)𝑥³+𝛾(𝐴)𝑥⁵ by the cubic polynomial 𝛿(𝐴)𝑥+𝜖(𝐴)𝑥³ that satisfies 𝐹(𝛿,𝜖)=0^𝜎𝛼𝑥+𝛽𝑥³+𝛾𝑥⁵−𝛿𝑥−𝜖𝑥³²𝑑𝑥⟶min,𝜕𝐹𝜕𝛿=0,𝜕𝐹𝜕𝜖=0. (11) By following this procedure, the coefficients 𝛿(𝐴) and 𝜖(𝐴) can be determined by using the expressions 𝛿(𝐴)=1,𝜖(𝐴)=12121𝛼−5𝛾𝜎99𝛽+10𝛾𝜎², (12) where the value of 𝜎 is fitted to satisfy (11). Thus, the equivalent representation form of (4) is given by 𝑑²𝑥𝑑𝑡²+𝛿(𝐴)𝑥+𝜖(𝐴)𝑥³≈0 (13) whose exact angular frequency-amplitude relationship is given by Ω_𝐶=𝜋𝜔_𝐶2𝐾𝑘²_𝐶, (14) where 𝜔_𝐶=𝛿(𝐴)(𝐴),𝑘+𝜖_𝐶=𝜖(𝐴)2(𝛿(𝐴)(𝐴.+𝜖)) (15)(16) Here, the coefficients 𝛼(𝐴), 𝛽(𝐴), and 𝛾(𝐴) are obtained by using (6)–(8). It is easy to show that these coefficients are given as 𝛼(𝐴)=6𝑎𝐴⁴𝐴𝐴²−80𝐽₁(𝐴)−16𝐴²−20𝐽₂(𝐴)𝛽(𝐴)=−𝑏,32𝑎𝐴⁴−𝐴𝐴²−60𝐽₁(𝐴)+27𝐴²−120𝐽₂(𝐴),𝛾(𝐴)=32𝑎𝐴⁴𝐴𝐴²−48𝐽₁(𝐴)−12𝐴²−16𝐽₂(𝐴), (17) where 𝐽₁(𝐴) and 𝐽₂(𝐴) are the first and second order Bessel functions of the first kind, and 𝐾(−𝐴²) and 𝐸(−𝐴²) are the complete elliptic integral of the first and second kinds for the modulus 𝐴, 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^∘≤𝐴≤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% for 0∘≤𝐴≤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 𝑓𝑦=ℎsgn𝑦𝑦^𝛼. (18) 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 ℎ is a constant parameter whose value depends on the physical system under consideration. We next use the transformation 𝑦=𝑥/𝐴 and write the corresponding equation of motion in the form 𝑑²𝑥𝑑𝑡²(𝑥)𝑓(𝑥)+𝑓=0,=𝑐₁sgn(𝑥)|𝑥|^𝛼;𝑐₁=ℎ𝐴^(𝛼−1);𝑥(0)=1,𝑥̇(0)=0. (19) Then, we use Chebyshev polynomial expansion and write the restoring force as 𝑓(𝑥)≈𝛼(𝐴)𝑥+𝛽(𝐴)𝑥³+𝛾(𝐴)𝑥⁵, (20) where 𝛼(𝐴)=3𝑐₁(𝑚−5)(𝑚−3)Γ[𝑚/2+1]2√)],𝛽(𝐴)𝜋Γ[(7+𝑚/2=−8𝑐₁(𝑚−5)(𝑚−1)Γ[𝑚/2+1]√)],𝛾(𝐴)=𝜋Γ[(7+𝑚/28𝑐₁(𝑚−3)(𝑚−1)Γ[𝑚/2+1]√)].𝜋Γ[(7+𝑚/2 (21) Here, Γ[𝑤] 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 1𝐸−06 for the interval 0≤𝐴≤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, ℎ=1, 𝐴=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 1𝐸−06.

Table 1 

Fitting parameter 𝜎 for selected values of the exponent 𝛼.

Figure 2 

Amplitude-time response curve by assuming that 𝛼=√2003, ℎ=1, 𝐴=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 𝑓𝑦=𝜆₁𝑦+𝜆₂sgn𝑦𝑦^𝛼11+𝛼₃𝑦^2𝛼2, (22) where 𝛼₁, 𝛼₂, 𝛼₃, 𝜆₁, and 𝜆₂ are system parameter values. If we introduce the transformation 𝑥=𝑦/𝐴, then the general oscillator differential equation becomes 𝑑²𝑥𝑑𝑡²+𝑓(𝑥)=0,𝑓(𝑥)=𝜆₁𝑥+𝜆₂𝐴^(𝛼1−1)sgn𝑦𝑦^𝛼11+𝛼₃𝐴²𝑦^2𝛼2, (23) with initial conditions 𝑥(0)=1 and 𝑥̇(0)=0. We next use Chebyshev polynomials to find the equivalent representation of 𝑓(𝑥) as 𝑓(𝑥)≈𝛼(𝐴)𝑥+𝛽(𝐴)𝑥³+𝛾(𝐴)𝑥⁵, (24) where 𝛼(𝐴)=𝜆₁+−1+(−1)^𝛼1𝐴^𝛼1−5𝜆₂Γ𝛼1+¹2×−802+𝛼₁3+𝛼₁+16𝐴²2+𝛼₁30+7𝛼₁−14𝛼₂𝛼₃+35𝐴⁴3+𝛼₁−2𝛼₂5+𝛼₁−2𝛼₂𝛼²₃×₂𝐹₁2+𝛼₁/2,𝛼₂,3+𝛼₁/2,−𝐴²𝛼₃Γ3+𝛼₁/2+162+𝛼₁1+𝐴²𝛼₃×53+𝛼₁+𝐴²15+2𝛼₁−4𝛼₂𝛼₃×₂𝐹₁4+𝛼₁/2,𝛼₂,3+𝛼₁/2,−𝐴²𝛼₃Γ3+𝛼₁×√/2𝜋3+𝛼₁−2𝛼₂5+𝛼₁−2𝛼₂𝛼²₃⁻¹,𝛽(𝐴)=√𝜋Γ7+𝛼₁24−1+(−1)^𝛼1𝐴^𝛼1−1×𝜆₂Γ𝛼1+¹2×5+𝛼₁73+𝛼₁×₂𝐹₁2+𝛼₁2,𝛼₂,3+𝛼₁2,−𝐴²𝛼₃−262+𝛼₁×₂𝐹₁4+𝛼₁2,𝛼₂,5+𝛼₁2,−𝐴²𝛼₃+202+𝛼₁4+𝛼₁×₂𝐹₁6+𝛼₁2,𝛼₂,7+𝛼₁2,−𝐴²𝛼₃⁻¹,𝛾(𝐴)=−1√𝜋Γ7+𝛼₁24−1+(−1)^𝛼1×𝐴^𝛼1−1𝜆₂Γ𝛼1+¹2×53+𝛼₁×5+𝛼₁×₂𝐹₁2+𝛼₁2,𝛼₂,3+𝛼₁2,−𝐴²𝛼₃+42+𝛼₁×−55+𝛼₁₂𝐹₁4+𝛼₁2,𝛼₂,5+𝛼₁2,−𝐴²𝛼₃+54+𝛼₁₂𝐹₁6+𝛼₁2,𝛼₂,7+𝛼₁2,−𝐴²𝛼₃⁻¹. (25) Here, ₂𝐹₁[𝑎,𝑏,𝑐,𝑧]=^∞_𝑛=0(𝑎)_𝑛(𝑏)_𝑛(𝑐)_𝑛𝑧𝑛!^𝑛 (26) 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, 𝜆₁=1, and 𝜆₂=−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 𝐴=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 𝐴. 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 𝐴, that is, for 𝐴<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 amplitude 𝐴. Here, 0.1≤𝐴≤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, 𝛼₃=−1, 𝜆₁=0, and 𝜆₂=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 𝐴=0.9953 with 𝜎=0.9985. This value is lower than the error value of −3.53% for 𝐴=0.9 reported in [14, 15]. Figure 7 shows the amplitude-time response curves with an initial oscillation amplitude value of 𝐴=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 𝐴 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 𝐴. Finally, we have plotted the angular frequency ratio versus the amplitude 𝐴. 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<𝐴<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.

Figure 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 𝐴=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 𝐴=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 amplitude 𝐴.

Of course, other equivalent restoring force representations can be studied by assuming different parameters values of 𝛼₁, 𝛼₂, 𝛼₃, 𝜆₁, and 𝜆₂ 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 𝑑²𝑥𝑑𝑡²(𝑥)𝑓(𝑥)=1+𝑓=0,𝐴²𝑥^𝑚,𝑥(0)=1,𝑥̇(0)=0, (27) where 𝑚 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 𝛼(𝐴)=3−−1/𝐴^2𝑚+1/𝐴^2𝑚(3+𝑚)(5+𝑚)Γ[1−𝑚/2]4√)],𝛽(𝐴)=4−𝜋Γ[(7−𝑚/2−1/𝐴^2𝑚−1/𝐴^2𝑚(1+𝑚)(5+𝑚)Γ[1−𝑚/2]√)],𝛾(𝐴)=4−𝜋Γ[(7−𝑚/2−1/𝐴^2𝑚+1/𝐴^2𝑚(1+𝑚)(3+𝑚)Γ[1−𝑚/2]√)].𝜋Γ[(7−𝑚/2 (28)

It was found in [8, 16, 17] that the exact angular frequency value of (27) when 𝑚=1 is given as 𝜔_ex=1.2533131𝐴. (29)

To find the approximate angular frequency value by using the enhanced cubication approach, we now substitute (28) into (14) and (16), to get that Ω_𝐶=𝜋(5/14)224𝜎²−48𝜎⁴−1893𝐴𝐾𝑘²,𝑘²=756−560𝜎²945−1120𝜎²+240𝜎⁴, (30)(31) 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.061800921 and 𝜎₂=1.3853631. In this way, we have found two possible approximate solutions to (27) since Ω_𝐶 is equal to (29), that is, Ω_𝐶=1.2533131/𝐴. Figure 9 shows the amplitude-time plots obtained by using the values of 𝜎₁, 𝜎₂, and the solution of (13) 𝑥(𝑡)=𝑐𝑛𝜔_𝐶𝑡,𝑘²_𝐶 (32) with initial conditions 𝑥(0)=𝐴=1, and 𝑥̇(0)=0. Notice from Figure 9(a) that when we use the root value of 𝜎₁, the estimated amplitude-time curve follows surprisingly well the exact numerical solution curve. In this case, the modulus of the elliptic function, 𝑘, has the value of 𝑘=3.13864𝑖, where 𝑖=√−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 𝑑=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 𝜎₂=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 Ω_𝐶, 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 𝜎₂ is not a valid solution to (27). Furthermore, one can numerically prove that this solution does not satisfy (27) at all times 𝑡. Thus, we can conclude that only the fitting parameter value of 𝜎₁=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 𝐴>0.001, our enhanced solution provided the same qualitative and quantitative behaviors than those of the numerical integration solutions of (27).

(a)

(b)

Figure 9 

Amplitude-time response curves of a singular oscillator with initial conditions 𝑥(0)=1,  𝑥̇(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.

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^∘≤𝐴≤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 1𝐸−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, 𝜆₁=1, and 𝜆₂=−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<𝐴≤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 𝜎₂=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.

Acknowledgment

This work was partially funded by the Tecnológico de Monterrey, Campus Monterrey through the Research Chairs in Nanotechnology and Intelligent Machines.