Approximate Solution for the Duffing-Harmonic Oscillator by the Enhanced Cubication Method

Abstract

The cubication and the equivalent nonlinearization methods are used to replace the original Duffing-harmonic oscillator by an approximate Duffing equation in which the coefficients for the linear and cubic terms depend on the initial oscillation amplitude. It is shown that this procedure leads to angular frequency values with a maximum relative error of 0.055%. This value is 21% lower than the relative errors attained by previously developed approximate solutions.

1. Introduction

Many techniques have been developed to obtain the approximate solution of second-order nonlinear differential equations such as multiple scales [1], harmonic balance [2], averaging [3], Lindstedt-Poincaré [4], to say a few. However, most of these methods could provide good approximate solutions if the nonlinear terms have small magnitude values and the system is subjected to small oscillation amplitudes. In an attempt to deal with strong nonlinearities and larger oscillation amplitudes for the Duffing-harmonic oscillator system, other types of solution techniques have been developed [5–14]. Among these techniques, there is the cubication approach introduced by Yuste and Sánchez in [15] which consists in replacing the system restoring force 𝑓(𝑥) by an equivalent cubic polynomial expression 𝑎₃𝑥³, where the value of 𝑎₃ is determined by using a weighted mean-square method or by using the principle of harmonic balance [16, 17]. Beléndez and coworkers used this idea and replaced the original second-order differential equation by the well-known Duffing equation which has an exact solution. During their solution process, they used Chebyshev polynomial expansions to replace the original restoring force by an equivalent form which provides approximate angular frequency expressions that are valid for the complete range of oscillations amplitudes [11–13].

Here in this paper, we combine the cubication and the nonlinearization methods to obtain the approximate angular frequency value of the Duffing-harmonic oscillator of the form: 𝑥̈+𝑓(𝑥)=0;𝑓(𝑥)=𝜀𝑥³𝐵+𝐴𝑥² (1.1) and show that it has the smallest relative error value when compare to the numerical integration and with others approximate solutions. In (1.1), 𝑥 represents system displacement, and 𝜀, 𝐴, and 𝐵 are constant parameters. This harmonic oscillator was first considered by Mickens to investigate the inclusion of higher harmonics in the method of harmonic balance [18]. Then, Radhakrishnan et al. studied (1.1) with 𝜀=1, 𝐴=𝜆, and 𝐵=1 and obtained the exact value of the circular frequency by numerical integration of its energy relation [19]. By applying the method of harmonic balance, Mickens in [2] derived an analytical approximation for the dimensionless form of (1.1) for which 𝜀=1, 𝐴=𝜆, and 𝐵=1 and found a good estimate value for the angular frequency expression. Tiwari et al. derived an approximate frequency-amplitude relation to (1.1) by assuming a single term solution and by following the Ritz procedure [5]. The first-order harmonic balance method via first Fourier coefficient was used by Hu in [6]. He found good agreement with the result obtained by the Ritz procedure. A new method to solve (1.1) by combining Newton's and the harmonic balance methods was derived by Lim et al. in [7] with results that are valid for the complete range of oscillation amplitudes. Öziş and Yildirim applied He's energy balance method to construct the frequency-amplitude response of an equivalent form of (1.1) with numerical results that agree well with the exact ones [8]. By applying the modified rational harmonic balance method, Beléndez and coworkers obtained the approximate solution of nonlinear oscillators in which the restoring force has a rational expression [9]. They obtained approximate solutions that agree well with the exact solutions for the whole range of values of oscillation amplitudes with a relative error value that is less than 0.40% [9]. The Duffing-harmonic oscillator has been also studied by Beléndez et al. using the second approximation of the modified homotopy perturbation method with a relative error of the frequency-amplitude value lower than 0.078% [10]. By using the iterative homotopy harmonic balancing approach, Guo and coworkers obtained an approximate solution to (1.1) with a discrepancy between the estimated angular frequency values and the exact ones as low as 0.094% [20]. Beléndez et al. in [12] developed a cubication method in which the restoring force is expanded in Chebyshev polynomials to obtain an equivalent cubic polynomial equation. They found that the maximum relative error between the approximate frequency-amplitude relation and the exact one does not exceed of 0.071%. As we may see, there are different methods that have been applied to solve (1.1) with the main purpose of finding the angular frequency values that are closer to the numerical ones. In the next section, we shall briefly review some of these solutions.

2. Review of Some Approximate Solutions

Here, we briefly review approximate solutions to (1.1) that have been derived by using different solution techniques. For instance, Mickens in [2] studied the dynamical response of the Duffing-harmonic oscillator: 𝑑²𝑦𝑑𝜏²+𝜀𝑦³𝐴𝑦²(0)+𝐵=0,𝑦=𝑦₁₀,𝑦̇(0)=0, (2.1) for which the parameters 𝜀, 𝐴, and 𝐵 are nonnegative. He used the following transformations: 𝑦=𝐵𝐴𝑥;𝜏=𝜀𝐵𝐴𝑡, (2.2) and wrote (2.1) in the dimensionless form: ̈𝑥𝑥+³𝑥²(0)+1=0,𝑥=𝑥₁₀,𝑥̇(0)=0, (2.3) then, he obtained the first approximate solution to (2.3) based on the method of harmonic balance and assumed that the exact angular frequency for the periodic solution of (2.3) can be determined from 𝜔²_𝑀𝑥₁₀=𝜙²𝑥²₁₀1+𝜙²𝑥²₁₀, (2.4) where 𝜋𝜙≡2𝐹𝜋/2,1/√2. (2.5) Here, 𝐹(𝜋/2,1/√2) is the complete elliptic integral of the first kind.

By following Ritz procedure, Tiwari et al. in [5] obtained the approximate angular frequency of (2.3): 𝜔²_𝑇𝑥₁₀2=1+𝑥²₁₀11+𝑥²₁₀−1, (2.6) which satisfies the limits: 𝑥₁₀issmall:𝜔𝑥₁₀²=34𝑥²₁₀+𝑂𝑥⁴₁₀;𝑥₁₀large:𝜔𝑥₁₀²1=1+𝑂𝑥²₁₀. (2.7)

Based on the homotopy method and by only considering the first approximation, He obtained the angular frequency of (2.3) which is given as [8]: 𝜔²_𝐻𝑥₁₀=34𝑥²₀31+4𝑥²₁₀⁻¹. (2.8)

Öziş and Yildirim used He's energy balance method to obtain the angular frequency of the Duffing-harmonic oscillator by writing (2.3) in its variational representation to construct its Hamiltonian form [8]. Then, they used a trail function and found that 𝜔²_𝑂𝑌𝑥₁₀2=1−𝑥²₁₀ln1+𝑥²₁₀1+𝑥²₁₀,/2 (2.9) which provides an approximate angular frequency expression to (2.3).

Furthermore, Lim et al. [7] introduced a new method for solving the Duffing-harmonic oscillator by combining Newton's method with the harmonic balance method and obtained a third-order approximation to the angular frequency value which is given by𝜔₃𝑥₁₀=Ω₃𝑥₁₀, (2.10) where Ω₃𝑥₁₀=𝐶𝑥₁₀𝐷𝑥₁₀. (2.11) The expressions for 𝐶(𝑥₁₀) and 𝐷(𝑥₁₀) are defined in [7].

On the other hand, Beléndez and coworkers in [12] used a cubication method and Chebyshev polynomials to obtain an equivalent expression for (2.3) in the form: 𝑥̈+𝛼𝑥₁₀𝑥+𝛽𝑥₁₀𝑥³≈0, (2.12) and solved (2.12) to obtain its displacement expression given by 𝑥(𝑡)=𝑥₁₀cn𝜔_𝐵𝑡,𝑘²_𝐵, (2.13) where the modulus 𝑘_𝐵 and the elliptic frequency of oscillation 𝜔_𝐵 are given by 𝑘_𝐵=𝑥₁₀𝛽2𝛼+2𝛽𝑥²₁₀;𝜔_𝐵=𝛼+𝛽𝑥²₁₀, (2.14) where 𝛼𝑥₁₀4=1+𝑥²₁₀51+𝑥²₁₀−2+24𝑥⁴₁₀11+𝑥²₁₀,𝛽𝑥−1₁₀=8𝑥⁴₁₀31−1+𝑥²₁₀+32𝑥⁶₁₀11−1+𝑥²₁₀. (2.15) Thus, the approximate circular frequency Ω_𝐵 of the Duffing-harmonic oscillator (2.3) is given by Ω_𝐵=𝜋𝜔_𝐵2𝐾𝑘²_𝐵. (2.16) Beléndez et al. showed that the maximum error attained by using (2.16) when compared to the exact one is not bigger than 0.071% [12].

By using Jacobi elliptic functions, Elías-Zùñiga et al. obtained an approximate expression to find the angular circular frequency of the nonhomogeneous representation of (2.1) by using the rational second-order Jacobi elliptic form solution [21]. Here, we have followed this approach to obtain the corresponding circular frequency Ω_𝐸𝐵 of (2.1) which has the form: Ω_𝐸𝐵=𝜋𝜔_𝐸𝐵2𝐾𝑘²_𝐸𝐵, (2.17) in which 𝑘_𝐸𝐵=𝐻₁𝐻₂;𝜔_𝐸𝐵=𝐻₃𝐻₄, (2.18) where 𝐻₁, 𝐻₂, 𝐻₃, and 𝐻₄ are given in the appendix.

In the next section, we shall enhance the cubication procedure used by Beléndez et al. in [12] by combining this approach with the equivalent nonlinearization method introduced by Cai and coworkers in [14], to derive an approximate angular frequency expression which is closer to its numerical integration solution value.

3. Solution Procedure

Here, we focus on the derivation of an approximate solution of the nonlinear Duffing-harmonic oscillator by combining the cubication and the nonlinearization methods. First, we write the restoring force in equivalent representation form that takes into account at least three terms of its Chebyshev polynomial expansion.

Let us consider that the nonlinear Duffing-harmonic oscillator has the form: 𝑑²𝑥𝑑𝑡²+𝑥³1+𝑥²=0, (3.1) with the initial conditions: 𝑥(0)=1+𝑥₁₀,𝑥̇(0)=0. (3.2) We next introduce the following change of variable into (3.1): 𝑥𝑦=𝑞,where𝑞=1+𝑥₁₀, (3.3) to get 𝑑²𝑦𝑑𝑡²𝑦𝑦=𝑞+𝑓=0,𝑓²𝑦³1+𝑞²𝑦². (3.4) Therefore, the initial conditions of (3.4) become 𝑦(0)=1,𝑑𝑦(0)𝑑𝑡=0. (3.5) Now, we follow the cubication procedure [12] and write the restoring force 𝑓(𝑦) in (3.4) as a function of the Chebyshev polynomial expansion: 𝑓𝑦=𝑞²𝑦³1+𝑞²𝑦²≈^𝑁_𝑛=0𝑏_2𝑛+1𝑥₁₀𝑇_2𝑛+1𝑦, (3.6) where the first three polynomials are 𝑇₁=𝑦;𝑇₃=4𝑦³−3𝑦;𝑇₅=16𝑦⁵−20𝑦³+5𝑦. (3.7) Notice that in [12], the restoring force was replaced by an equivalent form by using only two terms of the Chebyshev polynomial expansion to ensure a polynomial cubic equation. However, we consider three terms, that is, 𝑁=2 in (3.6), to replace the rational restoring force by a fifth-order polynomial equation: 𝑓𝑦=𝑞²𝑦³1+𝑞²𝑦²=𝑏₁𝑞𝑇₁𝑦+𝑏₃𝑞𝑇₃𝑦+𝑏₅𝑞𝑇₅𝑦=𝑏₁𝑞−3𝑏₃𝑞+5𝑏₅𝑞𝑦+4𝑏₃𝑞−20𝑏₅𝑞𝑦³+16𝑏₅𝑞𝑦⁵, (3.8) where 𝑏_2𝑛+1𝑞=2𝜋¹₋₁1−𝑦^2−1/2𝑓𝑦𝑇_2𝑛+1𝑦𝑑𝑦. (3.9) Then, from (3.9) we obtain 𝑏₁𝑞2=1−𝑞²+2𝑞²1+𝑞²,𝑏₃𝑞=2𝑞²+8𝑞⁴−2𝑞⁴1+𝑞²−6𝑞⁴1+𝑞²,𝑏₅𝑞2=−𝑞⁶16+𝑞⁴−11+𝑞²−151+𝑞²+𝑞²12−51+𝑞². (3.10) Therefore, we may write (3.4) as an equivalent cubic-quintic Duffing oscillator as follows: 𝑑²𝑦𝑑𝑡²+𝛼𝑞𝑦+𝛽𝑞𝑦³+𝛾𝑞𝑦⁵≈0, (3.11) in which 𝛼𝑞1=−𝑞⁶160+18𝑞⁴−𝑞⁶−61+𝑞²−1541+𝑞²+𝑞²144−701+𝑞²,𝛽𝑞=16𝑞⁶40+3𝑞⁴−21+𝑞²−381+𝑞²+𝑞²32−141+𝑞²,𝛾𝑞=−32𝑞⁶1+𝑞²𝑞⁴1+𝑞²−5+161+𝑞²−1+4𝑞²31+𝑞²−5. (3.12)

Since the cubication procedure requires a cubic polynomial representation form of the restoring force, we shall next transform the cubic-quintic restoring force of (3.11) into a cubic polynomial restoring force: 𝐹𝑦=𝛼𝑞𝑦+𝛽𝑞𝑦³+𝛾𝑞𝑦⁵≡𝛿𝑞𝑦+𝜖𝑞𝑦³ (3.13) and use the equivalent nonlinearization method [14] to find 𝛿(𝑞) and 𝜀(𝑞). In this method, we seek a polynomial of the form 𝛿(𝑞)𝑦+𝜖(𝑞)𝑦³ satisfying 𝐹(𝛿,𝜖)=^𝜎₀𝛼𝑞𝑦+𝛽𝑞𝑦³+𝛾𝑞𝑦⁵−𝛿𝑞𝑦−𝜖𝑞𝑦³²𝑑𝑦→min, (3.14) which requires 𝜕𝐹/𝜕𝛿(𝑞)=0 and 𝜕𝐹/𝜕𝜖(𝑞)=0. Notice that the value of 𝜎 is fitted to satisfy (3.14). Then, the equations that provide the coefficients 𝛿(𝑞) and 𝜖(𝑞) are 𝛿𝑞=1𝑞𝑞𝜎,𝜖𝑞=12121𝛼−5𝛾99𝛽𝑞+10𝛾𝑞𝜎², (3.15) with 𝜎=1.0436. Thus, the approximate equivalent cubic representation of (3.11) is given by 𝑑²𝑦𝑑𝑡²+𝛿𝑞𝑦+𝜖𝑞𝑦³≈0. (3.16) It is well known that (3.16) has an exact angular frequency-amplitude relationship given by Ω_𝐶=𝜋𝜔_𝐶2𝐾𝑘²_𝐶, (3.17) where 𝜔_𝐶=√𝑘𝛿+𝜖,_𝐶=𝜖2().𝛿+𝜖 (3.18)

In the next section, we shall compare our derived approximate circular frequency Ω_𝐶 given by (3.17) with respect to the numerical one and with some other approximate solutions.

4. Numerical Simulations

In this section, we compare the numerical integration solution of the angular frequency value of (3.1) [19] with the approximate solution given by (3.17) and with approximate solutions obtained by other methods.

First, let us begin by plotting the relative errors attain by using the angular frequency relations derived by Tiwari et al. [5], Mickens [2], Öziş and Yildirim [8], and He [8], since these are of the same order of magnitude. Figure 1 provides a comparison of the relative errors plotted against the oscillation amplitudes, 𝑥₀. As we may see from Figure 1, the maximum relative error values are 1.07%, 1.73%, 2.33%, and 2.81%, respectively.

Figure 1 

Relative errors for approximate frequency values obtained from Mickens [2], Tiwari et al. [5], He [8], and Öziş and Yildirim [8] derived solutions.

We next plot in Figure 2 the estimated relative error values obtained from the solutions derived by Lim et al. [7], Beléndez et al. [12], Elías-Zùñiga et al. [21], and the solution given by (3.17). In this Figure 2, we may see that the maximum relative error attained by Lim et al. approximate solution given by (2.11) is 0.1184%, which corresponds to an oscillation amplitude value of 𝑥₀=3.46. Besides, Figure 2 shows that for amplitude of oscillations 𝑥₀≥9, the relative errors attained by the combined Newton's and harmonic balance solution tend to approach to the relative error values of the elliptic balance solution. We may notice from Figure 2, that the elliptic balance solution relative error approaches to its maximum value of 0.082% at 𝑥₀=5.98. However, the maximum relative error attained by combining the cubication and nonlinearization methods is 0.055%, which is 21% lower than the relative error of 0.071% computed from (2.16) [12].

Figure 2 

Relative errors for approximate frequency values obtained from Lim et al. [7], Beléndez et al. [12], Elías-Zúñiga et al. [21] solutions and by our derived cubication-equivalent nonlinearization solution given by (3.17) with the parameter value of 𝜎=1.03915.

Figure 3 shows the corresponding relative error curves of the cubication, the nonlinearization, and the combined cubication and nonlinearization solutions plotted versus the amplitude of oscillations. We may see from Figure 3 that the nonlinearization solution gives the highest relative error value of 0.106% at 𝑥₀=0.8 with 𝜎=1.0457. However, it is clear from Figure 3 that our proposed combination procedure of the cubication and the nonlinearization methods, as described in Section 3, provides the smallest relative error value. Therefore, we may conclude that our approach enhances the precision of the cubication method [12] in which the restoring force is expanded as a function of Chebyshev polynomials.

Figure 3 

Relative errors for approximate frequency values obtained from the cubication, the nonlinearization, and the enhanced cubication methods.

5. Conclusions

In this paper, we have used the frequency-amplitude relationship Ω_𝐶 obtained by combining the cubication procedure [12] with the equivalent nonlinearization method for strongly nonlinear oscillators [14]. The combination of these methods, to solve strongly nonlinear oscillators of the Duffing-type, that provide a cubic-type polynomial expression for the restoring force is an expected result given the fact that nonlinear restoring forces can be equivalently represented by polynomial expressions of degrees one, two, or three [12, 15–17]. With our proposed approach, the original nonlinear equation of motion can also be replaced with the homogenous Duffing equation that has a well-known exact solution that depends on Jacobi elliptic functions. It is clear from the relative error values obtained by the expression given by (3.17), that our solution provides the best estimate to the angular frequency-amplitude value than those developed previously by other authors.

In a forthcoming paper, we will show the potential of combining the cubication and the equivalent nonlinearization methods to solve other strongly nonlinear oscillators of the Duffing-type.

Appendix

The expressions of 𝐻₁ through 𝐻₄ to compute 𝑘_𝐸𝐵 and 𝜔_𝐸𝐵 values in (2.18) are given by𝐻₁=16𝐴𝑎²₃𝑏₃−320+7𝑏₃−80+𝑏₃−28+5𝑏₃+𝐵−2048+𝑏₃−17408+𝑏₃−34688+𝑏₃−27648+7𝑏₃−1312+3𝑏₃−40+7𝑏₃,𝐻₂=28𝐴𝑎²₃1+𝑏₃−80+7𝑏₃−32+𝑏₃−13+5𝑏₃+𝐵−2048+𝑏₃−9344+𝑏₃−13952+𝑏₃−8280+7𝑏₃−184+21𝑏₃3+𝑏₃,𝐻₃=𝑎²₃8𝐴𝑎²₃80+𝑏₃304+7𝑏₃45+8−5𝑏₃𝑏₃+𝐵2048+𝑏₃9344+𝑏₃13952+𝑏₃8280+7𝑏₃−184+21𝑏₃3+𝑏₃𝜖,𝐻₄=8𝐴²𝑎⁴₃80+𝑏₃144+𝑏₃−61+7𝑏₃−18+5𝑏₃+4𝐴𝑎²₃𝐵640+𝑏₃1984+𝑏₃1088+𝑏₃−1262+7𝑏₃−141+𝑏₃15𝑏₃−8+𝐵²(2048+𝑏₃9728+𝑏₃14976+𝑏₃4480+𝑏₃−8416+7𝑏₃(−1020+𝑏₃−163+3𝑏₃22+7𝑏₃, (A.1) where 𝑎₃=1+𝑏₃𝑦₁₀, (A.2) and the constant 𝑏₃, that satisfies the condition |𝑎₃|>|𝑏₃| [2], is determined by solving the tenth-order polynomial equation as follows:−21𝐴𝐵𝑦²₁₀𝑏₃¹⁰+𝑏⁹₃63𝐵²+210𝐴𝐵𝑦²₁₀+𝑏⁸₃1134𝐵²+3906𝐴𝐵𝑦²₁₀+2688𝐴²𝑦⁴₁₀+𝑏⁷₃8559𝐵²+22596𝐴𝐵𝑦²₁₀+16128𝐴²𝑦⁴₁₀+𝑏⁶₃35076𝐵²+74779𝐴𝐵𝑦²₁₀+40320𝐴²𝑦⁴₁₀+𝑏⁵₃82296𝐵²+145274𝐴𝐵𝑦²₁₀+53760𝐴²𝑦⁴₁₀+𝑏⁴₃111168𝐵²+159592𝐴𝐵𝑦²₁₀+40320𝐴²𝑦⁴₁₀+𝑏³₃82560𝐵²+90528𝐴𝐵𝑦²₁₀+16128𝐴²𝑦⁴₁₀+𝑏²₃29184𝐵²+18560𝐴𝐵𝑦²₁₀+2688𝐴²𝑦⁴₁₀+𝑏₃3072𝐵²−2816𝐴𝐵𝑦²₁₀−1024𝐴𝐵𝑦²₁₀=0. (A.3)

Acknowledgment

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