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

Authors: Alex Elías-Zúñiga; Oscar Martínez-Romero; René K. Córdoba-Díaz
Journal: Mathematical Problems in Engineering
 (2012)

---

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


---

## Body

## 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 f(x) by an equivalent cubic polynomial expression a3x3, where the value of a3 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:(1.1)x¨+f(x)=0;f(x)=εx3(B+Ax2)
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. Ö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:
(2.1)d2ydτ2+εy3(Ay2+B)=0,y(0)=y10,y˙(0)=0,
for which the parameters ε, A, and B are nonnegative. He used the following transformations:
(2.2)y=BAx;τ=(εBA)t,
and wrote (2.1) in the dimensionless form:
(2.3)x¨+x3(x2+1)=0,x(0)=x10,x˙(0)=0,
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
(2.4)ωM2(x10)=ϕ2x1021+ϕ2x102,
where
(2.5)ϕ≡π2F(π/2,1/2).
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):
(2.6)ωT2(x10)=1+(2x102){11+x102-1},
which satisfies the limits:
(2.7)x10issmall:ω(x10)2=34x102+O(x104);x10large:ω(x10)2=1+O(1x102).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)ωH2(x10)=34x02(1+34x102)-1.Ö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.9)ωOY2(x10)=1-2x102ln(1+x1021+x102/2),
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)ω3(x10)=Ω3(x10),
where
(2.11)Ω3(x10)=C(x10)D(x10).
The expressions for C(x10) and D(x10) 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:
(2.12)x¨+α(x10)x+β(x10)x3≈0,
and solved (2.12) to obtain its displacement expression given by
(2.13)x(t)=x10cn(ωBt,kB2),
where the modulus kB and the elliptic frequency of oscillation ωB are given by
(2.14)kB=x10β2α+2βx102;ωB=α+βx102,
where
(2.15)α(x10)=1+4x102(51+x102-2)+24x104(11+x102-1),β(x10)=8x104(1-31+x102)+32x106(1-11+x102).
Thus, the approximate circular frequency ΩB of the Duffing-harmonic oscillator (2.3) is given by
(2.16)ΩB=πωB2K(kB2).
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, Eli´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:
(2.17)ΩEB=πωEB2K(kEB2),
in which
(2.18)kEB=H1H2;ωEB=H3H4,
where H1, H2, H3, and H4 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:(3.1)d2xdt2+x3(1+x2)=0,
with the initial conditions:
(3.2)x(0)=1+x10,x˙(0)=0.
We next introduce the following change of variable into (3.1):
(3.3)y=xq,whereq=1+x10,
to get
(3.4)d2ydt2+f(y)=0,f(y)=(q2y31+q2y2).
Therefore, the initial conditions of (3.4) become
(3.5)y(0)=1,dydt(0)=0.
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:
(3.6)f(y)=(q2y31+q2y2)≈∑n=0Nb2n+1(x10)T2n+1(y),
where the first three polynomials are
(3.7)T1=y;T3=4y3-3y;T5=16y5-20y3+5y.
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:
(3.8)f(y)=q2y31+q2y2=b1(q)T1(y)+b3(q)T3(y)+b5(q)T5(y)=[b1(q)-3b3(q)+5b5(q)]y+[4b3(q)-20b5(q)]y3+16b5(q)y5,
where
(3.9)b2n+1(q)=2π∫-11(1-y2)-1/2f(y)T2n+1(y)dy.
Then, from (3.9) we obtain
(3.10)b1(q)=1-2q2+2q21+q2,b3(q)=2q2+8q4-2q41+q2-6q41+q2,b5(q)=-2q6(16+q4-11+q2-151+q2+q2(12-51+q2)).
Therefore, we may write (3.4) as an equivalent cubic-quintic Duffing oscillator as follows:
(3.11)d2ydt2+α(q)y+β(q)y3+γ(q)y5≈0,
in which
(3.12)α(q)=-1q6(160+18q4-q6-61+q2-1541+q2+q2(144-701+q2)),β(q)=16q6(40+3q4-21+q2-381+q2+q2(32-141+q2)),γ(q)=-32q61+q2(q4(1+q2-5)+16(1+q2-1)+4q2(31+q2-5)).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)F(y)=α(q)y+β(q)y3+γ(q)y5≡δ(q)y+ϵ(q)y3
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)y3 satisfying
(3.14)F(δ,ϵ)=∫0σ(α(q)y+β(q)y3+γ(q)y5-δ(q)y-ϵ(q)y3)2dy→min,
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
(3.15)δ(q)=121(21α(q)-5γ(q)σ),ϵ(q)=19(9β(q)+10γ(q)σ2),
with σ=1.0436. Thus, the approximate equivalent cubic representation of (3.11) is given by
(3.16)d2ydt2+δ(q)y+ϵ(q)y3≈0.
It is well known that (3.16) has an exact angular frequency-amplitude relationship given by
(3.17)ΩC=πωC2K(kC2),
where
(3.18)ωC=δ+ϵ,kC=ϵ2(δ+ϵ).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.


## 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, x0. 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 Figure2 the estimated relative error values obtained from the solutions derived by Lim et al. [7], Beléndez et al. [12], Eli´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 x0=3.46. Besides, Figure 2 shows that for amplitude of oscillations x0≥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 x0=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.Figure3 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 x0=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Ω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.



---

*Source: 618750-2012-09-24.xml*