# Novel Analytical and Numerical Approximations to the Forced Damped Parametric Driven Pendulum Oscillator: Chebyshev Collocation Method

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

---

## Abstract

In this work, some novel approximate analytical and numerical solutions to the forced damped driven nonlinear (FDDN) pendulum equation and some relation equations of motion on the pivot vertically for arbitrary angles are obtained. The analytical approximation is derived in terms of the Jacobi elliptic functions with arbitrary elliptic modulus. For the numerical approximations, the Chebyshev collocation numerical method is introduced for analyzing the equation of motion. Moreover, the analytical approximation and numerical approximation using the Chebyshev collocation numerical method and the MATHEMATICA command Fit are compared with the Runge–Kutta (RK) numerical solution. Also, the maximum distance error to all obtained approximations is estimated with respect to the RK numerical solution. The obtained results help many authors to understand the mechanism of many phenomena related to the plasma physics, classical mechanics, quantum mechanics, optical fiber, and electronic circuits.


---

## Body

## 1. Introduction

The pendulum oscillator and some related equation have been used as a physical model to solve several natural problems related to bifurcations, oscillations, and chaos such as nonlinear plasma oscillations [1–9], Duffing oscillators [10–14], and Helmholtz oscillations [12], and many other applications can be found in [15–24]. There are few attempts for analyzing the equation of motion of the nonlinear damped pendulum taking the friction forces into account [25]. The approximate solution was obtained in form of the Jacobi elliptic functions. However, there are many others forces in addition to the friction force that affect the motion of the pendulum such as perturbed and periodic forces. These forces appear in different dynamic systems and cannot be neglected due to their great impact on the behavior of the oscillator. For instance, the unforced damped driven nonlinear pendulum equation/or the unforced damped parametric driven pendulum equation(1)φ..+2βφ.+ϕtsinφ=0,has been derived in detail in [26], where ϕt=ω02−εω2cosγt, ω02=g/l, β=μ/2ml, ω1=γ/ml, ω2=γ2/l, ε≪1 is a small parameter, and φ≡φt denotes the angular displacement. In (1), ω0 indicates the eigenfrequency of the system and β represents the damping coefficient. Here, the pendulum is modeled by a sphere of mass m, hanging at the end of a massless wire with length l and fixed to a supporting point “O,” swinging to and from in a vertical plane under the gravity acceleration “g.”. For β,ω2,F=0,0,0, the unforced undamped nonlinear pendulum oscillation/the unforced undamped Duffing oscillator is recovered [5]. (1) has only been analyzed numerically via the midpoint scheme, and based on our comprehensive survey, we did not find any attempt to find a semi-analytical solution to this equation. Motivated by the potential applications of the nonlinear oscillators, the forced damped driven nonlinear (FDDN) pendulum equation or sometimes called the forced damped parametric driven pendulum equation will be studied:(2)φ..+2βφ.+ϕtsinφ=FcosΩt.Also, some analytical approximations to (2) and some related equations will be derived for the first time and will be compared with the Runge–Kutta (RK) numerical solution. Moreover, the Chebyshev collocation numerical method [27–29] is introduced for analyzing both (1) and (2). Furthermore, the MATHEMATICA command Fit is devoted for analyzing the equation of motion. We graphically make a comparison between the analytical and numerical approximations, and the maximum distance error in the whole time domain is estimated.


## 2. Analytical Approximations to the FDDN Pendulum Equation

Let us now write the evolution equation in the form of the initial value problem (i.v.p.):(3)cφ..+2βφ.+ϕtsinφ=FcosΩt,φ0=φ0andφ′0=φ.0.0≤t≤T,where φt=0=φ0 indicates the oscillation amplitude.Using Chebyshev polynomial approximation, we can approximatesinφ as(4)sinφ≈φ−λφ3for−M≤φ≤M,where(5)λ≡λM=16+M5569−M2112+M31883.The errorEM of this approximation may be estimated via the following formula:(6)EM=3298M3−5378M2+1203M−12837for2π180≤M≤π2.For example, at the angleM=30°, the exact error equals Eex=0.0000608542 while the error according to formula (6) equals EM=0.0000455313 and the difference between them is given by E=Eex−EM=0.0000153229. The respective approximation for −30°≤φ≤30° reads as(7)sinφ≈φ−0.164389φ3≈φ−955φ3.Also, forM=75°, we obtain λ=2/13 which will be used as the default value in the present study. Consequently, i.v.p. (3) can be reduced to the following variable coefficient forced damped Duffing i.v.p.(8)cℚ≡φ..+2βφ.+ϕtφ−213φ3−FcosΩt=0,φ0=φ0andφ′0=φ.0.Suppose the solution of this problem is given by(9)cφ=θ+c1cosΩt+c2sinΩt,θ0≡θ0=φ0−c1andφ′0≡θ1=φ.0−c2Ω.The functionθ≡θt is a solution to the following ode:(10)θ..+2βθ.+ϕtθ−213θ3=0.Accordingly, we get(11)ℚ=−126cosΩtA1−126sinΩtA2−12612θ2A3+6θA4+A5ϕt,where the coefficients A1−A5 are given in Appendix 1.Now, for smallγ, we can define ϕt=ω02−εω2cosγt≈ω02−εω2=κ which leads to(12)ℚ≈−126cosΩtB1−126sinΩtB2−κ2612θ2A3+6θA4+A5,where the coefficients A1−A3 have the same values given in Appendix 1 while the values of coefficients of B1 and B2 are given in Appendix 2.The constantsc1 and c2 could be determined from the following system:(13)c3c13+3c22c1−26c1κ−52c2βΩ+26c1Ω2+26F=0,3c23+3c12c2−26c2κ+52c1βΩ+26c2Ω2=0.Eliminatingc2 from system (13), we have the following cubic equation:(14)−26c13F2κ2−3F2κΩ2−416β4Ω4−104β2κ2Ω2+208β2κΩ4−104β2Ω6+9F2κ2c13−624Fβ2Ω2c12+78F3−2704Fβ2Ω2κ+2704Fβ2Ω4=0.Also, by eliminatingc1 from system (13), we get(15)10816β4Ω4+2704β2κ2Ω2−5408β2κΩ4+2704β2Ω6c2+9F2κ2c23+312FβκΩ3−312Fβκ2Ωc22−5408Fβ3Ω3=0.We choose the least in magnitude pair of real rootsc1,c2 in (14) and (15). Accordingly, the final form of the analytical approximation to i.v.p. (3) is given by(16)φapproxt=θt+c1costΩ+c2sintΩ,with(17)θt=e−βt1+b2snftω|m2b1dnftω|msnftω|m+θ0cnftω|m,ft=2−330β2−329κ/κ%329γEtγ2|658εω2330β2−329κ,where(18)ω=−p2m−1,m=121−pp+qφ022+2θ12q,b1=δθ11−2mp,b2=p+qφ02−ω2ω,p=κ,q=−213κ.


## 3. Chebyshev Collocation Numerical Scheme for Analyzing FDDN Pendulum Equation

Now, the Chebyshev interpolation collocation method is introduced for analyzing i.v.p. (3) on the time interval 0,T. To do that, we first solve numerically the ode. Let φ^ be the RK numerical solution to the i.v.p. described by (3). Then, we assume that the solution is given in terms of Chebyshev polynomials:(19)φt=∑k=0nckTk2tT−1,where Tkt stands for the Chebyshev polynomial of the first kind and n denotes the highest degree of the Chebyshev polynomials involved in the linear combination.The collocation pointstk are defined as(20)tk=12T1+cos4k+12n+1π,for k=1,2,3,…,n.The additional equations are required for determining the values ofck. Thus, the values of the coefficients ck are found from the following linear system:(21)φ0=φ0,φ′0=φ.0,φ′′0=φ^′′0.φT=φ^T,φ′T=φ^′T,φ′′T=φ^′′T.φtk=φ^tkfork=5,6,…,n.In general, increasing the value ofn will not guarantee good approximations. Thus, we must choose some optimal value for n to our approximations. To this end, we define a range for possible n values, say(22)7≤nmin≤n≤nmax.We then find the optimal value forn within this range. Let φnt be the solution using formula (19) and let φRKt be the RK numerical solution to i.v.p. (3) on the interval 0≤t≤T. The following maximum distance error with respect to the RK numerical solution is defined:(23)ET,n=max0≤t≤Tφnt−φRKt.The optimal value forn on the range nmin≤n≤nmax will then be that for which the error ET,n is as small as possible. o verifies the validity of this claim. Let us use the following data: β,γ,ω0,ω2,ε,Ω,F,φ0,φ.0=0.2,0.2,1,1,0.2,1,0.1,0,0, as an example in i.v.p. (3). By solving this problem via both RK and Chebyshev collocation numerical methods in the interval 0≤t≤50 and estimating the error ET,n based on relation (23) for 7≤n≤60, we finally get the error ET,n associated with each number n as shown in Table 1. The results in Table 1 illustrate that the optimal value of n based on the mention data for β,γ,ω0,ω2,ε,Ω,F,φ0,φ.0 equals n=45 and the error value corresponding to n=45 equals E50,45=0.000964248. Also, the optimal polynomial according to the mentioned data reads as(24)φ45t=0.0511403t2−0.0177057t3+0.0340205t4−0.0823837t5+0.109124t6−0.0998239t7+0.0677353t8−0.0353909t9+0.0146451t10−0.00490664t11+0.00135445t12−0.000312451t13+0.0000609427t14−0.0000101493t15+1.45518×10−6t16−1.80894×10−7t17+1.96141×10−8t18−1.86462×10−9t19+1.561×10−10t20−1.15515×10−11t21+7.57989×10−13t22−4.42193×10−14t23+2.29822×10−15t24−1.06583×10−16t25+4.41541×10−18t26−1.63491×10−19t27+5.41107×10−21t28−1.6×10−22t29+4.22221×10−24t30−9.92634×10−26t31+2.07398×10−27t32−3.83854×10−29t33+6.26651×10−31t34−8.97497×10−33t35+1.11996×10−34t36−1.2071×10−36t37+1.1112×10−38t38−8.60994×10−41t39+5.50647×10−43t40−2.82885×10−45t41+1.12172×10−47t42−3.22145×10−50t43+5.96106×10−53t44−5.33464×10−56t45.Polynomial (24) allows us to estimate the cuts with the horizontal axis as well as the maxima and minima to the crest and the trough, respectively, as shown in Figure 1 and Table 2. Using the MATHEMATICA command Fit gives the following solution φMatht≡φMathematicat for β,γ,ω0,ω2,ε,Ω,F,φ0,φ.0=0.2,0.2,1,1,0.2,1,0.1,0,0:(25)φMatht=e−t/5t6365−2t523+125t494−229t332−21045t2619+33029t74−6344576sint+e−t/57873214642193096t6−17115684t5−426465793t4+10891280252t3−83720431296t2+147163503646t+393199198728cost+e−t/51544400−3575t7+93600t6−1480050t5+12725856t4−38075400t3−118006200t2+1213821180t+294883875−1173617,with error E=0.000213725. Figure 2 demonstrates the comparison between the approximations of i.v.p. (3) using RK numerical solution with MATHEMATICA command Fit (here, polynomial (25)) and Chebyshev collocation numerical solution (24). It is noted that all used techniques give highly accurate approximations with low errors as compared to the RK numerical solutions.Table 1NET;n712:731888:1870995:34902103:73001112:47806125:991511315:64721419:17891511:42471612:97251744:5793NET;n1856:02771922:74082045:86662186:01022241:44782354:48:002481:8262259:29082650:97742735:2609287:88678NET;n2924:00963011:86333111:92673213:06993316:415534160:44:00351021:51:0036455:30:0037293:27:003862:37:0039594:34:00NET;n401011:59:00411016:07:0042193:44:0043978:08:0044154:12:0045779:04:0046584:21:0047655:31:00481042:34:004904:05:00501060:25:00NET;n51135:10:0052626:08:0053401:51:0054666:50:0055999:51:0056893:24:0057622:23:0058199:58:00597:19746030:729Figure 1
The Chebyshev collocation and RK numerical approximations to i.v.p. (3) for β,γ,ω0,ω2,ε,Ω,φ0,φ˙0=0.2,0.2,1,1,0.2,1,0.1,0,0 and with n=45 is plotted in the φ,t plane and the cuts with the horizontal axis as well as the maxima and minima to the crest and the trough is determined.Table 2Zeros of the polynomial solution '45(t)Zeros of the derivative '045(t)KTkktkktkktk13:19931134:960911:999781126:674126:395941238:032824:930361230:43:0039:462881341:149337:982591330:43:00412:45831443:9479410:98441433:3364515:45621546:912513:96411536:4941618:51341649:9492613:96451636:4941721:6812716:981739:5245824:8139820:09031842:4811928:3547923:32911945:42261031:71651026:6742048:4281Figure 2
The comparison between the approximations of i.v.p. (3) using RK numerical solution with MATHEMATICA command Fit (here, polynomial (25)) and Chebyshev collocation numerical solution (24) for β,γ,ω0,ω2,ε,Ω,φ0,φ˙0=0.2,0.2,1,1,0.2,1,0.1,0,0 and n=45.
(a)(b)The semi-analytical solution (16) to i.v.p. (3) could be recovered as follows.Case (1). For β,ω2,F=0,0,0.1, the different approximations to the forced undamped Duffing oscillator with constant coefficients are introduced in Figure 3 with γ,ω0,ε,Ω,φ0,φ.0=0.1,1,0.1,2,0,0. The comparison between the RK method and the analytical approximation (15) is presented in Figure 3(a). The approximate solutions using the MATHEMATICA command Fit and RK method are displayed in Figure 3(b). In Figure 3(c), both RK and Chebyshev collocation numerical approximations are presented. Also, the maximum error for the analytical approximation (15) and MATHEMATICA command Fit and Chebyshev collocation numerical solutions as compared to the RK numerical approximation is estimated based on the following relation:(26)E∞|Type−solution=max0≤t≤30φType−solution−φRK.Accordingly, the maximum error of the three approximations for the present case is estimated as(27)E|semi−analy=0.000214826,E|Mathematica=0.00856153,En=47|Chebyshev=0.0000389093.Case (2). For β,ω2,F=0,1,0.1, the comparison between the analytical approximation (16) and the numerical approximations using the RK, MATHEMATICA command Fit, and Chebyshev collocation numerical methods to the forced undamped Duffing equation with variable coefficients is considered as shown in Figure 4 with γ,ω0,ε,Ω,φ0,φ.0=0.1,1,0.1,2,0,0. The maximum error of the three approximations to the present case is calculated as(28)E|semi−analy=0.00344896,E|Mathematica=0.0100608,En=47|Chebyshev=0.0000237867.Case (3). For β,ω2,F=0.1,1,0, the unforced damped Duffing equation with variable coefficients is recovered and its semi-analytical solution (16) is compared with the numerical approximations using RK, MATHEMATICA command Fit, and Chebyshev collocation numerical methods as demonstrated in Figure 5 with γ,ω0,ε,Ω,φ0,φ.0=0.1,1,0.1,2,0,0.1. In addition, the maximum error to the three approximations as compared to RK numerical approximation is estimated as(29)E|semi−analy=0.000447172,E|Mathematica=0.00133459,En=47|Chebyshev=2.05447×10−6.The MATHEMATICA code for the RK and Chebyshev collocation numerical approximations with the maximum error is given in Appendix 3.Case (4). For β,ω2,F=0.1,1,0.1, the general analytical approximation (16) to the forced damped parametric driven pendulum i.v.p. (3) is compared with the RK, MATHEMATICA command Fit, and Chebyshev collocation numerical solutions as elucidated in Figure 6. The influence of β on the amplitude of the semi-analytical solution (16) is investigated as shown in Figure 7. It is noted that the amplitude of the analytical approximation (15) decreases with the increase of β. Also, the maximum error according to relation (23) to the analytical approximation (16) and MATHEMATICA command Fit, and Chebyshev collocation numerical solutions as compared to the RK numerical solution is estimated as follows:(30)E|semi−analy=0.00326255,E|Mathematica=0.0128206,En=47|Chebyshev=0.000155082.Figure 3
The comparison between the semi-analytical solution (analytical approximation) (16) and the numerical approximations using the Chebyshev collocation method and RK numerical method as well as the MATHEMATICA command Fit to i.v.p. (3) for case (1):β,ω2,F=0,0,0.1.
(a)(b)(c)Figure 4
The comparison between the analytical approximation (16) and the numerical approximations using the Chebyshev collocation method and RK numerical method as well as the MATHEMATICA command Fit to i.v.p. (3) for case (2):β,ω2,F=0,1,0.1.
(a)(b)(c)Figure 5
The comparison between the analytical approximation (16) and the numerical approximations using the Chebyshev collocation method and RK numerical method as well as the MATHEMATICA command Fit to i.v.p. (3) for case (3):β,ω2,F=0.1,1,0.
(a)(b)(c)Figure 6
The comparison between the analytical approximation (16) and the numerical approximations using the Chebyshev collocation method and RK numerical method as well as the MATHEMATICA command Fit to i.v.p. (3) for case (4):β,ω2,F=0.1,1,0.1.
(a)(b)(c)Figure 7
The analytical approximation (16) to i.v.p. (3) plotted in the φ,tplane for different values to the coefficient of the damping term β.In all mentioned cases, it is observed that analytical approximation (16) to i.v.p. (3) and its related equations (here we mean the four mentioned cases) give highly accurate results as compared to the numerical approximations. It is observed that analytical approximation (16) is better than the MATHEMATICA approximation but less than the Chebyshev collocation numerical solution for all mentioned cases. However, in the fourth case, i.e., the forced damped parametric driven pendulum i.v.p. (3), the MATHEMATICA approximation is better than the analytical approximation (16). In general, all obtained approximations are characterized by their high accuracy. However, semi-analytical solution (16) is more stable than the Chebyshev collocation numerical solution against all relevant physical variables.


## 4. Conclusions

In this work, some effective and accurate analytical and numerical approximations to the forced damped parametric driven pendulum equation have been derived and investigated. The mentioned equation of motion has been reduced to the forced damped Duffing equation with variable coefficients in order to find its analytical solution. In terms of the Jacobi elliptic functions, the analytical approximation has been derived. For the numerical approximations, the Chebyshev collocation method has been used for analyzing the equation of motion and some related equations. It was noted that the analytical approximation could recover some special cases to the nonlinear pendulum oscillators. For instance, for undamping case, i.e., forβ=0, the solution to the forced undamped Duffing equation with variable coefficients has been recovered and examined. Also, for β,ω2=0,0, the solution to the forced undamped Duffing equation with constant coefficients has been recovered and discussed. The obtained approximations were compared with the RK numerical approximation and the MATHEMATICA command Fit approximation. Also, the maximum distance error has been estimated for all approximations as compared to the RK numerical approximation. It was found that the analytical approximation gives good results with high accuracy as compared to the numerical approximations. Furthermore, it was observed that the analytical approximation is better than the MATHEMATICA approximation but less than the Chebyshev collocation numerical solution for all mentioned cases except the case of the forced damped parametric driven pendulum i.v.p. (3), the MATHEMATICA approximation is slowly better than the analytical approximation (16). The methods used in this study could be extended to solve many nonlinear equations that control the different cases of pendulum oscillations [30–33]. In addition, the obtained results/solutions are useful for investigating several physical problems related to the oscillations in plasma physics, fluid mechanics, field theory, engineering science, solid state physics, and quantum mechanics.



---

*Source: 5454685-2022-06-22.xml*