# Analytical Approximate Solutions for the Cubic-Quintic Duffing Oscillator in Terms of Elementary Functions

Authors: A. Beléndez; M. L. Alvarez; J. Francés; S. Bleda; T. Beléndez; A. Nájera; E. Arribas
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

Accurate approximate closed-form solutions for the cubic-quintic Duffing oscillator are obtained in terms of elementary functions. To do this, we use the previous results obtained using a cubication method in which the restoring force is expanded in Chebyshev polynomials and the original nonlinear differential equation is approximated by a cubic Duffing equation. Explicit approximate solutions are then expressed as a function of the complete elliptic integral of the first kind and the Jacobi elliptic function cn. Then we obtain other approximate expressions for these solutions, which are expressed in terms of elementary functions. To do this, the relationship between the complete elliptic integral of the first kind and the arithmetic-geometric mean is used and the rational harmonic balance method is applied to obtain the periodic solution of the original nonlinear oscillator.


---

## Body

## 1. Introduction

The cubic-quintic Duffing equation is a differential equation with third- and/or fifth-power nonlinearity [1]. Due to the presence of this fifth-power nonlinearity added to the third nonlinearity of the common Duffing equation, this oscillator is difficult to handle. As Lai et al. pointed out [1], the cubic-quintic Duffing equation can be found in the modelling of different systems and we advise readers to consult this paper to learn more about these applications and a good set of references in which this strongly nonlinear equation can be found. Among these systems modelled by the cubic-quintic Duffing equation we would like to mention the nonlinear dynamics of a slender elastica, the compound KdV equation in nonlinear wave systems, or the propagation of a short electromagnetic pulse in a nonlinear medium [1–3].The purpose of this paper is to present analytical approximations to the periodic solutions to the cubic-quintic Duffing oscillator [1–3] in terms of elementary functions. To do this, we use the results previously obtained [4] using the Chebyshev series expansion of the restoring force [5, 6] for a quintic Duffing oscillator. Using this approximation, the original nonlinear differential equation is replaced by a cubic Duffing equation, which can be exactly solved. The replacement of the original nonlinear equation by an “approximate cubic Duffing equation” allows us to obtain an approximate frequency-amplitude relation as a function of the complete elliptic integral of the first kind and the solution in terms of the Jacobi elliptic function cn [4]. Using these first approximate equations, closed-form expressions for the approximate frequency and the solution in terms of elementary functions are obtained using the arithmetic-geometric mean combined with the (second-order) rational harmonic balance method.In relation to cubication methods, Bravo-Yuste [7] introduced a cubication procedure based on the harmonic balance method and Bravo-Yuste and Martín-Sánchez [8] also developed a weighted mean square method of cubication. Recently, Beléndez et al. [4, 9, 10] and Elías-Zùñiga et al. [11] employed similar procedures based on the use of Chebyshev polynomials. In particular, Elías-Zùñiga et al. combined the cubication method with the equivalent nonlinearization method developed by Cai et al. [12].In the early 1980’s and 1990’s, much effort was put on renormalization techniques and the use of the well-known arithmetic-geometric mean to analyze nonlinear oscillators and, in particular, the old famous simple pendulum [13]. Recently, and also analyzing the nonlinear pendulum, it has been shown that it is possible to obtain an approximate expression for the complete elliptic integral of the first kind by using a Taylor series expansion [14], which is similar to the second approximation derived using the arithmetic-geometric mean. Finally, the rational harmonic balance method (RHBM) [15] is more difficult to apply than the common harmonic balance method (HBM) for oscillators with complex nonlinearities. However, we will see that the RHBM is easy to apply to the cubic Duffing equation obtained when the Chebyshev polynomials expansion is taken into account. The RHBM presents the advantage versus the usual HBM that the RHBM includes approximations for all the harmonics contributing to the periodic solution [16].


## 2. Solution Procedure Using a Chebyshev Polynomials Expansion

A cubic-quintic Duffing oscillator is a conservative autonomous oscillatory system, which can be described by the following second-order differential equation with cubic-quintic nonlinearities [1]:
(2.1)d2xdt2+λ1x+λ3x3+λ5x5=0
with initial conditions
(2.2)x(0)=A,dxdt(0)=0,
where λ1x+λ3x3+λ5x5 is an odd function; x and t are dimensionless displacement and time variables, respectively; λ1, λ3, and λ5 are positive constant parameters [1]. The angular frequency of these oscillations is denoted by ω, which is a function of the initial amplitude A.The cubic-quintic Duffing equation is difficult to handle because of the presence of strong nonlinearity and it has no known closed-form solution. However, the quintic Duffing equation (λ1=1, λ3=0, and λ5=ε≠0 in (2.1)) has been recently solved approximately by means of a cubication procedure based on the expansion of the nonlinear restoring force in terms of the Chebyshev polynomials [4]. To apply this procedure, firstly a reduced variable y=x/A is introduced in (2.1) and (2.2) and the nonlinear function g(y)=λ1y+λ3A2y3+λ5A4y5 is expanded in terms of Chebyshev polynomials of the first kind Tn(y). Retaining only the first two terms in this expansion, we obtain
(2.3)g(y)=∑n=0∞b2n+1T2n+1(y)≈∑n=01b2n+1T2n+1(y)=b1T1(y)+b3T3(y)=(b1-3b3)y+4b3y3=αy+βy3,
where
(2.4)b2n+1=2π∫-1+111-y2g(y)T2n+1(y)dy,b1=λ1+34λ3A2+58λ5A4,b3=14λ3A2+516λ5A4.
Then the nonlinear differential equation in (2.1) can be then approximated by the (cubic-) Duffing equation
(2.5)d2ydt2+α(A)y+β(A)y3=0,
where
(2.6)α=λ1-516λ5A4,β=λ3A2+54λ5A4.Equation (2.5) can be exactly solved and its solution can be written in terms of the Jacobian elliptic function cn [15]. Then the approximate frequency and solution for the cubic-quintic equation in (2.1) will be the exact frequency and solution for the cubic equation in (2.5). Taking into account the well-known expressions for the frequency and the solution for the Duffing oscillator in (2.5), we can easily obtain the following analytical approximate expressions for the frequency and solution for the cubic-quintic Duffing oscillator as follows:
(2.7)ωa(A)=πλ1+λ3A2+(15/16)λ5A42K((8λ3A2+10λ5A4)/(16λ1+16λ3A2+15λ5A4)),(2.8)y(t)=cn(tλ1+λ3A2+1516λ5A4;8λ3A2+10λ5A416λ1+16λ3A2+15λ5A4),
where K(m) is the complete elliptic integral of the first kind
(2.9)K(m)=∫0π/2dθ1-msin2θ.


## 3. Closed-Form Expression for the Approximate Frequency in Terms of Elementary Functions Using the Arithmetic-Geometric Mean

As we can see from (2.7), the approximate frequency for the cubic-quintic Duffing oscillator is expressed in terms of the complete elliptic integral of the first kind, K(m). However, it is also possible to obtain a closed-form expression for the approximate frequency in terms of elementary functions. The complete elliptic integral of the first kind (2.9) cannot be expressed in terms of elementary functions, but can be numerically evaluated with high precision by a simple procedure based on the arithmetic-geometric mean. This is due to the fact that the arithmetic-geometric mean is the basis of Gauss’ method for the calculation of elliptic integrals [16–18]. The arithmetic-geometric mean of two numbers p and q, M(p,q) can be expressed in closed-form in terms of the complete elliptic integral of the first kind as [16]
(3.1)M(p,q)=(p+q)π4K((p-q)/(p+q))
which allows us to write K(m) in (2.9) as follows [16]:
(3.2)K(m)=π2M(1,1-m).The Legendre form of the arithmetic-geometric mean is given by [16]
(3.3)M(1,z)=∏n=0∞12(1+kn),
where z≡k0 and kn+1=2kn/(1+kn). From (3.3) it is possible to approximately calculate the arithmetic-geometric mean by choosing a finite number of terms in this equation as follows:
(3.4)M(1,z)≈∏n=0N12(1+kn).
From (2.7), (3.2), and (3.4) we obtain the following expression for the approximate frequency of the cubic-quintic Duffing oscillator:
(3.5)ωaagm(A)=116[(16λ1+8λ3A2+5λ5A4)1/4+(16λ1+16λ3A2+15λ5A4)1/4]2,
where only two terms (N=1) have been considered in (3.4)
(3.6)M(1,z)≈∏n=0112(1+kn)=14(1+k0)(1+k1)=14(1+k0)(1+2k01+k0)=14(1+z+2z)=14(1+z)2,
where z≡k0=1-m and m=(8λ3A2+10λ5A4)/(16λ1+16λ3A2+15λ5A4) (2.7). The maximum value of m reached for the cubic-quintic Duffing oscillator is limλ5A4→∞m=10/15<2/3 because for the case in which λ3≠0 and λ5=0 we obtain limλ3A2→∞m=8/16=1/2<2/3. We can see that (3.5) is a very simple and easily computable expression.From (3.2) and (3.6) it can be seen that it is possible to take into account the following approximation for the complete elliptic integral of the first kind:
(3.7)Ka(m)=2π[1+(1-m)1/4]2.
In Table 1 we present the percentage errors for comparison of the approximate expression for the complete elliptic integral of the first kind, Ka(m) (3.7), with the exact one, K(m), for the interval of values of m corresponding to the cubic-quintic Duffing oscillator (0≤m≤2/3). As we can see from this table, Ka(m) is a very good approximation for K(m) and the relative error is less than 0.087% for m≤2/3.Table 1
Percentage errors for comparison of the approximate expression for the complete elliptic integral of the first kind,Ka(m) (3.7), with the exact one, K(m).
m
K
a
(
m
)
K
(
m
)
(% Error)
0
1.570796
1.570796
0
0.1
1.612441
1.612441
0.00000075
0.2
1.659623
1.659624
0.000015
0.3
1.713888
1.713889
0.00010
0.4
1.777512
1.777519
0.00041
0.5
1.854049
1.854075
0.0014
0.6
1.949485
1.949568
0.0042
2/3
2.028783
2.028959
0.0087To verify if it is sufficient to consider only a two-term approximation in the Legendre form of the arithmetic-geometric mean, we compareωa(A) (2.7) and ωaagm(A) (3.5) for λ5≠0 and for λ3≠0 and λ5=0.Forλ5≠0 we obtain the following power series expansions for these equations for very large values of the oscillation amplitude A(3.8)ωa(A)≈π15λ58K(2/3)A2=0.749605λ5A2,(3.9)ωaagm(A)≈(1+31/4)25λ516A2=0.749670λ5A2.Furthermore, we have the following:(3.10)limλ5A4→0ωaagm(A)ωa(A)=1,limλ5A4→∞ωaagm(A)ωa(A)=1.000087.
The maximum relative error for the frequency using the arithmetic-geometric mean approximation in (3.5) instead of (2.7) is 0.0087%, which allows us to conclude that (3.6) is a suitable approximation to (3.3).Forλ3≠0 and λ5=0 we obtain the following power series expansions for these equations for very large values of the oscillation amplitude A:
(3.11)ωa(A)≈2πΓ(3/4)λ3Γ(1/4)A=0.847213λ3A,ωaagm(λ)≈(2+21/2+27/4)2λ38A=0.847225λ3A,
where Γ(z) is the Euler gamma function. Furthermore, we have the following equations:
(3.12)limλ3A2→0ωaagm(A)ωa(A)=1,limλ3A2→∞ωaagm(A)ωa(A)=1.000014.
Now the maximum relative error for the frequency using the arithmetic-geometric mean approximation in (3.5) instead of (2.7) is as low as 0.0014%, which allows us to conclude that (3.6) is a suitable approximation to (3.3).


## 4. Closed-Form Expression for the Solution in Terms of Elementary Functions Using the Rational Harmonic Balance Method

In order to obtain an expression for the approximate solution in terms of elementary functions we consider that the approximate frequency is given by (3.5) and we apply the rational harmonic balance method (RHBM). By defining a new independent variable replacing the time variable, τ=ωaagmt, (2.5) can be expressed as follows:
(4.1)ω2d2ydτ2+αy+βy3=0,
where we have written ω=ωaagm for simplicity. We propose the following approximation to the periodic solution of (4.1)
(4.2)y(τ)=(1-q)cosτ1-qcos2τ,0≤|q|<1.
Note that this expression satisfies the prescribed initial conditions independently of parameter q, which is to be determined. It is easy to verify that substituting q=-2B2/(1-B2) in (4.2), we obtain the rational representation
(4.3)y(τ)=(1+B2)cosτ1+B2cos(2τ),
which is the expression usually considered in the bibliography for the second-order rational harmonic balance method [15]. However, (4.2) has the advantage that it allows us to simplify the procedure to obtain the coefficients of the Fourier series expansion of y and y3. These series can be written as follows:
(4.4)yσ(τ)=∑n=1∞c2n+1(σ)cos[(2n+1)τ],σ=1,3,
where
(4.5)c2n+1(σ)=4π∫0π/2yσ(τ)cos[(2n+1)τ]dτ.
Substituting (4.4) into (4.1) the following result is obtained:
(4.6)(-ω2c1(1)+αc1(1)+βc1(3))cosτ+HOH=0,
where HOH stands for higher-order harmonics. Setting the coefficient of the lowest harmonic to zero gives
(4.7)(-ω2+α)c1(1)+βc1(3)=0,
which can be solved to obtain q as a function of A. Before doing this, it is necessary to obtain the values of the first Fourier coefficients c1(1) and c1(3) as a function of the unknown parameter q. To do this we firstly write the expression for yσ(τ) as follows
(4.8)yσ(τ)=(1-q)σcosστ(1-qcos2τ)-σ=(1-q)σcosστ∑m=1∞(-σm)(-1)mqmcos2mτ=(1-q)σ∑m=0∞(σ)mm!qmcos2m+στ=(1-q)σ∑m=0∞Γ(σ+m)Γ(σ)Γ(m+1)qmcos2m+στ
and then
(4.9)y(τ)=(1-q)∑m=0∞qmcos2m+1τ,y3(τ)=12(1-q)3∑m=0∞(m+1)(m+2)qmcos2m+3τ.
Then Fourier coefficients c1(1) and c1(3) take the expressions
(4.10)c1(1)=4π∫0π/2y(τ)cosτdτ=4π(1-q)∑m=0∞qm∫0π/2cos2m+2τdτ=4π(1-q)∑m=0∞πΓ(m+3/2)qm2Γ(m+2)=2q(1-q-1+q),c1(3)=4π∫0π/2y3(τ)cosτdτ=2π(1-q)3∑m=0∞(m+1)(m+2)qm∫0π/2cos2m+4τdτ=2π(1-q)3∑m=0∞π(m+1)(m+2)Γ(m+5/2)qm2Γ(m+3)=341-q.
Substitution of (4.10) into (4.7), yields the following expression
(4.11)2q(-ω2+α)(1-q-1+q)+34β1-q=0.
Solving for q gives
(4.12)q=-169β2[(ωaagm)2-α][4(ωaagm)2-4α-3β],
where α, β and ωaagm are given in (2.6) and (3.5). From (4.12) we obtain for λ5≠0(4.13)limλ5A4→∞q(A)=0.25072
and for λ3≠0 and λ5=0(4.14)limλ3A2→∞q(A)=0.16441.
In Figure 1 we have plotted parameter q as a function of A for different values of parameters λ1, λ3, and λ5. Equations (4.13) and (4.14) allow us to know the maximum value of parameter qfor λ5≠0 and λ3≠0 (λ5=0), respectively. In fact, for λ3≠0 and λ5=0 we have obtained an accurate approximate solution for the cubic Duffing equation.Figure 1
qas a function of A for cubic-quintic oscillators with the following parameters: (a) λ1=1, λ3=1, λ5=1 (red), (b) λ1=5, λ3=3, λ5=1 (green), and (c) λ1=1, λ3=10, λ5=1 (blue).


## 5. Comparison with the Exact Solution

The exact frequency for the cubic-quintic Duffing oscillator can be obtained by direct integration of (2.1) using the initial conditions in (2.2). The result is (see Appendix B in [1])
(5.1)ωe(A)=πk12∫0π/2(1+k2sin2t+k3sin4t)-1/2dt,
where parameters k1, k2, and k3 are given as follows:
(5.2)k1=λ1+12λ3A2+13λ5A4,k2=3λ3A2+2λ5A46λ1+3λ3A2+2λ5A4,k3=2λ5A46λ1+3λ3A2+2λ5A4.
For λ1≠0 it is possible to obtain the following expressions for very small values of A:
(5.3)ωe(A)≈λ1+38λ3λ1A2+80λ1λ5-21λ32256λ13/2A4+81λ33-304λ1λ3λ52048λ15/2A6+⋯,ωaagm(A)≈λ1+38λ3λ1A2+80λ1λ5-21λ32256λ13/2A4+81λ33-300λ1λ3λ52048λ15/2A6+⋯.
For λ5≠0 it is possible to obtain the following expression for very large values of λ5A4:
(5.4)ωe(A)≈0.746834λ5A2.
Based on (3.9) and (5.4), the approximation compared with respect to the exact frequency is
(5.5)limλ5A4→∞ωaagm(A)ωe(A)=1.0038.
As can be seen, the relative error of the analytical approximate frequency ωaagm is 0.38% in the limit as λ5A4→∞, which implies that the method proposed is suitable for solving (2.1) with λ5A4 much larger than unity. In Figure 2 we have plotted the comparison between the approximate and exact frequencies for different values of parameters λ1, λ3, and λ5.Figure 2
Comparison of approximate frequencies with exact frequencies for the cubic-quintic oscillator for (a)λ1=1, λ3=1, λ5=1 (red), (b) λ1=5, λ3=3, λ5=1 (green), and (c) λ1=1, λ3=10, λ5=1 (blue).For this cubic-quintic Duffing oscillator, Lai et al. obtained the following limits applying the Newton-harmonic balance method (NHBM) [1]:
(5.6)limλ5A4→∞ω1(A)ωe(A)=1.05856,limλ5A4→∞ω2(A)ωe(A)=1.01078,limλ5A4→∞ω3(A)ωe(A)=1.0023,
for the first-, second-, and third-order analytical approximations, respectively. As we can see, the approximate frequency obtained in this paper, combining a cubication procedure (based on Chebyshev polynomials expansion of the restoring force), the arithmetic geometric mean (for approximating the exact frequency of the cubic Duffing oscillator), and the second-order rational harmonic balance method, is more accurate than the first-order and second-order frequencies obtained using the NHBM and slightly less accurate (0.38% versus 0.23%) than the third-order frequency obtained using the NHBM.In Table2 we have presented the percentage errors for comparison of the approximate frequency with exact frequency for λ1=λ3=λ5=1. For these parameters (3.5) can be written as follows:
(5.7)ωaagm(A)=14[(1+12A2+52A4)1/4+(1+A2+152A4)1/4]2.
Taking into account (3.9), (4.2), and (4.13) it is easy to verify that the analytical approximate solution for the cubic-quintic Duffing oscillator has the following asymptotic expression for λ5A4→∞:
(5.8)yλ5A2→∞(t)=0.74928cos(0.749670λ51/2A2t)1-0.25072cos2(0.749670λ51/2A2t).
We analyze this example because this situation corresponds to the worst case, since it has the approximate frequency whose relative error is the greatest.Table 2
Percentage errors for comparison of the approximate frequency with exact frequency forλ1=λ3=λ5=1.
A
ω
a
agm
ω
e
(% error)
0.1
1.00377
1.00377
0.0000
0.3
1.03554
1.03554
0.0000
0.5
1.10657
1.10654
0.0027
1
1.52456
1.52359
0.064
3
7.29201
7.26863
0.32
5
19.2503
19.1815
0.36
8
48.4739
48.2946
0.37
10
75.4588
75.1774
0.37
20
300.355
299.223
0.38
50
1874.66
1867.57
0.38
70
3673.87
3659.98
0.38
100
7497.18
7468.83
0.38
300
67470.75
67215.57
0.38
500
187418.88
186709.04
0.38
700
367339.58
365949.25
0.38
1000
749670.06
746834.69
0.38Approximate and exact solutions for a few examples of the cubic-quintic Duffing oscillators are presented in Figures3, 4, 5, and 6 for one complete cycle. The exact periodic solutions (ye=xe/A) have been obtained by numerically integrating (2.1), and the proposed approximate periodic solutions (y=x/A) have been computed using (4.2). In these figures parameter h is defined as
(5.9)h=ωet2π.
which tends to 0.746834λ51/2A2t/2π as λ5A4→∞. The approximate solution, expressed as a function of h, which has been plotted in Figure 6 is
(5.10)yλ5A2→∞(h)=0.74928cos(6.3070448h)1-0.25072cos2(6.3070448h).
Parameters for oscillators analyzed in Figures 3–6 are presented in Table 3 in which the values for the L2 norm errors have been also included in order to produce a global estimation of the accuracy of the approximate solutions [18]
(5.11)L2=∫01|ye(h)-y(h)|2dh.
The maximum value of L2 norm is 0.0100931 and it is obtained for λ5A4→∞. In Figure 7 we have plotted the differences Δ=ye-y for the oscillators presented in Figures 3–6. Figures 3–7 show that (3.9) and (4.2) provide a good approximation to the exact periodic solution to the cubic-quintic Duffing oscillator.Table 3
Parameters for the cubic-quintic oscillators analyzed in Figures3–6.
λ
1
λ
3
λ
5
A
ω
e
ω
a
agm (% error)
q
L
2
1
10
1
1
2.96448
2.96504 (0.019)
0.15187
0.00055478
5
3
1
3
8.37877
8.39842 (0.23)
0.20980
0.0063554
1
1
1
3
7.26863
7.29201 (0.32)
0.23580
0.0085975
λ
5
A
4
→
∞
0.749670
λ
5
A
2
0.746834
λ
5
A
2 (0.38)
0.25072
0.0100931Figure 3
Comparison of the approximate solution (dashed line and triangles) with the numerical solution of (2.1) (continuous line and squares) for λ1=1, λ3=10, λ5=1, and A=1 (L2=0.00055478).Figure 4
Comparison of the approximate solution (dashed line and triangles) with the numerical solution of (2.1) (continuous line and squares) for λ1=5, λ3=3, λ5=1 and A=3 (L2=0.0063554).Figure 5
Comparison of the approximate solution (dashed line and triangles) with the numerical solution of (2.1) (continuous line and squares) for λ1=1, λ3=1, λ5=1 and A=3 (L2=0.0085975).Figure 6
Comparison of the approximate solution (dashed line and triangles) with the numerical solution of (2.1) (continuous line and squares) for λ5A4→∞ (L2=0.010093).Figure 7
Difference between normalized numerical and approximate solutions for: (orange)λ1=1, λ3=10, λ5=1, and A=1, (blue) λ1=5, λ3=3, λ5=1, and A=3; (green) λ1=1, λ3=1, λ5=1, and A=3; (red) λ5A4→∞.


## 6. Conclusions

Using the results obtained by applying a cubication method for the cubic-quintic Duffing oscillator, very simple and accurate expression for the frequency and the solution are obtained for the cubic-quintic Duffing oscillator. In this procedure the original nonlinear equation is replaced by a related “cubic Duffing equation” that approximates the original “cubic-quintic Duffing equation” closely enough to provide useful solutions. The cubic Duffing equation is exactly solved in terms of the complete elliptic integral of the first kind and the Jacobi elliptic function cn and these exact solutions are the approximate solutions to the cubic-quintic Duffing equation. Using the relationship between the complete elliptical integral of the first kind and the arithmetic-geometric mean and the second-order rational harmonic balance method, closed-form expressions for the approximate frequency and the solution are obtained in terms of elementary functions with a relative error for the approximate frequency less than 0.38% for all values of the oscillation amplitude. The results of the method presented in this paper have been found to be in good agreement with those of the exact solutions obtained numerically.



---

*Source: 286290-2012-10-04.xml*