# Solution to a Damped Duffing Equation Using He’s Frequency Approach

Authors: Alvaro H. S. Salas; Gilder-Cieza Altamirano; Manuel Sánchez-Chero
Journal: The Scientific World Journal
 (2022)

---

## Abstract

In this paper, we generalize He’s frequency approach for solving the damped Duffing equation by introducing a time varying amplitude. We also solve this equation by means of the homotopy method and the Lindstedt–Poincaré method. High accurate formulas for approximating the Jacobi elliptic function cn are formally derived using Chebyshev and Pade approximation techniques.


---

## Body

## 1. Introduction

The deep understanding of the mechanism of nonlinear oscillations has an effective role in interpreting the ambiguities of many natural, physical, and engineering phenomena in various fields of science. Accordingly, many researchers have been able to give correct scientific explanations about their scientific experiences based on a deep understanding of the characteristics of these phenomena after the clarity of the ambiguity about the phenomenon under study. In the framework of nonlinear dynamics, there is no doubt that the scenario of dynamic mechanism of the pendulum motion is one of the objects that have deserved more attention in modeling all kind of phenomena related to oscillations, bifurcations, and chaos. The simple pendulum has been used as a physical model to solve problems such as nonlinear plasma oscillations, Duffing oscillators, Helmholtz oscillations, rigid plates that satisfy the Johanessen performance criteria, transverse vibrations nonlinear of a plate carrying a concentrated mass, a beam supported by a double periodic axial oscillating mount, cracks subjected to concentrated forces, surface waves in a plasma column, coupled modes of nonlinear bending vibrations of a circular ring, double spin spacecraft, motion of spacecraft over slowly rotating asteroids, nonlinear vibration of clasped beams, the nonlinear equation of wave, and nonlinear mathematical models of DNA.It is known that the main objective of the numerical approaches is to find some numerical solutions to various realistic physical, engineering, and natural problems, especially when exact solutions are unavailable or extremely difficult to determine. There are many numerical approaches that were used for analyzing the family of the Duffing oscillator and Duffing–Helmholtz oscillator with constant coefficients. It is known that this family is integrable, i.e., its exact solution is available in the absence of the damping effect. On the other hand, if the damping effect and some other friction forces are taken into account, we get a nonintegral differential equation, i.e., its exact solution is not available.The nonlinear oscillators have many applications in science and engineering. One of such oscillators is the Duffing equation. George Duffing, a German engineer, wrote a comprehensive book about this in 1918. Since then, there has been a tremendous amount of work done on this equation, including the development of solution methods (both analytical and numerical), and the use of these methods to investigate the dynamic behavior of physical systems that are described by the various forms of the Duffing equation.Solution methods for nonlinear oscillators include the G/G method, modified mapping method and the extended mapping method, elliptic expansion method, modified (G/G)-expansion method, dynamical systems approach, the modified trigonometric function series method, generalized (G/G)-expansion method, tanh method, and the sn-ns method, among others [1–10].In this paper, we consider the damped and unforced Duffing equation. We solve it using an extended version of He’s frequency approach for the damped case. This oscillator was solved in [11] using a generalized elliptic functions. A simplification of the solution in [11] may be obtained by approximating the elliptic functions by means of trigonometric functions [12].


## 2. The Damped Duffing Oscillator

Let us consider the i.v.p.(1)x¨+2εx˙+αx+βx3=0,x0=x0andx′0=x˙0.The damped oscillator (1) is integrable only when α=8/9ε2. Equation (1) represents a damped Duffing oscillator. In the case when ε=0, we have an undamped Duffing oscillator which has an exact solution for any given arbitrary initial conditions. More precisely, the exact solution to the i.v.p.(2)x¨+αx+βx3=0,x0=x0andx′0=x˙0,is given by(3)xt=x0cnωt|m+b1snωt|mdnωt|m1+b2sn2ωt|m,where(4)b1=x˙0α+βx022+2βx˙02,b2=α+βx022α+βx022+2βx˙02−12,ω=α+βx022+2βx˙02,andm=12−αα+βx022+2βx˙02,α+βx022+2βx˙02≠0.The solution is periodic with period(5)T=4Kmω,where(6)Km=∫0π/2dθ1−msin2θ.We may use the following approximation formulas for evaluating the elliptic integralKm (see Table 1).Table 1
Approximations for the elliptic integralKm=∫0π/2dθ/1−msin2θ.
KmError on−1≤m≤0.2530m/59+19/110.16935711m2/26+23m/36+11/70.07184277m3/15+11m2/19+15m/43+35/230.050337319m4/32+19m3/28+2m2/29+5m/22+11/70.040271223m5/28+9m4/10−7m3/31−5m2/27+17m/41+51/320.0232432π/2−5πm/32/1−9m/160.000714295−19πm2/1152−2πm/9+π/2/1−25m/360.0000730467409πm2/9728−249πm/608+π/2/1025m2/4864−325m/304+13.85×10−6See Table2 for the approximation of 1/Km.Table 2
Approximate trigonometric solution to the i.v.p. (1) may be obtained in many ways.
1/KmError on−1≤m≤0.2520/33−5m/280.0322711−m2/13−8m/41+7/110.0104322−2m3/31−2m2/21−2m/13+9/140.0063857−3m4/44−3m3/35−m2/29−10m/71+7/110.00418393−m5/12−2m4/21+m3/124−m2/95−5m/31+19/300.003314872/π−9m/8π/1−5m/160.00025146419m2/160π−8m/5π+2/π/1−11m/200.00001933831025m2/2432π−325m/152π+2/π/409m2/4864−249m/304+11.36×10−6For example,(7)xt=c0cosωt+cos−1x0c0,where(8)c0=±x02+x˙02ω2,ω=2πT,T=4Kmα+βx02andm=βc022α+βc02.Let us consider the i.v.p.(9)x¨+αx+βx3=0,x0=Aandx′0=0.A very good approximate analytical solution is given by(10)xt=A1+λcosωt1+λcos2ωt,where(11)ω=πα+βA22Km,λ=π24Km2−1,m=A2β2α+A2β.Now, observe that the functionxt=cnt,m obeys the Duffing equation(12)x¨+1−2mx+2mx3=0,x0=1andx′0=0.Letx=xt be a continuous on 0,4Km function. Define(13)Emx=max0≤t≤4Kmcnt,m−xt.We have the following approximations for0≤m≤0.5(14)cnt,m≈x1t:=cosπ2Kmt,Emx1≈13m149+6m265≤0.0674for0≤m≤0.5,cnt,m≈x2t:=cos1−5m21−3m223t,Emx2≈14m2153+17m186≤0.06834for0≤m≤0.5,cnt,m≈x3t:=1−4m/7cos−3m2/23−5m/21+1t1−4/7mcos2−3m2/23−5m/21+1t,Emx3≈m100−m2210≤0.0089for0≤m≤0.5,cnt,m≈x4t:=1−2143m/3564cos−10m2/77−11m/46+1t1−53/81mcos2−10m2/77−11m/46+1t+7/132mcos4−10m2/77−11m/46+1t,Emx4≈2m91−4m289≤0.0045for0≤m≤0.5,cnt,m≈x5t:=1−m10cos1−5m21−3m223t+m10cos31−5m21−3m223t,Ex5≈5m69−11m2109≤0.0153for0≤m≤0.5,cnt,m≈x6t:=1−131m1500cos−3m223−5m21+1t+m12cos3−3m223−5m21+1t+m250cos5−3m223−5m21+1t,Emx6≈m23−5m262≤0.00781for0≤m≤0.5.The approximations above may be used to give trigonometric solution to the i.v.p. (9) whose exact solution reads(15)xt=Acnα+βA2t,βA22α+βA2.For example,(16)xt=Acos1−1033βA22α+βA2α+βA2t.The frequency-amplitude formulation for this solution is given by(17)ω2=33α+28A2β21089α+A2β.We obtain several frequency-amplitude formulations using formulas. The He’s frequency-amplitude formulation for the Duffing equationx¨+αx+βx3=0 establishes that(18)ω02=w′t2=α+34βA2,wherewt=α+34βA2t.


## 3. He’s Approach

Using He’s frequency approach, an approximate analytical solution in the absence of dampingε=0 may be obtained using(19)xapproxt=Acosωt+cos−1x0A,A=±−4α+3βx02±4α+3βx022+48βx˙026βandω=α+3β4A2.Let(20)wt=ωt.Observe thatw0=0. He’s idea for the undamped case is based on the following fact:(21)w′t2=ω02=fxx=α+3β4A2forx=32Aandfx=αx+βx3.Following this idea, for the damped case, we will replaceA with Aexp−ρt so that(22)w′t2=α+3β4A2exp−2ρtandw0=0.The approximate analytical solution for the damped Duffing equation (1) will then be(23)xapproxt=Aexp−ρtcoswt+cos−1x0A.From (22) it follows that(24)wt=∫0tα+3β4A2exp−2ρτdτ=Wt−W0,where(25)Wt=4α+3A2βe−2ρt2Aρ2αAeρt3A2β+4αe2ρtsinh−12αeρt3Aβ−A.The numberρ is a free parameter that is chosen in order to get as small residual error as possible. The default value is ρ=ε.


## 4. Homotopy Perturbation Method

We seek a solution in the ansatz form:(26)xt=v0t+pv1t+p2v2t+⋯.The homotopy is defined as follows:(27)Hx,p=x¨+2εx˙+αx+pβx3.We have(28)Hx,p=αv0t+2εv0′t+v0″t+αv1t+βv0t3+2εv1′t+v1″tp+αv2t+3βv1tv0t2+2εv2′t+v2″tp2+⋯.We now equate to zero the coefficients ofpjj=0,1,2,3,… and then we obtain an ode system. We solve these odes so that the functions vjtj=1,2,3,… do not contain secularity terms. The solutions are(29)v0t=c0e−εc1+tcosα−ε2c1+t,v1t=βc03e−3εc1+t16αα−ε23ε2−4α3εε2−αsin3α−ε2c1+t+α−ε23ε2−2αcos3α−ε2c1+t,v2t=3β2c05e−5εc1+t512α23ε2−4αα+3ε25ε2−9α−45ε2−9αα2−9αε2+9ε4cos3α−ε2c1+t+2α+3ε23α2−17αε2+15ε4cos5α−ε2c1+t+εα−ε2182ε2−α5ε2−9αsin3α−ε2c1+t−α+3ε230ε2−19αsin5α−ε2c1+tThe numbersc0 and c1 are obtained from the initial conditions x0=x0 and x′0=x˙0.


## 5. Improved Analytical Solution

### 5.1. First Approach

Assume the ansatz(30)xt=Aexp−ρtcosωt+cos−1x0A,f0=0.Then,(31)x¨+2εx˙+αx+βx3=14A3βcos3θe−3ρt+14Acosθe−3ρt3A2β−4e2ρtω′t2+4e2ρtα−2ερ+ρ2−Asinθeρ−tω″t+2ε−ρω′t.We will choose the functionω=ωt so that(32)3A2β−4e2ρtω′t2+4e2ρtα−2ερ+ρ2=0andω0=0,then,(33)ωt=∫0tα−2ερ+ρ2+34A2βe−2ρτdτ=Ωt−Ω0,where(34)Ωt=1ρα−2ερ+ρ2tanh−1α−2ερ+ρ2+3/4A2e−2tρβα−2ερ+ρ2−α−2ερ+ρ2+34A2e−2tρβ.Observe that(35)ω′t2=α−2ερ+ρ2+34A2βe−2ρt=fxxforx=32Aexp−ρt,fx=α−2ερ+ρ2x+βx3.In other words, whenε=ρ=0, this corresponds to the He’s frequency formulation for the Duffing equation:(36)x¨+α−2ερ+ρ2x+βx3=0.The numberρ is chosen in order to minimize the residual error. A default value for ρ is obtained by eliminating A from the system:(37)ω″0+2ε−ρω′0=0andx′0=x˙0.This last condition gives the sextic:(38)12βε2x˙02−2ε4α2+3αβx02−12βεx0x˙0+18βx˙02ρ+8α2+32αε2+9αβx02+24βε2x02−72βεx0x˙0+27βx˙02ρ2−224αε+16ε3+30βεx02−27βx0x˙0ρ3+44α+16ε2+9βx02ρ4−40ερ5+8ρ6=0.A root to this sextic nearρ=ε may be evaluated using the following approximation:(39)ρ≈ε−6βεαx02+2εx˙0x0+x˙024α−ε22+6βx02α+ε2+21βεx0x˙0+9βx˙0264α−ε24+3β3βx0413α2+12αε2+16ε4+4εx˙0x052α−ε22+x0257αβ+54βε2+4x˙0222α−ε22+3βx029α+41ε2+8x02α−ε227α+8ε2+360βεx˙03x0+81βx˙04.The value ofA is given by(40)A=±3βx02+8ερ−4α−4ρ2±4α+ρρ−2ε−3βx022+48βx02α+2ρρ−ε+2ρx˙0x0+x˙026β.


### 5.2. Second Approach

The exact solution to the undamped Duffing equationx¨+αx+βx3=0 is given by(41)xt=c0cnwt+c1,m,where(42)wt=α+βc02tandm=βc022α+βc02.Then,(43)w′t2=α+βc02andm=βc022w′t2.Let us replace α with α−2ρε+ρ2 and c0 with c0exp−ρt so that(44)w′t2=α−2ρε+ρ2+βc02exp−2ρtandm=mt=βc02exp−2ρt2w′t2=1/21+α−2ρε+ρ2/βc02exp2ρt.For the damped Duffing equation we will define the solution in the form:(45)xt=c0exp−ρtcnλ∫0tα−2ρε+ρ2+βc02exp−2ρtdτ+c1,1/21+α−2ρε+ρ2/βc02exp2ρt.The constantsc0 and c1 are determined from the initial conditions. On the other hand,(46)∫0tα−2ρε+ρ2+βc02exp−2ρtdτ=Wt−W0,where(47)Wt=1ρα−2ερ+ρ2tanh−1βc02e−2ρtα−2ερ+ρ2+1−α+βc02e−2ρt−2ερ+ρ2.The numbersρ and λ are free parameters that we choose in order to get as small residual error as possible. The default ρ value is ρ=2/3ε and the default value for λ is λ=1. Observe that this ansatz will give the exact solution for the integrable case, i.e., when ρ=2/3ε and α=8/9ε2.


## 5.1. First Approach

Assume the ansatz(30)xt=Aexp−ρtcosωt+cos−1x0A,f0=0.Then,(31)x¨+2εx˙+αx+βx3=14A3βcos3θe−3ρt+14Acosθe−3ρt3A2β−4e2ρtω′t2+4e2ρtα−2ερ+ρ2−Asinθeρ−tω″t+2ε−ρω′t.We will choose the functionω=ωt so that(32)3A2β−4e2ρtω′t2+4e2ρtα−2ερ+ρ2=0andω0=0,then,(33)ωt=∫0tα−2ερ+ρ2+34A2βe−2ρτdτ=Ωt−Ω0,where(34)Ωt=1ρα−2ερ+ρ2tanh−1α−2ερ+ρ2+3/4A2e−2tρβα−2ερ+ρ2−α−2ερ+ρ2+34A2e−2tρβ.Observe that(35)ω′t2=α−2ερ+ρ2+34A2βe−2ρt=fxxforx=32Aexp−ρt,fx=α−2ερ+ρ2x+βx3.In other words, whenε=ρ=0, this corresponds to the He’s frequency formulation for the Duffing equation:(36)x¨+α−2ερ+ρ2x+βx3=0.The numberρ is chosen in order to minimize the residual error. A default value for ρ is obtained by eliminating A from the system:(37)ω″0+2ε−ρω′0=0andx′0=x˙0.This last condition gives the sextic:(38)12βε2x˙02−2ε4α2+3αβx02−12βεx0x˙0+18βx˙02ρ+8α2+32αε2+9αβx02+24βε2x02−72βεx0x˙0+27βx˙02ρ2−224αε+16ε3+30βεx02−27βx0x˙0ρ3+44α+16ε2+9βx02ρ4−40ερ5+8ρ6=0.A root to this sextic nearρ=ε may be evaluated using the following approximation:(39)ρ≈ε−6βεαx02+2εx˙0x0+x˙024α−ε22+6βx02α+ε2+21βεx0x˙0+9βx˙0264α−ε24+3β3βx0413α2+12αε2+16ε4+4εx˙0x052α−ε22+x0257αβ+54βε2+4x˙0222α−ε22+3βx029α+41ε2+8x02α−ε227α+8ε2+360βεx˙03x0+81βx˙04.The value ofA is given by(40)A=±3βx02+8ερ−4α−4ρ2±4α+ρρ−2ε−3βx022+48βx02α+2ρρ−ε+2ρx˙0x0+x˙026β.


## 5.2. Second Approach

The exact solution to the undamped Duffing equationx¨+αx+βx3=0 is given by(41)xt=c0cnwt+c1,m,where(42)wt=α+βc02tandm=βc022α+βc02.Then,(43)w′t2=α+βc02andm=βc022w′t2.Let us replace α with α−2ρε+ρ2 and c0 with c0exp−ρt so that(44)w′t2=α−2ρε+ρ2+βc02exp−2ρtandm=mt=βc02exp−2ρt2w′t2=1/21+α−2ρε+ρ2/βc02exp2ρt.For the damped Duffing equation we will define the solution in the form:(45)xt=c0exp−ρtcnλ∫0tα−2ρε+ρ2+βc02exp−2ρtdτ+c1,1/21+α−2ρε+ρ2/βc02exp2ρt.The constantsc0 and c1 are determined from the initial conditions. On the other hand,(46)∫0tα−2ρε+ρ2+βc02exp−2ρtdτ=Wt−W0,where(47)Wt=1ρα−2ερ+ρ2tanh−1βc02e−2ρtα−2ερ+ρ2+1−α+βc02e−2ρt−2ερ+ρ2.The numbersρ and λ are free parameters that we choose in order to get as small residual error as possible. The default ρ value is ρ=2/3ε and the default value for λ is λ=1. Observe that this ansatz will give the exact solution for the integrable case, i.e., when ρ=2/3ε and α=8/9ε2.


## 6. Lindstedt–Poincaré Method

We seek a solution in the ansatz form:(48)xt=u0t+βu1t+β2u2t+⋯,then.(49)x¨++2εx˙+αx+βx3=αu0t+2εu0′t+u0″t+αu1t+2εu1′t+u1″t+u0t3β+αu2t+2εu2′t+u2″t+3u1tu0t2β2+⋯.We now equate to zero the coefficients ofβjj=0,1,2,3,… and then we obtain an ode system. We solve these odes so that the functions ujtj=1,2,3,… do not contain secularity terms. The solutions are(50)u0t=c0e−εc1+tcosα−ε2c1+t,u1t=c03e−3εc1+t16α3ε2−4α3ε2−2αcos3α−ε2c1+t−3εα−ε2sin3α−ε2c1+t,u2t=3c05e−5εc1+t−45ε2−9αα2−9αε2+9ε4cos3α−ε2c1+t+2α+3ε23α2−17αε2+15ε4cos5α−ε2c1+t+εα−ε2182ε2−α5ε2−9αsin3α−ε2c1+t−α+3ε230ε2−19αsin5α−ε2c1+t512α23ε2−4αα+3ε25ε2−9α.The constantsc0 and c1 are obtained from the initial conditions.


## 7. Numerical Solution

We make use of the following backward finite differences formulas for the first and second derivatives:(51)x˙ti≈3xi−4−16xi−3+36xi−2−48xi−1+25xi12Δtfori≥4,x¨ti≈11xi−4−56xi−3+114xi−2−104xi−1+35xi12Δt2fori≥4.The discretized ode reads(52)11xi−4−56xi−3+114xi−2−104xi−1+35xi12Δt2+2ε3xi−4−16xi−3+36xi−2−48xi−1+25xi12Δt+αxi+βxi3=0i=4,5,6,…,n,n=⌊TΔt⌋,ti=i⋅Δt.The valuesx1=xt1, x2=xt2, and x3=xt3 are obtained from some reasonable analytical or numerical solution. The value x4 is found from the cubic:(53)βx43+α+510εΔt+712Δt2x4+x06εΔt+11−8x14εΔt+7+6x212εΔt+19−8x312εΔt+1312Δt2=0.We choose the closest tox3 real root to the cubic in equation (53). Suppose we already found the values xk for k=4,5,…,i−1. Then, the value of xi is obtained by choosing the closest to xi−1 real root to the cubic in equation (52). Thus, we solve the ode recursively using Tartaglia’s formula for the cubic. For small Δt, we may use the following practical formula:(54)xi=xi−1−p+3xi−12xi−13+pxi−1+qp2−3qxi−13pxi−12+6xi−14,where(55)p=6εΔt+11xi−4−84εΔt+7xi−3+612εΔt+19xi−2−812εΔt+13xi−112βΔt2,q=12αΔt2+50εΔt+3512βΔt2.Remark 1.
We also may use the following formulas based on the Runge–Kutta method.(56)xi=xi−1+Δt26fΔt2+xi−1,12Δtxi−1+xi−112Δtfxi−1,xi−1,xi−1+xi−1+fΔt2+xi−1,14ΔtΔtfxi−1,xi−1,xi−1+2xi−1+xi−1,12ΔtfΔt2+xi−1,12Δtxi−1+xi−1,12Δtfxi−1,xi−1,xi−1+xi−1+xi−1+fxi−1,xi−1,xi−1+Δtxi−1,(57)xi=xi−1+Δt62fΔt2+xi−1,12Δtxi−1+xi−1,12Δtfxi−1,xi−1,xi−1+xi−1+fΔt2+xi−1,14ΔtΔtfxi−1,xi−1,xi−1+xi−1+2xi−1+xi−1,12ΔtfΔt2+xi−1,12Δtxi−1+xi−1,12Δtfxi−1,xi−1,xi−1+xi−1+xi−1+fΔt+xi−1,Δt12ΔtfΔt2+xi−1,12Δtfxi−1,xi−1,xi−1+xi−1+xi−1+xi−1,ΔtfΔt2+xi−1,12Δtxi−1+xi−1,12Δtfxi−1,xi−1,xi−1+xi−1+xi−1+xi−1+fxi−1,xi−1,xi−1.
The above formulas are used to solve any i.v.p. having the form:(58)x¨=ft,x,x˙,x0=x0andx′0=x˙0.
The initial values for solving the recurrences (56)–(57) are x0 and x˙0.


## 8. Analysis and Discussion

Let us examine the accuracy of the obtained solutions. Letε=0.1, α=β=1, x0=1, x˙0=0, and 0≤t≤35. The problem to be solved reads(59)x¨+0.2ε+x+x3=0,x0=1andx′0=0.The trigonometric solution is(60)xtrigot=1.00251exp−0.0938627tcos5.340453.96015+3.015082.71828−0.187725t+0.995006log0.5025093.96015+3.015082.71828−0.187725t−1.−0.995006log0.5025093.96015+3.015082.71828−0.187725t+1.−0.68891+0.070781.The Jacobi solution is(61)xJacobit=1.00227exp−0.0932826tcnωt|12+1.97113exp0.186565,where(62)ωt=−10.50791.00455e−0.186565t+0.990045+10.4554tanh−10.99501e0.0932826t0.990045e0.186565t+1.00455+5.61143.See Figure1 for comparison with numerical solution.Figure 1
Approximate analytical solutions. Comparison between the numerical solution and the analytical approximate solutions-elliptic and trigonometric.



---

*Source: 5009722-2022-07-11.xml*