# An Exact Solution to the Quadratic Damping Strong Nonlinearity Duffing Oscillator

Authors: Alvaro H. Salas; S. A. El-Tantawy; Noufe H. Aljahdaly
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The nonlinear equations of motion such as the Duffing oscillator equation and its family are seldom addressed in intermediate instruction in classical dynamics; this one is problematic because it cannot be solved in terms of elementary functions before. Thus, in this work, the stability analysis of quadratic damping higher-order nonlinearity Duffing oscillator is investigated. Hereinafter, some new analytical solutions to the undamped higher-order nonlinearity Duffing oscillator in the form of Weierstrass elliptic function are obtained. Posteriorly, a novel exact analytical solution to the quadratic damping higher-order nonlinearity Duffing equation under a certain condition (not arbitrary initial conditions) and in the form of Weierstrass elliptic function is derived in detail for the first time. Furthermore, the obtained solutions are camped to the Runge–Kutta fourth-order (RK4) numerical solution.


---

## Body

## 1. Introduction

In the last few decades, considerable work has been invested in developing new methods for analytical and numerical solutions for several differential equations [1–4] such as strongly nonlinear oscillators, but it is still difficult to obtain convergent results in cases of strong nonlinearity. The Duffing oscillator model is one of the most important and most popular models in dynamic systems due to its importance in explaining many nonlinear phenomena in science and engineering. For instance, the nonlinear Duffing equation was devoted for studying the dynamical behavior of oscillations in plasma physics [5], plasma physics and rigid rotator [6], magneto-elastic mechanical systems [7], nonlinear vibration of beams and plates [8], fluid flow induced vibration [9], large amplitude oscillation of centrifugal governor systems [10], etc. Since the early century so far, a lot of researchers have been devoted carefully to study the Duffing oscillator model in several physical situations [11]. So far the Duffing oscillator model has been an effective and applicable example to nonlinear oscillations in dynamic systems [12]. In addition to its mentioned applications, it also has various applications in the mechanical engineering [9, 10] and electrical engineering [13, 14]. It has been solved numerically and analytically via using various methods in order to recognize its physical characters. Some methods which consider the conservative Duffing osciullators are the harmonic balance technique [15], the frequency-amplitude formulation [16], the homotopy perturbation method [10, 17], the coupled homotopy-variational approach [18], and the differential transform [19]. It is known that if the damping effect is taken into account in the conservative Duffing oscillators, the model becomes a nonconservative system and the oscillation amplitude reduces over time. There are many methods that have been devoted to solve the nonconservative Duffing oscillators, whether they are semianalytical or numerical solutions [20–26]. However, most of the analytical methods are unable to deal with the Duffing equation in the presence of a damping effect. Therefore, all published solutions are either approximate numerical solutions, semianalytical or approximate analytical solutions, or analytical solutions under certain restrictions [20–31].The most popular formula for conservative Duffing oscillator for arbitrary initial conditions is given by [32, 33](1)x¨t+pxt+qx3t=0,x0=x0,x′0=x˙0.This equation is completely integrable, and its exact analytical solutions in the Jacobi elliptic function forms have been obtained [32, 33]. However, equation (1) without damping effect is not realistic physically, and if the damping effect is taken into account, the following nonconservative Duffing oscillator for arbitrary initial conditions may be obtained as [5, 6, 34–37](2)x¨+εx˙+px+qx3=0,x0=x0&x′0=x˙0.In generally, equation (2) is not integrable differential equation due to the presence damping term εx˙, and many methods have been devoted to solve this equation, as was previously mentioned in all the upper references such as [5, 6, 34–37].The generalized form of the Duffing oscillator equation which has been modelled in an arbitrary form of the nonlinear differential equation may be written as follows:(3)x¨+εx˙x˙p−1+fx=0.In equation (3), the second term x˙x˙p−1 refers to the nonlinear damping term which is proportional to the p th power of x˙ and ε is a positive or negative damping parameter where its sign depends on the type of the physical application. The third term in equation (3), i.e., fx=∑j=1nanxn involves the linear term a1x and the nonlinear terms of x which is known as the restoring force of the dynamical system. As we mentioned earlier, many authors studied the linear damping p=1 for linear and cubic restoring forces x and x3 [5, 6, 34–37]. On the contrary, there are a few attempts to solve and investigate the damping Duffing oscillator for quadratic restoring force p=2 [38–40]. The elliptic-harmonic balance method has been applied to find an approximate analytical solution in the Jacobi elliptic function form to a quadratic damping system with linear restoring force: x¨+εx˙2+αx=0 [38]. In [40], some solutions to the quadratic damping Helmholtz–Duffing equation x¨+εx˙2+αx+βx2+γx3=0 in the form of classical Jacobi elliptic functions have been discussed in details. Given the difficulty of solving quadratic damping Duffing equation with higher-order nonlinearity, this problem remains one of the most important hot topics that require a more careful study to obtain more accurate solutions.In this paper, the quadratic damping Duffing oscillator with linear and cubic restoring forcesfx=px+qx3 is considered under the initial conditions x0=x0 and x˙0=x˙0. Thus, the mathematical model of this physical situation may be written as the initial value problem (IVP):(4)x¨+εx˙2+px+qx3=0,x0=x0&x′0=x˙0.The IVP (4) is considered a nonintegrable differential equation but may admit exact analytical solution under a certain condition. In this study, we will derive this condition in detail in order to find an exact solution to the IVP (4) in terms of the Weierstrass elliptic function for the first time. Moreover, before finding the exact solution to the IVP (4), we will study its stability in order to analyze its characteristics to recognize its dynamical behavior.The rest of the paper is presented in the following points. The stability analysis of the quadratic damping higher-order nonlinearity Duffing equation will be investigated in Section2. Some possible analytical solutions to the Duffing equation and its family in the form of Jacobi elliptic function and Weierstrass elliptic function will be discussed in detail. Firstly, we will discuss some possible solutions to the IVP (2) in terms of Jacobi elliptic function depending on the discriminant of the IVP (2). Thereafter, the exact analytical solutions to both the IVP (2) and IVP (4) in terms of Weierstrass elliptic function will be derived and obtained in details. The summary of our investigations will be introduced in Section 4.


## 2. Stability Analysis to the Quadratic Damping Duffing Equation

The first equation in the IVP (4) may be written as a dynamical system as follows:(5)x˙=fx,y=y,y˙=gx,y=−εy2−px−qx3.The first integral reads(6)e2εx4ε4−2ε2p+4ε3px+4ε3qx3−6ε2qx2+6εqx−3q+4ε4y2=h,System (5) becomes Hamilton if the following condition is fulfilled(7)∂xfx,y+∂ygx,y=0,with(8)fx,y=∂yhx,y,gx,y=−∂xhx,y,It is noted that condition (7) is fulfilled; thus, system (5) becomes Hamilton with the following Hamiltonian:(9)Hx,y=e2εx4ε4−2ε2p+4ε3px+4ε3qx3−6ε2qx2+6εqx−3q+4ε4y2.Now, ifpq<0, then system (5) has three equilibrium points:(10)x1,y1=0,0,x2,y2=−pq,0,x3,y3=−−pq,0,and, for pq>0, there is only one equilibrium point:(11)x1,y1=0,0.The Jacobian matrix for system (5) at equilibrium point xe,ye reads(12)Axe,ye=0−p−3qxe21−2yeε.Thus, the determinant and trace of the Jacobian matrixAxe,ye may be estimated as(13)DetAxe,ye=p+3qxe2,TrAxe,ye=−2εye.The behavior of the dynamical system near the equilibria is described by the Poincaré Diagram as given in [41].Let us analyze the equilibrium points.Forxe,ye=0,0, we obtain(14)DetA0,0=p,TrA0,0=0.The origin becomes a central point ifp>0 and a saddle if p<0, as shown in Figures 1(a) and 1(b), respectively.Forpq<0 and xe,ye=−p/q,0 or xe,ye=−−p/q,0, we obtain(15)DetAxe,ye=−2p,TrAxe,ye=0.It is clear that ifp<0, the equilibrium point becomes a central point and a saddle point is covered if p>0, as illustrated in Figures 2(a) and 2(b), respectively.Figure 1
Phase Portrait of system (5) for the equilibrium point xe,ye=0,0: (a) for central point p>0 where p=1, q=−2, and ε=0.2 and (b) for saddle point p<0 where p=−1, q=−2, and ε=0.2.
(a)(b)Figure 2
Phase Portrait of system (5) for the equilibrium points xe,ye=−p/q,0 or xe,ye=−−p/q,0 for pq<0: (a) for central point p<0 where p=−2, q=0.4, and ε=0.2 and (b) for saddle point p>0 where p=2, q=−0.4, and ε=0.2.
(a)(b)


## 3. Analytical Solutions to the Quadratic Damping Duffing Equation

Before embarking on solving the IVP (4), we will discuss some possible solutions to both the IVP (1) and the IVP (2) in terms of Weierstrass elliptic function.


### 3.1. The Analytical Solution of Undamping Duffing Equation in Terms of Weierstrass Elliptic Function

Let us solve the following undamped IVP:(16)x″t+pxt+qx3t=0,x0=x0&x′0=x˙0,where x02+x˙02q≠0, and the primes represent the differentiation with respect to the independent variable say “t.” Here, it is assumed that p, q, u0, and u˙0 are real numbers. To find an analytical solution to the IVP (16) in the form of Weierstrass elliptic function ℘, the following ansatz is assumed:(17)xt=A+B1+C℘t+t0;g2,g3,where A, B, C, g2, g3, and t0are some constants to be determined later and the condition p+qx022+2qx˙02>0 must be fulfilled. Note that the Weierstrass elliptic function ℘≡℘t+t0;g2,g3≡℘z;g2,g3 satisfies the following relations:(18)℘′2=4℘3−g2℘−g3,℘″=−g22+6℘2.Substituting both equations (17) and (18) into equation (16) and after some simplifications, we finally obtain(19)121+C℘∑j=03ℝj℘j=0,where the values of the coefficients ℝj are given in Appendix A. Equating the coefficients ℝj to zero, we get a system of algebraic equations, and by solving this system, we obtain(20)B=−A2+4p3A2q+p,C=123A2q+p,g2=112−3A4q2−6A2pq+p2,g3=1216p9A4q2+18A2pq+p2.Thus, the general exact solution to the undamping strong nonlinear Duffing equationx″t+pxt+qx3t=0 in terms of Weierstrass elliptic function is given by(21)xt=A−6AA2q+p3A2q+p+12℘t+t0;1/12−3q2A4−6pqA2+p2,1/216p9q2A4+18pqA2+p2.Inserting the initial conditionsx0=x0 and x′0=x˙0 into solution (21), we finally get the value of t0 as follows:(22)t0=±℘−1x0−A−BCA−x0;g2,g3,and the value of A may be determined from the condition(23)x″0+px0+qx30=0,and this condition is equivalent to(24)qA4+2pA2−2px02+qx04+2x˙02=0.The roots of equation (24) give the value of A as follows:(25)A=±−p±Δq,where Δ=p+qx022+2qx˙02.Remark 1.
Solution (21) is periodic with period:(26)T=2∫e1∞dx4x3−g2x−g3,where e1 is the greatest real root to the cubic(27)4z3−g2z−g3=0.


### 3.2. An Analytical Solution of Quadratic Damping Duffing Equation in Terms of Weierstrass Elliptic Function

Here, to find an exact solution to the IVP (4) in terms of Weierstrass elliptic function, the following ansatz is introduced:(28)xt=A+B1+C℘t,where ℘t≡℘t−t0;g2,g3.Inserting ansatz (28) and equation (18) into x¨+εx˙2+px+qx3=0, we get Rt=0, with(29)Rt=∑j=04Lj℘jt,where the values of the coefficients Lj are given in Appendix B. Equating the coefficients Lj to zero gives an algebraic equation system, and by solving this system, we obtain(30)B=6AA2q+p2A3qε−3A2q+2Apε−p,(31)C=−122A3qε−3A2q+2Apε−p,(32)g2=112−4q2ε2A6−q8pε2+3qA4−2p2pε2+3qA2+p2,(33)g3=W1×W2,with(34)W1=A3qε+Apε+p2165A3qε−3A2q+5Apε−p,W2=8q3ε3A9−12q3ε2A8+6q2ε4pε2+3qA7−9q24pε2+3qA6+6pqε4pε2+11qA5−9pq4pε2+7qA4+2p2ε4pε2+27qA3−3p24pε2+7qA2+6p3εA−p3.The value ofA may be determined from the following cubic equation:(35)P≡−3q−8pε2+4ε3q+8pε2A−24qε2A2+32qε3A3=0.We choose the first root of equation (35) as the default value to A or any real root. Now, we must determine the value of t0 so that(36)x0=x0,x′0=x˙0.Let(37)c0=℘t0;g2,g3,c1=℘′t0;g2,g3,t0=−℘−1c0;g2,g3.Making use of the addition formula for the Weierstrass function gives(38)xt=A+B1+C1/4c1−℘t;g2,g3/c0−℘t;g2,g32−c0−℘t;g2,g3,where the values of A, B, and C are defined in equations (30), (31), and (35).Using the initial conditions given in equation (36) will give us(39)A+Bc0C+1=x0,−Bc1Cc0C+12=x˙0.By solving system (39), we can get the values of c0 and c1 as(40)c0=x0−A−BCA−x0,c1=−Bx˙0CA−x02.Note that the value ofA must also satisfy the condition(41)R0=0=2qε2A6−3qε2x0ε−1A5+2pε2+6qx02ε2−3qx0ε+3qA4−ε6px0ε−5p+2qx03ε+3qx02A3+32px02ε2−3px0ε+2p+qx03εA2−ε2px03ε−3px02−3x˙02A+px03ε−6px02−3qx04−3x˙02x0ε−6x˙02.The resultant of the polynomials given in equations (35) and (41) with respect to A reads(42)Resultant=−432q5ε64x0p+qx02ε3+6p+18qx02ε2+30qx0ε+17q×4x˙02ε4+4x0p+qx02ε3−2p+3qx02ε2+6qx0ε−3q3.This resultant must vanish in order to get exact solutions to the problem under consideration.We may consider the condition(43)4x˙02ε4+4x0p+qx02ε3−2p+3qx02ε2+6qx0ε−3q=0.From this condition, we can calculate the value ofq as(44)q=−2ε2p2εx0−1+2ε2x˙024ε3x03−6ε2x02+6εx0−3.Thus, the IVP (4) has been solved analytically under restriction (43) or (44).A comparison between the analytical solution (38) and the approximate numerical RK4 solution is introduced in Figure 3. It is observed that our solution is in agreement with the numerical simulation solution.Figure 3
A comparison between the analytical solution (38) (dashed curve) and numerical solution using RK4 (dotted curve). Here, p=−0.7, ε=0.1, x0=1, and x˙0=0.25.In the limiting case, i.e., forε⟶0, the exact solution to the undamping Duffing problem (1) for arbitrary initial conditions may be covered(45)xt=A−6Ap+A2qp+3A2q+12℘t−t0;1/12−3q2A4−6pqA2+p2,1/216p9q2A4+18pqA2+p2,where(46)t0=℘−13A3q+3A2qx0+5pA+px012A−x0;112−3q2A4−6pqA2+p2,1216p9q2A4+18pqA2+p2,A=±−p±p+qx022+2qx˙02q.Note that both solutions (21) and (45) are equivalent.


## 3.1. The Analytical Solution of Undamping Duffing Equation in Terms of Weierstrass Elliptic Function

Let us solve the following undamped IVP:(16)x″t+pxt+qx3t=0,x0=x0&x′0=x˙0,where x02+x˙02q≠0, and the primes represent the differentiation with respect to the independent variable say “t.” Here, it is assumed that p, q, u0, and u˙0 are real numbers. To find an analytical solution to the IVP (16) in the form of Weierstrass elliptic function ℘, the following ansatz is assumed:(17)xt=A+B1+C℘t+t0;g2,g3,where A, B, C, g2, g3, and t0are some constants to be determined later and the condition p+qx022+2qx˙02>0 must be fulfilled. Note that the Weierstrass elliptic function ℘≡℘t+t0;g2,g3≡℘z;g2,g3 satisfies the following relations:(18)℘′2=4℘3−g2℘−g3,℘″=−g22+6℘2.Substituting both equations (17) and (18) into equation (16) and after some simplifications, we finally obtain(19)121+C℘∑j=03ℝj℘j=0,where the values of the coefficients ℝj are given in Appendix A. Equating the coefficients ℝj to zero, we get a system of algebraic equations, and by solving this system, we obtain(20)B=−A2+4p3A2q+p,C=123A2q+p,g2=112−3A4q2−6A2pq+p2,g3=1216p9A4q2+18A2pq+p2.Thus, the general exact solution to the undamping strong nonlinear Duffing equationx″t+pxt+qx3t=0 in terms of Weierstrass elliptic function is given by(21)xt=A−6AA2q+p3A2q+p+12℘t+t0;1/12−3q2A4−6pqA2+p2,1/216p9q2A4+18pqA2+p2.Inserting the initial conditionsx0=x0 and x′0=x˙0 into solution (21), we finally get the value of t0 as follows:(22)t0=±℘−1x0−A−BCA−x0;g2,g3,and the value of A may be determined from the condition(23)x″0+px0+qx30=0,and this condition is equivalent to(24)qA4+2pA2−2px02+qx04+2x˙02=0.The roots of equation (24) give the value of A as follows:(25)A=±−p±Δq,where Δ=p+qx022+2qx˙02.Remark 1.
Solution (21) is periodic with period:(26)T=2∫e1∞dx4x3−g2x−g3,where e1 is the greatest real root to the cubic(27)4z3−g2z−g3=0.


## 3.2. An Analytical Solution of Quadratic Damping Duffing Equation in Terms of Weierstrass Elliptic Function

Here, to find an exact solution to the IVP (4) in terms of Weierstrass elliptic function, the following ansatz is introduced:(28)xt=A+B1+C℘t,where ℘t≡℘t−t0;g2,g3.Inserting ansatz (28) and equation (18) into x¨+εx˙2+px+qx3=0, we get Rt=0, with(29)Rt=∑j=04Lj℘jt,where the values of the coefficients Lj are given in Appendix B. Equating the coefficients Lj to zero gives an algebraic equation system, and by solving this system, we obtain(30)B=6AA2q+p2A3qε−3A2q+2Apε−p,(31)C=−122A3qε−3A2q+2Apε−p,(32)g2=112−4q2ε2A6−q8pε2+3qA4−2p2pε2+3qA2+p2,(33)g3=W1×W2,with(34)W1=A3qε+Apε+p2165A3qε−3A2q+5Apε−p,W2=8q3ε3A9−12q3ε2A8+6q2ε4pε2+3qA7−9q24pε2+3qA6+6pqε4pε2+11qA5−9pq4pε2+7qA4+2p2ε4pε2+27qA3−3p24pε2+7qA2+6p3εA−p3.The value ofA may be determined from the following cubic equation:(35)P≡−3q−8pε2+4ε3q+8pε2A−24qε2A2+32qε3A3=0.We choose the first root of equation (35) as the default value to A or any real root. Now, we must determine the value of t0 so that(36)x0=x0,x′0=x˙0.Let(37)c0=℘t0;g2,g3,c1=℘′t0;g2,g3,t0=−℘−1c0;g2,g3.Making use of the addition formula for the Weierstrass function gives(38)xt=A+B1+C1/4c1−℘t;g2,g3/c0−℘t;g2,g32−c0−℘t;g2,g3,where the values of A, B, and C are defined in equations (30), (31), and (35).Using the initial conditions given in equation (36) will give us(39)A+Bc0C+1=x0,−Bc1Cc0C+12=x˙0.By solving system (39), we can get the values of c0 and c1 as(40)c0=x0−A−BCA−x0,c1=−Bx˙0CA−x02.Note that the value ofA must also satisfy the condition(41)R0=0=2qε2A6−3qε2x0ε−1A5+2pε2+6qx02ε2−3qx0ε+3qA4−ε6px0ε−5p+2qx03ε+3qx02A3+32px02ε2−3px0ε+2p+qx03εA2−ε2px03ε−3px02−3x˙02A+px03ε−6px02−3qx04−3x˙02x0ε−6x˙02.The resultant of the polynomials given in equations (35) and (41) with respect to A reads(42)Resultant=−432q5ε64x0p+qx02ε3+6p+18qx02ε2+30qx0ε+17q×4x˙02ε4+4x0p+qx02ε3−2p+3qx02ε2+6qx0ε−3q3.This resultant must vanish in order to get exact solutions to the problem under consideration.We may consider the condition(43)4x˙02ε4+4x0p+qx02ε3−2p+3qx02ε2+6qx0ε−3q=0.From this condition, we can calculate the value ofq as(44)q=−2ε2p2εx0−1+2ε2x˙024ε3x03−6ε2x02+6εx0−3.Thus, the IVP (4) has been solved analytically under restriction (43) or (44).A comparison between the analytical solution (38) and the approximate numerical RK4 solution is introduced in Figure 3. It is observed that our solution is in agreement with the numerical simulation solution.Figure 3
A comparison between the analytical solution (38) (dashed curve) and numerical solution using RK4 (dotted curve). Here, p=−0.7, ε=0.1, x0=1, and x˙0=0.25.In the limiting case, i.e., forε⟶0, the exact solution to the undamping Duffing problem (1) for arbitrary initial conditions may be covered(45)xt=A−6Ap+A2qp+3A2q+12℘t−t0;1/12−3q2A4−6pqA2+p2,1/216p9q2A4+18pqA2+p2,where(46)t0=℘−13A3q+3A2qx0+5pA+px012A−x0;112−3q2A4−6pqA2+p2,1216p9q2A4+18pqA2+p2,A=±−p±p+qx022+2qx˙02q.Note that both solutions (21) and (45) are equivalent.


## 4. Conclusion

The equation of motion with higher-order nonlinear restoring force and higher-order nonlinear damping force in which the damping is proportional to the square of the velocity has been carefully studied. As an example of this type of equation of motion, the stability analysis of the quadratic damping higher-order nonlinear Duffing oscillator equation using the Bifurcation method has been investigated. Then comes the main goal, which is to find some exact solutions to the Duffing oscillator equation and its family. In the form of Weierstrass elliptic function, the exact solution of the undamping higher-order nonlinear Duffing equation has been discussed. As for the second and most important goal of this study, it is to obtain an exact solution to the quadratic damping higher-order nonlinear Duffing oscillator equation. In general, the damping higher-order nonlinear Duffing oscillator equation is nonintegrable differential equation but under a certain condition it becomes integrable. Accordingly, the condition that makes this equation integrable has been deduced, and then, this solution has been obtained in detail in terms of the Weierstrass elliptic function. Furthermore, a comparison between the obtained solutions and the numerical RK4 solution to the Duffing equation and its family has been carried out. Also, in the limiting case, i.e., forε⟶0, the exact solution of the quadratic damping higher-order nonlinear Duffing oscillator equation has been reduced to the exact solution of the undamping higher-order nonlinear Duffing equation.In future work, the Weierstrass elliptic functions will be devoted to find some exact and approximate analytical solutions to the forced higher-order damping cubic-quintic Duffing oscillator with external constant and periodic forces.



---

*Source: 8875589-2021-01-19.xml*