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

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.


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 f x by an equivalent cubic polynomial expression a 3 x 3 , where the value of a 3 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 .

Mathematical Problems in Engineering
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: and show that it has the smallest relative error value when compare to the numerical integration and with others approximate solutions.In 1.1 , x represents system displacement, and ε, A, and B 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, A λ, and B 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, A λ, and B 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.Özis ¸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.

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: for which the parameters ε, A, and B are nonnegative.He used the following transformations: and wrote 2.1 in the dimensionless form: 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 where Here, F π/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 : which satisfies the limits: x 10 is small: ω x 10 .

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 :

2.8
Özis ¸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 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 where The expressions for C x 10 and D x 10 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: ẍ α x 10 x β x 10 x 3 ≈ 0, 2.12 and solved 2.12 to obtain its displacement expression given by

2.15
Thus, the approximate circular frequency Ω B of the Duffing-harmonic oscillator 2.3 is given by

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 Ω EB of 2.1 which has the form: where H 1 , H 2 , H 3 , and H 4 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.

Solution Procedure
Here, we focus on the derivation of an approximate solution of the nonlinear Duffingharmonic 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: with the initial conditions: x 0 1 x 10 , ẋ 0 0.

3.2
We next introduce the following change of variable into 3.1 : , where q 1 x 10 , 3.3 to get 3.4 Therefore, the initial conditions of 3.4 become y 0 1, dy dt 0 0.

3.5
Now, we follow the cubication procedure 12 and write the restoring force f y in 3.4 as a function of the Chebyshev polynomial expansion: where the first three polynomials are T 1 y; T 3 4y 3 − 3y; T 5 16y 5 − 20y 3 5y.

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, N 2 in 3.6 , to replace the rational restoring force by a fifth-order polynomial equation: f y q 2 y 3 1 q 2 y 2 b 1 q T 1 y b 3 q T 3 y b 5 q T 5 y b 1 q − 3b 3 q 5b 5 q y 4b 3 q − 20b 5 q y 3 16b 5 q y 5 ,

3.8
where f y T 2n 1 y dy.

3.9
Mathematical Problems in Engineering 7 Then, from 3.9 we obtain

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: F y α q y β q y 3 γ q y 5 ≡ δ q y q y 3 3.13 and use the equivalent nonlinearization method 14 to find δ q and ε q .In this method, we seek a polynomial of the form δ q y q y 3 satisfying F δ, σ 0 α q y β q y 3 γ q y 5 − δ q y − q y 3 2 dy −→ min, 3.14 Mathematical Problems in Engineering which requires ∂F/∂δ q 0 and ∂F/∂ q 0. Notice that the value of σ is fitted to satisfy 3.14 .Then, the equations that provide the coefficients δ q and q are δ q 1 21 21α q − 5γ q σ , q 1 9 9β q 10γ q σ 2 ,

3.15
with σ 1.0436.Thus, the approximate equivalent cubic representation of 3.11 is given by d 2 y dt 2 δ q y q y 3 ≈ 0.

3.16
It is well known that 3.16 has an exact angular frequency-amplitude relationship given by 3.17 where

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

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 , Özis ¸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, x 0 .As we may see from Figure 1, the maximum relative error values are 1.07%, 1.73%, 2.33%, and 2.81%, respectively.
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 x 0 3.46.Besides, Figure 2 shows that for amplitude of oscillations x 0 ≥ 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 x 0 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 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 x 0 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.

Conclusions
In this paper, we have used the frequency-amplitude relationship Ω C 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.A.3

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