# Some Novel Solutions to a Quadratically Damped Pendulum Oscillator: Analytical and Numerical Approximations

Authors: Alvaro H. Salas; Wedad Albalawi; M. R. Alharthi; S. A. El-Tantawy
Journal: Complexity
 (2022)

---

## Abstract

In this paper, some novel analytical and numerical techniques are introduced for solving and analyzing nonlinear second-order ordinary differential equations (ODEs) that are associated to some strongly nonlinear oscillators such as a quadratically damped pendulum equation. Two different analytical approximations are obtained: for the first approximation, the ansatz method with the help of Chebyshev approximate polynomial is employed to derive an approximation in the form of trigonometric functions. For the second analytical approximation, a novel hybrid homotopy with Krylov–Bogoliubov–Mitropolsky method (HKBMM) is introduced for the first time for analyzing the evolution equation. For the numerical approximation, both the finite difference method (FDM) and Galerkin method (GM) are presented for analyzing the strong nonlinear quadratically damped pendulum equation that arises in real life, such as nonlinear phenomena in plasma physics, engineering, and so on. Several examples are discussed and compared to the Runge–Kutta (RK) numerical approximation to investigate and examine the accuracy of the obtained approximations. Moreover, the accuracy of all obtained approximations is checked by estimating the maximum residual and distance errors.


---

## Body

## 1. Introduction

Duffing-type equation is one of the most important second-degree differential equations that is used to describe many different phenomena [1–6]. The Duffing equation can be used for describing a nonlinear oscillator with a cubic nonlinearity, and the standard form of this equation reads as x¨t+fx=0, with fx=∑i=1∞αixi being the only odd polynomial where i=1,3,5,…. George Duffing, a German engineer, is the first person who did arrive at this equation and used it in the study of many different oscillators [3]. He also prepared a book in this regard and explained in it many applications that use this equation in the interpretation of many natural phenomena. Since then, there has been a tremendous amount of research works done about this equation of motion and some related equations, including (un)damped Duffing oscillator x¨t+βx˙t+fx=0, forced Duffing oscillator x¨t+fx=Ft, forced damped Duffing oscillator x¨t+βx˙t+fx=Ft, and many other oscillators with odd polynomials and complicated damping term [7–10]. Moreover, there is another type of oscillator that combines both odd and even polynomials, which is called the Helmholtz–Duffing (HD) oscillator x¨t+fx+gx=0 (here, fx=∑i=1∞αixi is only odd polynomial where i=1,3,5,… and gx=∑i=2∞γixi is only even polynomial where i=2,4,6,…) and some related oscillators such as (un)damped HD oscillator x¨t+βx˙t+fx+gx=0, forced HD oscillator x¨t+fx+gx=Ft, forced damped HD oscillator x¨t+βx˙t+fx+gx=Ft and many other HD oscillators with complicated damping term and complicated polynomials [11–15]. All these oscillations have several applications in various fields of science, e.g., oscillations in electronic circuits, oscillations in different plasma models, pendulum oscillator, etc. Due to the importance of these equations, many studies have been conducted to find some analytical and numerical solutions to accurately describe the engineering and physical systems associated with these oscillations [17–19].As a contribution to the literature, in this article, we present some novel analytical and numerical solutions to the complicated damped HD-type oscillator for a given arbitrary initial conditions by means of both elliptic (exact solution) and trigonometric functions (approximate solution). First, we follow the work of Sugie [20] where the author obtained the equation of motion of underwater pendulum and studied the stability of this oscillator. This equation is called the quadratically damped pendulum equation [20]:(1)θ..+2εθ.θ.+ω02sinθ=0,where ε represents the coefficient of the damping term and ω0 indicates restoring coefficient per unit of the moment of inertia. For small θ, equation (1) can be approximated as follows:(2)θ..+2εθ.θ.+ω02θ=0.Numerous oscillators with quadrature damping have been investigated over a wide range of different fields [21–26]. There are many methods for solving nonlinear differential equations. There are many analytical and numerical methods that dealt with solving different differential equations, and some of these methods can be found in Refs. [28–34]. In this paper, we will consider four different methods for solving and analyzing equation (1). First, we will solve this equation using the effective ansatz method in order to find some analytical approximations. In the second method, the hybrid homotopy Krylov–Bogoliubov–Mitropolsky method (HKBMM) will be employed to find an approximate solution with high accuracy. On the other hand, two highly accurate numerical schemes which are called the finite difference method (FDM) and Galerkin Hats method (GHM) will be introduced for analyzing evolution equation (1).


## 2. Analytical Approximations

In this section, two different approximations will be obtained. For the first approximation, the ansatz method with the help of Chebyshev approximate polynomial is employed to obtain an approximation in trigonometric form. For the second approximation, the new HKBMM is introduced.


### 2.1. First Approach: Ansatz Method and Trigonometric Solution

Let us rewrite evolution equation (1) in the form of the initial value problem (i.v.p.):(3)θ..+2εθ.θ.+ω02sinθ=0,θ0=θ0andθ′0=θ.0.Based on Chebyshev polynomial approximation, the value ofsinθ can be expanded as(4)sinθ≈θ−λθ3,θ.θ.≈r0θ.+r1θ.3,forθ.≤M,where(5)r0=M21−12,r1=1M2−2,λ=213.Other possible choices forr0,r1,λ can be considered as(6)r0=5M16,r1=3548M,λ=16.Forω0,θ.0=1,0, the following approximation is obtained:(7)M=−0.900775εθ0+0.240105ε+0.905583θ0.Next, we replace the original i.v.p. (3) by the following approximate i.v.p.(8)ℝ≡θ..+2εr0θ.+r1θ.3+ω02θ−λθ3=0,θ0=θ0θ′0=θ.0.Assume that the solutions to the i.v.p. (8) have the following formulas:(9)θTrigon=c0exp−ρtcosft+cos−1θ0c0,θTrigon=c0exp−ρtsinft+sin−1θ0c0,with the initial conditions (ICs)(10)f0=0,θ0=θ0,θ′0=θ.0,where the number ρ and the function f≡ft are chosen in order to get the smallest possible residual error Rt:(11)Rt=θ″+2εθ′θ′+ω02θ−λθ3.Now by substituting ansatz (9) into the i.v.p. (8), we get(12)ℝ=Ϝ1sinϕ+Ϝ2cosϕ+Ϝ3sin3ϕ+Ϝ4cos3ϕ,b2−4ac,with(13)Ϝ1=12c0e−3ερt−3c02εf′3r1−3c02ε3ρ2r1+4εr0e2ερt−4ερe2ερtf′−2f″e2ερt,Ϝ2=14c0e−3ερt−6c02ε2ρr1+4e2ερtf′2−3c02λω02−6c02ε4ρ3r1−8ε2ρr0e2ερt+4ε2ρ2e2ερt+4ω02e2ερt,Ϝ3=−12c03εr1e−3ερt3ε2ρ2−f′2f′,Ϝ4=−14c03e−3ερt−6ε2ρr1f′2+λω02+2ε4ρ3r1,where ϕ=ft+cos−1θ0/c0.ForϜ2=0, we get(14)f′2=4e2ερtε2ρρ−2r0+ω02−3c02λω02+2ε4ρ3r16c02ε2ρr1+4e2ερt=A+B1+Cexp2ρεt,with(15)A=ε2ρρ−2r0+ω02,B=−2ε2ρ2−ω02λ2ε2ρr1+1+2ε2ρr0,C=23c02ε2ρr1.Integrating equation (14) leads to(16)f=Ft−F0,with(17)Ft=1ερAtanh−1A+B/Ce2ερt+1A−A+Btanh−1A+B/Ce2ερt+1A+B.Inserting the value off given in equation (16) into ansatz (9) and applying the ICs θ′0=θ.0, the value of c0 can be determined from the following quartic equation:(18)42ε2ρ2θ02−2ε2ρr0θ02+θ02ω02+2ερθ0θ.0+θ.02+−4ε2ρ2+8ε2ρr0−4ω02−3λθ02ω02+12ε3ρ2r1θ0θ.0+6ε2ρr1θ.02c02+32ε4ρ3r1+λω02c04=0.Solution (9) θTrigon is presented in Figures 1(a) and 1(b) for θ0=0 and θ0=π/6 at ε,ω0,θ.0=0.2,1,0.1. Moreover, in the same figure, solution (9) is compared to the RK numerical approximation and the maximum distance error according to the following relation is calculated:(19)Error≡L∞=max0≤t≤TRK−θTrigon.Figure 1
Both trigonometric solution (9) and RK numerical approximation are plotted in θ,t−plane.
(a)(b)The maximum distance error according to relation (19) for θ0=0 and θ0=π/6 at ε,ω0,θ.0=0.2,1,0.1 is, respectively, estimated as(20)L∞|θ0=0=0.00515945,L∞|θ0=π6=0.0511852.It is noted that the accuracy of solution (9) becomes good and acceptable for small θ0, but for large value of θ0, the accuracy of solution (9) reduces as shown in Figure 1.Also, the maximum residual error is defined as(21)ETθ=max0≤t≤Tθ″t+2εθ′tθ′t+ω02sinθt.This is another form for the error to check the accuracy of the obtained approximations.


### 2.2. Second Approach: HKBMM

Let us consider the i.v.p.(22)x¨+ω02x=Ft,x,x˙,x0=x0andx′0=x˙0for0≤t≤T.Suppose that the physical problem described by (22) involves some small parameters ε1, ε2,…., εr. Let x≡xt be the solution to the i.v.p. (22) and assume that(23)Ft,x,x˙≡0whenε1=ε2=…=εr=0.The solutionx≡xt depends not only on t but also on the parameters ε1, ε2,…,εr, so that we can rewrite x≡xt as(24)x=xt;ε1,ε2,…,εr.Let us multiply each parameter by some other parameterp and consider the following p−parametric solution:(25)xp=xpt;pε1,pε2,…,pεr.Accordingly, the functionxp may be written in a power series as follows:(26)xp=u0+pu1+p2u2+⋯,where uk depends on t only, say uk=ukt.Based on Krylov–Bogoliubov–Mitropolsky method (KBMM), the solution of equation (22) is assumed to be(27)xp=acosψ+∑n=1Npnuna,ψ+OpN+1,where each un is a periodic function of ψ and a and ψ are assumed to vary with time according to(28)dadt≡a˙=∑n=1NpnAna+OpN+1,(29)dψdt≡ψ˙=ω0+∑n=1Npnψna+OpN+1,where a≡at and ψ≡ψt.Moreover, the hybrid homotopy KBMM (HKBMM) is suggested to be(30)Hxp,t=x¨p+ω02xp−pFt,xp,x˙p.The next step is to write the residualHpx,t as a power series in p:(31)Hxp,t=x¨p+ω02xp+pΥ1+p2Υ2+p3Υ3+⋯.For the determination of the unknown functionsun, ψn, An, and a, we equate to zero the coefficients Υn in equation (22) and then we can get a system of ODEs. To avoid the so-called secularity, we choose only the solutions that do not contain cosψ nor sinψ. For N=2 (the first approximation), we may use the following formulas (we neglected all terms containing pj for j≥2):(32)x˙=pω0u1,ψ−aψ1sinψ+A1cosψ−aω0sinψ,x¨=pω02u1,ψ,ψ−2aψ1ω0cosψ−2A1ω0sinψ−aω02cosψ,x2=12a2cos2ψ+a22+2apu1cosψ,x3=34a3cosψ+14a3cos3ψ+p32a2u1cos2ψ+3a2u12,x4=12a4cos2ψ+18a4cos4ψ+3a48+p3a3u1cosψ+a3u1cos3ψ,x5=58a5cosψ+516a5cos3ψ+116a5cos5ψ+p52a4u1cos2ψ+58a4u1cos4ψ+15a4u18,and(33)xx˙=paω0cosψu1,ψ+a2ψ1sinψ−cosψ+12aA1cos2ψ+aA12−au1ω0sinψ,−a2ω0sinψcosψ,x˙2=p−2aω02sinψu1,ψ+a2ψ1ω0+a2ψ1ω0−cos2ψ−aA1ω0sin2ψ−12a2ω02cos2ψ+12a2ω02,x˙3=p32a2ω03u1,ψ−32a2ω03cos2ψu1,ψ−94a3ψ1ω02sinψ+34a3ψ1ω02sin3ψ+34a2A1ω02cosψ−34a2A1ω02cos3ψ−34a3ω03sinψ+14a3ω03sin3ψ,x2x˙=p12a2ω0u1,ψ+12a2ω0cos2ψu1,ψ−14a3ψ1sinψ−14a3ψ1sin3ψ+34a2A1cosψ+14a2A1cos3ψ−a2u1ω0sin2ψ−14a3ω0sinψ−14a3ω0sin3ψ,with(34)u1,ψ,ψ=∂ψ2u1,u1,ψ=∂ψu1.The approximate analytical solution is obtained by puttingp=1. However, we may keep the parameter p and then we may use it as a residual minimization parameter. The optimal value to p will be near p=1.Now, the proposed method can be applied for investigating the i.v.p. (8):(35)θ..+ω02θ=Ft,θ,θ.,θ0=θ0andθ′0=θ.0,where in our case, x=θ and(36)Ft,θ,θ.=−2εr0θ.+r1θ.3+λω02θ3.Observe that whenε1=ε⟶0 and ε2=λ⟶0, we get F=0.In equations (27)–(29), for N=1 and λ=2/13, we have(37)xp=acosψ+pu1a,ψ,a˙=A1a,ψ˙=ω0+ψ1a.The homotopy to equation (35) is written as(38)Hpθ,t=θ..p+ω02θp−p−2εr0θ.+r1θ.3+213ω02θ3.The substitution of equation (37) into equation (38) leads to(39)Hpθ,t=144ω02u1,ψ,ψ+4u1ω02+14−3a3λω02−8aψ1ω0cosψ+14−6a3εr1ω03−8aεr0ω0−8A1ω0sinψ−14a3λω02cos3ψ+12a3εr1ω03sin3ψp+⋯.We must have(40)144ω02u1,ψ,ψ+4u1ω02+14−3a3λω02−8aψ1ω0cosψ+14−6a3εr1ω03−8aεr0ω0−8A1ω0sinψ−14a3λω02cos3ψ+12a3εr1ω03sin3ψ=0.The coefficients ofcosψ and sinψ must be vanished to eliminate the secularity. Accordingly, we have(41)ψ1a=−352a2ω0,A1a=−34a3εr1ω02−aεr0.Thus, equation (40) reduces to(42)144ω02u1,ψ,ψ+4u1ω02+14a3λω02cos3ψ+12a3εr1ω03sin3ψ=0.Solving equation (42), the following particular solution without any secularity terms is obtained:(43)u1a,ψ=−1208a3cos3ψ−13εr1ω0sin3ψ.From equations (37) and (41), the functions a≡at and ψ≡ψt can be determined:(44)a˙=−34εr1ω02a3−εr0a,ψ˙=ω0−352ω0a2.By solving system (44), we get(45)a=2r0/3r11+4r0/3A2r1ω02e2tεr0−1⋅ω0,and(46)ψ=ω0+r013ω0r1t+B−126εr1ω0loge2ερr0t+3A2r1ω024r0e2tεr0t−1.We finally get the analytical approximation in its first approximationp=1.(47)θ=acosψ−1208a3cos3ψ−13εr1ω0sin3ψ.We now introduce three optimization parameters by replacingε,ψ=ρε,κψ and then the approximation (47) can be modified to be(48)θp,ρ,κ=aρcosψρ,κ−p208aρ3cos3ψρ,κ−13εr1ω0sin3ψρ,κ,with(49)aρ=2r0/3r11+4r0/3A2r1ω02e2tερr0−1ω0,and(50)ψρ,κ=κω0t+B−loge2ερr0t+3A2r1ω02/4r0e2tερr0t−1−2ερr0t26ερr1ω0.The numbersρ, κ, and p are free parameters that we choose in order to minimize the residual error(51)Rt=θ″t+2εθ′tθ′t+ω02sinθt.The default parameter values areρ=κ=p=1. The constants A and B are determined from the initial conditions (ICs) θ0=θ0 and θ′0=θ0.Following the same procedure above, we can get some higher-order approximations. For example, forN=3, the following solution is obtained:(52)θ=acosψ+S1a3+S2a5+S3a7,where the coefficients S1,2,3 are defined in Appendix. The values of a,ψ associated to this solution can be determined from the following equations:(53)a˙=−εr0a−352εr0+13r1ω02a3−3ε5r0−16r1ω02+429ε2r0r12ω021664a5−27εr1ω0221+2197ε2r12ω02173056a7,and(54)ψ˙=ω0−ε2r022ω0+3ε2r02−2ω02104ω0a2+3ω0494ε2r0r1+507ε2r12ω02−510816a4+3ω037687ε2r12ω02−412249728a6.The approximate solution of the i.v.p. (35) using the HKBMM is introduced in Figures 2(a) and 2(b) for θ0=0 and θ0=π/6. Also, the maximum distance error L∞ is estimated for the two cases as follows:(55)L∞|θ0=0=0.00145602,L∞|θ0=π6=0.0209098.Figure 2
The approximate solutions using both HKBMM and RK numerical method are plotted inθ,t− plane.
(a)(b)It is observed that approximate solution of the i.v.p. (35) using the HKBMM is characterized by high accuracy and more stability for long time and for arbitrary values of θ0. Also, this approximation is better than the trigonometric solution (9) θTrigon as shown from Figures 1 and 2 as well as from the values of the errors.


## 2.1. First Approach: Ansatz Method and Trigonometric Solution

Let us rewrite evolution equation (1) in the form of the initial value problem (i.v.p.):(3)θ..+2εθ.θ.+ω02sinθ=0,θ0=θ0andθ′0=θ.0.Based on Chebyshev polynomial approximation, the value ofsinθ can be expanded as(4)sinθ≈θ−λθ3,θ.θ.≈r0θ.+r1θ.3,forθ.≤M,where(5)r0=M21−12,r1=1M2−2,λ=213.Other possible choices forr0,r1,λ can be considered as(6)r0=5M16,r1=3548M,λ=16.Forω0,θ.0=1,0, the following approximation is obtained:(7)M=−0.900775εθ0+0.240105ε+0.905583θ0.Next, we replace the original i.v.p. (3) by the following approximate i.v.p.(8)ℝ≡θ..+2εr0θ.+r1θ.3+ω02θ−λθ3=0,θ0=θ0θ′0=θ.0.Assume that the solutions to the i.v.p. (8) have the following formulas:(9)θTrigon=c0exp−ρtcosft+cos−1θ0c0,θTrigon=c0exp−ρtsinft+sin−1θ0c0,with the initial conditions (ICs)(10)f0=0,θ0=θ0,θ′0=θ.0,where the number ρ and the function f≡ft are chosen in order to get the smallest possible residual error Rt:(11)Rt=θ″+2εθ′θ′+ω02θ−λθ3.Now by substituting ansatz (9) into the i.v.p. (8), we get(12)ℝ=Ϝ1sinϕ+Ϝ2cosϕ+Ϝ3sin3ϕ+Ϝ4cos3ϕ,b2−4ac,with(13)Ϝ1=12c0e−3ερt−3c02εf′3r1−3c02ε3ρ2r1+4εr0e2ερt−4ερe2ερtf′−2f″e2ερt,Ϝ2=14c0e−3ερt−6c02ε2ρr1+4e2ερtf′2−3c02λω02−6c02ε4ρ3r1−8ε2ρr0e2ερt+4ε2ρ2e2ερt+4ω02e2ερt,Ϝ3=−12c03εr1e−3ερt3ε2ρ2−f′2f′,Ϝ4=−14c03e−3ερt−6ε2ρr1f′2+λω02+2ε4ρ3r1,where ϕ=ft+cos−1θ0/c0.ForϜ2=0, we get(14)f′2=4e2ερtε2ρρ−2r0+ω02−3c02λω02+2ε4ρ3r16c02ε2ρr1+4e2ερt=A+B1+Cexp2ρεt,with(15)A=ε2ρρ−2r0+ω02,B=−2ε2ρ2−ω02λ2ε2ρr1+1+2ε2ρr0,C=23c02ε2ρr1.Integrating equation (14) leads to(16)f=Ft−F0,with(17)Ft=1ερAtanh−1A+B/Ce2ερt+1A−A+Btanh−1A+B/Ce2ερt+1A+B.Inserting the value off given in equation (16) into ansatz (9) and applying the ICs θ′0=θ.0, the value of c0 can be determined from the following quartic equation:(18)42ε2ρ2θ02−2ε2ρr0θ02+θ02ω02+2ερθ0θ.0+θ.02+−4ε2ρ2+8ε2ρr0−4ω02−3λθ02ω02+12ε3ρ2r1θ0θ.0+6ε2ρr1θ.02c02+32ε4ρ3r1+λω02c04=0.Solution (9) θTrigon is presented in Figures 1(a) and 1(b) for θ0=0 and θ0=π/6 at ε,ω0,θ.0=0.2,1,0.1. Moreover, in the same figure, solution (9) is compared to the RK numerical approximation and the maximum distance error according to the following relation is calculated:(19)Error≡L∞=max0≤t≤TRK−θTrigon.Figure 1
Both trigonometric solution (9) and RK numerical approximation are plotted in θ,t−plane.
(a)(b)The maximum distance error according to relation (19) for θ0=0 and θ0=π/6 at ε,ω0,θ.0=0.2,1,0.1 is, respectively, estimated as(20)L∞|θ0=0=0.00515945,L∞|θ0=π6=0.0511852.It is noted that the accuracy of solution (9) becomes good and acceptable for small θ0, but for large value of θ0, the accuracy of solution (9) reduces as shown in Figure 1.Also, the maximum residual error is defined as(21)ETθ=max0≤t≤Tθ″t+2εθ′tθ′t+ω02sinθt.This is another form for the error to check the accuracy of the obtained approximations.


## 2.2. Second Approach: HKBMM

Let us consider the i.v.p.(22)x¨+ω02x=Ft,x,x˙,x0=x0andx′0=x˙0for0≤t≤T.Suppose that the physical problem described by (22) involves some small parameters ε1, ε2,…., εr. Let x≡xt be the solution to the i.v.p. (22) and assume that(23)Ft,x,x˙≡0whenε1=ε2=…=εr=0.The solutionx≡xt depends not only on t but also on the parameters ε1, ε2,…,εr, so that we can rewrite x≡xt as(24)x=xt;ε1,ε2,…,εr.Let us multiply each parameter by some other parameterp and consider the following p−parametric solution:(25)xp=xpt;pε1,pε2,…,pεr.Accordingly, the functionxp may be written in a power series as follows:(26)xp=u0+pu1+p2u2+⋯,where uk depends on t only, say uk=ukt.Based on Krylov–Bogoliubov–Mitropolsky method (KBMM), the solution of equation (22) is assumed to be(27)xp=acosψ+∑n=1Npnuna,ψ+OpN+1,where each un is a periodic function of ψ and a and ψ are assumed to vary with time according to(28)dadt≡a˙=∑n=1NpnAna+OpN+1,(29)dψdt≡ψ˙=ω0+∑n=1Npnψna+OpN+1,where a≡at and ψ≡ψt.Moreover, the hybrid homotopy KBMM (HKBMM) is suggested to be(30)Hxp,t=x¨p+ω02xp−pFt,xp,x˙p.The next step is to write the residualHpx,t as a power series in p:(31)Hxp,t=x¨p+ω02xp+pΥ1+p2Υ2+p3Υ3+⋯.For the determination of the unknown functionsun, ψn, An, and a, we equate to zero the coefficients Υn in equation (22) and then we can get a system of ODEs. To avoid the so-called secularity, we choose only the solutions that do not contain cosψ nor sinψ. For N=2 (the first approximation), we may use the following formulas (we neglected all terms containing pj for j≥2):(32)x˙=pω0u1,ψ−aψ1sinψ+A1cosψ−aω0sinψ,x¨=pω02u1,ψ,ψ−2aψ1ω0cosψ−2A1ω0sinψ−aω02cosψ,x2=12a2cos2ψ+a22+2apu1cosψ,x3=34a3cosψ+14a3cos3ψ+p32a2u1cos2ψ+3a2u12,x4=12a4cos2ψ+18a4cos4ψ+3a48+p3a3u1cosψ+a3u1cos3ψ,x5=58a5cosψ+516a5cos3ψ+116a5cos5ψ+p52a4u1cos2ψ+58a4u1cos4ψ+15a4u18,and(33)xx˙=paω0cosψu1,ψ+a2ψ1sinψ−cosψ+12aA1cos2ψ+aA12−au1ω0sinψ,−a2ω0sinψcosψ,x˙2=p−2aω02sinψu1,ψ+a2ψ1ω0+a2ψ1ω0−cos2ψ−aA1ω0sin2ψ−12a2ω02cos2ψ+12a2ω02,x˙3=p32a2ω03u1,ψ−32a2ω03cos2ψu1,ψ−94a3ψ1ω02sinψ+34a3ψ1ω02sin3ψ+34a2A1ω02cosψ−34a2A1ω02cos3ψ−34a3ω03sinψ+14a3ω03sin3ψ,x2x˙=p12a2ω0u1,ψ+12a2ω0cos2ψu1,ψ−14a3ψ1sinψ−14a3ψ1sin3ψ+34a2A1cosψ+14a2A1cos3ψ−a2u1ω0sin2ψ−14a3ω0sinψ−14a3ω0sin3ψ,with(34)u1,ψ,ψ=∂ψ2u1,u1,ψ=∂ψu1.The approximate analytical solution is obtained by puttingp=1. However, we may keep the parameter p and then we may use it as a residual minimization parameter. The optimal value to p will be near p=1.Now, the proposed method can be applied for investigating the i.v.p. (8):(35)θ..+ω02θ=Ft,θ,θ.,θ0=θ0andθ′0=θ.0,where in our case, x=θ and(36)Ft,θ,θ.=−2εr0θ.+r1θ.3+λω02θ3.Observe that whenε1=ε⟶0 and ε2=λ⟶0, we get F=0.In equations (27)–(29), for N=1 and λ=2/13, we have(37)xp=acosψ+pu1a,ψ,a˙=A1a,ψ˙=ω0+ψ1a.The homotopy to equation (35) is written as(38)Hpθ,t=θ..p+ω02θp−p−2εr0θ.+r1θ.3+213ω02θ3.The substitution of equation (37) into equation (38) leads to(39)Hpθ,t=144ω02u1,ψ,ψ+4u1ω02+14−3a3λω02−8aψ1ω0cosψ+14−6a3εr1ω03−8aεr0ω0−8A1ω0sinψ−14a3λω02cos3ψ+12a3εr1ω03sin3ψp+⋯.We must have(40)144ω02u1,ψ,ψ+4u1ω02+14−3a3λω02−8aψ1ω0cosψ+14−6a3εr1ω03−8aεr0ω0−8A1ω0sinψ−14a3λω02cos3ψ+12a3εr1ω03sin3ψ=0.The coefficients ofcosψ and sinψ must be vanished to eliminate the secularity. Accordingly, we have(41)ψ1a=−352a2ω0,A1a=−34a3εr1ω02−aεr0.Thus, equation (40) reduces to(42)144ω02u1,ψ,ψ+4u1ω02+14a3λω02cos3ψ+12a3εr1ω03sin3ψ=0.Solving equation (42), the following particular solution without any secularity terms is obtained:(43)u1a,ψ=−1208a3cos3ψ−13εr1ω0sin3ψ.From equations (37) and (41), the functions a≡at and ψ≡ψt can be determined:(44)a˙=−34εr1ω02a3−εr0a,ψ˙=ω0−352ω0a2.By solving system (44), we get(45)a=2r0/3r11+4r0/3A2r1ω02e2tεr0−1⋅ω0,and(46)ψ=ω0+r013ω0r1t+B−126εr1ω0loge2ερr0t+3A2r1ω024r0e2tεr0t−1.We finally get the analytical approximation in its first approximationp=1.(47)θ=acosψ−1208a3cos3ψ−13εr1ω0sin3ψ.We now introduce three optimization parameters by replacingε,ψ=ρε,κψ and then the approximation (47) can be modified to be(48)θp,ρ,κ=aρcosψρ,κ−p208aρ3cos3ψρ,κ−13εr1ω0sin3ψρ,κ,with(49)aρ=2r0/3r11+4r0/3A2r1ω02e2tερr0−1ω0,and(50)ψρ,κ=κω0t+B−loge2ερr0t+3A2r1ω02/4r0e2tερr0t−1−2ερr0t26ερr1ω0.The numbersρ, κ, and p are free parameters that we choose in order to minimize the residual error(51)Rt=θ″t+2εθ′tθ′t+ω02sinθt.The default parameter values areρ=κ=p=1. The constants A and B are determined from the initial conditions (ICs) θ0=θ0 and θ′0=θ0.Following the same procedure above, we can get some higher-order approximations. For example, forN=3, the following solution is obtained:(52)θ=acosψ+S1a3+S2a5+S3a7,where the coefficients S1,2,3 are defined in Appendix. The values of a,ψ associated to this solution can be determined from the following equations:(53)a˙=−εr0a−352εr0+13r1ω02a3−3ε5r0−16r1ω02+429ε2r0r12ω021664a5−27εr1ω0221+2197ε2r12ω02173056a7,and(54)ψ˙=ω0−ε2r022ω0+3ε2r02−2ω02104ω0a2+3ω0494ε2r0r1+507ε2r12ω02−510816a4+3ω037687ε2r12ω02−412249728a6.The approximate solution of the i.v.p. (35) using the HKBMM is introduced in Figures 2(a) and 2(b) for θ0=0 and θ0=π/6. Also, the maximum distance error L∞ is estimated for the two cases as follows:(55)L∞|θ0=0=0.00145602,L∞|θ0=π6=0.0209098.Figure 2
The approximate solutions using both HKBMM and RK numerical method are plotted inθ,t− plane.
(a)(b)It is observed that approximate solution of the i.v.p. (35) using the HKBMM is characterized by high accuracy and more stability for long time and for arbitrary values of θ0. Also, this approximation is better than the trigonometric solution (9) θTrigon as shown from Figures 1 and 2 as well as from the values of the errors.


## 3. Numerical Solution

In this section, some effective and highly accurate numerical schemes will be introduced for analyzing evolution equation (3). Both the finite difference method (FDM) and Galerkin Hats method (GHM) are presented below.


### 3.1. Numerical Approximation via FDM

First, let us discuss and apply the FDM on the general second-order ODE. Thus, the following general second-order is introduced:(56)x¨=Ft,x,x˙,x0=x0,x′0=x˙0and0≤t≤T.Choose some positive integern≥6 and divide the interval 0,T into n-subintervals by means of the knots ti=ih, where h=T/n/i=0,1,2,…,n. Then, the first and second-order derivatives can be approximated as follows:(57)x′ti≈−12xi−5+75xi−4−200xi−3+300xi−2−300xi−1+137xi60h,x″ti≈−10xi−5+61xi−4−156xi−3+214xi−2−154xi−1+45xi12h2.Consequently, the following discrete version to ODE (56) for i=5,6,… is obtained:(58)−10xi−5+61xi−4−156xi−3+214xi−2−154xi−1+45xi12h2=Fti,xi,−12xi−5+75xi−4−200xi−3+300xi−2−300xi−1+137xi60h.The values ofx1, x2, x3, and x4 are obtained from some numerical or approximate analytical solution to the i.v.p. (56). System (58) may be solved recursively.The above algorithm can be applied for analyzing the i.v.p. (note hereθt≡xt without loss of generality):(59)θ..=Ft,θ,θ.,θ0=θ0&θ′0=θ.0,with(60)Ft,θ,θ.=−2εθ.θ.−ω02sinθ.The numerical approximation using FDM is plotted in Figures3 and 4 for different values of n,θ0. In Figures 3(a) and 3(b), the FDM numerical approximation is plotted against n=150 and n=300, respectively. Moreover, the effect of θ0 on the numerical approximation is illustrated in Figures 4(a) and 4(b) for θ0=0 and θ0=π/6, respectively. Furthermore, the maximum distance error is calculated for all mentioned cases as follows:(61)L∞=0.0279013,forn,θ0=150,0,L∞=0.00173893,forn,θ0=300,0,L∞=0.0429752,forn,θ0=150,π6,L∞=0.00284604,forn,θ0=300,π6.Figure 3
The approximate solutions using both FDM and RK numerical method are plotted inθ,t−plane for (a)n,θ0=150,0 and (b) n,θ0=300,0.
(a)(b)Figure 4
The approximate solutions using both FDM and RK numerical method are plotted inθ,t−plane for (a)n,θ0=150,π/6 and (b) n,θ0=300,π/6.
(a)(b)It is clear that the accuracy of the FDM numerical approximation increases with increasingn. Also, this approximation is stable against the long time intervals and arbitrary angle θ0. Moreover, this approximation is better than the trigonometric solution (9) θTrigon.


### 3.2. Numerical Approximation via Galerkin Hats Method

First, let us consider a polynomial second-order forced damped i.v.p.(62)x¨+2εx˙+Px=0,x0=0andx′0=x˙0,where x≡xt.Let us consider the i.v.p.(63)y¨+2εy˙+Py=0,y0=y0andy′0=y˙0.Suppose thatx=xt;c1 is the solution to the i.v.p.(64)x¨+2εx˙+Px=0,x0=0andx′0=c1.Define(65)yt=xt+c0;c1.The constantsc0 and c1 are determined from the system(66)x0+c0;c1=y0andx′0+c0;c1=y˙0.Thus, problemproblem (61) reduces to (60) by problem (63) reduces to (62).Some particular cases to the i.v.p. (62) are defined as(67)x¨+2εx˙+nt+ptx=ft,x¨+2εx˙+nt+ptx+qtx2=ft,x¨+2εx˙+nt+ptx+qtx2+rtx3=ft.We will use the same idea for the linear case.Here, we start to discuss the linear case in system (67):(68)x¨+2εx˙+ptx=0,x0=0andx′0=x˙0.An approximate solution to the i.v.p. (68) is assumed to be in the following ansatz form:(69)x=∑k=1nckφkt,where the functions φkt are the so-called linear Galerkin hats.Let us investigate the present problem in the interval0≤t≤T and by choosing some positive integer n≥2 and define the step h=T/n and let ξj=jh=jT/n for j=0,1,2,⋯. The functions φkt for k=1,2,…,n are defined as(70)φkt=ct−ξk−1h,ifξk−1≤t≤ξk,−t−ξk+1h,ifξk≤t≤ξk+1.For an illustration, see Figure1.Figure1 Galerkin hats for n=7 and T=5.Some properties of these functions can be illustrated as follows:(71)φjtφktdt=0,for j−k≥2 and t∈0,T, and(72)∫0Tφjpt=2Tp+1n,for j≥1 and p=1,2,3,….In general, forj−k=1 and r,s=0,1,2,3,⋯, we have(73)∫0Tφjrtφkst=Tr!s!nr+s+1!.Using the formula(74)∫0TφjtxtNdt=hN+1N+2∑k=0N−1k+1cj−1N−k+cj+1N−kcjk+2N+1cjN,for any N≥0 and c0=cn+1=0, and assuming that ajt≡aj= const, we may evaluate easily the following integration:(75)∫0TPxφjtdt=∑N=0maN∫0TφjtxtNdtfor anyj,∫0Tx″tφjtdt=−∫0Tx′tφj′tdt=cj−1−2cj+cj+1hfor anyj,∫0Tx′tx′tφjtdt=n2Tcj−1−cj−22+cj−cj−12.Now, let us return to the original i.v.p. (3):(76)θ..+2εθ.θ.+ω02sinθ=0,θ0=0andθ′0=θ.0.Using the Chebyshev approximation,(77)sinθ≈θ−40241θ3+θ5131forθ≤π2.Another approximation may be obtained by minimizing the integral(78)∫−π2π2aθ+bθ3+cθ5−sinθ2dθ⟶min.The minimization procedure yields the values(79)min=π2−336011975040−2661120π2+171720π4−2664π6+13π8π11≈5.5×10−9,with(80)a=10523760−2592π2+19π4π7=0.99977,b=−252018480+17π2π2−120π9=−0.16583,c=1663201008−112π2+π4π11=7.5742×10−3,sinθ≈θ−16θ3+θ5132.Then, i.v.p. (76) can be reduced to the following approximate i.v.p.(81)θ..+2εθ.θ.+ω02θ−40241θ3+θ5131=0,θ0=0andθ′0=θ.0.Assume that the solution to i.v.p. (81) is given by(82)θ=∑k=1nckφkt.The Galerkin method solves the following algebraic system of nonlinear equations:(83)nTcj−2−2cj−1+cj+nTcj−2−cj−12+cj−cj−12ε+ω02Sj=0,with(84)Sj=Tω025544ncj5+Tω02cj−12772ncj4+Tω025cj−12−779240ncj3+Tω02cj−110cj−12−23113860ncj2+Tω0225cj−14−693cj−12+462027720ncj+Tω0227720n5cj−25+10cj−1cj−24+35cj−12−77cj−23+20cj−13−462cj−1cj−22+25cj−14−693cj−12+4620cj−2+12cj−15cj−14−154cj−12+1540,for j=2,3,…, where c0=0 and c1=T/nθ.0. System (83) can be solved recursively.Forj=2, the value of c2 can be determined from the following quintic equation:(85)2310n6εθ.0−4620n6θ.0∓2310n6εθ.0+1540n4T2ω02θ.0−154n2T4ω02θ.03+5T6ω02θ.052310n6+27720n6±27720n6ε+4620n4T2ω02−693n2T4ω02θ.02+25T6ω02θ.0427720n5Tc2−T2ω02θ.0231n2−10T2θ.0213860n4c22−Tω0277n2−5T2θ.029240n3c23+T2ω02θ.02772n2c24+Tω025544nc25=0.We choose the real root to equation (86) that is closest to c1. Suppose we already found the values c2, c3,…, ck−1 for some k≥2. Then, the value of ck is found by solving the quintic in equation (82) by letting j=k. However, since ck−ck−12=±ck−ck−1, we must solve two quintics in ck. We choose the real root to the two quintics that is closest to ck−1. Let us introduce the quintic(86)λ5ck5+λ4ck4+λ3ck3+λ2ck2+λ1ck+λ0=0,Any of the two quintics (82) for j=k. Then, for sufficiently large n, the value of ck may be estimated by means of the formula(87)z≈z0−λ1+5λ5z04+4λ4z03+3λ3z02+2λ2z0λ0+λ5z05+λ4z04+λ3z03+λ2z02+λ1z0λ12−λ0λ2+15λ52z08+24λ4λ5z07+10λ42z06+17λ3λ5z06+15λ3λ4z05+9λ2λ5z05+6λ32z04+9λ2λ4z04+8λ2λ3z03+2λ1λ4z03−10λ0λ5z03+3λ22z02+3λ1λ3z02−6λ0λ4z02+3λ1λ2z0−3λ0λ3z0where z0=ck−1.For arbitrary ICs, i.e., for any values toθ0,θ.0, the following ansatz is assumed:(88)θt=θ0φ0t+Tnθ.0φ2t+∑k=2nckφkt.Then, we get(89)∫0Tφjtθ″tdt=ncj−2−2cj−1+cjT,and(90)2ε∫0Tφjtθ′tθ′tdt=εnTcj−cj−12sgncj−cj−1−cj−2−cj−12sgncj−2−cj−1.Using the approximation(91)sinθ≈θ−16θ3+1132θ5,gives(92)∫0Tφjtsinθtdt≈Tω0227720n5cj−25−231cj−23+4620cj−2+60cj−15+5cj5−231cj3+4620cj+25cj−14cj−2+cj+4cj−135cj−22+5cj2−462+3cj−125cj−23−231cj−2+5cj3−231cj+2cj−15cj−24−231cj−22+5cj4−231cj2+9240.The following definitions are used in our analysis:(93)c0=θ0,c1=θ0+Tnθ.0,cj=0forj<0orj>n.The algebraic system to be solved forj=2,3,…,n reads as(94)ncj−2−2cj−1+cjT+εnTcj−cj−12sgncj−cj−1−cj−2−cj−12sgncj−2−cj−1+Tω0227720n5cj−25−231cj−23+4620cj−2+60cj−15+5cj5−231cj3+4620cj+25cj−14cj−2+cj+4cj−135cj−22+5cj2−462+3cj−125cj−23−231cj−2+5cj3−231cj+2cj−15cj−24−231cj−22+5cj4−231cj2+9240=0.The Galerkin numerical approximation to i.v.p. (3) is introduced in Figures 5 and 6 for different values of n,θ0,θ.0. Also, the maximum distance error is estimated for different values of n,θ0,θ.0 as follows:(95)L∞=0.0198616,forn,θ0=150,0,L∞=0.000434117,forn,θ0=1000,0,L∞=0.105879,forn,θ0=150,π6,L∞=0.0146958,forn,θ0=1000,π6.Figure 5
The approximate solutions using both GM and RK numerical method are plotted inθ,t−plane for (a)n,θ0=150,0 and (b)n,θ0=1000,0.
(a)(b)Figure 6
The approximate solutions using both GM and RK numerical method are plotted inθ,t−plane for (a)n,θ0=150,π/6 and (b)n,θ0=1000,π/6.
(a)(b)It is clear that increasing the number of hatsn leads to the increase of approximation accuracy, i.e., the error shrinks with the enhancement of the number of hats n. Moreover, it is observed that the Galerkin numerical approximation is characterized by high accuracy compared to RK numerical approximation.


## 3.1. Numerical Approximation via FDM

First, let us discuss and apply the FDM on the general second-order ODE. Thus, the following general second-order is introduced:(56)x¨=Ft,x,x˙,x0=x0,x′0=x˙0and0≤t≤T.Choose some positive integern≥6 and divide the interval 0,T into n-subintervals by means of the knots ti=ih, where h=T/n/i=0,1,2,…,n. Then, the first and second-order derivatives can be approximated as follows:(57)x′ti≈−12xi−5+75xi−4−200xi−3+300xi−2−300xi−1+137xi60h,x″ti≈−10xi−5+61xi−4−156xi−3+214xi−2−154xi−1+45xi12h2.Consequently, the following discrete version to ODE (56) for i=5,6,… is obtained:(58)−10xi−5+61xi−4−156xi−3+214xi−2−154xi−1+45xi12h2=Fti,xi,−12xi−5+75xi−4−200xi−3+300xi−2−300xi−1+137xi60h.The values ofx1, x2, x3, and x4 are obtained from some numerical or approximate analytical solution to the i.v.p. (56). System (58) may be solved recursively.The above algorithm can be applied for analyzing the i.v.p. (note hereθt≡xt without loss of generality):(59)θ..=Ft,θ,θ.,θ0=θ0&θ′0=θ.0,with(60)Ft,θ,θ.=−2εθ.θ.−ω02sinθ.The numerical approximation using FDM is plotted in Figures3 and 4 for different values of n,θ0. In Figures 3(a) and 3(b), the FDM numerical approximation is plotted against n=150 and n=300, respectively. Moreover, the effect of θ0 on the numerical approximation is illustrated in Figures 4(a) and 4(b) for θ0=0 and θ0=π/6, respectively. Furthermore, the maximum distance error is calculated for all mentioned cases as follows:(61)L∞=0.0279013,forn,θ0=150,0,L∞=0.00173893,forn,θ0=300,0,L∞=0.0429752,forn,θ0=150,π6,L∞=0.00284604,forn,θ0=300,π6.Figure 3
The approximate solutions using both FDM and RK numerical method are plotted inθ,t−plane for (a)n,θ0=150,0 and (b) n,θ0=300,0.
(a)(b)Figure 4
The approximate solutions using both FDM and RK numerical method are plotted inθ,t−plane for (a)n,θ0=150,π/6 and (b) n,θ0=300,π/6.
(a)(b)It is clear that the accuracy of the FDM numerical approximation increases with increasingn. Also, this approximation is stable against the long time intervals and arbitrary angle θ0. Moreover, this approximation is better than the trigonometric solution (9) θTrigon.


## 3.2. Numerical Approximation via Galerkin Hats Method

First, let us consider a polynomial second-order forced damped i.v.p.(62)x¨+2εx˙+Px=0,x0=0andx′0=x˙0,where x≡xt.Let us consider the i.v.p.(63)y¨+2εy˙+Py=0,y0=y0andy′0=y˙0.Suppose thatx=xt;c1 is the solution to the i.v.p.(64)x¨+2εx˙+Px=0,x0=0andx′0=c1.Define(65)yt=xt+c0;c1.The constantsc0 and c1 are determined from the system(66)x0+c0;c1=y0andx′0+c0;c1=y˙0.Thus, problemproblem (61) reduces to (60) by problem (63) reduces to (62).Some particular cases to the i.v.p. (62) are defined as(67)x¨+2εx˙+nt+ptx=ft,x¨+2εx˙+nt+ptx+qtx2=ft,x¨+2εx˙+nt+ptx+qtx2+rtx3=ft.We will use the same idea for the linear case.Here, we start to discuss the linear case in system (67):(68)x¨+2εx˙+ptx=0,x0=0andx′0=x˙0.An approximate solution to the i.v.p. (68) is assumed to be in the following ansatz form:(69)x=∑k=1nckφkt,where the functions φkt are the so-called linear Galerkin hats.Let us investigate the present problem in the interval0≤t≤T and by choosing some positive integer n≥2 and define the step h=T/n and let ξj=jh=jT/n for j=0,1,2,⋯. The functions φkt for k=1,2,…,n are defined as(70)φkt=ct−ξk−1h,ifξk−1≤t≤ξk,−t−ξk+1h,ifξk≤t≤ξk+1.For an illustration, see Figure1.Figure1 Galerkin hats for n=7 and T=5.Some properties of these functions can be illustrated as follows:(71)φjtφktdt=0,for j−k≥2 and t∈0,T, and(72)∫0Tφjpt=2Tp+1n,for j≥1 and p=1,2,3,….In general, forj−k=1 and r,s=0,1,2,3,⋯, we have(73)∫0Tφjrtφkst=Tr!s!nr+s+1!.Using the formula(74)∫0TφjtxtNdt=hN+1N+2∑k=0N−1k+1cj−1N−k+cj+1N−kcjk+2N+1cjN,for any N≥0 and c0=cn+1=0, and assuming that ajt≡aj= const, we may evaluate easily the following integration:(75)∫0TPxφjtdt=∑N=0maN∫0TφjtxtNdtfor anyj,∫0Tx″tφjtdt=−∫0Tx′tφj′tdt=cj−1−2cj+cj+1hfor anyj,∫0Tx′tx′tφjtdt=n2Tcj−1−cj−22+cj−cj−12.Now, let us return to the original i.v.p. (3):(76)θ..+2εθ.θ.+ω02sinθ=0,θ0=0andθ′0=θ.0.Using the Chebyshev approximation,(77)sinθ≈θ−40241θ3+θ5131forθ≤π2.Another approximation may be obtained by minimizing the integral(78)∫−π2π2aθ+bθ3+cθ5−sinθ2dθ⟶min.The minimization procedure yields the values(79)min=π2−336011975040−2661120π2+171720π4−2664π6+13π8π11≈5.5×10−9,with(80)a=10523760−2592π2+19π4π7=0.99977,b=−252018480+17π2π2−120π9=−0.16583,c=1663201008−112π2+π4π11=7.5742×10−3,sinθ≈θ−16θ3+θ5132.Then, i.v.p. (76) can be reduced to the following approximate i.v.p.(81)θ..+2εθ.θ.+ω02θ−40241θ3+θ5131=0,θ0=0andθ′0=θ.0.Assume that the solution to i.v.p. (81) is given by(82)θ=∑k=1nckφkt.The Galerkin method solves the following algebraic system of nonlinear equations:(83)nTcj−2−2cj−1+cj+nTcj−2−cj−12+cj−cj−12ε+ω02Sj=0,with(84)Sj=Tω025544ncj5+Tω02cj−12772ncj4+Tω025cj−12−779240ncj3+Tω02cj−110cj−12−23113860ncj2+Tω0225cj−14−693cj−12+462027720ncj+Tω0227720n5cj−25+10cj−1cj−24+35cj−12−77cj−23+20cj−13−462cj−1cj−22+25cj−14−693cj−12+4620cj−2+12cj−15cj−14−154cj−12+1540,for j=2,3,…, where c0=0 and c1=T/nθ.0. System (83) can be solved recursively.Forj=2, the value of c2 can be determined from the following quintic equation:(85)2310n6εθ.0−4620n6θ.0∓2310n6εθ.0+1540n4T2ω02θ.0−154n2T4ω02θ.03+5T6ω02θ.052310n6+27720n6±27720n6ε+4620n4T2ω02−693n2T4ω02θ.02+25T6ω02θ.0427720n5Tc2−T2ω02θ.0231n2−10T2θ.0213860n4c22−Tω0277n2−5T2θ.029240n3c23+T2ω02θ.02772n2c24+Tω025544nc25=0.We choose the real root to equation (86) that is closest to c1. Suppose we already found the values c2, c3,…, ck−1 for some k≥2. Then, the value of ck is found by solving the quintic in equation (82) by letting j=k. However, since ck−ck−12=±ck−ck−1, we must solve two quintics in ck. We choose the real root to the two quintics that is closest to ck−1. Let us introduce the quintic(86)λ5ck5+λ4ck4+λ3ck3+λ2ck2+λ1ck+λ0=0,Any of the two quintics (82) for j=k. Then, for sufficiently large n, the value of ck may be estimated by means of the formula(87)z≈z0−λ1+5λ5z04+4λ4z03+3λ3z02+2λ2z0λ0+λ5z05+λ4z04+λ3z03+λ2z02+λ1z0λ12−λ0λ2+15λ52z08+24λ4λ5z07+10λ42z06+17λ3λ5z06+15λ3λ4z05+9λ2λ5z05+6λ32z04+9λ2λ4z04+8λ2λ3z03+2λ1λ4z03−10λ0λ5z03+3λ22z02+3λ1λ3z02−6λ0λ4z02+3λ1λ2z0−3λ0λ3z0where z0=ck−1.For arbitrary ICs, i.e., for any values toθ0,θ.0, the following ansatz is assumed:(88)θt=θ0φ0t+Tnθ.0φ2t+∑k=2nckφkt.Then, we get(89)∫0Tφjtθ″tdt=ncj−2−2cj−1+cjT,and(90)2ε∫0Tφjtθ′tθ′tdt=εnTcj−cj−12sgncj−cj−1−cj−2−cj−12sgncj−2−cj−1.Using the approximation(91)sinθ≈θ−16θ3+1132θ5,gives(92)∫0Tφjtsinθtdt≈Tω0227720n5cj−25−231cj−23+4620cj−2+60cj−15+5cj5−231cj3+4620cj+25cj−14cj−2+cj+4cj−135cj−22+5cj2−462+3cj−125cj−23−231cj−2+5cj3−231cj+2cj−15cj−24−231cj−22+5cj4−231cj2+9240.The following definitions are used in our analysis:(93)c0=θ0,c1=θ0+Tnθ.0,cj=0forj<0orj>n.The algebraic system to be solved forj=2,3,…,n reads as(94)ncj−2−2cj−1+cjT+εnTcj−cj−12sgncj−cj−1−cj−2−cj−12sgncj−2−cj−1+Tω0227720n5cj−25−231cj−23+4620cj−2+60cj−15+5cj5−231cj3+4620cj+25cj−14cj−2+cj+4cj−135cj−22+5cj2−462+3cj−125cj−23−231cj−2+5cj3−231cj+2cj−15cj−24−231cj−22+5cj4−231cj2+9240=0.The Galerkin numerical approximation to i.v.p. (3) is introduced in Figures 5 and 6 for different values of n,θ0,θ.0. Also, the maximum distance error is estimated for different values of n,θ0,θ.0 as follows:(95)L∞=0.0198616,forn,θ0=150,0,L∞=0.000434117,forn,θ0=1000,0,L∞=0.105879,forn,θ0=150,π6,L∞=0.0146958,forn,θ0=1000,π6.Figure 5
The approximate solutions using both GM and RK numerical method are plotted inθ,t−plane for (a)n,θ0=150,0 and (b)n,θ0=1000,0.
(a)(b)Figure 6
The approximate solutions using both GM and RK numerical method are plotted inθ,t−plane for (a)n,θ0=150,π/6 and (b)n,θ0=1000,π/6.
(a)(b)It is clear that increasing the number of hatsn leads to the increase of approximation accuracy, i.e., the error shrinks with the enhancement of the number of hats n. Moreover, it is observed that the Galerkin numerical approximation is characterized by high accuracy compared to RK numerical approximation.


## 4. Conclusion

In this paper, the quadratically damped pendulum equation with strong nonlinearity has been solved and analyzed using some novel and effective analytical and numerical techniques. In the beginning, the ansatz method was devoted to find an analytical approximation to the mentioned equation in the form of trigonometric functions. Also, in this study, for the first time, a new hybrid homotopy with Krylov–Bogoliubov–Mitropolsky method (HKBMM) was applied for analyzing the evolution equation and deriving an analytical approximation with high accuracy. Moreover, both the finite difference method (FDM) and Galerkin method (GM) were employed for analyzing the present evolution equation as well as some related oscillators. The obtained approximations were graphically compared with each other. Furthermore, the maximum residual error for each approximation was estimated. In the GM, we derived the iterative schemes for finding the coefficients that appear in the linear Galerkin hat combination in the ansatz form solution for the evolution equation. These coefficients may be found iteratively. It was found that the numerical approximations are more accurate than analytical ones, but both give good accuracy. Also, it was observed that the obtained results become reasonably good for small initial speed. One of the most important features of Galerkin method is that it gives more stable approximations for any values to the physical parameters and for long time. Thus, this method can be devoted for studying and investigating different pendulum oscillators for any nonlinearity [16, 27].


### 4.1. Future Work

There are many oscillators that may be solved using the proposed method. Examples include(96)x¨+2εx˙+αx+βx2+γx3=FcosΩt,x¨+2εx˙+αx+βx+γx5=FcosΩt,x¨−ε1−x2x˙+αx+βx2+γx3=FcosΩt,x¨+2εx˙+α−Q0cosγtsinx=FcosΩt.


## 4.1. Future Work

There are many oscillators that may be solved using the proposed method. Examples include(96)x¨+2εx˙+αx+βx2+γx3=FcosΩt,x¨+2εx˙+αx+βx+γx5=FcosΩt,x¨−ε1−x2x˙+αx+βx2+γx3=FcosΩt,x¨+2εx˙+α−Q0cosγtsinx=FcosΩt.



---

*Source: 7803798-2022-05-28.xml*