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

In this paper, some novel analytical and numerical techniques are introduced for solving and analyzing nonlinear second-order ordinary diﬀerential equations (ODEs) that are associated to some strongly nonlinear oscillators such as a quadratically damped pendulum equation. Two diﬀerent analytical approximations are obtained: for the ﬁrst 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 ﬁrst time for analyzing the evolution equation. For the numerical approximation, both the ﬁnite diﬀerence 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.


Introduction
Duffing-type equation is one of the most important seconddegree differential equations that is used to describe many different phenomena [1][2][3][4][5][6]. e 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) + f(x) � 0, with f(x) � ∞ i�1 α i x i 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 € , and many other oscillators with odd polynomials and complicated damping term [7][8][9][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) + f(x) + g(x) � 0 (here, f(x) � ∞ i�1 α i x i is only odd polynomial where i � 1, 3, 5, . . . and g(x) � ∞ i�2 c i x i is only even polynomial where i � 2, 4, 6, . . .) and some related oscillators such as (un)damped HD oscillator € x(t) + β _ x(t) + f(x) + g(x) � 0, forced HD oscillator € x(t) + f(x) + g(x) � F(t), forced damped HD oscillator € x(t) + β _ x(t) + f(x) + g(x) � F(t) and many other HD oscillators with complicated damping term and complicated polynomials [11][12][13][14][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][18][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.
is equation is called the quadratically damped pendulum equation [20]: 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: Numerous oscillators with quadrature damping have been investigated over a wide range of different fields [21][22][23][24][25][26].
ere are many methods for solving nonlinear differential equations.
ere are many analytical and numerical methods that dealt with solving different differential equations, and some of these methods can be found in Refs. [28][29][30][31][32][33][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).

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.

First Approach: Ansatz Method and Trigonometric
Solution. Let us rewrite evolution equation (1) in the form of the initial value problem (i.v.p.): Based on Chebyshev polynomial approximation, the value of sin θ can be expanded as where Other possible choices for (r 0 , r 1 , λ) can be considered as For (ω 0 , θ . 0 ) � (1, 0), the following approximation is obtained: Next, we replace the original i.v.p. (3) by the following approximate i.v.p.
Also, the maximum residual error is defined as is is another form for the error to check the accuracy of the obtained approximations.

Second Approach: HKBMM.
Let us consider the i.v.p.
Suppose that the physical problem described by (22) involves some small parameters ε 1 , ε 2 ,. . .., ε r . Let x ≡ x(t) be the solution to the i.v.p. (22) and assume that e solution x ≡ x(t) depends not only on t but also on the parameters ε 1 , ε 2 , . . . , ε r , so that we can rewrite x ≡ x(t) as Let us multiply each parameter by some other parameter p and consider the following p−parametric solution: x p � x p t; pε 1 , pε 2 , . . . , pε r .
Accordingly, the function x p may be written in a power series as follows: Complexity where u k depends on t only, say u k � u k (t).
Based on Krylov-Bogoliubov-Mitropolsky method (KBMM), the solution of equation (22) is assumed to be where each u n is a periodic function of ψ and a and ψ are assumed to vary with time according to where a ≡ a(t) and ψ ≡ ψ(t). Moreover, the hybrid homotopy KBMM (HKBMM) is suggested to be e next step is to write the residual H p (x, t) as a power series in p: For the determination of the unknown functions u n , ψ n , A n , and a, we equate to zero the coefficients Υ n in equation (22) and then we can get a system of ODEs. To avoid the socalled 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 p j for j ≥ 2): and x _ x � p aω 0 cos(ψ)u 1,ψ + a 2 ψ 1 sin(ψ)(−cos(ψ))

Complexity
with e approximate analytical solution is obtained by putting p � 1. However, we may keep the parameter p and then we may use it as a residual minimization parameter. e optimal value to p will be near p � 1. Now, the proposed method can be applied for investigating the i.v.p. (8): where in our case, x � θ and Observe that when ε 1 � ε ⟶ 0 and ε 2 � λ ⟶ 0, we get F � 0.
In equations (27)- (29), for N � 1 and λ � 2/13, we have e homotopy to equation (35) is written as e substitution of equation (37) into equation (38) leads to We must have e coefficients of cos(ψ) and sin(ψ) must be vanished to eliminate the secularity. Accordingly, we have us, equation (40) reduces to Solving equation (42), the following particular solution without any secularity terms is obtained: From equations (37) and (41), the functions a ≡ a(t) and ψ ≡ ψ(t) can be determined: By solving system (44), we get and We finally get the analytical approximation in its first approximation (p � 1).
Following the same procedure above, we can get some higher-order approximations. For example, for N � 3, the following solution is obtained: where the coefficients S 1,2,3 are defined in Appendix. e values of (a, ψ) associated to this solution can be determined from the following equations: and e 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:

Complexity
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.

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.

Numerical Approximation via FDM.
First, let us discuss and apply the FDM on the general second-order ODE. us, the following general second-order is introduced: Choose some positive integer n ≥ 6 and divide the interval [0, T] into n-subintervals by means of the knots t i � ih, where h � T/n/(i � 0, 1, 2, . . . , n). en, the first and second-order derivatives can be approximated as follows: Consequently, the following discrete version to ODE (56) for i � 5, 6, . . . is obtained: e values of x 1 , x 2 , x 3 , and x 4 are obtained from some numerical or approximate analytical solution to the i.v.p. (56). System (58) may be solved recursively. e above algorithm can be applied for analyzing the i.v.p. (note here θ(t) ≡ x(t) without loss of generality): e numerical approximation using FDM is plotted in Figures 3 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) It is clear that the accuracy of the FDM numerical approximation increases with increasing n. 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 .

Numerical Approximation via Galerkin Hats Method.
First, let us consider a polynomial second-order forced damped i.v.p.
e constants c 0 and c 1 are determined from the system us, problemproblem (61) reduces to (60) by problem (63) reduces to (62).
Some particular cases to the i.v.p. (62) are defined as (67) We will use the same idea for the linear case.
Here, we start to discuss the linear case in system (67): An approximate solution to the i.v.p. (68) is assumed to be in the following ansatz form: where the functions φ k (t) are the so-called linear Galerkin hats. Let us investigate the present problem in the interval 0 ≤ 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, · · ·. e functions φ k (t) for k � 1, 2, . . . , n are defined as For an illustration, see Figure 1. Figure 1 Galerkin hats for n � 7 and T � 5. Some properties of these functions can be illustrated as follows: for |j − k| ≥ 2 and t ∈ [0, T], and for j ≥ 1 and p � 1, 2, 3, . . .. In general, for |j − k| � 1 and r, s � 0, 1, 2, 3, · · ·, we have Using the formula for any N ≥ 0 and c 0 � c n+1 � 0, and assuming that a j (t) ≡ a j � const, we may evaluate easily the following integration: Now, let us return to the original i.v.p. (3): Using the Chebyshev approximation, Another approximation may be obtained by minimizing the integral π 2 −π 2 e minimization procedure yields the values min � π 2 − 3360 11975040 − 2661120π 2 + 171720π 4 − 2664π 6 + 13π 8 Assume that the solution to i.v.p. (81) is given by e Galerkin method solves the following algebraic system of nonlinear equations: with for j � 2, 3, . . ., where c 0 � 0 and c 1 � T/nθ . 0 . System (83) can be solved recursively.
For j � 2, the value of c 2 can be determined from the following quintic equation: We choose the real root to equation (86) that is closest to c 1 . Suppose we already found the values c 2 , c 3 ,. . ., c k−1 for some k ≥ 2. en, the value of c k is found by solving the quintic in equation (82) by letting j � k. However, since ���������� (c k − c k− 1 ) 2 � ±(c k − c k−1 ), we must solve two quintics in c k . We choose the real root to the two quintics that is closest to c k−1 . Let us introduce the quintic Any of the two quintics (82) for j � k. en, for sufficiently large n, the value of c k may be estimated by means of the formula where z 0 � c k−1 . For arbitrary ICs, i.e., for any values to (θ 0 , θ . 0 ), the following ansatz is assumed: and Using the approximation e following definitions are used in our analysis: e algebraic system to be solved for j � 2, 3, . . . , n reads as Complexity 11 e 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 , θ It is clear that increasing the number of hats n 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.

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. e 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. ese 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. us, this method can be devoted for studying and investigating different pendulum oscillators for any nonlinearity [16,27].