# Some Novel Analytical Approximations to the (Un)damped Duffing–Mathieu Oscillators

Authors: Haifa A. Alyousef; Alvaro H. Salas; Sadah A. Alkhateeb; S. A. El-Tantawy
Journal: Journal of Mathematics
 (2022)

---

## Abstract

Some novel exact solutions and approximations to the damped Duffing–Mathieu-type oscillator with cubic nonlinearity are obtained. This work is divided into two parts: in the first part, some exact solutions to both damped and undamped Mathieu oscillators are obtained. These solutions are expressed in terms of the Mathieu functions of the first kind. In the second part, the equation of motion to the damped Duffing–Mathieu equation (dDME) is solved using some effective and highly accurate approaches. In the first approach, the nonintegrable dDME with cubic nonlinearity is reduced to the integrable dDME with linear term having undermined optimal parameter (maybe called reduced method). Using a suitable technique, we can determine the value of the optimal parameter and then an analytical approximation is obtained in terms of the Mathieu functions. In the second approach, the ansatz method is employed for deriving an analytical approximation in terms of trigonometric functions. In the third approach, the homotopy perturbation technique with the extended Krylov–Bogoliubov–Mitropolskii (HKBM) method is applied to find an analytical approximation to the dDME. Furthermore, the dDME is solved numerically using the Runge–Kutta (RK) numerical method. The comparison between the analytical and numerical approximations is carried out. All obtained approximations can help a large number of researchers interested in studying the nonlinear oscillations and waves in plasma physics and many other fields because many evolution equations related to the nonlinear waves and oscillations in a plasma can be reduced to the family of Mathieu-type equation, Duffing-type equation, etc.


---

## Body

## 1. Introduction

Several physical and natural phenomena related to biology, chemistry, physics, engineering problems, and so on can be modelled by both ordinary differential equations (ODEs) and partial differential equations (PDEs) for studying the nonlinear self-excited oscillators [1–12]. Also, in many real-life problems, some internal and external forces that can affect the system under consideration cannot be neglected. For example, the friction and collisional force and many others that affect on the motion of particles, whether in solid, liquid, gas or plasma physics, cannot be neglected. Therefore, these forces must be included in the mathematical models that will be used for studying the natural and physical problem, such as investigating the nonlinear oscillations in various plasma models [13–23]. Interest in the study of nonlinear oscillations in a plasma is due to its many potential applications. Nowadays, plasma processing is seen as an important and effective technology which has been able to enter into all modern industries. In addition, plasma had a great credit for the modern technology in electronics, medicine, agriculture, biomedicine, automobiles, optics, aerospace, telecommunications, solar energy, polymers, papers, textiles, etc. [16, 24]. Therefore, we focus our attention on the applications of the family of Duffing-type equation for modelling the nonlinear oscillations in a plasma. The Duffing-type oscillator is one of the most popular differential equations that has spread widely due to its ability to explain many nonlinear phenomena in various fields of science and in mechanical systems and engineering. This equation is a mathematical model described by a second order of ODEs with a nonlinear spring force. It is used for describing the motion of (un)damped oscillator rather than simple harmonic motion. Motivated by potential engineering and plasma physics applications in addition to many others applications in electronic circuits and microcontroller, the family of Duffing equation such as (un)damped Duffing equation with lower and higher-order nonlinearities, (un)damped Duffing–Helmholtz equation, and (un)damped Mathieu–Duffing equation have received wide attention due to their ability for investigating the mechanism of a rigid pendulum oscillator, oscillations in different plasma models, and so on. This family of second-order differential equations has provided useful and successful models for investigating the nonlinear oscillations and chaotic nature. The biggest challenge in the study of the dynamics of nonlinear mechanical systems is to find some real solutions (including the analytical and numerical solutions) to the evolution equations that are used for describing the characteristics of nonlinear phenomena under consideration. Accordingly, studying the solutions of many equations of motion to various oscillators is one of the most difficult tasks facing many researchers.The solutions associated to mentioned evolution equations and many other related equations have been studied extensively due to the fact that such equations arise in a variety of realistic problems. For instance, the fluid equations of electron-ion unmagnetized cold plasma were reduced to Mathieu equation in order to investigate the electron waves [13]. During this study, the authors assumed that the density perturbations of the plasma species are only time-dependent functions and do not depend on space. Also, the basic set of fluid equations to a multicomponent complex plasma consisting of inertial two types of dust grains (both positive and negative charges) as well as inertialess Maxwellian species including electrons and ions were reduced to a Mathieu-type equation for studying the excitation of dust-acoustic oscillations [14]. Moreover, a nonlinear Van der Pol–Mathieu-type equation was derived for the dust grain density in order to investigate the dynamics of dust-acoustic oscillations in a dusty plasma consisting of inertial Boltzmann distributed species (electrons and ions) and inertialess dust grains [15]. A modified Van der Pol–Duffing oscillator with forced term was used in the study and was used in modelling the dynamics of nonlinear oscillations in different plasma models [16]. More recently, the multistage method was used for solving the damping Duffing equation with forced term in order to model the oscillations in a complex unmagnetized plasma [17].Due to the importance of the family of Duffing-type equation and motivated by the mentioned studies, we focus our attention on the analysis of the so-called (un)damped Mathieu–Duffing oscillator with a twin-well potential [25].(1)ℝ≡x¨+α−Q0cosΩtx+βx3=0,and(2)ℝ1≡x¨+εx˙+α−Q0cosΩtx+βx3=0,for studying the vibrating/oscillating behavior of systems described by (1) and (2), where α,β>0 are, respectively, the stiffness coefficients of the linear and nonlinear terms and ε represents the coefficients of damping term. The longitudinal loading is periodic, Ω and Q0 are the frequency and excitation strength of the periodic loading, respectively. The total energy of the undamped Mathieu–Duffing oscillator (DMO) according to (1) is defined by H=H0+H1, where H0=1/2x˙2−1/2αx2+1/4βx4 and H1=−1/2x2Q0cosΩt are, respectively, the unperturbed and perturbation in the Hamiltonian of (1). There is another form to the equation of periodic motion which is called (un)damped Mathieu–Helmholtz oscillator.(3)x¨+α−Q0cosΩtx+βx2=0,x¨+εx˙+α−Q0cosΩtx+βx2=0.The analytical solution to the damped Mathieu–Helmholtz oscillator (3) was obtained using the finite Fourier series expansion [26]. Note that in evolution equation (3), the nonlinear term βx2 is different from the nonlinear term in equations (1) and (2).The objectives of our study are to find some novel solutions to the (un)damped Duffing–Mathieu-type oscillator, under the initial conditionsx0=0 and x′0=x˙0. Two cases for Duffing–Mathieu-type oscillator will be discussed. In the first case, we will get some exact solutions to (un)damped Mathieu equation in terms of the Mathieu functions of the first kind. In the second case, the damped Duffing–Mathieu oscillator (dDMO) will be solved analytically and numerically using some different approaches. In the first approach, the cubic nonlinear term in equation (1) βx3 is replaced by the linear term βκx, where the constant κ≥0 represents an optimal parameter. Then, the nonintegrable dDMO reduces to an integrable one which has an exact solution but with undermined parameter κ. Using a suitable technique, we can determine the value of the optimal parameter κ. Thus, we can get an analytical approximation to the dDMO (2) in terms of the Mathieu functions. For the second approach, the ansatz method with the help of the solution to the undamped Duffing oscillator is employed to derive an analytical approximation to the dDMO (2) in the form of trigonometric functions. Furthermore, the homotopy perturbation technique with the extended Krylov–Bogoliubov–Mitropolskii (KBM) which is called HKBM method is also devoted for solving the dDME (2) for arbitrary physical parameters [27, 28].


## 2. Mathematical Analysis

Here, we proceed to find some approximate solutions to both undamped and damped Duffing–Mathieu oscillators (1) and (2), respectively. Below we discuss the different approaches for solving the mentioned equations.


### 2.1. An Exact Solution to (Un)damped Mathieu Equation

Both undamped Mathieu equation, i.e., (1) for β=0, and the damped Mathieu equation, i.e., (2) for β=0, have exact solutions in the form of Mathieu functions. First, let us find the solution of the following damped Mathieu equation:(4)x¨+εx˙+α−Q0cosΩtx=0,with subjected to the initial conditions (ICs): x0=x0 and x′0=x˙0.Using the following MATHEMATICA command(5)gt_≔x′′t+2εx′t−α−Q0CosΩtxtNDSolvegt==0&&x0==x0&&x′0==x˙0,xt,t1,1,2//FullSimplify1,we can get the exact solution to (4) as follows:(6)xt=e−εtx0MC1t,εMC10,ε+2εx0+x˙0ΩMS1t,εMathieuSPrime4α−ε2/Ω2,2Q0/Ω2,0,with(7)MC1t,ε=MathieuC4α−ε2Ω2,2Q0Ω2,Ω2t,MS1t,ε=MathieuS4α−ε2Ω2,2Q0Ω2,Ω2t,where MathieuS and MathieuC are the Mathieu functions of the first kind or sometimes called sine-elliptic and cosine-elliptic, respectively. For ε=0, the damped Mathieu equation (4) reduces to the following undamped Mathieu equation:(8)x¨−α+Q0cosΩtx=0,and solution (6) reduces to the following one:(9)xt=x0MC1t,0MC10,0+2x˙0MS1t,0ΩMathieuSPrime4α/Ω2,2Q0/Ω2,0.The exact solution (6) is compared with RK numerical solution as shown in Figures 1(a) and 1(b) for different values of x0.Figure 1
Both exact solution (6) and RK numerical solution to the damped Mathieu equation (4) are compared with each other for different values to the initial amplitude x0.
(a)(b)


### 2.2. Some Analytical Approximations to the Damped Duffing–Mathieu Oscillator

Here, we proceed to discuss two techniques (the hybridp−expansion method and the ansatz method) for finding some analytical approximations to dDMO (2). For studying dDMO (2), first we rewrite equation (2) in the following new i.v.p.(10)ℝ1=0,x0=x0&x′0=x˙0.


#### 2.2.1. First Approach: Reduced Method

Forβ≠0 and α>0, we may obtain simple approximation to the i.v.p. (10) by replacing the cubic term βx3 by the linear term βκx, where the constant κ≥0 which is used as an optimal parameter to reduce the residual error. Accordingly, dDMO (2) of cubic nonlinearity reduces to the following dDMO with linear term. Thus, we can replace the i.v.p. (10) by the following new i.v.p.(11)x¨+2εx˙+α−Q0cosΩtx+βκx=0,x0=x0andx′0=x˙0.Thus, the exact solution to the i.v.p. (11) is expressed by(12)xκ≡xκt=e−εtx0MC2t,ε,κMC20,ε,κ+2εx0+x˙0ΩMS2t,ε,κMathieuSPrime4βκ−ε2+α/Ω2,2Q0/Ω2,0,with(13)MC2t,ε,κ=MathieuC4βκ−ε2+αΩ2,2Q0Ω2,Ω2t,MS2t,ε,κ=MathieuS4βκ−ε2+αΩ2,2Q0Ω2,Ω2t.The residual is defined as(14)Rκt=x¨κ+εx˙κ+α−Q0cosΩtxκ+βxκ3.A suitable value ofκ can be obtained by solving the equation Rκt0=0 for some t0>0, say t0=1. Making use of the Padé approximate technique, for 0<t0≤1 and x0≈0, we get(15)κsuitable=Y1Y2,with(16)Y1=3x02Q0t0x02−αt0x02−2εt0x˙0x0−4t0x˙02−2x˙0x0×4Q0t0x02−4αt0x02+3βt0x04−8εt0x˙0x0−6t0x˙02−8x˙0x0,Y2=−26αQ0t02x04+12βQ0t02x06−52εQ0t02x˙0x03+13Q02t02x04−52Q0t0x˙0x03−30Q0t02x˙02x02+13α2t02x04−12αβt02x06+52αεt02x˙0x03+52αt0x˙0x03+30αt02x˙02x02+9β2t02x08−24βεt02x˙0x05−24βt0x˙0x05+12βt02x˙02x04+52ε2t02x˙02x02+104εt0x˙02x02+60εt02x˙03x0+60t0x˙03x0+12t02x˙04+52x˙02x02.For a givenκ, the residual error LRκ of the approximation (12) to the i.v.p. (10) is defined as(17)LRκ=max0≤t≤TRκt=max0≤t≤Tx¨κ+2εx˙κ+α−Q0cosΩtxκ+βxκ3.The optimal value ofκoptimal for the parameter κ on the interval 0≤t≤T is defined as(18)κoptimal=minκ≥0LRκ.Let us apply the obtained approximation (12) for investigating the properties of the damping oscillations to the dDMO (2) at different values to the physical parameters α,β,ε,Ω,Q0,x0,x˙0. The profile of the approximation (12) using the values of κsuitable (given in equation (15)) and κoptimal (given in (18)) is displayed in Figures 2 and 3 at α,β,ε,Ω,Q0,x0,x˙0=4,1,0.1,0.5,0.1,0,0.2 and α,β,ε,Ω,Q0,x0,x˙0=4,1,0.1,0.5,0.1,π/6,0.2, respectively. The obtained results showed that this approximation gives results with good and acceptable accuracy.Figure 2
The profile of the approximation (12) using the values of κsuitable (given in equation (15)) and κoptimal (given in (18)) is plotted against t for α,β,ε,Ω,Q0,x0,x˙0=4,1,0.1,0.5,0.1,0,0.2.
(a)(b)Figure 3
The profile of the approximation (12) using the values of κsuitable (given in equation (15)) and κoptimal (given in (18)) is plotted against t for α,β,ε,Ω,Q0,x0,x˙0=4,1,0.1,0.5,0.1,π/6,0.2.
(a)(b)


#### 2.2.2. The Ansatz Method for Solving dDMO

Now, we can summarize the main points to get some approximations to the i.v.p. (10) in the following steps.Step 1.
Let us assume the following ansatz:(19)xt=yft,where the time-dependent function f≡ft can be determined later and the function y≡yt represents the solution of the following i.v.p.(20)y¨+2εy˙+α−Q0y+βy3=0,y0=x0&y′0=x˙0.Step 2.
Inserting ansatz (19) into i.v.p. (10), we get(21)ℝ1=y′f−2εf′2+2εf′+f′′+yfα−Q0cosΩt+f′2Q0−α+βy3f1−f′2.Step 3.
For vanishing the coefficient ofyf in equation (21), we get(22)f′2=−α+Q0cosΩtQ0−α.Step 4.
Integrating equation (22) and with the help of f0=0, we have(23)ft=2ΩEΩ2t,2Q0Q0−α,where E stands for the EllipticE function.
Now, the value offt has been determined but the solution of the i.v.p. (20) needs to be determined. Thus, we are faced with two things: either we use one of the solutions found in the literature [29] or try to find another solution in the form of trigonometric functions.Step 5.
In this step, we proceed to find a solution to i.v.p. (20) in the form of trigonometric functions. Without loss of generality, i.v.p. (20) can be redefined as(24)ℝ2=v¨+2εv˙+pv+qv3=0,v0=v0&v′0=v˙0,where i.v.p. (24) is the same as i.v.p. (20) with p=α−Q0, q=β, and vt=yt.Step 6.
Our objective here is to derive another solution that does not involve elliptic functions but an elementary solution. To do that, we assumeε>0, and for limt⟶∞vt=0, we get(25)vt=c0e−ρtcosht+arccosx0c0,where h0=0 and ht is undermined function.Step 7.
By substituting solution (25) into i.v.p. (24), we have(26)ℝ2=c0sinθe−ρt−2εh′+2ρh′−h′′+14c03qcos3θe−3ρt+14c0cosθe−3ρt3c02q+4e2ρt−2ερ−h′2+p+ρ2.Step 8.
For vanishing the coefficient ofcosθ in equation (26): 3c02q+4e2ρt−2ερ−h′2+p+ρ2=0, we have(27)h′=±124p−8ερ+4ρ2+3c02qe−2ρt,and by solving (27) with h0=0, we get(28)ht=Ht−H0,with(29)Ht=1ρΠtanh−11+3c02qe−2ρt4Π−34c02qe−2ρt+Π,H0=1ρΠtanh−11+3c02q4Π−34c02q+Π,where Π=p−2ερ+ρ2.Step 9.
The numberc0 is obtained from the condition v′0=v˙0 and it is a solution to the quartic(30)3qc04+4p−8ερ+4ρ2−3qv02c02−4pv02−2ερv02+2ρ2v02+2ρv0v˙0+v˙02=0,where the number ρ is a free/optimal parameter that is chosen in order to minimize the residual error. Its default value is ρ=ε.Step 10.
Finally, the trigonometric approximation to i.v.p. (10) is obtained:(31)xt=yft=c0e−ρftcoshft+arccosx0c0.Step 11.
Also, we can solve i.v.p. (20) using RK numerical method and then replacing t⟶ft (given in equation (23)). The following MATHEMATICA command is introduced for this purpose:(32)gt_≔y′′t+2εy′t+α−Q0yt+βyt3;RKt_≔NDSolvegt==0&&y0==x0&&y′0==x˙0,yt,t1,1,2,(33)xt_≔RKft.
Both analytical and numerical approximations (31) and (33) to i.v.p. (10) are, respectively, plotted against the RK numerical solution as illustrated in Figures 4 and 5. Also, at α,β,ε,Ω,Q0,x˙0=4,1,0.1,0.5,0.1,0.2, the maximum global distance of both approximations (31) and (33) is estimated for different values to x0 as(34)Ldx0=0=max0≤t≤30RK−xtApprox.26=0.000506901,Ldx0=0=max0≤t≤30RK−xtApprox.28=0.000482946,Ldx0=π6=max0≤t≤30RK−xtApprox.26=0.00415468,Ldx0=π6=max0≤t≤30RK−xtApprox.28=0.00266014.
The obtained results show the high accuracy and efficiency of the obtained approximations (31) and (33). Moreover, these approximations are stable against long time and for arbitrary values of the physical parameters.Figure 4
Both trigonometric approximation (31) and RK numerical simulation to the damped Duffing–Mathieu problem (10) are compared with each other for different values to the initial amplitude x0.
(a)(b)Figure 5
Both RK numerical simulation (32) and RK numerical simulation (33) using the definition of t⟶ft to the damped Duffing–Mathieu problem (10) are compared with each other for different values to the initial amplitude x0.
(a)(b)


### 2.3. The Homotopy Extended Krylov–Bogoliubov–Mitropolskii Method

The homotopy extended Krylov–Bogoliubov–Mitropolskii (HKBM) method may be used for solving both conservative and nonconservative oscillators. Based on this method (more details can be found in [27, 28]), i.v.p. (10) can be redefined as(35)x¨+α−Q0x+pεx˙+Q01−cosΩtx+βx3=0,x0=x0andx′0=x˙0,0≤t≤T,where xp≡xpt indicates the solution of i.v.p. (35) while the solution of the dDMO (2) is obtained for p=1. For ω0=α−Q0, ϕt=Q01−cosΩt and α−Q0>0, i.v.p. (35) can be written in the following reduced form:(36)x¨+ω02x+pεx˙+ϕtx+βx3=0,x0=x0andx′0=x˙0,0≤t≤T.According to the HKBM method, the following ansatz solution is introduced:(37)xp=acosψ+∑n=1Npnuna,ψ+OpN+1,where each un≡una,ψ is a periodic function in ψ, and both amplitude a and phase ψ are assumed to vary with time and subject to the conditions(38)dadt≡a˙=∑n=1NpnAna+OpN+1,dψdt≡ψ˙=ω0+∑n=1Npnψna+OpN+1,where a≡at and ψ≡ψt.Inserting ansatz solution (37) and using (38) and after several tedious calculations, we can determine the unknown time-dependent functions un,ψn,An,a. To avoid the so-called secularity, we choose only the solution that does not contain cosψ nor sinψ. For N=1, we get(39)a˙=−εat,ψ˙=3βat2+4ϕt8ω0+ω0,u1a,ψ=a3βcos3ψ32ω02,and(40)xpt=acosψ+pβ32ω02a3cos3ψ.By solving system (39), we have(41)a=c0e−εt,ψ=e−2εt16εΩω08εe2εt2c1Ωω0−Q0Ωt+sinΩt+2αΩt+3βc02Ωe2εt−1.The first-order approximate solution is obtained forp=1:(42)xt=x1t=acosψ+β32ω02a3cos3ψ,where the values of a,ψ are defined in (41) while the constants c0 and c1 can be obtained from the initial conditions.The comparison between the HKBM first-order approximate solution (37) and the RK numerical simulation is reported as shown in Figures 6(a) and 6(b) for x0=0 and x0=π/6, respectively. Also, both HKBM first-order approximate solution (42) and RK numerical simulation are, respectively, compared with each other for weak ε=0.1 and strong ε=0.5 damping as illustrated in Figures 7(a) and 7(b). Furthermore, at α,β,Ω,Q0,x˙0=4,1,0.5,0.1,0.2 and for different values to x0,ε, the maximum global distance error to the HKBM first-order approximate solution (42) is estimated as(43)Ldx0=0=max0≤t≤30RK−xtHKBM37=0.00110165,Ldx0=π6=max0≤t≤30RK−xtHKBM37=0.0061929,Ldε=0.1=max0≤t≤30RK−xtHKBM37=0.00110165,Ldε=0.5=max0≤t≤30RK−xtHKBM37=0.00493896.Figure 6
Both HKBM first-order approximate solution (42) and RK numerical simulation are plotted against different values of the initial angle x0.
(a)(b)Figure 7
Both HKBM first-order approximate solution (42) and RK numerical simulation are plotted against different values of the damping parameter ε.
(a)(b)It is clear that the HKBM first-order approximate solution (42) is characterized by high accuracy and more stability at long time.


## 2.1. An Exact Solution to (Un)damped Mathieu Equation

Both undamped Mathieu equation, i.e., (1) for β=0, and the damped Mathieu equation, i.e., (2) for β=0, have exact solutions in the form of Mathieu functions. First, let us find the solution of the following damped Mathieu equation:(4)x¨+εx˙+α−Q0cosΩtx=0,with subjected to the initial conditions (ICs): x0=x0 and x′0=x˙0.Using the following MATHEMATICA command(5)gt_≔x′′t+2εx′t−α−Q0CosΩtxtNDSolvegt==0&&x0==x0&&x′0==x˙0,xt,t1,1,2//FullSimplify1,we can get the exact solution to (4) as follows:(6)xt=e−εtx0MC1t,εMC10,ε+2εx0+x˙0ΩMS1t,εMathieuSPrime4α−ε2/Ω2,2Q0/Ω2,0,with(7)MC1t,ε=MathieuC4α−ε2Ω2,2Q0Ω2,Ω2t,MS1t,ε=MathieuS4α−ε2Ω2,2Q0Ω2,Ω2t,where MathieuS and MathieuC are the Mathieu functions of the first kind or sometimes called sine-elliptic and cosine-elliptic, respectively. For ε=0, the damped Mathieu equation (4) reduces to the following undamped Mathieu equation:(8)x¨−α+Q0cosΩtx=0,and solution (6) reduces to the following one:(9)xt=x0MC1t,0MC10,0+2x˙0MS1t,0ΩMathieuSPrime4α/Ω2,2Q0/Ω2,0.The exact solution (6) is compared with RK numerical solution as shown in Figures 1(a) and 1(b) for different values of x0.Figure 1
Both exact solution (6) and RK numerical solution to the damped Mathieu equation (4) are compared with each other for different values to the initial amplitude x0.
(a)(b)


## 2.2. Some Analytical Approximations to the Damped Duffing–Mathieu Oscillator

Here, we proceed to discuss two techniques (the hybridp−expansion method and the ansatz method) for finding some analytical approximations to dDMO (2). For studying dDMO (2), first we rewrite equation (2) in the following new i.v.p.(10)ℝ1=0,x0=x0&x′0=x˙0.


### 2.2.1. First Approach: Reduced Method

Forβ≠0 and α>0, we may obtain simple approximation to the i.v.p. (10) by replacing the cubic term βx3 by the linear term βκx, where the constant κ≥0 which is used as an optimal parameter to reduce the residual error. Accordingly, dDMO (2) of cubic nonlinearity reduces to the following dDMO with linear term. Thus, we can replace the i.v.p. (10) by the following new i.v.p.(11)x¨+2εx˙+α−Q0cosΩtx+βκx=0,x0=x0andx′0=x˙0.Thus, the exact solution to the i.v.p. (11) is expressed by(12)xκ≡xκt=e−εtx0MC2t,ε,κMC20,ε,κ+2εx0+x˙0ΩMS2t,ε,κMathieuSPrime4βκ−ε2+α/Ω2,2Q0/Ω2,0,with(13)MC2t,ε,κ=MathieuC4βκ−ε2+αΩ2,2Q0Ω2,Ω2t,MS2t,ε,κ=MathieuS4βκ−ε2+αΩ2,2Q0Ω2,Ω2t.The residual is defined as(14)Rκt=x¨κ+εx˙κ+α−Q0cosΩtxκ+βxκ3.A suitable value ofκ can be obtained by solving the equation Rκt0=0 for some t0>0, say t0=1. Making use of the Padé approximate technique, for 0<t0≤1 and x0≈0, we get(15)κsuitable=Y1Y2,with(16)Y1=3x02Q0t0x02−αt0x02−2εt0x˙0x0−4t0x˙02−2x˙0x0×4Q0t0x02−4αt0x02+3βt0x04−8εt0x˙0x0−6t0x˙02−8x˙0x0,Y2=−26αQ0t02x04+12βQ0t02x06−52εQ0t02x˙0x03+13Q02t02x04−52Q0t0x˙0x03−30Q0t02x˙02x02+13α2t02x04−12αβt02x06+52αεt02x˙0x03+52αt0x˙0x03+30αt02x˙02x02+9β2t02x08−24βεt02x˙0x05−24βt0x˙0x05+12βt02x˙02x04+52ε2t02x˙02x02+104εt0x˙02x02+60εt02x˙03x0+60t0x˙03x0+12t02x˙04+52x˙02x02.For a givenκ, the residual error LRκ of the approximation (12) to the i.v.p. (10) is defined as(17)LRκ=max0≤t≤TRκt=max0≤t≤Tx¨κ+2εx˙κ+α−Q0cosΩtxκ+βxκ3.The optimal value ofκoptimal for the parameter κ on the interval 0≤t≤T is defined as(18)κoptimal=minκ≥0LRκ.Let us apply the obtained approximation (12) for investigating the properties of the damping oscillations to the dDMO (2) at different values to the physical parameters α,β,ε,Ω,Q0,x0,x˙0. The profile of the approximation (12) using the values of κsuitable (given in equation (15)) and κoptimal (given in (18)) is displayed in Figures 2 and 3 at α,β,ε,Ω,Q0,x0,x˙0=4,1,0.1,0.5,0.1,0,0.2 and α,β,ε,Ω,Q0,x0,x˙0=4,1,0.1,0.5,0.1,π/6,0.2, respectively. The obtained results showed that this approximation gives results with good and acceptable accuracy.Figure 2
The profile of the approximation (12) using the values of κsuitable (given in equation (15)) and κoptimal (given in (18)) is plotted against t for α,β,ε,Ω,Q0,x0,x˙0=4,1,0.1,0.5,0.1,0,0.2.
(a)(b)Figure 3
The profile of the approximation (12) using the values of κsuitable (given in equation (15)) and κoptimal (given in (18)) is plotted against t for α,β,ε,Ω,Q0,x0,x˙0=4,1,0.1,0.5,0.1,π/6,0.2.
(a)(b)


### 2.2.2. The Ansatz Method for Solving dDMO

Now, we can summarize the main points to get some approximations to the i.v.p. (10) in the following steps.Step 1.
Let us assume the following ansatz:(19)xt=yft,where the time-dependent function f≡ft can be determined later and the function y≡yt represents the solution of the following i.v.p.(20)y¨+2εy˙+α−Q0y+βy3=0,y0=x0&y′0=x˙0.Step 2.
Inserting ansatz (19) into i.v.p. (10), we get(21)ℝ1=y′f−2εf′2+2εf′+f′′+yfα−Q0cosΩt+f′2Q0−α+βy3f1−f′2.Step 3.
For vanishing the coefficient ofyf in equation (21), we get(22)f′2=−α+Q0cosΩtQ0−α.Step 4.
Integrating equation (22) and with the help of f0=0, we have(23)ft=2ΩEΩ2t,2Q0Q0−α,where E stands for the EllipticE function.
Now, the value offt has been determined but the solution of the i.v.p. (20) needs to be determined. Thus, we are faced with two things: either we use one of the solutions found in the literature [29] or try to find another solution in the form of trigonometric functions.Step 5.
In this step, we proceed to find a solution to i.v.p. (20) in the form of trigonometric functions. Without loss of generality, i.v.p. (20) can be redefined as(24)ℝ2=v¨+2εv˙+pv+qv3=0,v0=v0&v′0=v˙0,where i.v.p. (24) is the same as i.v.p. (20) with p=α−Q0, q=β, and vt=yt.Step 6.
Our objective here is to derive another solution that does not involve elliptic functions but an elementary solution. To do that, we assumeε>0, and for limt⟶∞vt=0, we get(25)vt=c0e−ρtcosht+arccosx0c0,where h0=0 and ht is undermined function.Step 7.
By substituting solution (25) into i.v.p. (24), we have(26)ℝ2=c0sinθe−ρt−2εh′+2ρh′−h′′+14c03qcos3θe−3ρt+14c0cosθe−3ρt3c02q+4e2ρt−2ερ−h′2+p+ρ2.Step 8.
For vanishing the coefficient ofcosθ in equation (26): 3c02q+4e2ρt−2ερ−h′2+p+ρ2=0, we have(27)h′=±124p−8ερ+4ρ2+3c02qe−2ρt,and by solving (27) with h0=0, we get(28)ht=Ht−H0,with(29)Ht=1ρΠtanh−11+3c02qe−2ρt4Π−34c02qe−2ρt+Π,H0=1ρΠtanh−11+3c02q4Π−34c02q+Π,where Π=p−2ερ+ρ2.Step 9.
The numberc0 is obtained from the condition v′0=v˙0 and it is a solution to the quartic(30)3qc04+4p−8ερ+4ρ2−3qv02c02−4pv02−2ερv02+2ρ2v02+2ρv0v˙0+v˙02=0,where the number ρ is a free/optimal parameter that is chosen in order to minimize the residual error. Its default value is ρ=ε.Step 10.
Finally, the trigonometric approximation to i.v.p. (10) is obtained:(31)xt=yft=c0e−ρftcoshft+arccosx0c0.Step 11.
Also, we can solve i.v.p. (20) using RK numerical method and then replacing t⟶ft (given in equation (23)). The following MATHEMATICA command is introduced for this purpose:(32)gt_≔y′′t+2εy′t+α−Q0yt+βyt3;RKt_≔NDSolvegt==0&&y0==x0&&y′0==x˙0,yt,t1,1,2,(33)xt_≔RKft.
Both analytical and numerical approximations (31) and (33) to i.v.p. (10) are, respectively, plotted against the RK numerical solution as illustrated in Figures 4 and 5. Also, at α,β,ε,Ω,Q0,x˙0=4,1,0.1,0.5,0.1,0.2, the maximum global distance of both approximations (31) and (33) is estimated for different values to x0 as(34)Ldx0=0=max0≤t≤30RK−xtApprox.26=0.000506901,Ldx0=0=max0≤t≤30RK−xtApprox.28=0.000482946,Ldx0=π6=max0≤t≤30RK−xtApprox.26=0.00415468,Ldx0=π6=max0≤t≤30RK−xtApprox.28=0.00266014.
The obtained results show the high accuracy and efficiency of the obtained approximations (31) and (33). Moreover, these approximations are stable against long time and for arbitrary values of the physical parameters.Figure 4
Both trigonometric approximation (31) and RK numerical simulation to the damped Duffing–Mathieu problem (10) are compared with each other for different values to the initial amplitude x0.
(a)(b)Figure 5
Both RK numerical simulation (32) and RK numerical simulation (33) using the definition of t⟶ft to the damped Duffing–Mathieu problem (10) are compared with each other for different values to the initial amplitude x0.
(a)(b)


## 2.2.1. First Approach: Reduced Method

Forβ≠0 and α>0, we may obtain simple approximation to the i.v.p. (10) by replacing the cubic term βx3 by the linear term βκx, where the constant κ≥0 which is used as an optimal parameter to reduce the residual error. Accordingly, dDMO (2) of cubic nonlinearity reduces to the following dDMO with linear term. Thus, we can replace the i.v.p. (10) by the following new i.v.p.(11)x¨+2εx˙+α−Q0cosΩtx+βκx=0,x0=x0andx′0=x˙0.Thus, the exact solution to the i.v.p. (11) is expressed by(12)xκ≡xκt=e−εtx0MC2t,ε,κMC20,ε,κ+2εx0+x˙0ΩMS2t,ε,κMathieuSPrime4βκ−ε2+α/Ω2,2Q0/Ω2,0,with(13)MC2t,ε,κ=MathieuC4βκ−ε2+αΩ2,2Q0Ω2,Ω2t,MS2t,ε,κ=MathieuS4βκ−ε2+αΩ2,2Q0Ω2,Ω2t.The residual is defined as(14)Rκt=x¨κ+εx˙κ+α−Q0cosΩtxκ+βxκ3.A suitable value ofκ can be obtained by solving the equation Rκt0=0 for some t0>0, say t0=1. Making use of the Padé approximate technique, for 0<t0≤1 and x0≈0, we get(15)κsuitable=Y1Y2,with(16)Y1=3x02Q0t0x02−αt0x02−2εt0x˙0x0−4t0x˙02−2x˙0x0×4Q0t0x02−4αt0x02+3βt0x04−8εt0x˙0x0−6t0x˙02−8x˙0x0,Y2=−26αQ0t02x04+12βQ0t02x06−52εQ0t02x˙0x03+13Q02t02x04−52Q0t0x˙0x03−30Q0t02x˙02x02+13α2t02x04−12αβt02x06+52αεt02x˙0x03+52αt0x˙0x03+30αt02x˙02x02+9β2t02x08−24βεt02x˙0x05−24βt0x˙0x05+12βt02x˙02x04+52ε2t02x˙02x02+104εt0x˙02x02+60εt02x˙03x0+60t0x˙03x0+12t02x˙04+52x˙02x02.For a givenκ, the residual error LRκ of the approximation (12) to the i.v.p. (10) is defined as(17)LRκ=max0≤t≤TRκt=max0≤t≤Tx¨κ+2εx˙κ+α−Q0cosΩtxκ+βxκ3.The optimal value ofκoptimal for the parameter κ on the interval 0≤t≤T is defined as(18)κoptimal=minκ≥0LRκ.Let us apply the obtained approximation (12) for investigating the properties of the damping oscillations to the dDMO (2) at different values to the physical parameters α,β,ε,Ω,Q0,x0,x˙0. The profile of the approximation (12) using the values of κsuitable (given in equation (15)) and κoptimal (given in (18)) is displayed in Figures 2 and 3 at α,β,ε,Ω,Q0,x0,x˙0=4,1,0.1,0.5,0.1,0,0.2 and α,β,ε,Ω,Q0,x0,x˙0=4,1,0.1,0.5,0.1,π/6,0.2, respectively. The obtained results showed that this approximation gives results with good and acceptable accuracy.Figure 2
The profile of the approximation (12) using the values of κsuitable (given in equation (15)) and κoptimal (given in (18)) is plotted against t for α,β,ε,Ω,Q0,x0,x˙0=4,1,0.1,0.5,0.1,0,0.2.
(a)(b)Figure 3
The profile of the approximation (12) using the values of κsuitable (given in equation (15)) and κoptimal (given in (18)) is plotted against t for α,β,ε,Ω,Q0,x0,x˙0=4,1,0.1,0.5,0.1,π/6,0.2.
(a)(b)


## 2.2.2. The Ansatz Method for Solving dDMO

Now, we can summarize the main points to get some approximations to the i.v.p. (10) in the following steps.Step 1.
Let us assume the following ansatz:(19)xt=yft,where the time-dependent function f≡ft can be determined later and the function y≡yt represents the solution of the following i.v.p.(20)y¨+2εy˙+α−Q0y+βy3=0,y0=x0&y′0=x˙0.Step 2.
Inserting ansatz (19) into i.v.p. (10), we get(21)ℝ1=y′f−2εf′2+2εf′+f′′+yfα−Q0cosΩt+f′2Q0−α+βy3f1−f′2.Step 3.
For vanishing the coefficient ofyf in equation (21), we get(22)f′2=−α+Q0cosΩtQ0−α.Step 4.
Integrating equation (22) and with the help of f0=0, we have(23)ft=2ΩEΩ2t,2Q0Q0−α,where E stands for the EllipticE function.
Now, the value offt has been determined but the solution of the i.v.p. (20) needs to be determined. Thus, we are faced with two things: either we use one of the solutions found in the literature [29] or try to find another solution in the form of trigonometric functions.Step 5.
In this step, we proceed to find a solution to i.v.p. (20) in the form of trigonometric functions. Without loss of generality, i.v.p. (20) can be redefined as(24)ℝ2=v¨+2εv˙+pv+qv3=0,v0=v0&v′0=v˙0,where i.v.p. (24) is the same as i.v.p. (20) with p=α−Q0, q=β, and vt=yt.Step 6.
Our objective here is to derive another solution that does not involve elliptic functions but an elementary solution. To do that, we assumeε>0, and for limt⟶∞vt=0, we get(25)vt=c0e−ρtcosht+arccosx0c0,where h0=0 and ht is undermined function.Step 7.
By substituting solution (25) into i.v.p. (24), we have(26)ℝ2=c0sinθe−ρt−2εh′+2ρh′−h′′+14c03qcos3θe−3ρt+14c0cosθe−3ρt3c02q+4e2ρt−2ερ−h′2+p+ρ2.Step 8.
For vanishing the coefficient ofcosθ in equation (26): 3c02q+4e2ρt−2ερ−h′2+p+ρ2=0, we have(27)h′=±124p−8ερ+4ρ2+3c02qe−2ρt,and by solving (27) with h0=0, we get(28)ht=Ht−H0,with(29)Ht=1ρΠtanh−11+3c02qe−2ρt4Π−34c02qe−2ρt+Π,H0=1ρΠtanh−11+3c02q4Π−34c02q+Π,where Π=p−2ερ+ρ2.Step 9.
The numberc0 is obtained from the condition v′0=v˙0 and it is a solution to the quartic(30)3qc04+4p−8ερ+4ρ2−3qv02c02−4pv02−2ερv02+2ρ2v02+2ρv0v˙0+v˙02=0,where the number ρ is a free/optimal parameter that is chosen in order to minimize the residual error. Its default value is ρ=ε.Step 10.
Finally, the trigonometric approximation to i.v.p. (10) is obtained:(31)xt=yft=c0e−ρftcoshft+arccosx0c0.Step 11.
Also, we can solve i.v.p. (20) using RK numerical method and then replacing t⟶ft (given in equation (23)). The following MATHEMATICA command is introduced for this purpose:(32)gt_≔y′′t+2εy′t+α−Q0yt+βyt3;RKt_≔NDSolvegt==0&&y0==x0&&y′0==x˙0,yt,t1,1,2,(33)xt_≔RKft.
Both analytical and numerical approximations (31) and (33) to i.v.p. (10) are, respectively, plotted against the RK numerical solution as illustrated in Figures 4 and 5. Also, at α,β,ε,Ω,Q0,x˙0=4,1,0.1,0.5,0.1,0.2, the maximum global distance of both approximations (31) and (33) is estimated for different values to x0 as(34)Ldx0=0=max0≤t≤30RK−xtApprox.26=0.000506901,Ldx0=0=max0≤t≤30RK−xtApprox.28=0.000482946,Ldx0=π6=max0≤t≤30RK−xtApprox.26=0.00415468,Ldx0=π6=max0≤t≤30RK−xtApprox.28=0.00266014.
The obtained results show the high accuracy and efficiency of the obtained approximations (31) and (33). Moreover, these approximations are stable against long time and for arbitrary values of the physical parameters.Figure 4
Both trigonometric approximation (31) and RK numerical simulation to the damped Duffing–Mathieu problem (10) are compared with each other for different values to the initial amplitude x0.
(a)(b)Figure 5
Both RK numerical simulation (32) and RK numerical simulation (33) using the definition of t⟶ft to the damped Duffing–Mathieu problem (10) are compared with each other for different values to the initial amplitude x0.
(a)(b)


## 2.3. The Homotopy Extended Krylov–Bogoliubov–Mitropolskii Method

The homotopy extended Krylov–Bogoliubov–Mitropolskii (HKBM) method may be used for solving both conservative and nonconservative oscillators. Based on this method (more details can be found in [27, 28]), i.v.p. (10) can be redefined as(35)x¨+α−Q0x+pεx˙+Q01−cosΩtx+βx3=0,x0=x0andx′0=x˙0,0≤t≤T,where xp≡xpt indicates the solution of i.v.p. (35) while the solution of the dDMO (2) is obtained for p=1. For ω0=α−Q0, ϕt=Q01−cosΩt and α−Q0>0, i.v.p. (35) can be written in the following reduced form:(36)x¨+ω02x+pεx˙+ϕtx+βx3=0,x0=x0andx′0=x˙0,0≤t≤T.According to the HKBM method, the following ansatz solution is introduced:(37)xp=acosψ+∑n=1Npnuna,ψ+OpN+1,where each un≡una,ψ is a periodic function in ψ, and both amplitude a and phase ψ are assumed to vary with time and subject to the conditions(38)dadt≡a˙=∑n=1NpnAna+OpN+1,dψdt≡ψ˙=ω0+∑n=1Npnψna+OpN+1,where a≡at and ψ≡ψt.Inserting ansatz solution (37) and using (38) and after several tedious calculations, we can determine the unknown time-dependent functions un,ψn,An,a. To avoid the so-called secularity, we choose only the solution that does not contain cosψ nor sinψ. For N=1, we get(39)a˙=−εat,ψ˙=3βat2+4ϕt8ω0+ω0,u1a,ψ=a3βcos3ψ32ω02,and(40)xpt=acosψ+pβ32ω02a3cos3ψ.By solving system (39), we have(41)a=c0e−εt,ψ=e−2εt16εΩω08εe2εt2c1Ωω0−Q0Ωt+sinΩt+2αΩt+3βc02Ωe2εt−1.The first-order approximate solution is obtained forp=1:(42)xt=x1t=acosψ+β32ω02a3cos3ψ,where the values of a,ψ are defined in (41) while the constants c0 and c1 can be obtained from the initial conditions.The comparison between the HKBM first-order approximate solution (37) and the RK numerical simulation is reported as shown in Figures 6(a) and 6(b) for x0=0 and x0=π/6, respectively. Also, both HKBM first-order approximate solution (42) and RK numerical simulation are, respectively, compared with each other for weak ε=0.1 and strong ε=0.5 damping as illustrated in Figures 7(a) and 7(b). Furthermore, at α,β,Ω,Q0,x˙0=4,1,0.5,0.1,0.2 and for different values to x0,ε, the maximum global distance error to the HKBM first-order approximate solution (42) is estimated as(43)Ldx0=0=max0≤t≤30RK−xtHKBM37=0.00110165,Ldx0=π6=max0≤t≤30RK−xtHKBM37=0.0061929,Ldε=0.1=max0≤t≤30RK−xtHKBM37=0.00110165,Ldε=0.5=max0≤t≤30RK−xtHKBM37=0.00493896.Figure 6
Both HKBM first-order approximate solution (42) and RK numerical simulation are plotted against different values of the initial angle x0.
(a)(b)Figure 7
Both HKBM first-order approximate solution (42) and RK numerical simulation are plotted against different values of the damping parameter ε.
(a)(b)It is clear that the HKBM first-order approximate solution (42) is characterized by high accuracy and more stability at long time.


## 3. Conclusion

Given the importance of nonlinear oscillations in plasma physics and engineering and their strong connection to the family of the Duffing-type oscillator, in this work, some exact solutions to the damped and undamped Mathieu equations as well as some analytical approximations to the damped Duffing–Mathieu oscillator (dDMO) using different approaches have been obtained. The exact solutions to both damped and undamped Mathieu equation have been obtained in the terms of Mathieu functions of the first kind. These solutions are numerically compared with the Runge–Kutta (RK) numerical simulation. It was observed that both exact and numerical solutions are completely matched with each other in the whole time interval. On the other hand, the dDMO has been solved using some different approaches. In the first one, the nonintegrable dDMO with cubic nonlinear termβx3 has been reduced to an integrable dDMO with linear term βκx in which κ is undermined optimal parameter. The kappa optimal parameter κ has been determined using a suitable technique as we discussed in the text above. After determining the kappa optimal parameter, a highly accurate analytical approximation has been obtained in terms of the Mathieu functions. In the second approach, a highly accurate analytical approximation has been derived in detail in terms of trigonometric functions using the ansatz method. In the third technique, the homotopy extended Krylov–Bogoliubov–Mitropolskii (HKBM) method was used for getting an effective analytical approximation. Furthermore, the dDMO has been analyzed numerically using the RK numerical method. The comparison between all obtained approximations and the RK numerical solutions has been carried out. Moreover, the maximum global distance error in the whole time interval to all obtained approximations has been estimated. All obtained approximations are characterized by the high accuracy and efficiency in addition to being more stable for a long time.


### 3.1. Future Work

We may solve the following oscillators using of the methods described in this paper:


#### 3.1.1. Future Idea I

Cubic-quintic Duffing–Mathieu equation:(44)x¨+ω02x+2εx˙+ϕtx+βx3+γx5=0,x0=x0&x′0=x0,0≤t≤T.


#### 3.1.2. Future Idea II

Forced damped Duffing–Mathieu equation:(45)x¨+ω02x+2εx˙+ϕtx+βx3=Ft,x0=x0&x′0=x0,0≤t≤T.


#### 3.1.3. Future Idea III

The forced Van der Pol–Duffing oscillator:(46)x¨−ε1−x2x˙+ω02x+βx3=Ft,x0=x0&x′0=x0,0≤t≤T,and many others oscillators.


## 3.1. Future Work

We may solve the following oscillators using of the methods described in this paper:


### 3.1.1. Future Idea I

Cubic-quintic Duffing–Mathieu equation:(44)x¨+ω02x+2εx˙+ϕtx+βx3+γx5=0,x0=x0&x′0=x0,0≤t≤T.


### 3.1.2. Future Idea II

Forced damped Duffing–Mathieu equation:(45)x¨+ω02x+2εx˙+ϕtx+βx3=Ft,x0=x0&x′0=x0,0≤t≤T.


### 3.1.3. Future Idea III

The forced Van der Pol–Duffing oscillator:(46)x¨−ε1−x2x˙+ω02x+βx3=Ft,x0=x0&x′0=x0,0≤t≤T,and many others oscillators.


## 3.1.1. Future Idea I

Cubic-quintic Duffing–Mathieu equation:(44)x¨+ω02x+2εx˙+ϕtx+βx3+γx5=0,x0=x0&x′0=x0,0≤t≤T.


## 3.1.2. Future Idea II

Forced damped Duffing–Mathieu equation:(45)x¨+ω02x+2εx˙+ϕtx+βx3=Ft,x0=x0&x′0=x0,0≤t≤T.


## 3.1.3. Future Idea III

The forced Van der Pol–Duffing oscillator:(46)x¨−ε1−x2x˙+ω02x+βx3=Ft,x0=x0&x′0=x0,0≤t≤T,and many others oscillators.



---

*Source: 2715767-2022-05-05.xml*