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.
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 , with being the only odd polynomial where . 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 , forced Duffing oscillator , forced damped Duffing oscillator , 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 (here, is only odd polynomial where and is only even polynomial where ) and some related oscillators such as (un)damped HD oscillator , forced HD oscillator , forced damped HD oscillator 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]:where represents the coefficient of the damping term and 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–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.):
Based on Chebyshev polynomial approximation, the value of can be expanded aswhere
Other possible choices for can be considered as
For , the following approximation is obtained:
Next, we replace the original i.v.p. (3) by the following approximate i.v.p.
Assume that the solutions to the i.v.p. (8) have the following formulas:with the initial conditions (ICs)where the number and the function are chosen in order to get the smallest possible residual error :
Now by substituting ansatz (9) into the i.v.p. (8), we getwithwhere .
For , we getwith
Integrating equation (14) leads towith
Inserting the value of given in equation (16) into ansatz (9) and applying the ICs , the value of can be determined from the following quartic equation:
Solution (9) is presented in Figures 1(a) and 1(b) for and at . 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:
(a)
(b)
The maximum distance error according to relation (19) for and at is, respectively, estimated as
It is noted that the accuracy of solution (9) becomes good and acceptable for small , but for large value of , the accuracy of solution (9) reduces as shown in Figure 1.
Also, the maximum residual error is defined as
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.
Suppose that the physical problem described by (22) involves some small parameters , ,…., . Let be the solution to the i.v.p. (22) and assume that
The solution depends not only on but also on the parameters , , so that we can rewrite as
Let us multiply each parameter by some other parameter and consider the following parametric solution:
Accordingly, the function may be written in a power series as follows:where depends on only, say .
Based on Krylov–Bogoliubov–Mitropolsky method (KBMM), the solution of equation (22) is assumed to bewhere each is a periodic function of and and are assumed to vary with time according towhere and .
Moreover, the hybrid homotopy KBMM (HKBMM) is suggested to be
The next step is to write the residual as a power series in :
For the determination of the unknown functions , , , and , we equate to zero the coefficients 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 nor . For (the first approximation), we may use the following formulas (we neglected all terms containing for ):andwith
The approximate analytical solution is obtained by putting . However, we may keep the parameter and then we may use it as a residual minimization parameter. The optimal value to will be near .
Now, the proposed method can be applied for investigating the i.v.p. (8):where in our case, and
Observe that when and , we get .
In equations (27)–(29), for and , we have
The homotopy to equation (35) is written as
The substitution of equation (37) into equation (38) leads to
We must have
The coefficients of and must be vanished to eliminate the secularity. Accordingly, we have
Thus, equation (40) reduces to
Solving equation (42), the following particular solution without any secularity terms is obtained:
From equations (37) and (41), the functions and can be determined:
By solving system (44), we getand
We finally get the analytical approximation in its first approximation .
We now introduce three optimization parameters by replacing and then the approximation (47) can be modified to bewithand
The numbers , , and are free parameters that we choose in order to minimize the residual error
The default parameter values are . The constants and are determined from the initial conditions (ICs) and .
Following the same procedure above, we can get some higher-order approximations. For example, for , the following solution is obtained:where the coefficients are defined in Appendix. The values of associated to this solution can be determined from the following equations:and
The approximate solution of the i.v.p. (35) using the HKBMM is introduced in Figures 2(a) and 2(b) for and . Also, the maximum distance error is estimated for the two cases as follows:
(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 . Also, this approximation is better than the trigonometric solution (9) 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:
Choose some positive integer and divide the interval into -subintervals by means of the knots , where . Then, the first and second-order derivatives can be approximated as follows:
Consequently, the following discrete version to ODE (56) for is obtained:
The values of , , , and 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 without loss of generality):with
The numerical approximation using FDM is plotted in Figures 3 and 4 for different values of . In Figures 3(a) and 3(b), the FDM numerical approximation is plotted against and , respectively. Moreover, the effect of on the numerical approximation is illustrated in Figures 4(a) and 4(b) for and , respectively. Furthermore, the maximum distance error is calculated for all mentioned cases as follows:
(a)
(b)
(a)
(b)
It is clear that the accuracy of the FDM numerical approximation increases with increasing . Also, this approximation is stable against the long time intervals and arbitrary angle . Moreover, this approximation is better than the trigonometric solution (9) .
3.2. Numerical Approximation via Galerkin Hats Method
First, let us consider a polynomial second-order forced damped i.v.p.where .
Let us consider the i.v.p.
Suppose that is the solution to the i.v.p.
Define
The constants and are determined from the system
Thus, problemproblem (61) reduces to (60) by problem (63) reduces to (62).
Some particular cases to the i.v.p. (62) are defined as
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 are the so-called linear Galerkin hats.
Let us investigate the present problem in the interval and by choosing some positive integer and define the step and let for . The functions for are defined as
For an illustration, see Figure 1.
Figure 1 Galerkin hats for and .
Some properties of these functions can be illustrated as follows:for and , andfor and .
In general, for and , we have
Using the formulafor any and , and assuming that 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
The minimization procedure yields the valueswith
Then, i.v.p. (76) can be reduced to the following approximate i.v.p.
Assume that the solution to i.v.p. (81) is given by
The Galerkin method solves the following algebraic system of nonlinear equations:withfor , where and . System (83) can be solved recursively.
For , the value of can be determined from the following quintic equation:
We choose the real root to equation (86) that is closest to . Suppose we already found the values , ,…, for some . Then, the value of is found by solving the quintic in equation (82) by letting . However, since , we must solve two quintics in . We choose the real root to the two quintics that is closest to . Let us introduce the quintic
Any of the two quintics (82) for . Then, for sufficiently large , the value of may be estimated by means of the formulawhere .
For arbitrary ICs, i.e., for any values to , the following ansatz is assumed:
Then, we getand
Using the approximationgives
The following definitions are used in our analysis:
The algebraic system to be solved for reads as