# Analytical Approximant to a Quadratically Damped Duffing Oscillator

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

---

## Abstract

The Duffing oscillator of a system with strong quadratic damping is considered. We give an 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. We also solve the oscillator by menas of He’s homotopy method and the famous Krylov–Bogoliubov–Mitropolsky method. The approximant allows estimating the points at which the solution crosses the horizontal axis.


---

## Body

## 1. Introduction

In the standard textbooks, usually the systems with linear damping are considered. Due to their simplicity and the existence of an exact analytical solution, the problem is discussed in details. Unfortunately, in reality, the systems and damping are usually not linear. In recent times, a number of articles have appeared in the literature which deal with the phenomenon of a linear oscillator subject to a quadratic damping force [1–7]. Most elementary textbooks deal with viscous damping for the obvious reason that it involves a linear dependance on the velocity of the oscillator and presents the simplest situation where an exact analytical treatment is possible. In general, this involves the analysis of a second-order ordinary differential equation (ODE) of the Liénard type, namely,(1)x¨+fxx˙+gx=0.Nonlinear equations of motion such as this are seldom addressed in intermediate instruction in classical dynamics; this one is problematic because it cannot be solved in terms of elementary functions. 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 Duffing oscillator (fx=εx˙ and gx=αx+βx3):(2)x¨+εx˙x˙+αx+βx3=0,x0=x0andx′0=x˙0.The quadratically damped oscillator (2) is never critically damped or overdamped, and that to first order in the damping constant, the oscillation frequency is identical to the natural frequency. In the abscence of damping, we obtain the Duffing equation(3)x¨+αx+βx3=0,x0=x0andx′0=x˙0.Equation (3) admits the exact analytical solution [8].(4)xt=x0cnωt|m+x˙0|ωsnωt|mdnωt|m1+bsnωt|m2,where(5)ω=Δ,m=121−αΔ,b=12α+βx02Δ−1,Δ=α+βx022+2βx02≠0.


## 2. Solution Procedure

In what follows, we will assume that ifx=xt is a solution to the i.v.p. (2), then(6)limt⟶∞xt=0.Our aim is to give an approximate analytical solution to the i.v.p. (2). The residual function R=Rt is defined as follows:(7)Rt=x¨+εx˙x˙+αx+βx3=x¨+δεx˙2+αx+βx3,δ=±1.


### 2.1. First Approach

The ansatz is assumed as(8)xt=c0e−ρtcosωt+cos−1x0c0.We have(9)Rt=14c0cosθe−ρt4α+3βc02e−2ρt+4ρ2−4ω++14βc03cos3θe−3ρt+12c02δεcos2θρ2−ωe−2ρt++c02δερωsin2θe−2ρt+12c02δερ2+ωe−2ρt+2c0ρωsinθe−ρt,δ=±1.Since for smallt,e−2ρt≈1, we will choose the value of ω so that(10)4α+3βc02+4ρ2−4ω=0.Then,(11)ω=α+34βc02+ρ2.The value ofc0 is found from the initial condition x′0=x˙0, and it reads(12)c0=±3βx02−4α−4ρ2±16α+ρ22+3β8x02α+3ρ2+3βx04+32ρx˙0x0+16x˙026β.In the case whenβ⟶0, we define(13)c0=±x02α+2ρ2+2ρx˙0x0+x˙02α+ρ2.The numberρ is a free parameter that is chosen in order to minimize the residual error.


### 2.2. Second Approach

#### 2.2.1. He’s Homotopy Method

We will approximate the expressionεx˙x˙ by means of the formula(14)εx˙x˙≈r0x˙+r1x˙3+r2x˙5,x˙≤M,where(15)r0=−3−3M36ε,r1=2/383−93Mε,r2=−42/323−33M3ε.The homotopy equation is defined as(16)Hx,p=x¨+αx+pr0x˙+r1x˙3+r2x˙5+βx3.Following He’s approach [9–16], we assume the solution in the ansatz form(17)xt=exp−μptr0cosωt+B+py1ωt+B,whereω=α+pω1.Letτ=ωt+B. Then,(18)Hx,p=116−2αr08r0+5α2r04r2+6αr02r1−16μsinτ+4r03r02β−4ω1cosτ+5α5/2r05r2+4α3/2r03r1sin3τ+4r03βcos3τ−α5/2r05r2sin5τ+16αy1τ+16αy1′′τ.p+⋯.We define(19)ω1=3r02β4,μ=1168r0+5α2r04r2+6αr02r1.Then,(20)Hx,p=5α5/2r05r2+4α3/2r03r1sin3τ+4r03βcos3τ−α5/2r05r2sin5τ+16αy1τ+16αy1′′τ.Solving the odeHx,p=0 for y1=y1τ gives(21)y1τ=12r03βcos3τ+3α3/2r035αr02r2+4r1sin3τ−α5/2r05r2sin5τ384α.The first-order homotopy approximation will then be(22)r0coswt+1384α12r03βcos3wt+3α3/2r035αr02r2+4r1sin3wt−α5/2r05r2sin5wt,wherewt=α+3β4r02t+B,xt=exp−1168r0+6αr1r02+5α2r2r04t.The constantsr0 and r1 are determined from the initial conditions(23)x0=x0,x′0=x˙0.


### 2.3. Third Approach

#### 2.3.1. The Krylov–Bogoliubov–Mitropolsky Method (KBM)

The Krylov–Bogoliubov–Mitropolsky Method (KBM) is a technique to give an approximate analytical solution to the weakly nonlinear second-order equation(24)d2udt2+ω02u=εfu,dudt.Whenε=0, the solution of (24) may be expressed as(25)u=acosω0+θ,where a and θ are constants. For the case when ε>0 is small, Krylov and Bogoliubov (1947) assumed that the solution is still given by (25) but with time-varying a and θ and subject to the condition(26)dudt=−aω0sinϕ,ϕ=ω0t+θ.In the general case, the solution is assumed in the ansatz form(27)u=acosψ+∑n=1Nεnuna,ψ+OεN+1,where each un is a periodic function of ψ with a period 2π and a and ψ are assumed to vary with time according to(28)dadt=∑n=1NεnAna+OεN+1,dψdt=ω0+∑n=1Nεnψna+OεN+1.In order to uniquely determineAn and ψn, we require that no un contains cosψ. Let N=3. Then,(29)dudt=−aω0sinψ+ω0u1,ψ−aψ1sinψ+A1cosψε++A1u1,a+ω0u2,ψ+ψ1u1,ψ−aψ2sinψ+A2cosψε2++A2u1,a+A1u2,a+ω0u3,ψ+ψ2u1,ψ+ψ1u2,ψ−aψ3sinψ+A3cosψε3+⋯,d2udt2=−aω02cosψ+ω02u1,ψψ−2aψ1ω0cosψ−2A1ω0sinψε++2A1ω0u1,aψ+2ψ1ω0u1,ψψ+ω02u2,ψψ+sinψ−aA1ψ˙1−2A1ψ1−2A2ω0+cosψA1A˙1−a2ψ2ω0+ψ12ε2++A˙1A1u1,a+A12u1,aa+2A1ψ1u1,aψ+2A1ω0u2,aψ+2A2ω0u1,aψ+A1ψ˙1u1,ψ+2ψ2ω0u1,ψψ+2ψ1ω0u2,ψψ+ψ12u1,ψψ+ω02u3,ψψ+sinψ−aA2ψ˙1−aA1ψ˙2−2A2ψ1−2A1ψ2−2A3ω0+cosψ−2aψ3ω0−2aψ1ψ2+A2A˙1+A1A˙2ε3+⋯.Here,(30)ψ˙1=ψ1′t,u2,ψ=∂u2∂ψ,u2,ψψ=∂2u2∂ψ2,u1,aψ=∂2u1∂ψ∂a, etc.Let us consider the i.v.p.(31)x¨+r0x˙+r1x˙3+r2x˙5+αx+βx3=0,x0=x0,x′0=x˙0for0≤t≤T.The values forr0, r1, and r2 are given by (15). Using KBM, we obtain the following odes for determining a=at and ψ=ψt:(32)a˙=−12εs0a+3εβs0−2α2s116αa3+ε−195β2s0−96α2βs1−992α3ε2s0s12−320α4s2+1602α3ε2s02s21024α2a5−−ε189β2s1+117α3ε2s13+720α2βs2+370α3ε2s0s1s28192αa7−εs21701β2+1782α3ε2s12+2170α3ε2s0s273728a9−−29875α4ε3s1s22196608a11−34625α5ε3s23589824a13.(33)ψ˙=8α−ε2s028α+3β8α−ε2s0264α3/2a2+−15β2−57αβε2s0s1+9α3ε2s12+40α3ε2s0s2256α3/2a4++341β3−223α3βε2s12−310α3βε2s0s2+400α5ε2s1s28192α5/2a6++α3/2ε2s2−471βs1+2200α2s224576a8+13091α5/2βε2s22196608a10.(34)s0=−3−3M36,s1=2/383−93M,s2=−42/323−33M3,whereM=max0≤t≤Tx˙t.The expression for the KBM solution is(35)x=xt=acosψ++a33072α2β−2016a2αβ2+1251a4β3−576αβε2s02+2304α3ε2s0s1−3456a2α2βε2s0s1+2592a2α4ε2s12−477a4α3βε2s12−4756a4α3βε2s0s2+4080a4α5ε2s1s2+321a6α4βε2s1s2+2000a6α6ε2s22+1300a8α5βε2s22cos3ψ98304α3++a3ε884736αβs0−1741824a2β2s0+1179648α3s1+884736a2α2βs1−943488a4αβ2s1−221184α2ε2s02s1+580608a2α3ε2s0s12+383616a4α4ε2s13+1474560a2α4s2+1787904a4α3βs2−257184a6α2β2s2−552960a2α3ε2s02s2+2572800a4α4ε2s0s1s2+1939680a6α5ε2s12s2+2790400a6α5ε2s0s22+2912400a8α6ε2s1s22+1226125a10α7ε2s23sin3ψ37748736α5/2−−a7−144β3+3024α3βε2s12+2672α3βε2s0s2−3840α5ε2s1s2+7548a2α4βε2s1s2−5200a2α6ε2s22+285a4α5βε2s22cos7ψ4718592α3−−a7ε−1440β2s1+1728α3ε2s13+3072α2βs2−1362a2αβ2s2+9040α3ε2s0s1s2+16530a2α4ε2s12s2+9200a2α4ε2s0s22+30065a4α5ε2s1s22+15900a6α6ε2s23sin7ψ4718592α3/2++a9αε2s22268βs1−400α2s2+2477a2αβs2cos9ψ7864320+a9εs2−450β2+954α3ε2s12+784α3ε2s0s2+3417a2α4ε2s1s2+2580a4α5ε2s22sin9ψ4718592α−−1177a11α2βε2s22cos11ψ47185920−a11α7/2ε3s223564s1+4915a2αs2sin11ψ113246208+1175a13α9/2ε3s23sin13ψ792723456.The odes (32) and (33) may be solved numerically. However, since we are interested in obtaining analytical solutions, we may limit ourselves to the following approximate solutions:(36)a=at=2Ar04r0er0t+3αr1A2er0t−1,ψ=ψt=αt+B−β2α3/2r1r0t+loger0t+3A2er0t−1αr14r0,where the values for r0, r1, and r2 are given by (15). The constants A and B are determined from the initial conditions.


### 2.4. Fourth Approach

We assume the ansatz(37)xt=c0e−ρtcosft+cos−1x0c0,f0=0.We have(38)Rt=14c0exp−3ρt4e2ρtα+3c02β+4e2ρtρ2−4e2ρtf′t2cosθ+2c0eρtδερ2+c02βcos3θ+2c0eρtδεf′t2+2c0eρtδεcos2θρ2−f′t2+4e2ρt2ρf′t−f″tsinθ+4c0eρtδερf′tsin2θ,δ=±1.We will choose the functionf=ft so that(39)4e2ρtα+3c02β+4e2ρtρ2−4e2ρtf′t2=0,f′t>0.Then,(40)ft=Ft−F0,Ft=−4α+ρ2+3βc02e−2ρt2βc0ρ4α+ρ2e2ρt/βc02+3βc04α+ρ2e2ρtβc02+3−2α+ρ2eρtsinh−12α+ρ2eρt3βc0.The value ofc0 is found from the initial condition x′0=x˙0, and it reads(41)c0=±3βx02−4α−4ρ2±16ρ4+8ρ24α+9βx02+4α+3βx022+96βρx0x˙0+48βx˙026β.In the case whenβ⟶0, we define(42)c0=±12x02+x02α+3ρ2+4ρx˙0x0+2x˙02α+ρ2.


## 2.1. First Approach

The ansatz is assumed as(8)xt=c0e−ρtcosωt+cos−1x0c0.We have(9)Rt=14c0cosθe−ρt4α+3βc02e−2ρt+4ρ2−4ω++14βc03cos3θe−3ρt+12c02δεcos2θρ2−ωe−2ρt++c02δερωsin2θe−2ρt+12c02δερ2+ωe−2ρt+2c0ρωsinθe−ρt,δ=±1.Since for smallt,e−2ρt≈1, we will choose the value of ω so that(10)4α+3βc02+4ρ2−4ω=0.Then,(11)ω=α+34βc02+ρ2.The value ofc0 is found from the initial condition x′0=x˙0, and it reads(12)c0=±3βx02−4α−4ρ2±16α+ρ22+3β8x02α+3ρ2+3βx04+32ρx˙0x0+16x˙026β.In the case whenβ⟶0, we define(13)c0=±x02α+2ρ2+2ρx˙0x0+x˙02α+ρ2.The numberρ is a free parameter that is chosen in order to minimize the residual error.


## 2.2. Second Approach

### 2.2.1. He’s Homotopy Method

We will approximate the expressionεx˙x˙ by means of the formula(14)εx˙x˙≈r0x˙+r1x˙3+r2x˙5,x˙≤M,where(15)r0=−3−3M36ε,r1=2/383−93Mε,r2=−42/323−33M3ε.The homotopy equation is defined as(16)Hx,p=x¨+αx+pr0x˙+r1x˙3+r2x˙5+βx3.Following He’s approach [9–16], we assume the solution in the ansatz form(17)xt=exp−μptr0cosωt+B+py1ωt+B,whereω=α+pω1.Letτ=ωt+B. Then,(18)Hx,p=116−2αr08r0+5α2r04r2+6αr02r1−16μsinτ+4r03r02β−4ω1cosτ+5α5/2r05r2+4α3/2r03r1sin3τ+4r03βcos3τ−α5/2r05r2sin5τ+16αy1τ+16αy1′′τ.p+⋯.We define(19)ω1=3r02β4,μ=1168r0+5α2r04r2+6αr02r1.Then,(20)Hx,p=5α5/2r05r2+4α3/2r03r1sin3τ+4r03βcos3τ−α5/2r05r2sin5τ+16αy1τ+16αy1′′τ.Solving the odeHx,p=0 for y1=y1τ gives(21)y1τ=12r03βcos3τ+3α3/2r035αr02r2+4r1sin3τ−α5/2r05r2sin5τ384α.The first-order homotopy approximation will then be(22)r0coswt+1384α12r03βcos3wt+3α3/2r035αr02r2+4r1sin3wt−α5/2r05r2sin5wt,wherewt=α+3β4r02t+B,xt=exp−1168r0+6αr1r02+5α2r2r04t.The constantsr0 and r1 are determined from the initial conditions(23)x0=x0,x′0=x˙0.


## 2.2.1. He’s Homotopy Method

We will approximate the expressionεx˙x˙ by means of the formula(14)εx˙x˙≈r0x˙+r1x˙3+r2x˙5,x˙≤M,where(15)r0=−3−3M36ε,r1=2/383−93Mε,r2=−42/323−33M3ε.The homotopy equation is defined as(16)Hx,p=x¨+αx+pr0x˙+r1x˙3+r2x˙5+βx3.Following He’s approach [9–16], we assume the solution in the ansatz form(17)xt=exp−μptr0cosωt+B+py1ωt+B,whereω=α+pω1.Letτ=ωt+B. Then,(18)Hx,p=116−2αr08r0+5α2r04r2+6αr02r1−16μsinτ+4r03r02β−4ω1cosτ+5α5/2r05r2+4α3/2r03r1sin3τ+4r03βcos3τ−α5/2r05r2sin5τ+16αy1τ+16αy1′′τ.p+⋯.We define(19)ω1=3r02β4,μ=1168r0+5α2r04r2+6αr02r1.Then,(20)Hx,p=5α5/2r05r2+4α3/2r03r1sin3τ+4r03βcos3τ−α5/2r05r2sin5τ+16αy1τ+16αy1′′τ.Solving the odeHx,p=0 for y1=y1τ gives(21)y1τ=12r03βcos3τ+3α3/2r035αr02r2+4r1sin3τ−α5/2r05r2sin5τ384α.The first-order homotopy approximation will then be(22)r0coswt+1384α12r03βcos3wt+3α3/2r035αr02r2+4r1sin3wt−α5/2r05r2sin5wt,wherewt=α+3β4r02t+B,xt=exp−1168r0+6αr1r02+5α2r2r04t.The constantsr0 and r1 are determined from the initial conditions(23)x0=x0,x′0=x˙0.


## 2.3. Third Approach

### 2.3.1. The Krylov–Bogoliubov–Mitropolsky Method (KBM)

The Krylov–Bogoliubov–Mitropolsky Method (KBM) is a technique to give an approximate analytical solution to the weakly nonlinear second-order equation(24)d2udt2+ω02u=εfu,dudt.Whenε=0, the solution of (24) may be expressed as(25)u=acosω0+θ,where a and θ are constants. For the case when ε>0 is small, Krylov and Bogoliubov (1947) assumed that the solution is still given by (25) but with time-varying a and θ and subject to the condition(26)dudt=−aω0sinϕ,ϕ=ω0t+θ.In the general case, the solution is assumed in the ansatz form(27)u=acosψ+∑n=1Nεnuna,ψ+OεN+1,where each un is a periodic function of ψ with a period 2π and a and ψ are assumed to vary with time according to(28)dadt=∑n=1NεnAna+OεN+1,dψdt=ω0+∑n=1Nεnψna+OεN+1.In order to uniquely determineAn and ψn, we require that no un contains cosψ. Let N=3. Then,(29)dudt=−aω0sinψ+ω0u1,ψ−aψ1sinψ+A1cosψε++A1u1,a+ω0u2,ψ+ψ1u1,ψ−aψ2sinψ+A2cosψε2++A2u1,a+A1u2,a+ω0u3,ψ+ψ2u1,ψ+ψ1u2,ψ−aψ3sinψ+A3cosψε3+⋯,d2udt2=−aω02cosψ+ω02u1,ψψ−2aψ1ω0cosψ−2A1ω0sinψε++2A1ω0u1,aψ+2ψ1ω0u1,ψψ+ω02u2,ψψ+sinψ−aA1ψ˙1−2A1ψ1−2A2ω0+cosψA1A˙1−a2ψ2ω0+ψ12ε2++A˙1A1u1,a+A12u1,aa+2A1ψ1u1,aψ+2A1ω0u2,aψ+2A2ω0u1,aψ+A1ψ˙1u1,ψ+2ψ2ω0u1,ψψ+2ψ1ω0u2,ψψ+ψ12u1,ψψ+ω02u3,ψψ+sinψ−aA2ψ˙1−aA1ψ˙2−2A2ψ1−2A1ψ2−2A3ω0+cosψ−2aψ3ω0−2aψ1ψ2+A2A˙1+A1A˙2ε3+⋯.Here,(30)ψ˙1=ψ1′t,u2,ψ=∂u2∂ψ,u2,ψψ=∂2u2∂ψ2,u1,aψ=∂2u1∂ψ∂a, etc.Let us consider the i.v.p.(31)x¨+r0x˙+r1x˙3+r2x˙5+αx+βx3=0,x0=x0,x′0=x˙0for0≤t≤T.The values forr0, r1, and r2 are given by (15). Using KBM, we obtain the following odes for determining a=at and ψ=ψt:(32)a˙=−12εs0a+3εβs0−2α2s116αa3+ε−195β2s0−96α2βs1−992α3ε2s0s12−320α4s2+1602α3ε2s02s21024α2a5−−ε189β2s1+117α3ε2s13+720α2βs2+370α3ε2s0s1s28192αa7−εs21701β2+1782α3ε2s12+2170α3ε2s0s273728a9−−29875α4ε3s1s22196608a11−34625α5ε3s23589824a13.(33)ψ˙=8α−ε2s028α+3β8α−ε2s0264α3/2a2+−15β2−57αβε2s0s1+9α3ε2s12+40α3ε2s0s2256α3/2a4++341β3−223α3βε2s12−310α3βε2s0s2+400α5ε2s1s28192α5/2a6++α3/2ε2s2−471βs1+2200α2s224576a8+13091α5/2βε2s22196608a10.(34)s0=−3−3M36,s1=2/383−93M,s2=−42/323−33M3,whereM=max0≤t≤Tx˙t.The expression for the KBM solution is(35)x=xt=acosψ++a33072α2β−2016a2αβ2+1251a4β3−576αβε2s02+2304α3ε2s0s1−3456a2α2βε2s0s1+2592a2α4ε2s12−477a4α3βε2s12−4756a4α3βε2s0s2+4080a4α5ε2s1s2+321a6α4βε2s1s2+2000a6α6ε2s22+1300a8α5βε2s22cos3ψ98304α3++a3ε884736αβs0−1741824a2β2s0+1179648α3s1+884736a2α2βs1−943488a4αβ2s1−221184α2ε2s02s1+580608a2α3ε2s0s12+383616a4α4ε2s13+1474560a2α4s2+1787904a4α3βs2−257184a6α2β2s2−552960a2α3ε2s02s2+2572800a4α4ε2s0s1s2+1939680a6α5ε2s12s2+2790400a6α5ε2s0s22+2912400a8α6ε2s1s22+1226125a10α7ε2s23sin3ψ37748736α5/2−−a7−144β3+3024α3βε2s12+2672α3βε2s0s2−3840α5ε2s1s2+7548a2α4βε2s1s2−5200a2α6ε2s22+285a4α5βε2s22cos7ψ4718592α3−−a7ε−1440β2s1+1728α3ε2s13+3072α2βs2−1362a2αβ2s2+9040α3ε2s0s1s2+16530a2α4ε2s12s2+9200a2α4ε2s0s22+30065a4α5ε2s1s22+15900a6α6ε2s23sin7ψ4718592α3/2++a9αε2s22268βs1−400α2s2+2477a2αβs2cos9ψ7864320+a9εs2−450β2+954α3ε2s12+784α3ε2s0s2+3417a2α4ε2s1s2+2580a4α5ε2s22sin9ψ4718592α−−1177a11α2βε2s22cos11ψ47185920−a11α7/2ε3s223564s1+4915a2αs2sin11ψ113246208+1175a13α9/2ε3s23sin13ψ792723456.The odes (32) and (33) may be solved numerically. However, since we are interested in obtaining analytical solutions, we may limit ourselves to the following approximate solutions:(36)a=at=2Ar04r0er0t+3αr1A2er0t−1,ψ=ψt=αt+B−β2α3/2r1r0t+loger0t+3A2er0t−1αr14r0,where the values for r0, r1, and r2 are given by (15). The constants A and B are determined from the initial conditions.


## 2.3.1. The Krylov–Bogoliubov–Mitropolsky Method (KBM)

The Krylov–Bogoliubov–Mitropolsky Method (KBM) is a technique to give an approximate analytical solution to the weakly nonlinear second-order equation(24)d2udt2+ω02u=εfu,dudt.Whenε=0, the solution of (24) may be expressed as(25)u=acosω0+θ,where a and θ are constants. For the case when ε>0 is small, Krylov and Bogoliubov (1947) assumed that the solution is still given by (25) but with time-varying a and θ and subject to the condition(26)dudt=−aω0sinϕ,ϕ=ω0t+θ.In the general case, the solution is assumed in the ansatz form(27)u=acosψ+∑n=1Nεnuna,ψ+OεN+1,where each un is a periodic function of ψ with a period 2π and a and ψ are assumed to vary with time according to(28)dadt=∑n=1NεnAna+OεN+1,dψdt=ω0+∑n=1Nεnψna+OεN+1.In order to uniquely determineAn and ψn, we require that no un contains cosψ. Let N=3. Then,(29)dudt=−aω0sinψ+ω0u1,ψ−aψ1sinψ+A1cosψε++A1u1,a+ω0u2,ψ+ψ1u1,ψ−aψ2sinψ+A2cosψε2++A2u1,a+A1u2,a+ω0u3,ψ+ψ2u1,ψ+ψ1u2,ψ−aψ3sinψ+A3cosψε3+⋯,d2udt2=−aω02cosψ+ω02u1,ψψ−2aψ1ω0cosψ−2A1ω0sinψε++2A1ω0u1,aψ+2ψ1ω0u1,ψψ+ω02u2,ψψ+sinψ−aA1ψ˙1−2A1ψ1−2A2ω0+cosψA1A˙1−a2ψ2ω0+ψ12ε2++A˙1A1u1,a+A12u1,aa+2A1ψ1u1,aψ+2A1ω0u2,aψ+2A2ω0u1,aψ+A1ψ˙1u1,ψ+2ψ2ω0u1,ψψ+2ψ1ω0u2,ψψ+ψ12u1,ψψ+ω02u3,ψψ+sinψ−aA2ψ˙1−aA1ψ˙2−2A2ψ1−2A1ψ2−2A3ω0+cosψ−2aψ3ω0−2aψ1ψ2+A2A˙1+A1A˙2ε3+⋯.Here,(30)ψ˙1=ψ1′t,u2,ψ=∂u2∂ψ,u2,ψψ=∂2u2∂ψ2,u1,aψ=∂2u1∂ψ∂a, etc.Let us consider the i.v.p.(31)x¨+r0x˙+r1x˙3+r2x˙5+αx+βx3=0,x0=x0,x′0=x˙0for0≤t≤T.The values forr0, r1, and r2 are given by (15). Using KBM, we obtain the following odes for determining a=at and ψ=ψt:(32)a˙=−12εs0a+3εβs0−2α2s116αa3+ε−195β2s0−96α2βs1−992α3ε2s0s12−320α4s2+1602α3ε2s02s21024α2a5−−ε189β2s1+117α3ε2s13+720α2βs2+370α3ε2s0s1s28192αa7−εs21701β2+1782α3ε2s12+2170α3ε2s0s273728a9−−29875α4ε3s1s22196608a11−34625α5ε3s23589824a13.(33)ψ˙=8α−ε2s028α+3β8α−ε2s0264α3/2a2+−15β2−57αβε2s0s1+9α3ε2s12+40α3ε2s0s2256α3/2a4++341β3−223α3βε2s12−310α3βε2s0s2+400α5ε2s1s28192α5/2a6++α3/2ε2s2−471βs1+2200α2s224576a8+13091α5/2βε2s22196608a10.(34)s0=−3−3M36,s1=2/383−93M,s2=−42/323−33M3,whereM=max0≤t≤Tx˙t.The expression for the KBM solution is(35)x=xt=acosψ++a33072α2β−2016a2αβ2+1251a4β3−576αβε2s02+2304α3ε2s0s1−3456a2α2βε2s0s1+2592a2α4ε2s12−477a4α3βε2s12−4756a4α3βε2s0s2+4080a4α5ε2s1s2+321a6α4βε2s1s2+2000a6α6ε2s22+1300a8α5βε2s22cos3ψ98304α3++a3ε884736αβs0−1741824a2β2s0+1179648α3s1+884736a2α2βs1−943488a4αβ2s1−221184α2ε2s02s1+580608a2α3ε2s0s12+383616a4α4ε2s13+1474560a2α4s2+1787904a4α3βs2−257184a6α2β2s2−552960a2α3ε2s02s2+2572800a4α4ε2s0s1s2+1939680a6α5ε2s12s2+2790400a6α5ε2s0s22+2912400a8α6ε2s1s22+1226125a10α7ε2s23sin3ψ37748736α5/2−−a7−144β3+3024α3βε2s12+2672α3βε2s0s2−3840α5ε2s1s2+7548a2α4βε2s1s2−5200a2α6ε2s22+285a4α5βε2s22cos7ψ4718592α3−−a7ε−1440β2s1+1728α3ε2s13+3072α2βs2−1362a2αβ2s2+9040α3ε2s0s1s2+16530a2α4ε2s12s2+9200a2α4ε2s0s22+30065a4α5ε2s1s22+15900a6α6ε2s23sin7ψ4718592α3/2++a9αε2s22268βs1−400α2s2+2477a2αβs2cos9ψ7864320+a9εs2−450β2+954α3ε2s12+784α3ε2s0s2+3417a2α4ε2s1s2+2580a4α5ε2s22sin9ψ4718592α−−1177a11α2βε2s22cos11ψ47185920−a11α7/2ε3s223564s1+4915a2αs2sin11ψ113246208+1175a13α9/2ε3s23sin13ψ792723456.The odes (32) and (33) may be solved numerically. However, since we are interested in obtaining analytical solutions, we may limit ourselves to the following approximate solutions:(36)a=at=2Ar04r0er0t+3αr1A2er0t−1,ψ=ψt=αt+B−β2α3/2r1r0t+loger0t+3A2er0t−1αr14r0,where the values for r0, r1, and r2 are given by (15). The constants A and B are determined from the initial conditions.


## 2.4. Fourth Approach

We assume the ansatz(37)xt=c0e−ρtcosft+cos−1x0c0,f0=0.We have(38)Rt=14c0exp−3ρt4e2ρtα+3c02β+4e2ρtρ2−4e2ρtf′t2cosθ+2c0eρtδερ2+c02βcos3θ+2c0eρtδεf′t2+2c0eρtδεcos2θρ2−f′t2+4e2ρt2ρf′t−f″tsinθ+4c0eρtδερf′tsin2θ,δ=±1.We will choose the functionf=ft so that(39)4e2ρtα+3c02β+4e2ρtρ2−4e2ρtf′t2=0,f′t>0.Then,(40)ft=Ft−F0,Ft=−4α+ρ2+3βc02e−2ρt2βc0ρ4α+ρ2e2ρt/βc02+3βc04α+ρ2e2ρtβc02+3−2α+ρ2eρtsinh−12α+ρ2eρt3βc0.The value ofc0 is found from the initial condition x′0=x˙0, and it reads(41)c0=±3βx02−4α−4ρ2±16ρ4+8ρ24α+9βx02+4α+3βx022+96βρx0x˙0+48βx˙026β.In the case whenβ⟶0, we define(42)c0=±12x02+x02α+3ρ2+4ρx˙0x0+2x˙02α+ρ2.


## 3. Analysis and Discussion

In this section, we will compare the accuracy of the solution methods using the different approaches described in the previous section.Example 1.
Let us consider the i.v.p.(43)x¨+0.2x˙x˙+x+2x3=0,x0=0,x′0=0.1.
The approximate analytical solution using the first approach (see formula (8)) for ρ=0.007191 is(44)xapproxt=0.0992665e−0.007191tsin1.00739t.
It is shown in Figure1.
The solution obtained by means of He’s homotopy method (see Figure2) equals(45)xHet=e−0.0089515t−0.0000427888sin3.00025t+7.85294−3.49e−6sin5.00041t+7.85294−0.0999944cos1.00008t+7.85294−1.03e−6cos3.00025t+7.85294−6.9e−8cos5.00041t+7.85294.
Using the KBM, we obtained the following solution (Figure3):(46)xKBMt=6.736e−9sin−1.89144log0.850343−1.e0.00345092t+4.2634+5.94212i−4.99347t+1.32e−6−1.43e−6e0.00345092t−sin−1.13486log0.850343−1.e0.00345092t+2.55804+3.56527i−2.99608t1.1e−6e0.00345092t−9.2e−7cos−1.13486log0.850343−1.e0.00345092t+2.55804+3.56527i−2.99608t+0.00844854−0.0198709e0.00345092t+0.0116841e0.00690184t,cos−0.378287log0.850343−1.e0.00345092t+0.85268+1.18842i−0.998695t0.0922354e0.00345092t−0.07843170.850343−e0.00345092t2.
Now, using the fourth approach, we get the solution (see Figure4)(47)xtrigt=−0.0992665e−0.00719tsin69.5414.00021+0.0591231e−0.01438t−4.8811e0.00719t202.977e0.01438t+34.00021+0.0591231e−0.01438tsinh−18.22551e0.00719t0.01478+1.e0.01438t+249.895.Figure 1
Comparison between the analytical approximant (44) (formula (8)) and the Runge–Kutta numerical solution for the i.v.p (43).Figure 2
He’s homotopy method-error comparison with the Runge–Kutta numerical solution.Figure 3
KBM-error comparison with the Runge–Kutta numerical solution.Figure 4
Method using the fourth approach-error comparison with the Runge–Kutta numerical solution.Example 2.
Let us consider the i.v.p.(48)x¨+0.2x˙x˙+x+10x3=0,x0=0,x′0=0.1.
The approximate analytical solution forρ=0.0084 using the fourth approach is(49)xapproxt=0.097e−0.0084tsin59.52312.067−4.0003+0.280334e−0.017t−0.6113956.672−e0.0084t42.81e0.017t+34.0003+0.28e−0.017tsinh−13.778e0.0084t0.28+4.0003e0.017t.
It is shown in Figure5.
The obtained results may be applied to solve the pendulum equation with quadratic damping(50)x¨+εx˙x˙+ω2sinx=0,x0=x0,x′0=x˙0.
Indeed, we may use the approximation(51)sinx≈x−319x3forx≤π3,and then, we replace i.v.p. (43) with the i.v.p.(52)x¨+εx˙x˙+ω2x−3ω219x3=0,x0=x0,x′0=x˙0.Figure 5
Comparison between the analytical approximant (49) and the Runge–Kutta numerical solution for the i.v.p (48).


## 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 oscillator x¨+εx˙x˙+αx=0. A similar approach may be employed to study the quadratically damped cubic-quintic oscillator x¨+εx˙x˙+αx+βx3+γx5=0. Also, a more general quadratically damped oscillator x¨+εx˙x˙+hx=0 may be solved for any odd parity function hx. In future work, we will study quadratically damped forced oscillators having the form x¨+εx˙x˙+hx=Ft for any continuous functions hx and Ft.



---

*Source: 3131253-2022-02-07.xml*