# Analytical Approximant to a Quadratically Damped Forced Cubic-Quintic Duffing Oscillator

Authors: Alvaro H. S. Salas
Journal: The Scientific World Journal
 (2022)

---

## Abstract

The cubic-quintic Duffing oscillator of a system with strong quadratic damping and forcing is considered. We give elementary approximate analytical solution to this oscillator in terms of exponential and trigonometric functions. We compare the analytical approximant with the Runge–Kutta numerical solution. The approximant allows us to estimate the points at which the solution crosses the horizontal axis.


---

## Body

## 1. Introduction

In this paper, some novel analytical and numerical techniques are introduced for analyzing and solving nonlinear ordinary differential equations (NODEs) that are associated to some strongly nonlinear oscillators such as a quadratically damped cubic-quintic Duffing equation. There are many numerical and analytical approaches that were applied for solving the second-order nonlinear oscillator equations. For instance, both the homotopy perturbation method (HBM) and MTS technique were applied for analyzing a forced Van der Pol (VdP) generalized oscillator to obtain the amplitudes of the forced harmonic and super and subharmonic oscillatory states[1]. Also, Melnikov’s method was employed for analyzing a mVdPD equation and deriving analytical criteria for the appearance of horseshoe chaos in chemical oscillations [2]. He et al. [3] used the Poincare–Lindstedt technique (PLT) for solving and analyzing the hybrid Rayleigh–Van der Pol–Duffing equation. Moreover, the homotopy analysis method (HAM) was employed for analyzing the DVdP oscillator [4]. Both methods of differentiable dynamics and Lie symmetry reduction method were devoted for analyzing the DVdP-type oscillator [5]. The principal feature associated with quadratic damping is a discontinuous jump of the damping force in the equation of motion whenever the velocity vanishes such that the frictional force always opposes the motion.In this paper, we will consider the following quadratically damped and forced cubic-quintic Duffing oscillator:(1)x¨+εx˙x˙+αx+βx3+γx5=Ft,x0=x0andx′0=x˙0.The quadratically damped oscillator (1) is never critically damped or overdamped. In the absence of damping and forcing, we obtain the cubic-quintic Duffing equation(2)x¨+αx+βx3+γx5=0,x0=x0andx′0=x˙0.Equation (2) admits exact analytical solution that is expressed in terms of the Jacobian elliptic functions. Other solution methods may be found in [1–8].


## 2. Solution Procedure

### 2.1. First Case: Undamped and Unforced Cubic-Quintic Duffing Equation

Let us consider the i.v.p.(3)x¨+αx+βx3+γx5=0,x0=x0andx′0=x˙0.Assume the ansatz(4)xt=vt1+λv2t+μv4t,where the function v=vt is the solution to some Duffing equation(5)v¨+pv+qv3=0,v0=v0andv′0=v˙0.The numbersv0 and v˙0 are determined from the initial conditions. Observe that(6)v˙0=−x˙0λv02+μv04+13/2μv04−1.From (5), it follows that(7)v˙2=v˙02+pv02+qv042−pv2−qv42.Observe that(8)12x˙2+12αx2+14βx4+16γx6−αx022−βx044−16γx06−x˙022=0.We have(9)0=12x˙2+12αx2+14βx4+16γx6−αx022−βx044−16γx06−x˙022=1121+λv2+μv4∑j=06Sjv2j,where(10)S0=6pv02+3qv04+6v˙02−6αx02−3βx04−2γx06−6x˙02,S1=−3−2α+2p+6αλx02+3βλx04+2γλx06+6λx˙02,S2=−3−4αλ−β+4μpv02+2μqv04+q+4μv˙02+6αλ2x02+6αμx02+3βλ2x04+3βμx04+2γλ2x06+2γμx06+6λ2x˙02+6μx˙02,S3=6αλ2+12αμ+3βλ+2γ+12μp−6αλ3x02−36αλμx02−3βλ3x04−18βλμx04−2γλ3x06−12γλμx06−6λ3x˙02−36λμx˙02,S4=3μ4αλ+β+2μpv02+μqv04+2q+2μv˙02−6αλ2x02−6αμx02−3βλ2x04−3βμx04−2γλ2x06−2γμx06−6λ2x˙02−6μx˙02,S5=−3μ2−2α+2p+6αλx02+3βλx04+2γλx06+6λx˙02,S6=−μ23q+6αμx02+3βμx04+2γμx06+6μx˙02.Eliminatingp, q, μ, v0, and v˙0 from the system S0=S1=S2=S3=S4=S5=S6=0 gives(11)6αx02+3βx04+2γx06+6x˙02λ3−6αλ2−3βλ−2γ46αx02+3βx04+2γx06+6x˙022λ3−24α6αx02+3βx04+2γx06+6x˙02λ2+312α2−6αβx02−3β2x04−2βγx06−6βx˙02λ+9αβ+6αγx02+3βγx04+2γ2x06+6γx˙02=0.The values forp, q, and μ are(12)p=122α−6αλx02−3βλx04−2γλx06−6λx˙02,q=1312αλ+3β−18αλ2x02−30αμx02−9βλ2x04−15βμx04−6γλ2x06−10γμx06−18λ2x˙02−30μx˙02,μ=6αλ2+3βλ+2γ−6αλ3x02−3βλ3x04−2γλ3x06−6λ3x˙0212−2α+6αλx02+3βλx04+2γλx06+6λx˙02.The value ofv0 is found from the condition x0=x0:(13)v0=±1−λx02±λx02−12−4μx042μx02.On the other hand, the solution to i.v.p. (5) is given by(14)vt=ccnqc2+pt+cn−1v0c,c2q2qc2+p,c2q2qc2+p,where(15)c=±−p±p+qv022+2qv˙02q.The number(16)Δ=p+qv022+2qv˙02,is called the discriminant to Duffing equation (5). In the case when Δ>0, this solution may also be written in the form(17)vt=v0cntω|m+v˙0/ωdntω|msntω|m1+p+qv02−ω/2ωsntω|m2,where(18)ω=−p2m−1,m=121±pp+qv022+2qv˙02.Suppose thatΔ<0. Define(19)δ=2pv02+qv04+2v˙02.SinceΔ<0, necessarily q<0. From the equality(20)δ=p2−Δ−q,it is evident that δ>0.Then, the solution to i.v.p. (5) reads(21)vt=δ−q4−2δ/−q41+u0cnω¯x|m¯+u˙0/ωsnω¯x|m¯dnω¯x|m¯/1+b¯snω¯x|m¯2,where(22)ω¯=Δ¯,m¯=12±α¯2Δ¯,b¯=α¯+β¯u022Δ¯−12,u0=A+v0A−v0andu˙0=2Av˙02A−v02,Δ¯=α¯+u02β¯2+2β¯u˙02,α¯=123A2q−p,β¯=12A2q+p,A=δ−q4,δ=2pv02+qv04+2v˙02.Observe that(23)Δ¯=α¯+u02β¯2+2β¯u˙02=2p2−qδ+−qδ>0.Example 1.
Let us consider the i.v.p.(24)x¨−1.25819x−2.40487x3+0.299493x5=0,x0=−0.261431andx′0=0.770472.
The exact solution is given by(25)xexactt=vt1−1.2914v2t+0.231351v4t,where(26)vt=−0.250491cn0.916322t|0.0909919+0.738957dn0.916322t|0.0909919sn0.916322t|0.09099191−0.0823122sn20.916322t|0.0909919.
See Figure1.Figure 1
Comparison between the exact solution and the Runge–Kutta numerical solution for i.v.p. (24).


### 2.2. Solution by Means of He’s Frequency Method

Let(27)fx=αx+βx3+γx5.He’s method assumes the solution in the ansatz form(28)xt=Acosωt+cos−1x0A,A≠0.The frequency is evaluated by means of the formula(29)ω=fxxatx=32A.Then,(30)ω=ωHe=α+3β4A2+9γ16A4.The numberA is found from the initial condition x′0=x˙0 so that(31)A=±x˙02α+9A4γ/16+3A2β/4+x02.


### 2.3. Solution by Means of a Simple Trigonometric Ansatz

As in He’s approach, we assume the solution in ansatz form (28) so that(32)x¨+αx+βx3+γx5=116A5γcos5θ=+116cos3θ5A5γ+4A3β+18Acosθ8α+5A4γ+6A2β−8ω.We choose the frequencyω so that(33)8α+5A4γ+6A2β−8ω=0.Then,(34)ω=ωtrigo=α+3β4A2+10γ16A4.This last formula looks like He’s formula (30). The difference is ωtrigo−ωHe=γ/16A4. This suggests to consider the following κ-parameter solution:(35)xκ,t=xt=Acosα+3β4A2+9γ16κA4t+cos−1x0A,where the number A is a solution to the sextic(36)116γλA6+3β4−116γλx02A4+α−3βx024A2−αx02−x˙02=0.The numberλ is chosen in order to get as small residual error as possible.


### 2.4. Solution by Means of an Improved Trigonometric Ansatz

#### 2.4.1. First Improved Ansatz

Let us consider the i.v.p.(37)x¨+αx+βx3+γx5=0,x0=Aandx′0=0.Assume the ansatz(38)xt=A1+λ+μcosωt1+λcos2ωt+μcos4ωt.Let(39)Rt=x″t+αxt+βx3t+γx5t.The numbersλ,μ, and ω are found from the conditions(40)R0=R″0=R″″0=0,λ=−24α2−αω+2A8γ2+5A6βγ+6αA4γ+3A4β2−A4γω+7αA2β−A2βω−3ω2α+A4γ+A2β4α+2A4γ+3A2β+2ω,μ=4α+2A4γ+3A2β−4ω4α+2A4γ+3A2β+2ω.The frequencyω is found from the quadratic equation(41)34α2+A42A2γ+3β5A2γ+6β+α40A4γ+51A2β−54α+2A4γ+3A2βω−14ω2=0.


#### 2.4.2. Second Improved Ansatz

Let us consider the i.v.p.(42)x¨+αx+βx3+γx5=0,x0=Aandx′0=0.Assume the ansatz(43)xt=A1+λ+μ+υcosωt1+λcos2ωt+μcos4ωt+υcos6ωt.The numbersλ,μ,υ, and ω are found from the conditions(44)R0=R″0=R40=R60=0,λ=−334α3+20α2ω−24αω2+10A12γ3+37A10βγ2+50αA8γ2+45A8β2γ+10A8γ2ω+118αA6βγ+18A6β3+25A6βγω+74α2A4γ+69αA4β2+30αA4γω+15A4β2ω−24A4γω2+85α2A2β+35αA2βω−24A2βω2−30ω3α+A4γ+A2β34α2+40αω+10A8γ2+27A6βγ+40αA4γ+18A4β2+20A4γω+51αA2β+30A2βω+16ω2,μ=334α2+10A8γ2+27A6βγ+40αA4γ+18A4β2+51αA2β−34ω234α2+40αω+10A8γ2+27A6βγ+40αA4γ+18A4β2+20A4γω+51αA2β+30A2βω+16ω2,υ=−34α2−20αω+10A8γ2+27A6βγ+40αA4γ+18A4β2−10A4γω+51αA2β−15A2βω−14ω234α2+40αω+10A8γ2+27A6βγ+40αA4γ+18A4β2+20A4γω+51αA2β+30A2βω+16ω2.The frequencyω is found from the cubic equation(45)−4α+3A2β+2A4γ124α2+186A2αβ+63A4β2+160A4αγ+102A6βγ+40A8γ2+844α+3A2β+2A4γω+160ω3=0.


### 2.5. Homotopy Method

Consider the homotopy(46)Hx,p=x¨+αx+pβx3+γx5,and assume the solution in then ansatz form(47)xt=y0ωt+py1ωt+p2y2ωt+p3y3ωt+p4y4ωt+p5y5ωt+⋯,where ω=α+pω1+p2ω2+p3ω2+p4ω2+p5ω2+⋯.Plugging the expression forxt into Hx,p and equating the coefficients of pjj=0,1,2,3,… will give an ode system. Solving this system so that no secular terms appear will give the following expressions: y.


### 2.6. Second Case: Quadratically Damped and Unforced Oscillator

Our aim is to giveapproximate analytical solution to i.v.p. (1). Define the residual function R=Rt as follows:(48)Rt=x¨+εx˙x˙+αx+βx3+γx5=x¨±εx˙2+αx+βx3+γx5.


#### 2.6.1. First Approach

Assume the ansatz(49)xt=c0e−ρtcosft+cos−1x0c0.Letθ=ft+cos−1x0/c0. We have(50)Rt=c08e4tρα+8e4tρρ2+6e2tρβc02+5γc04−8e4tρf′t2cosθ+116e−5tργcos5θf′tc05+116e−5tρc034e2tρβ+5γc02cos3θ+e−2tρερc02f′tsin2θ+12e−2tρεc02ρ2−f′t2cos2θ+12e−2tρεc02ρ2+f′t2+18e−5tρ+e−tρc02ρf′t−f″tsinθ.We will choose the functionf=ft so that(51)8e4tρα+8e4tρρ2−8e4tρf′t2+6e2tρβc02+5γc04=0,andf0=0andf′t>0.Then,(52)ft=122∫0t8α+8ρ2+6e−2ρτβc02+5e−4ρτγc04dτ=Ft−F0,where(53)Ft=8α+ρ2+6βc02e−2ρt+5γc04e−4ρt40γρ5γc04+6βc02e2ρt+8α+ρ2e4ρt,−52γ5γc04+6βc02e2ρt+8α+ρ2e4ρt+20γα+ρ2e2ρttanh−13βc02+8α+ρ2e2ρt22α+ρ25γc04+6βc02e2ρt+8α+ρ2e4ρt−310βe2ρttanh−15γc02+3βe2ρt5γ5γc04+6βc02e2ρt+8α+ρ2e4ρt.The value ofc0 is found from the initial condition x′0=x˙0, and it is a solution to the sextic(54)−8αx02+2ρ2x02+2ρx˙0x0+x˙02+24α+4ρ2−3βx02z2+6β−5γx02z4+5γz6=0.The numberρ is a free parameter that is chosen in order to minimize the residual error. In particular, when ε⟶0, we obtain approximate trigonometric solution to the undamped cubic-quintic Duffing equation(55)x¨+αx+βx3+γx5=0,x0=x0andx′0=x˙0.Example 2.
Let us consider the i.v.p.(56)x¨+0.25x˙x˙+x+5x3+10x5=0,x0=0andx′0=0.1.
The approximate analytical solution forρ=0.0084 is(57)uapproxt=xt=−0.0982083e−0.0084tcosft+π2,where the function ft is given by (52) with ε=0.25 , α=1, β=5,γ=10, c0=−0.0983, and ρ=0.0084 (see Figure 2).
The obtained results may be applied to solve the pendulum equation with quadratic damping(58)θ..+εθ.θ.+ω2sinθ=0,θ0=θ0andθ′0=θ.0.
Indeed, we may use the approximation(59)sinθ≈θ−θ36+θ5131forx≤π2,and then we replace i.v.p. (58) with the i.v.p.(60)θ..+εθ.θ.+ω2θ−θ36+θ5131=0,θ0=θ0andθ′0=θ.0.Figure 2
Comparison between analytical approximant (57) and the Runge–Kutta numerical solution for i.v.p. (56).Example 3.
Let us consider the i.v.p.(61)θ..+0.2θ.θ.+sinθ=0,θ=0=20°andθ′0=0.1.
The optimal value forρ is ρ=0.0225. The value of c0 is c0=−0.365609 (see Figure3).Figure 3
Comparison between the analytical approximant and the Runge–Kutta numerical solution for i.v.p. (61) with ρ=0.0225andc0=−0.365609.


#### 2.6.2. Second Approach

Let us consider the i.v.p.(62)x¨+εx˙x˙+αx+βx3+γx5=0,x0=x0andx′0=x˙0.Assume the ansatz(63)xt=xρ,κt=e−ρtyτ,τ=τt=1−exp−εκtεκ,where the function y=yτ is the exact solution to the i.v.p.(64)y¨τ+αyτ+βy3τ+γy5τ=0,y0=x0andx′0=x˙0+ρx0.The numbersρ and κ are free parameters that are chosen in order to minimize the residual error(65)Rt=x¨t+εx˙tx˙t+αxt+βx3t+γx5t.Observe that whenε⟶0, 1−exp−εκt/εκ⟶t and then we obtain the exact solution to the undamped and unforced cubic-quintic oscillator (3). So, we expect accurate approximate analytical solution for small ε. This approach is more accurate, but here the solution involves elliptic functions and the solution is not elementary.


#### 2.6.3. Third Approach

Let us consider the i.v.p.(66)x¨+εx˙x˙+αx+βx3+γx5=0,x0=x0andx′0=x˙0.Assume the ansatz(67)xt=xρ,λt=ce−ρtcosft+d1+λcos2ft+d.Proceeding in a similar way as in the first approach, we may choose(68)ft=∫0t2λ+2α+ρ2+βc2e−2ρτ10λ+4dτ.Here we have two free parametersρ and λ that are chosen in order to get as less residual error as possible. The numbers c and d are determined from the initial conditions as follows:(69)d=±cos−1±x0c2−λx02.The numberc is a solution to the octic(70)−2λ2λ+1λ+2x06α+ρ2+λx046αλ2+16αλ+8α+6λ2ρ2+16λρ2+8ρ2−βλ2x02−βλx02c2+−6αλ2x02−14αλx02−4αx02+3βλ2x04+2βλx04−6λ2ρ2x02−24λρ2x02−20λρx˙0x0−10λx˙02−8ρ2x02−8ρx˙0x0−4x˙02c4+2αλ+4α+2λρ2+4ρ2−3βλx02−βx02c6+βc8=0.


### 2.7. Third Case: Quadratically Damped and Forced Oscillator

Let us consider the i.v.p.(71)x¨+εx˙x˙+αx+βx3+γx5=F0cosωt+F1sinωt,x0=x0andx′0=x˙0.Assume the ansatz(72)xt=ut+d0cosωt+d1sinωt,where the function u=ut is a solution to some i.v.p.(73)u¨+εu˙u˙+pu+qu3+γu5=0,u0=x0−d0andu′0=x˙0−d1ω.The suitable constantsp, q, r, d0, and d1 are to be determined. Define the residual function(74)Rt=x¨t+εx˙tx˙t+αxt+βx3t+γx5t−F0cosωt−F1sinωt.The expressionRt contains many terms. Equating some coefficients of vt, v′t, sinjωt, and cosjωt (j=1,2,3,…), we obtain the following algebraic system:(75)8αd1+6βd13+6βd02d1+5γd15+10γd02d13+5γd04d1−8d1ω2−8F1=0,(76)8αd0+6βd03+6βd12d0+5γd05+10γd12d03+5γd14d0−8d0ω2−8F0=0,(77)8α+12βd02+12βd12+15γd04+30γd12d02+15γd14−8p=0,(78)β+5γd02+5γd12−q=0,(79)γ−r=0.From equations (77)–(79), we obtain(80).p=α+38d02+d124β+5γd02+d12,q=β+5γd02+d12,r=γ.Decupling equations (75) and (76) by means of their eliminants gives(81)5γF04+10γF12F02+5γF14d05+6βF04+6βF12F02d03+8αF04−8F04ω2d0−8F05=0,(82)5γF04+10γF12F02+5γF14d15+6βF14+6βF02F12d13+8αF14−8F14ω2d1−8F15=0.We choose the least in magnitude real rootd0 to quintic (81) and the least in magnitude real root d1 to quintic (82). Assuming that the forces F0 and F1 are small in magnitude, we have the following approximations for these roots:(83)d0=4F0α−ω224α−ω23+3βF02+F12,d1=4F1α−ω224α−ω23+3βF02+F12.More precise approximations are(84)d0=2F0α−ω2α−ω23/2±α−ω23+3βF02+F12,d1=2F1α−ω2α−ω23/2±α−ω23+3βF02+F12.


## 2.1. First Case: Undamped and Unforced Cubic-Quintic Duffing Equation

Let us consider the i.v.p.(3)x¨+αx+βx3+γx5=0,x0=x0andx′0=x˙0.Assume the ansatz(4)xt=vt1+λv2t+μv4t,where the function v=vt is the solution to some Duffing equation(5)v¨+pv+qv3=0,v0=v0andv′0=v˙0.The numbersv0 and v˙0 are determined from the initial conditions. Observe that(6)v˙0=−x˙0λv02+μv04+13/2μv04−1.From (5), it follows that(7)v˙2=v˙02+pv02+qv042−pv2−qv42.Observe that(8)12x˙2+12αx2+14βx4+16γx6−αx022−βx044−16γx06−x˙022=0.We have(9)0=12x˙2+12αx2+14βx4+16γx6−αx022−βx044−16γx06−x˙022=1121+λv2+μv4∑j=06Sjv2j,where(10)S0=6pv02+3qv04+6v˙02−6αx02−3βx04−2γx06−6x˙02,S1=−3−2α+2p+6αλx02+3βλx04+2γλx06+6λx˙02,S2=−3−4αλ−β+4μpv02+2μqv04+q+4μv˙02+6αλ2x02+6αμx02+3βλ2x04+3βμx04+2γλ2x06+2γμx06+6λ2x˙02+6μx˙02,S3=6αλ2+12αμ+3βλ+2γ+12μp−6αλ3x02−36αλμx02−3βλ3x04−18βλμx04−2γλ3x06−12γλμx06−6λ3x˙02−36λμx˙02,S4=3μ4αλ+β+2μpv02+μqv04+2q+2μv˙02−6αλ2x02−6αμx02−3βλ2x04−3βμx04−2γλ2x06−2γμx06−6λ2x˙02−6μx˙02,S5=−3μ2−2α+2p+6αλx02+3βλx04+2γλx06+6λx˙02,S6=−μ23q+6αμx02+3βμx04+2γμx06+6μx˙02.Eliminatingp, q, μ, v0, and v˙0 from the system S0=S1=S2=S3=S4=S5=S6=0 gives(11)6αx02+3βx04+2γx06+6x˙02λ3−6αλ2−3βλ−2γ46αx02+3βx04+2γx06+6x˙022λ3−24α6αx02+3βx04+2γx06+6x˙02λ2+312α2−6αβx02−3β2x04−2βγx06−6βx˙02λ+9αβ+6αγx02+3βγx04+2γ2x06+6γx˙02=0.The values forp, q, and μ are(12)p=122α−6αλx02−3βλx04−2γλx06−6λx˙02,q=1312αλ+3β−18αλ2x02−30αμx02−9βλ2x04−15βμx04−6γλ2x06−10γμx06−18λ2x˙02−30μx˙02,μ=6αλ2+3βλ+2γ−6αλ3x02−3βλ3x04−2γλ3x06−6λ3x˙0212−2α+6αλx02+3βλx04+2γλx06+6λx˙02.The value ofv0 is found from the condition x0=x0:(13)v0=±1−λx02±λx02−12−4μx042μx02.On the other hand, the solution to i.v.p. (5) is given by(14)vt=ccnqc2+pt+cn−1v0c,c2q2qc2+p,c2q2qc2+p,where(15)c=±−p±p+qv022+2qv˙02q.The number(16)Δ=p+qv022+2qv˙02,is called the discriminant to Duffing equation (5). In the case when Δ>0, this solution may also be written in the form(17)vt=v0cntω|m+v˙0/ωdntω|msntω|m1+p+qv02−ω/2ωsntω|m2,where(18)ω=−p2m−1,m=121±pp+qv022+2qv˙02.Suppose thatΔ<0. Define(19)δ=2pv02+qv04+2v˙02.SinceΔ<0, necessarily q<0. From the equality(20)δ=p2−Δ−q,it is evident that δ>0.Then, the solution to i.v.p. (5) reads(21)vt=δ−q4−2δ/−q41+u0cnω¯x|m¯+u˙0/ωsnω¯x|m¯dnω¯x|m¯/1+b¯snω¯x|m¯2,where(22)ω¯=Δ¯,m¯=12±α¯2Δ¯,b¯=α¯+β¯u022Δ¯−12,u0=A+v0A−v0andu˙0=2Av˙02A−v02,Δ¯=α¯+u02β¯2+2β¯u˙02,α¯=123A2q−p,β¯=12A2q+p,A=δ−q4,δ=2pv02+qv04+2v˙02.Observe that(23)Δ¯=α¯+u02β¯2+2β¯u˙02=2p2−qδ+−qδ>0.Example 1.
Let us consider the i.v.p.(24)x¨−1.25819x−2.40487x3+0.299493x5=0,x0=−0.261431andx′0=0.770472.
The exact solution is given by(25)xexactt=vt1−1.2914v2t+0.231351v4t,where(26)vt=−0.250491cn0.916322t|0.0909919+0.738957dn0.916322t|0.0909919sn0.916322t|0.09099191−0.0823122sn20.916322t|0.0909919.
See Figure1.Figure 1
Comparison between the exact solution and the Runge–Kutta numerical solution for i.v.p. (24).


## 2.2. Solution by Means of He’s Frequency Method

Let(27)fx=αx+βx3+γx5.He’s method assumes the solution in the ansatz form(28)xt=Acosωt+cos−1x0A,A≠0.The frequency is evaluated by means of the formula(29)ω=fxxatx=32A.Then,(30)ω=ωHe=α+3β4A2+9γ16A4.The numberA is found from the initial condition x′0=x˙0 so that(31)A=±x˙02α+9A4γ/16+3A2β/4+x02.


## 2.3. Solution by Means of a Simple Trigonometric Ansatz

As in He’s approach, we assume the solution in ansatz form (28) so that(32)x¨+αx+βx3+γx5=116A5γcos5θ=+116cos3θ5A5γ+4A3β+18Acosθ8α+5A4γ+6A2β−8ω.We choose the frequencyω so that(33)8α+5A4γ+6A2β−8ω=0.Then,(34)ω=ωtrigo=α+3β4A2+10γ16A4.This last formula looks like He’s formula (30). The difference is ωtrigo−ωHe=γ/16A4. This suggests to consider the following κ-parameter solution:(35)xκ,t=xt=Acosα+3β4A2+9γ16κA4t+cos−1x0A,where the number A is a solution to the sextic(36)116γλA6+3β4−116γλx02A4+α−3βx024A2−αx02−x˙02=0.The numberλ is chosen in order to get as small residual error as possible.


## 2.4. Solution by Means of an Improved Trigonometric Ansatz

### 2.4.1. First Improved Ansatz

Let us consider the i.v.p.(37)x¨+αx+βx3+γx5=0,x0=Aandx′0=0.Assume the ansatz(38)xt=A1+λ+μcosωt1+λcos2ωt+μcos4ωt.Let(39)Rt=x″t+αxt+βx3t+γx5t.The numbersλ,μ, and ω are found from the conditions(40)R0=R″0=R″″0=0,λ=−24α2−αω+2A8γ2+5A6βγ+6αA4γ+3A4β2−A4γω+7αA2β−A2βω−3ω2α+A4γ+A2β4α+2A4γ+3A2β+2ω,μ=4α+2A4γ+3A2β−4ω4α+2A4γ+3A2β+2ω.The frequencyω is found from the quadratic equation(41)34α2+A42A2γ+3β5A2γ+6β+α40A4γ+51A2β−54α+2A4γ+3A2βω−14ω2=0.


### 2.4.2. Second Improved Ansatz

Let us consider the i.v.p.(42)x¨+αx+βx3+γx5=0,x0=Aandx′0=0.Assume the ansatz(43)xt=A1+λ+μ+υcosωt1+λcos2ωt+μcos4ωt+υcos6ωt.The numbersλ,μ,υ, and ω are found from the conditions(44)R0=R″0=R40=R60=0,λ=−334α3+20α2ω−24αω2+10A12γ3+37A10βγ2+50αA8γ2+45A8β2γ+10A8γ2ω+118αA6βγ+18A6β3+25A6βγω+74α2A4γ+69αA4β2+30αA4γω+15A4β2ω−24A4γω2+85α2A2β+35αA2βω−24A2βω2−30ω3α+A4γ+A2β34α2+40αω+10A8γ2+27A6βγ+40αA4γ+18A4β2+20A4γω+51αA2β+30A2βω+16ω2,μ=334α2+10A8γ2+27A6βγ+40αA4γ+18A4β2+51αA2β−34ω234α2+40αω+10A8γ2+27A6βγ+40αA4γ+18A4β2+20A4γω+51αA2β+30A2βω+16ω2,υ=−34α2−20αω+10A8γ2+27A6βγ+40αA4γ+18A4β2−10A4γω+51αA2β−15A2βω−14ω234α2+40αω+10A8γ2+27A6βγ+40αA4γ+18A4β2+20A4γω+51αA2β+30A2βω+16ω2.The frequencyω is found from the cubic equation(45)−4α+3A2β+2A4γ124α2+186A2αβ+63A4β2+160A4αγ+102A6βγ+40A8γ2+844α+3A2β+2A4γω+160ω3=0.


## 2.4.1. First Improved Ansatz

Let us consider the i.v.p.(37)x¨+αx+βx3+γx5=0,x0=Aandx′0=0.Assume the ansatz(38)xt=A1+λ+μcosωt1+λcos2ωt+μcos4ωt.Let(39)Rt=x″t+αxt+βx3t+γx5t.The numbersλ,μ, and ω are found from the conditions(40)R0=R″0=R″″0=0,λ=−24α2−αω+2A8γ2+5A6βγ+6αA4γ+3A4β2−A4γω+7αA2β−A2βω−3ω2α+A4γ+A2β4α+2A4γ+3A2β+2ω,μ=4α+2A4γ+3A2β−4ω4α+2A4γ+3A2β+2ω.The frequencyω is found from the quadratic equation(41)34α2+A42A2γ+3β5A2γ+6β+α40A4γ+51A2β−54α+2A4γ+3A2βω−14ω2=0.


## 2.4.2. Second Improved Ansatz

Let us consider the i.v.p.(42)x¨+αx+βx3+γx5=0,x0=Aandx′0=0.Assume the ansatz(43)xt=A1+λ+μ+υcosωt1+λcos2ωt+μcos4ωt+υcos6ωt.The numbersλ,μ,υ, and ω are found from the conditions(44)R0=R″0=R40=R60=0,λ=−334α3+20α2ω−24αω2+10A12γ3+37A10βγ2+50αA8γ2+45A8β2γ+10A8γ2ω+118αA6βγ+18A6β3+25A6βγω+74α2A4γ+69αA4β2+30αA4γω+15A4β2ω−24A4γω2+85α2A2β+35αA2βω−24A2βω2−30ω3α+A4γ+A2β34α2+40αω+10A8γ2+27A6βγ+40αA4γ+18A4β2+20A4γω+51αA2β+30A2βω+16ω2,μ=334α2+10A8γ2+27A6βγ+40αA4γ+18A4β2+51αA2β−34ω234α2+40αω+10A8γ2+27A6βγ+40αA4γ+18A4β2+20A4γω+51αA2β+30A2βω+16ω2,υ=−34α2−20αω+10A8γ2+27A6βγ+40αA4γ+18A4β2−10A4γω+51αA2β−15A2βω−14ω234α2+40αω+10A8γ2+27A6βγ+40αA4γ+18A4β2+20A4γω+51αA2β+30A2βω+16ω2.The frequencyω is found from the cubic equation(45)−4α+3A2β+2A4γ124α2+186A2αβ+63A4β2+160A4αγ+102A6βγ+40A8γ2+844α+3A2β+2A4γω+160ω3=0.


## 2.5. Homotopy Method

Consider the homotopy(46)Hx,p=x¨+αx+pβx3+γx5,and assume the solution in then ansatz form(47)xt=y0ωt+py1ωt+p2y2ωt+p3y3ωt+p4y4ωt+p5y5ωt+⋯,where ω=α+pω1+p2ω2+p3ω2+p4ω2+p5ω2+⋯.Plugging the expression forxt into Hx,p and equating the coefficients of pjj=0,1,2,3,… will give an ode system. Solving this system so that no secular terms appear will give the following expressions: y.


## 2.6. Second Case: Quadratically Damped and Unforced Oscillator

Our aim is to giveapproximate analytical solution to i.v.p. (1). Define the residual function R=Rt as follows:(48)Rt=x¨+εx˙x˙+αx+βx3+γx5=x¨±εx˙2+αx+βx3+γx5.


### 2.6.1. First Approach

Assume the ansatz(49)xt=c0e−ρtcosft+cos−1x0c0.Letθ=ft+cos−1x0/c0. We have(50)Rt=c08e4tρα+8e4tρρ2+6e2tρβc02+5γc04−8e4tρf′t2cosθ+116e−5tργcos5θf′tc05+116e−5tρc034e2tρβ+5γc02cos3θ+e−2tρερc02f′tsin2θ+12e−2tρεc02ρ2−f′t2cos2θ+12e−2tρεc02ρ2+f′t2+18e−5tρ+e−tρc02ρf′t−f″tsinθ.We will choose the functionf=ft so that(51)8e4tρα+8e4tρρ2−8e4tρf′t2+6e2tρβc02+5γc04=0,andf0=0andf′t>0.Then,(52)ft=122∫0t8α+8ρ2+6e−2ρτβc02+5e−4ρτγc04dτ=Ft−F0,where(53)Ft=8α+ρ2+6βc02e−2ρt+5γc04e−4ρt40γρ5γc04+6βc02e2ρt+8α+ρ2e4ρt,−52γ5γc04+6βc02e2ρt+8α+ρ2e4ρt+20γα+ρ2e2ρttanh−13βc02+8α+ρ2e2ρt22α+ρ25γc04+6βc02e2ρt+8α+ρ2e4ρt−310βe2ρttanh−15γc02+3βe2ρt5γ5γc04+6βc02e2ρt+8α+ρ2e4ρt.The value ofc0 is found from the initial condition x′0=x˙0, and it is a solution to the sextic(54)−8αx02+2ρ2x02+2ρx˙0x0+x˙02+24α+4ρ2−3βx02z2+6β−5γx02z4+5γz6=0.The numberρ is a free parameter that is chosen in order to minimize the residual error. In particular, when ε⟶0, we obtain approximate trigonometric solution to the undamped cubic-quintic Duffing equation(55)x¨+αx+βx3+γx5=0,x0=x0andx′0=x˙0.Example 2.
Let us consider the i.v.p.(56)x¨+0.25x˙x˙+x+5x3+10x5=0,x0=0andx′0=0.1.
The approximate analytical solution forρ=0.0084 is(57)uapproxt=xt=−0.0982083e−0.0084tcosft+π2,where the function ft is given by (52) with ε=0.25 , α=1, β=5,γ=10, c0=−0.0983, and ρ=0.0084 (see Figure 2).
The obtained results may be applied to solve the pendulum equation with quadratic damping(58)θ..+εθ.θ.+ω2sinθ=0,θ0=θ0andθ′0=θ.0.
Indeed, we may use the approximation(59)sinθ≈θ−θ36+θ5131forx≤π2,and then we replace i.v.p. (58) with the i.v.p.(60)θ..+εθ.θ.+ω2θ−θ36+θ5131=0,θ0=θ0andθ′0=θ.0.Figure 2
Comparison between analytical approximant (57) and the Runge–Kutta numerical solution for i.v.p. (56).Example 3.
Let us consider the i.v.p.(61)θ..+0.2θ.θ.+sinθ=0,θ=0=20°andθ′0=0.1.
The optimal value forρ is ρ=0.0225. The value of c0 is c0=−0.365609 (see Figure3).Figure 3
Comparison between the analytical approximant and the Runge–Kutta numerical solution for i.v.p. (61) with ρ=0.0225andc0=−0.365609.


### 2.6.2. Second Approach

Let us consider the i.v.p.(62)x¨+εx˙x˙+αx+βx3+γx5=0,x0=x0andx′0=x˙0.Assume the ansatz(63)xt=xρ,κt=e−ρtyτ,τ=τt=1−exp−εκtεκ,where the function y=yτ is the exact solution to the i.v.p.(64)y¨τ+αyτ+βy3τ+γy5τ=0,y0=x0andx′0=x˙0+ρx0.The numbersρ and κ are free parameters that are chosen in order to minimize the residual error(65)Rt=x¨t+εx˙tx˙t+αxt+βx3t+γx5t.Observe that whenε⟶0, 1−exp−εκt/εκ⟶t and then we obtain the exact solution to the undamped and unforced cubic-quintic oscillator (3). So, we expect accurate approximate analytical solution for small ε. This approach is more accurate, but here the solution involves elliptic functions and the solution is not elementary.


### 2.6.3. Third Approach

Let us consider the i.v.p.(66)x¨+εx˙x˙+αx+βx3+γx5=0,x0=x0andx′0=x˙0.Assume the ansatz(67)xt=xρ,λt=ce−ρtcosft+d1+λcos2ft+d.Proceeding in a similar way as in the first approach, we may choose(68)ft=∫0t2λ+2α+ρ2+βc2e−2ρτ10λ+4dτ.Here we have two free parametersρ and λ that are chosen in order to get as less residual error as possible. The numbers c and d are determined from the initial conditions as follows:(69)d=±cos−1±x0c2−λx02.The numberc is a solution to the octic(70)−2λ2λ+1λ+2x06α+ρ2+λx046αλ2+16αλ+8α+6λ2ρ2+16λρ2+8ρ2−βλ2x02−βλx02c2+−6αλ2x02−14αλx02−4αx02+3βλ2x04+2βλx04−6λ2ρ2x02−24λρ2x02−20λρx˙0x0−10λx˙02−8ρ2x02−8ρx˙0x0−4x˙02c4+2αλ+4α+2λρ2+4ρ2−3βλx02−βx02c6+βc8=0.


## 2.6.1. First Approach

Assume the ansatz(49)xt=c0e−ρtcosft+cos−1x0c0.Letθ=ft+cos−1x0/c0. We have(50)Rt=c08e4tρα+8e4tρρ2+6e2tρβc02+5γc04−8e4tρf′t2cosθ+116e−5tργcos5θf′tc05+116e−5tρc034e2tρβ+5γc02cos3θ+e−2tρερc02f′tsin2θ+12e−2tρεc02ρ2−f′t2cos2θ+12e−2tρεc02ρ2+f′t2+18e−5tρ+e−tρc02ρf′t−f″tsinθ.We will choose the functionf=ft so that(51)8e4tρα+8e4tρρ2−8e4tρf′t2+6e2tρβc02+5γc04=0,andf0=0andf′t>0.Then,(52)ft=122∫0t8α+8ρ2+6e−2ρτβc02+5e−4ρτγc04dτ=Ft−F0,where(53)Ft=8α+ρ2+6βc02e−2ρt+5γc04e−4ρt40γρ5γc04+6βc02e2ρt+8α+ρ2e4ρt,−52γ5γc04+6βc02e2ρt+8α+ρ2e4ρt+20γα+ρ2e2ρttanh−13βc02+8α+ρ2e2ρt22α+ρ25γc04+6βc02e2ρt+8α+ρ2e4ρt−310βe2ρttanh−15γc02+3βe2ρt5γ5γc04+6βc02e2ρt+8α+ρ2e4ρt.The value ofc0 is found from the initial condition x′0=x˙0, and it is a solution to the sextic(54)−8αx02+2ρ2x02+2ρx˙0x0+x˙02+24α+4ρ2−3βx02z2+6β−5γx02z4+5γz6=0.The numberρ is a free parameter that is chosen in order to minimize the residual error. In particular, when ε⟶0, we obtain approximate trigonometric solution to the undamped cubic-quintic Duffing equation(55)x¨+αx+βx3+γx5=0,x0=x0andx′0=x˙0.Example 2.
Let us consider the i.v.p.(56)x¨+0.25x˙x˙+x+5x3+10x5=0,x0=0andx′0=0.1.
The approximate analytical solution forρ=0.0084 is(57)uapproxt=xt=−0.0982083e−0.0084tcosft+π2,where the function ft is given by (52) with ε=0.25 , α=1, β=5,γ=10, c0=−0.0983, and ρ=0.0084 (see Figure 2).
The obtained results may be applied to solve the pendulum equation with quadratic damping(58)θ..+εθ.θ.+ω2sinθ=0,θ0=θ0andθ′0=θ.0.
Indeed, we may use the approximation(59)sinθ≈θ−θ36+θ5131forx≤π2,and then we replace i.v.p. (58) with the i.v.p.(60)θ..+εθ.θ.+ω2θ−θ36+θ5131=0,θ0=θ0andθ′0=θ.0.Figure 2
Comparison between analytical approximant (57) and the Runge–Kutta numerical solution for i.v.p. (56).Example 3.
Let us consider the i.v.p.(61)θ..+0.2θ.θ.+sinθ=0,θ=0=20°andθ′0=0.1.
The optimal value forρ is ρ=0.0225. The value of c0 is c0=−0.365609 (see Figure3).Figure 3
Comparison between the analytical approximant and the Runge–Kutta numerical solution for i.v.p. (61) with ρ=0.0225andc0=−0.365609.


## 2.6.2. Second Approach

Let us consider the i.v.p.(62)x¨+εx˙x˙+αx+βx3+γx5=0,x0=x0andx′0=x˙0.Assume the ansatz(63)xt=xρ,κt=e−ρtyτ,τ=τt=1−exp−εκtεκ,where the function y=yτ is the exact solution to the i.v.p.(64)y¨τ+αyτ+βy3τ+γy5τ=0,y0=x0andx′0=x˙0+ρx0.The numbersρ and κ are free parameters that are chosen in order to minimize the residual error(65)Rt=x¨t+εx˙tx˙t+αxt+βx3t+γx5t.Observe that whenε⟶0, 1−exp−εκt/εκ⟶t and then we obtain the exact solution to the undamped and unforced cubic-quintic oscillator (3). So, we expect accurate approximate analytical solution for small ε. This approach is more accurate, but here the solution involves elliptic functions and the solution is not elementary.


## 2.6.3. Third Approach

Let us consider the i.v.p.(66)x¨+εx˙x˙+αx+βx3+γx5=0,x0=x0andx′0=x˙0.Assume the ansatz(67)xt=xρ,λt=ce−ρtcosft+d1+λcos2ft+d.Proceeding in a similar way as in the first approach, we may choose(68)ft=∫0t2λ+2α+ρ2+βc2e−2ρτ10λ+4dτ.Here we have two free parametersρ and λ that are chosen in order to get as less residual error as possible. The numbers c and d are determined from the initial conditions as follows:(69)d=±cos−1±x0c2−λx02.The numberc is a solution to the octic(70)−2λ2λ+1λ+2x06α+ρ2+λx046αλ2+16αλ+8α+6λ2ρ2+16λρ2+8ρ2−βλ2x02−βλx02c2+−6αλ2x02−14αλx02−4αx02+3βλ2x04+2βλx04−6λ2ρ2x02−24λρ2x02−20λρx˙0x0−10λx˙02−8ρ2x02−8ρx˙0x0−4x˙02c4+2αλ+4α+2λρ2+4ρ2−3βλx02−βx02c6+βc8=0.


## 2.7. Third Case: Quadratically Damped and Forced Oscillator

Let us consider the i.v.p.(71)x¨+εx˙x˙+αx+βx3+γx5=F0cosωt+F1sinωt,x0=x0andx′0=x˙0.Assume the ansatz(72)xt=ut+d0cosωt+d1sinωt,where the function u=ut is a solution to some i.v.p.(73)u¨+εu˙u˙+pu+qu3+γu5=0,u0=x0−d0andu′0=x˙0−d1ω.The suitable constantsp, q, r, d0, and d1 are to be determined. Define the residual function(74)Rt=x¨t+εx˙tx˙t+αxt+βx3t+γx5t−F0cosωt−F1sinωt.The expressionRt contains many terms. Equating some coefficients of vt, v′t, sinjωt, and cosjωt (j=1,2,3,…), we obtain the following algebraic system:(75)8αd1+6βd13+6βd02d1+5γd15+10γd02d13+5γd04d1−8d1ω2−8F1=0,(76)8αd0+6βd03+6βd12d0+5γd05+10γd12d03+5γd14d0−8d0ω2−8F0=0,(77)8α+12βd02+12βd12+15γd04+30γd12d02+15γd14−8p=0,(78)β+5γd02+5γd12−q=0,(79)γ−r=0.From equations (77)–(79), we obtain(80).p=α+38d02+d124β+5γd02+d12,q=β+5γd02+d12,r=γ.Decupling equations (75) and (76) by means of their eliminants gives(81)5γF04+10γF12F02+5γF14d05+6βF04+6βF12F02d03+8αF04−8F04ω2d0−8F05=0,(82)5γF04+10γF12F02+5γF14d15+6βF14+6βF02F12d13+8αF14−8F14ω2d1−8F15=0.We choose the least in magnitude real rootd0 to quintic (81) and the least in magnitude real root d1 to quintic (82). Assuming that the forces F0 and F1 are small in magnitude, we have the following approximations for these roots:(83)d0=4F0α−ω224α−ω23+3βF02+F12,d1=4F1α−ω224α−ω23+3βF02+F12.More precise approximations are(84)d0=2F0α−ω2α−ω23/2±α−ω23+3βF02+F12,d1=2F1α−ω2α−ω23/2±α−ω23+3βF02+F12.


## 3. Further Applications

Suppose we are given a quadratically damped oscillator(85)x¨+εx˙x˙+hx=0,x0=x0andx′0=x˙0,where h−x=−hx is an odd continuous function. Suppose that x≤M. Then, we may approximate the function hx by means of the following quintic polynomial:(86)hx≈αx+βx3+γx5,where(87)α=h−M/2−hM/2−5+33h−1/22−3M+5+33h1/22−3M−33−5h−1/22+3M−h1/22+3M32M,(88)β=2−8h−M/2+8hM/2+7+53h−1/22−3M−7+53h1/22−3M+53−7h−1/22+3M−h1/22+3M3M3,(89)γ=422h−M/2−2hM/2−1+3h−1/22−3M+1+3h1/22−3M−3−1h−1/22+3M−h1/22+3M3M5.Then, the i.v.p. (85) is reduced to the i.v.p. (1).Example 4.
Consider the motion of a satellite along a path that is equidistant from two identical massive stars with mutually interacting gravitational fields. If the distance between the two stars is 2d and the coordinate of the satellite motion isx, then the equation of motion of the satellite is given as(90)x¨+εx˙x˙+2mxd2+x23/2=0,x0=x0andx′0=x˙0,where m is the mass of a star and the restoring force is(91)hx=2mxd2+x23/2.
The nonlinear restoring force is an irrational force because of the bottom square root. The restoring force spikes near the origin. The spikes indicate the point when the satellite is most influenced by the mutual gravitational field of the stars. Away from the origin, the restoring force decreases gradually and approaches the horizontal axis asymptotically. This means that the satellite is far away from the stars and experiences a much smaller gravitational force.


## 4. Conclusions

We have obtained approximate analytical solutions to the quadratically damped Duffing oscillator equation by means of an elementary approach. We introduced a parameter technique that allowed us to optimize the obtained solution. The results are also valid for the linear quadratically damped oscillatorx¨+εx˙x˙+αx=0. Also, a more general quadratically damped oscillator x¨+εx˙x˙+hx=0 may be solved for any odd parity function hx. We also show the way to solve quadratically damped forced oscillators having the form x¨+εx˙x˙+hx=Ft for any continuous functions hx and Ft with h−x=−hx.The quadratically damped cubic-quintic oscillator having both forcing term and quadratic damping term has been analyzed analytically using some highly accurate approaches. The proposed analytical techniques may be applied to solve other strongly nonlinear oscillators.



---

*Source: 8125305-2022-09-13.xml*