# Some Novel Approaches for Analyzing the Unforced and Forced Duffing–Van der Pol Oscillators

Authors: Alvaro H. Salas; Wedad Albalawi; S. A. El-Tantawy; L. S. El-Sherif
Journal: Journal of Mathematics
 (2022)

---

## Abstract

In the current investigation, both unforced and forced Duffing–Van der Pol oscillator (DVdPV) oscillators with a strong nonlinearity and external periodic excitations are analyzed and investigated analytically and numerically using some new and improved approaches. The new approach is constructed based on Krylov–Bogoliubov–Metroolsky method (KBMM). One of the most important features of this approach is that we do not need to solve a system of differential equations, but only solve a system of algebraic equations. Moreover, the ease and faster of applying this method gives high-accurate results and this approach is better than many approaches in the literature. This approach is applied for analyzing (un)forced DVdP oscillators. Also, some improvements are made to He’s frequency-amplitude formulation in order to solve unforced DVdP oscillator to obtain high-accurate results. Furthermore, the He’s homotopy perturbation method (He’s HPM) is employed for analyzing unforced DVdP oscillator. The comparison between all mentioned approaches is carried out. The application of our approach is not limited to (un)forced DVdPV oscillators only but can be applied to analyze many higher-order nonlinearity oscillators for any odd power and it gives more accurate results than other approaches. Both used methods and obtained approximations will help many researchers in general and plasma physicists in particular in the interpretation of their results.


---

## Body

## 1. Introduction

The study of the dynamics of nonlinear oscillators is one of the topics of great importance for many researchers due to its many important applications in various areas of physics, applied mathematics, and practical engineering [1–17]. Differential equations are one of the best and most successful models in modeling many nonlinear dynamical systems. For instance, Duffing-type equation is one of the most famous and successful equations that has been used for modeling and interpreting many nonlinear oscillations in many different dynamical systems such as electrical circuit, optical stability, the buckled beam, and different oscillations in a plasma [16–19]. In plasma physics, there are many evolution equations that can be reduced to Duffing-type equation, Helmholtz-type equation, Duffing–Helmhlotz equation, and Mathieu equation in order to investigate the various oscillations that occur within complicated plasma systems [20–23]. There is another type of equation of motion that was used for modeling the nonlinear oscillations in biology, electronics, engineering, plasma physics, and chemistry which is called Van der Pol–Duffing (VdPD) (sometimes called Duffing–Van der Pol (DVdP)) equation and its family [24, 25]. For example, a forced modified VdPD (mVdPD) equation was adopted for investigating the strong nonlinear oscillations in plasma [26]. Also, a mVdPD equation with asymmetric potential was used for modeling the nonlinear chemical dynamics [27]. Many numerical and analytical approaches were applied for solving the second-order nonlinear oscillator equations. For example, both HBM and MTSs techniques were devoted for analyzing a forced Van der Pol (VdP) generalized oscillator to obtain the amplitudes of the forced harmonic, superharmonic, and subharmonic oscillatory states [26]. Melnikov’s method was used for analyzing a mVdPD equation to derive analytical criteria for the appearance of horseshoe chaos in chemical oscillations [27]. He et al. [16] used the Poincare ´–Lindstedt technique (PLT) for solving and analyzing the Hybrid Rayleigh–van der pol–Duffing equation. Also, the homotopy analysis method (HAM) was used for analyzing DVdP oscillator [28]. Both methods of differentiable dynamics and Lie symmetry reduction method were devoted for analyzing the DVdP-type oscillator [29]. Moreover, DVdP oscillator was solved numerically via Adomian’s decomposition method (ADM) [30]. Based on this approach, the equation of motion is converted to a system of first-order differential equations and then was solved to obtain a numerical approximation. Moreover, the authors made a comparison with Lindsted’s method (LM) approximation. They found that the obtained approximation using ADM is better than LM. However, in the two approaches, the approximations become convergence and more accurate only in the short time domain but these approximations become dis-convergence and not accurate for long time domain. Most methods in the literature lead to complicated formulas for the obtained approximations and the analysis of such solutions are much difficult or not convergence for a long time. However, the Krylov–Bogoliubov–Mitropolsky method (KBMM) was adopted for deriving the periodic steady-state solutions to the following DVdP driven oscillator [6].(1)x¨−ε1−x2x˙+ω02x+βx3=fcosΩt,x0=x0,x′0=x˙0,where the overdot indicates the derivative with respect to “t”. Recently, Salas et al. [11] applied the ansatz method, HBM, PLT, and KBMM for analyzing the forced VdP oscillator and found that KBMM gives more accurate approximations. Motivated by the investigations in Ref. [6, 11], we proceed to analyze the DVdP oscillator using a new effective and simplification technique based on KBMM. In our approach, we will prove that the new approach does not demand to solve any ordinary differential equations (odes). Moreover, we will prove that the new suggested approach gives high-accurate and convergence approximations in the whole time domain. Note that for ε=0, i.v.p. (1) reduces to the forced Duffing oscillator whose general solution is well known [11]. Moreover, in this investigation, we try to improve He’s frequency-amplitude formulation to be suitable for analyzing the DVdP oscillator. Also, the He’s homotopy perturbation method (He’s HPM) will apply for analyzing and investigating the DVdP oscillator.The rest of this paper is introduced in the following fashion: in Section2, the new suggested approach is introduced. The analytical approximations to the unforced DVdP oscillator is reported in Section 3 using the new mentioned approach, the He’s HPM, and improved He’s frequency-amplitude formulation. Moreover, in Section 4, the new mentioned approach is devoted for getting an analytical approximation to the forced DVdP oscillator. The obtained results are summarized in Section 5.


## 2. New Approach Based on KBM for Solving Strongly Nonlinear Oscillators

Let us consider the following general form to the second-order i.v.p.:(2)x¨+ω02x+Fx,x˙=0,x0=x0,x′0=x˙0,where the expression Fx,0 is an odd polynomial of x.To introduce ap-parameter solution to i.v.p. (2), we rewrite this problem in the following new form:(3)Rpx≡x¨+ω02x+pFx,x˙=0,x0=x0,x′0=x˙0,where xp≡xpt≡xp,t indicates the solution to i.v.p. (3), we can call a p-parameter solution. Then, the solution to the original i.v.p. (2) can be obtained for p=1.Assuming the solution to i.v.p. (3) is given by the following ansatz form:(4)xp=acosψ+∑k=1Npk∑j=1ka2j+1∑i=1jr2i+1,j,kcos2i+1ψ+s2i+1,j,ksin2i+1ψ,where the functions a≡at and ψ≡ψt are assumed to vary with time according to the following ordinary differential equations (odes):(5)a˙=dadt≡a′t=∑k=1Npk∑j=kk+1Bk,ja2j−1s,(6)ψ˙=dψdt≡ψ′t=ω0+pC1,0a2+∑k=2Npk∑j=0kCk,ja2j.We plug the ansatz (4) with the relations (5) and (6) into Rpx given in equation (3) and equating the coefficients of pn (n=1,2,3,…,) a2j+1t, cos 2i+1ψt, and sin 2i+1ψt (i, j=1,2,3,…,N) to zero, then we can get a system of algebraic equations. By solving this system, we can determine the unknown coefficients r2i+1,j,k, s2i+1,j,k, Bk,j, and Ck,j.Remark 1.
We can obtain another method by replacing (6) with(7)ψ˙=ω02+pC1,0a2+∑k=2Npk∑j=0kCk,ja2j.
In the below section, we will use this method for analyzing both the unforced DVdP oscillator, i.e., i.v.p. (1) for f=0 and the forced DVdP oscillator (1).


## 3. Analytical Approximations to the Unforced Duffing–Van Der Pol Oscillator

Here, we can analyze the unforced DVdP oscillator, i.e., i.v.p. (1) for f=0 using three different approaches, namely:(i)
Our new approach mentioned in the above section(ii)
The He’s HPM(iii)
Improved He’s frequency-amplitude formulation (He’s FAF)


### 3.1. Our New Approach

By usingN=2 in solution (4), then the solution to the following unforced problem (3)(8)x¨+ω02x+p−ε1−x2x˙+βx3=0,can be introduced in following form:(9)xp=acosψ+∑k=12a2k+1pk∑j=1kr2j+1,kcos2j+1ψ+s2j+1,ksin2j+1ψ.Accordingly, we get(10)xp=acosψ+pa3A1t+p2a3A2t+a5A3t,with(11)A1t=r3,1,1cos3ψ+s3,1,1sin3ψ,A2t=r3,1,2cos3ψ+s3,1,2sin3ψ,A3t=r3,2,2cos3ψ+r5,2,2cos5ψ+s3,2,2sin3ψ+s5,2,2sin5ψ,(12)a˙=pB1,2at3+B1,1at+p2B2,3at5+B2,2at3,ψ˙=ω0+pC1,0at2+p2C2,2at4+C2,1at2+C2,0.On the other hand(13)Rp=D1p+D2p2+Op3,with(14)D1=S1asinψ+a3S2sinψ+S3cosψ+S4cos3ψ+S5sin3ψ,D2=a5S6sinψ+S7cosψ+S8sin3ψ+S9cos3ψ+S10sin5ψ+S11cos5ψ+a3S12sinψ+S13cosψ+S14sin3ψ+S15cos3ψ+aS16cosψ,where the values of coefficients Sii=1,2,…,16 are defined in Appendix (i).Equating the coefficients ofp and p2 to zero in(15)Rp=x¨+ω02x+p−ε1−x2x˙+βx3,also, equating to zero the coefficients of a2j+1, cos2i+1ψ, and sin2i+1ψ, where (i, j=1,2,3,…,N), i.e., Si=0, we get an algebraic system. The solution of this system yields(16)B1,1=ε2,B1,2=−ε8,B2,2=−3βε16ω02,B2,3=βε32ω02,C1,0=3β8ω0,C2,0=−ε28ω0,C2,1=ε28ω0,C2,2=−15β2+7ε2ω02256ω03r3,1,1=β32ω02,r3,1,2=−ε2128ω02,r3,2,2=−21β2+ε2ω021024ω04,r5,2,2=3β2−5ε2ω023072ω04,s3,1,1=−ε32ω0,s3,1,2=−3βε128ω03,s3,2,2=3βε256ω03,s5,2,2=−βε384ω03.From the above values in equation (11), we have(17)A1t=132ω02βcos3ψ−εω0sin3ψ,A2t=−ε128ω033βsin3ψ+εω0cos3ψ,A3t=−13072ω04321β2+ε2ω02cos3ψ+5ε2ω02−3β2cos5ψ+4βεω02sin5ψ−9sin3ψ.Accordingly, the odes for determining the functionsa,ψ read(18)a˙=p12εat−18εat3+p2−3βε16ω02at3+βε32ω02at5,ψ˙=ω0+p3β8ω0at2+p2−ε28ω0+ε28ω0at2−15β2+7ε2ω02256ω03at4.Forp=1, the value of a˙ given in equation (18) reduces to(19)a˙|p=1=12εa−18εa3+−3βε16ω02a3+βε32ω02a5.The amplitude for the limit cycle is obtained from the conditiona˙|p=1=0:(20)12εa−18εa3+−3βε16ω02a3+βε32ω02a5=0.Solving equation (20) gives(21)A=ω022−β9β−4ω02/ω04+4β+3.Observe that(22)limβ⟶0A=2,this is called the cycle amplitude for the VdP oscillator.We can use the following Chebyshev approximation in order to facilitate the solution to the ode system (18):(23)a˙|p=1=12εa−18εa3+−3βε16ω02a3+βε32ω02a5≈116ε8−βω02a−116ε2+βω02a3,=κa−μa3,with(24)κ=116ε8−βω02,μ=116ε2+βω02.By solving equation (23), we get(25)a=c0κeκtc02μe2κt−1+κ.Also, the expression for determiningψ can be obtained from the second equation in (18) for p=1 whose solution reads(26)ψ=∑i=14Wi+c1,where the values of Wii=1,2,3,4 are defined in Appendix (ii). The constants c0 and c1 are determined from the initial conditions (ICs) x0=x0 and x′0=x˙0.In all above expressions, forp=1, the approximation to the following i.v.p. can be obtained:(27)x¨−ε1−x2x˙+ω02x+βx3=0,x0=x0andx′0=x˙0,


### 3.2. He’s Homotopy Perturbation Method

Moreover, the approximate solution to the DVdP i.v.p. (1) using He’s HPM is obtained. Briefly, He’s HPM can be used for a series of nonlinear oscillators differential equations which many classical perturbation methods failed to solve them or to give some accurate solutions. This method suggests the solution in the following ansatz:(28)xLP=∑i=0∞pixiωt,with(29)ω2=ω02+∑i=1∞piωi.Substitute equations (28) and (29) into i.v.p. (1), and by collecting the coefficients of same powers of p, we finally obtain some reduced equations. We have(30)ℝF=x0τ+x0″τω0+βx0τ3+ω02x1τ+12−2εω0+τω1x0′τ+εω0x0τ2x0′τ+ω1x0″τ+ω02x1″τ+12τω1x03τ−FcosτΩω0p+Op2,where τ=tω.Equating to zero the coefficients ofpj and solving the resulting odes gives(31)x0τ=Acosτ,x1τ=A332ω02βcos3τ−εω0sin3τ+−−A3εω0+16c2ω02+8Aω1−6A3βτ16ω02sinτ+8Aω1+16c1ω02−5A3β+8Aεω0−2A3εω0τ16ω02cosτ+Fω02−Ω2cosΩω0τ,where A represents the amplitude of the oscillator.Secularity terms in the last expression are not allowed so that the coefficients ofcost and sint must be equal to zero which lead to(32)−5A3β+8Aω1+16c1ω02=0,−2AA2−4εω0=0,ω016c2ω0−A3ε=0,8Aω1−6A3β=0.Solving the last system gives(33)A=2,c1=−β2ω02,c2=ε2ω0,ω1=3β.Then, the following solution toOp2 and for p=1 is obtained:(34)xLP=2cosωt+14ω02βcos3ωt−εω0sin3ωt+Fω02−Ω2cosΩωω0t,where ω=3β+ω02 and ω02≠Ω2. This solution is obtained under the initial conditions (ICs)(35)x0=2+β4ω02+Fω02−Ω2,x˙0=−3εΘ4ω0.Solution (34) recovers the unforced case for F=0.


### 3.3. Improved He’s Frequency-Amplitude Formulation

To demonstrate the general idea of He’s frequency-amplitude formulation [31–35], let us consider the following oscillator:(36)x¨+fx=0,where fx indicates the nonlinear restoring force. The following conditions are hold: f0=0 and fx/x>0 .He considered Duffing oscillator(37)x¨+x+εx3=0,x0=A,x˙0=0,where A represents the amplitude of the oscillator. Based on He’s principle, we have(38)fx=x+εx3ψ˙2=ω2=fxx|x=3/2A=1+3A22,where ω denotes the frequency of oscillator.Now, by considering the following DVdP oscillator:(39)x¨−ε1−x2x˙+ω02x+βx3=0,in this case, the function fx reads(40)fx=ω02x+βx3,which leads to(41)ψ˙2=ω2=fxx|x=3/2A=ω02+34βA2.Since the amplitude now depends on time, we will reason heuristically to determine it(42)A=at,a˙=εκ2a−ε8a3,with(43)κ=ω022−β9β−4ω02/ω04+4+3β4β.Note that κ⟶1 as β⟶0.Then, the improved He’s solution becomes(44)xt=atcosωt,with(45)at=2κeεκt/24κ/c02+eεκt−1,ωt=c1+∫0tω02+34βa2τdτ=Wt−W0+c1,where(46)Wt=1εκ23βκ+ω02tanh−13βκeεκt/4κ/c02+eεκt−1+ω023βκ+ω02−2ω0coth−1ω03βκeεκt/4κ/c02+eεκt−1+ω02.The constantsc0 and c1 are determined from the ICs x0=x0 and x′0=x˙0. The amplitude for the limit cycle reads(47)rβ=limt⟶+∞at=9β2−4βω02+4ω04+3β+2ω02ω0⟶2asβ⟶0.As a numerical example, we can use the same model and data that were given in Ref. [36], which lead to the following unforced DVdP i.v.p. (27):(48)x¨−0.11−x2x˙+x+0.01x3=0,x0=2,x′0=0.Solution (10) and RK numerical approximation to i.v.p. (48) are graphically mapped as shown in Figure 1. Moreover, the approximation (34) using He’s HMP and the approximation (44) using the improved He’s FAF are compared with the obtained analytical approximation (10) and RK numerical approximation as illustrated in Figure 1. In addition, the maximum distance error in the whole time domain 0≤t≤50 with respect to RK numerical approximation is estimated(49)LR=max0≤t≤50RK−xApprox=0.0108535,LR=max0≤t≤50RK−xHe′sHPM=0.125709,LR=max0≤t≤50RK−xHe′sFAF=0.164318.Figure 1
The approximate solutions to the unforced i.v.p. (48) using the analytical approximation (10) and RK numerical approximation as well as He’s FAF, He’s HPM, and PLT approximations are plotted in the x,t− plane.
(a)(b)(c)It is clear that the analytical approximation (10) and RK numerical approximation are very compatible with each other. Also, they are more accurate than He’s FAF and He’s HPM approximations.


## 3.1. Our New Approach

By usingN=2 in solution (4), then the solution to the following unforced problem (3)(8)x¨+ω02x+p−ε1−x2x˙+βx3=0,can be introduced in following form:(9)xp=acosψ+∑k=12a2k+1pk∑j=1kr2j+1,kcos2j+1ψ+s2j+1,ksin2j+1ψ.Accordingly, we get(10)xp=acosψ+pa3A1t+p2a3A2t+a5A3t,with(11)A1t=r3,1,1cos3ψ+s3,1,1sin3ψ,A2t=r3,1,2cos3ψ+s3,1,2sin3ψ,A3t=r3,2,2cos3ψ+r5,2,2cos5ψ+s3,2,2sin3ψ+s5,2,2sin5ψ,(12)a˙=pB1,2at3+B1,1at+p2B2,3at5+B2,2at3,ψ˙=ω0+pC1,0at2+p2C2,2at4+C2,1at2+C2,0.On the other hand(13)Rp=D1p+D2p2+Op3,with(14)D1=S1asinψ+a3S2sinψ+S3cosψ+S4cos3ψ+S5sin3ψ,D2=a5S6sinψ+S7cosψ+S8sin3ψ+S9cos3ψ+S10sin5ψ+S11cos5ψ+a3S12sinψ+S13cosψ+S14sin3ψ+S15cos3ψ+aS16cosψ,where the values of coefficients Sii=1,2,…,16 are defined in Appendix (i).Equating the coefficients ofp and p2 to zero in(15)Rp=x¨+ω02x+p−ε1−x2x˙+βx3,also, equating to zero the coefficients of a2j+1, cos2i+1ψ, and sin2i+1ψ, where (i, j=1,2,3,…,N), i.e., Si=0, we get an algebraic system. The solution of this system yields(16)B1,1=ε2,B1,2=−ε8,B2,2=−3βε16ω02,B2,3=βε32ω02,C1,0=3β8ω0,C2,0=−ε28ω0,C2,1=ε28ω0,C2,2=−15β2+7ε2ω02256ω03r3,1,1=β32ω02,r3,1,2=−ε2128ω02,r3,2,2=−21β2+ε2ω021024ω04,r5,2,2=3β2−5ε2ω023072ω04,s3,1,1=−ε32ω0,s3,1,2=−3βε128ω03,s3,2,2=3βε256ω03,s5,2,2=−βε384ω03.From the above values in equation (11), we have(17)A1t=132ω02βcos3ψ−εω0sin3ψ,A2t=−ε128ω033βsin3ψ+εω0cos3ψ,A3t=−13072ω04321β2+ε2ω02cos3ψ+5ε2ω02−3β2cos5ψ+4βεω02sin5ψ−9sin3ψ.Accordingly, the odes for determining the functionsa,ψ read(18)a˙=p12εat−18εat3+p2−3βε16ω02at3+βε32ω02at5,ψ˙=ω0+p3β8ω0at2+p2−ε28ω0+ε28ω0at2−15β2+7ε2ω02256ω03at4.Forp=1, the value of a˙ given in equation (18) reduces to(19)a˙|p=1=12εa−18εa3+−3βε16ω02a3+βε32ω02a5.The amplitude for the limit cycle is obtained from the conditiona˙|p=1=0:(20)12εa−18εa3+−3βε16ω02a3+βε32ω02a5=0.Solving equation (20) gives(21)A=ω022−β9β−4ω02/ω04+4β+3.Observe that(22)limβ⟶0A=2,this is called the cycle amplitude for the VdP oscillator.We can use the following Chebyshev approximation in order to facilitate the solution to the ode system (18):(23)a˙|p=1=12εa−18εa3+−3βε16ω02a3+βε32ω02a5≈116ε8−βω02a−116ε2+βω02a3,=κa−μa3,with(24)κ=116ε8−βω02,μ=116ε2+βω02.By solving equation (23), we get(25)a=c0κeκtc02μe2κt−1+κ.Also, the expression for determiningψ can be obtained from the second equation in (18) for p=1 whose solution reads(26)ψ=∑i=14Wi+c1,where the values of Wii=1,2,3,4 are defined in Appendix (ii). The constants c0 and c1 are determined from the initial conditions (ICs) x0=x0 and x′0=x˙0.In all above expressions, forp=1, the approximation to the following i.v.p. can be obtained:(27)x¨−ε1−x2x˙+ω02x+βx3=0,x0=x0andx′0=x˙0,


## 3.2. He’s Homotopy Perturbation Method

Moreover, the approximate solution to the DVdP i.v.p. (1) using He’s HPM is obtained. Briefly, He’s HPM can be used for a series of nonlinear oscillators differential equations which many classical perturbation methods failed to solve them or to give some accurate solutions. This method suggests the solution in the following ansatz:(28)xLP=∑i=0∞pixiωt,with(29)ω2=ω02+∑i=1∞piωi.Substitute equations (28) and (29) into i.v.p. (1), and by collecting the coefficients of same powers of p, we finally obtain some reduced equations. We have(30)ℝF=x0τ+x0″τω0+βx0τ3+ω02x1τ+12−2εω0+τω1x0′τ+εω0x0τ2x0′τ+ω1x0″τ+ω02x1″τ+12τω1x03τ−FcosτΩω0p+Op2,where τ=tω.Equating to zero the coefficients ofpj and solving the resulting odes gives(31)x0τ=Acosτ,x1τ=A332ω02βcos3τ−εω0sin3τ+−−A3εω0+16c2ω02+8Aω1−6A3βτ16ω02sinτ+8Aω1+16c1ω02−5A3β+8Aεω0−2A3εω0τ16ω02cosτ+Fω02−Ω2cosΩω0τ,where A represents the amplitude of the oscillator.Secularity terms in the last expression are not allowed so that the coefficients ofcost and sint must be equal to zero which lead to(32)−5A3β+8Aω1+16c1ω02=0,−2AA2−4εω0=0,ω016c2ω0−A3ε=0,8Aω1−6A3β=0.Solving the last system gives(33)A=2,c1=−β2ω02,c2=ε2ω0,ω1=3β.Then, the following solution toOp2 and for p=1 is obtained:(34)xLP=2cosωt+14ω02βcos3ωt−εω0sin3ωt+Fω02−Ω2cosΩωω0t,where ω=3β+ω02 and ω02≠Ω2. This solution is obtained under the initial conditions (ICs)(35)x0=2+β4ω02+Fω02−Ω2,x˙0=−3εΘ4ω0.Solution (34) recovers the unforced case for F=0.


## 3.3. Improved He’s Frequency-Amplitude Formulation

To demonstrate the general idea of He’s frequency-amplitude formulation [31–35], let us consider the following oscillator:(36)x¨+fx=0,where fx indicates the nonlinear restoring force. The following conditions are hold: f0=0 and fx/x>0 .He considered Duffing oscillator(37)x¨+x+εx3=0,x0=A,x˙0=0,where A represents the amplitude of the oscillator. Based on He’s principle, we have(38)fx=x+εx3ψ˙2=ω2=fxx|x=3/2A=1+3A22,where ω denotes the frequency of oscillator.Now, by considering the following DVdP oscillator:(39)x¨−ε1−x2x˙+ω02x+βx3=0,in this case, the function fx reads(40)fx=ω02x+βx3,which leads to(41)ψ˙2=ω2=fxx|x=3/2A=ω02+34βA2.Since the amplitude now depends on time, we will reason heuristically to determine it(42)A=at,a˙=εκ2a−ε8a3,with(43)κ=ω022−β9β−4ω02/ω04+4+3β4β.Note that κ⟶1 as β⟶0.Then, the improved He’s solution becomes(44)xt=atcosωt,with(45)at=2κeεκt/24κ/c02+eεκt−1,ωt=c1+∫0tω02+34βa2τdτ=Wt−W0+c1,where(46)Wt=1εκ23βκ+ω02tanh−13βκeεκt/4κ/c02+eεκt−1+ω023βκ+ω02−2ω0coth−1ω03βκeεκt/4κ/c02+eεκt−1+ω02.The constantsc0 and c1 are determined from the ICs x0=x0 and x′0=x˙0. The amplitude for the limit cycle reads(47)rβ=limt⟶+∞at=9β2−4βω02+4ω04+3β+2ω02ω0⟶2asβ⟶0.As a numerical example, we can use the same model and data that were given in Ref. [36], which lead to the following unforced DVdP i.v.p. (27):(48)x¨−0.11−x2x˙+x+0.01x3=0,x0=2,x′0=0.Solution (10) and RK numerical approximation to i.v.p. (48) are graphically mapped as shown in Figure 1. Moreover, the approximation (34) using He’s HMP and the approximation (44) using the improved He’s FAF are compared with the obtained analytical approximation (10) and RK numerical approximation as illustrated in Figure 1. In addition, the maximum distance error in the whole time domain 0≤t≤50 with respect to RK numerical approximation is estimated(49)LR=max0≤t≤50RK−xApprox=0.0108535,LR=max0≤t≤50RK−xHe′sHPM=0.125709,LR=max0≤t≤50RK−xHe′sFAF=0.164318.Figure 1
The approximate solutions to the unforced i.v.p. (48) using the analytical approximation (10) and RK numerical approximation as well as He’s FAF, He’s HPM, and PLT approximations are plotted in the x,t− plane.
(a)(b)(c)It is clear that the analytical approximation (10) and RK numerical approximation are very compatible with each other. Also, they are more accurate than He’s FAF and He’s HPM approximations.


## 4. Analytical Approximation to the Forced Duffing–Van Der Pol Oscillator

Let us consider the following forced DVdP i.v.p.:(50)ℝF=x¨−ε1−x2x˙+ω02x+βx3−FcosΩt=0,x0=x0,x′0=x˙0.Assume that the solution to i.v.p. (50) is given by the following ansatz:(51)xt=yt+d0cosΩt+d1sinΩt,where y≡yt is a solution to the unforced DVdP oscillator(52)y¨−ε1−x2y˙+ω02y+βy3=0,y0=x0−d0,y′0=x˙0−d1Ω.Putting solution (51) into (50), we have(53)ℝF=RcosΩt+SsinΩt+h.o.t,with(54)R=3βd034+34βd12d0+14d1d02εΩ+14d13εΩ−d1εΩ+d0ω02−d0Ω2−F+yt23βd0+d1εΩ+2d0εy′tyt,S=34βd1d02+3βd134−14d03εΩ−14d12d0εΩ+d0εΩ+d1ω02−d1Ω2−yt23βd1+d0εΩ+2d1εy′tyt,where h.o.t. represents higher-order terms. By neglecting yt and y′t from system (54) at R,S=0,0, then the constants d0 and d1 can be determined from the system(55)3βd034+34βd12d0+14d1d02εΩ+14d13εΩ−d1εΩ+d0ω02−d0Ω2−F=0,34βd1d02+3βd134−14d03εΩ−14d12d0εΩ+d0εΩ+d1ω02−d1Ω2=0.From this system, we get(56)Y0+Y1d0+Y2d02+Y3d03=0,(57)Z0+Z1d1+Z2d12+Z3d13=0,where the values of Yii=0,1,2,3 and Zii=0,1,2,3 are defined in Appendix (iii). We choose the least in magnitude real roots to cubics (56) and (57).As a numerical example, the two approximations (51) and (34) according to ICs (35) for i.v.p. (50) are displayed in Figure 2 for ε,ω0,β,F,Ω=0.1,1,0.01,1,5.2. Also, the maximum distance error for the two approximations is estimated as follows:(58)LR=max0≤t≤50RK−xApprox=0.0227101,LR=max0≤t≤50RK−xHe′sHMP=0.106029.Figure 2
The two approximations (51) and (34) according to ICs (35) for the forced i.v.p. (50) are plotted in the x,t− plane.
(a)(b)On the other side, for arbitrary ICs, the analytical approximation (51) versus the RK numerical approximation is presented in Figure 3 for ω0,β,F,Ω,x0,x˙0=1,0.01,1,5.2,0,0.183 and different values to ε. Also, the maximum distance error at the same values of the physical parameters and ICs mentioned in Figure 3 is estimated as follows:(59)LR|ε=0.1=max0≤t≤40RK−xApprox=0.00178818,LR|ε,=−0.1=max0≤t≤40RK−xApprox=0.00005143.Figure 3
Both analytical approximation (51) and RK numerical approximation to the forced i.v.p. (50) for different values of ε are plotted in the x,t− plane.
(a)(b)We can conclude that in all cases, both two analytical approximations (10) and (51) for the unforced and forced DVdP oscillators are more accurate and convergence as compared to He’s FAF and He’s HPM.


## 5. Conclusion

Both higher-order nonlinearity unforced Duffing–Van der Pol (DVdP) oscillator and forced DVdP oscillator having linear and cubic nonlinear terms have been analyzed using some effectiveness and more accurate approaches. The new approach was constructed based on the Krylov–Bogoliubov–Metroolsky method (KBMM). The new approach was discussed in detail for the two issues. In our analysis, we only stopped at the first approximation because it is sufficient in all cases. Also, this new approach can be used for analyzing many strongly nonlinear oscillators. Moreover, the new approach does not demand to solve any ordinary differential equations (odes) because we only simply equate the coefficients of the trigonometric functionsat and p to zero in order to get a simple system of algebraic equations. Then, this system becomes very easy to solve it to determine the undetermined coefficients. Also, this new method is also characterized by being direct and fast. Moreover, it is characterized by high-accuracy if it is compared with other methods in the literature. Also, we improved He’s frequency-amplitude formulation technique in order to solve unforced DVdP oscillator to obtain high-accurate results. Furthermore, the unforced DVdP oscillator was analyzed via He’s homotopy perturbation method. The maximum distance error in the whole time domain with respect to Runge–Kutta numerical approach has been estimated. It was found that our new approach is better than all method approaches. Moreover, the new approach can be devoted for analyzing many strong nonlinearity oscillators with any odd power. Also, the new approach can be applied for arbitrary initial conditions.Future work: the ansatz that has been used in this paper is called the KBM first-Variant. This approach cannot recover He’s amplitude formula. On the other side, we will use another new ansatz which maybe called the KBM second-Variant, in this case, He’s amplitude formula can be recovered.



---

*Source: 2174192-2022-07-31.xml*