The Vibrational Motion of a Dynamical System Using Homotopy Perturbation Technique

T. S. Amer¹, A. A. Galal², Shimaa Elnaggar^2,3

¹Mathematics Department, Faculty of Science, Tanta University, Tanta, Egypt.
²Department of Physics and Engineering Mathematics, Faculty of Engineering, Tanta University, Tanta, Egypt.
³Kafrelsheikh Higher Institute of Engineering and Technology, Kafrelsheikh, Egypt.

DOI: 10.4236/am.2020.1111073                

Abstract

This paper outlines the vibrational motion of a nonlinear system with a spring of linear stiffness. Homotopy perturbation technique (HPT) is used to obtain the asymptotic solution of the governing equation of motion. The numerical solution of this equation is obtained using the fourth order Runge-Kutta method (RKM). The comparison between both solutions reveals high consistency between them which confirms that, the accuracy of the obtained solution using aforementioned perturbation technique. The time history of the attained solution is represented through some plots to reveal the good effect of the different parameters of the considered system on the motion at any instant. The conditions of the stability of the attained solution are presented and discussed.

Keywords

Homotopy Technique, Nonlinear Vibration, Lagrange’s Equation, Stability

1. Introduction

Many problems related to mathematicians, physicists, biologists, chemists and engineers are formulated in differential equations whether linear or nonlinear. The solutions of a linear one can be obtained easily using some of well-established methods on the contrary with nonlinear differential Equations (NDE) that we often refuge to approximate solutions. Nonlinear oscillations had shed the interest of many scientists due to that most of the problems dealing with vibrations are nonlinear, see [1] [2].

Since it is difficult to find the exact solutions of such equations, many researchers have turned their attention to obtain the approximate solutions of these problems using perturbation techniques [3] [4] [5] through inserting a small parameter like Krylov-Bogolubov-Mitropolski technique (KBMT), Lindstedt Poincar? technique (LPT), the multiple scales technique (MST), averaging technique (AT) and others. These solutions are given by expansions in powers of the small parameter. On the other side, if this parameter can’t be introduced into the solution procedure of these equations, the researchers have faced a difficulty to obtain the solutions of the classical perturbation problems. Therefore, a great difficulty in solving these equations will be arisen by using the traditional perturbation methods which are not convenient for solving the strongly vibrations problems.

Over the past two decades, many mathematicians and physicists have done their great efforts to find new mathematical tools to deal with the dynamical systems that mathematically described by nonlinear differential equations [5] [6] [7]. In [5], the authors have overcome this difficulty through establishing effective ways to deal with the nonlinearity that arises in these problems. A combination between homotopy perturbation technique (HPT) and Frobenius approach was stated in [6] to get the exact solution of NDE. However, for non-linear conservative systems, the generalization of some perturbation methods has overcome this restriction. Moreover, the generalization of the HPT of He’s [7] and LPT has fructified the desired results for strongly nonlinear vibrations [8] - [14]. The analytic and numerical solutions of three different problems are obtained in [15] using HPT and RKM respectively. Another analytic method called Energy Balance Method (EBM) is used to compare the attained results with each other.

A max-min method is presented in [9] to get the approximate solutions of NDE and it is applied to some examples. In [10], the authors have applied a modified HPT to obtain more accurate approximate solution of a coupled two strongly NDE. Another version of HPT is applied in [12] for excited nonlinear problems and the Mathieu equation is studied as the simplest example. The same author suggested a modified version of HPT by absorbing MST in [13] for nonlinear oscillators systems. Recently, he outlines his work on asymptotic solution for delayed self-feedback of a nonlinear dynamical problem in [14] and the stability configuration is presented.

In [16], the author compared the results obtained by HPT with the results obtained by homotopy analysis technique to reveal that HPT results are more accurate than the second one. HPT was applied in [17] to obtain the asymptotic solutions of the non-linear equation of Fredholm integral of second kind. In [18], the authors suggested a modification of HPT which is considered an adjustment tool for obtaining the periodic solutions of nonlinear oscillatory systems. Thereby, the difficult problem can be transformed into an easier one to be handled. A new technique is introduced in [19] through a combination between HPT and modification of LPT in order to get the solutions of certain non-smooth oscillators while the work in [20] can be considered as a standard HPT in which it can be used as an application of several nonlinear differential equations.

In [21], the authors presented a modification of HPT to get the solution of a dynamical model consists of the motion of a rod in a circular surface without slipping. The obtained results are in good agreement with the numerical ones. This problem was treated in [22] using the modified harmonic balance technique [23]. A combination of MPT and Laplace transformation to achieve the asymptotic solution of the governing equation of motion of the same problem is studied in [24], in which the stability of the obtained solutions is examined. The approximate solutions of some tested vibrating systems are obtained in [25] using a modification of HPT and Amplitude frequency formulation (AFF).

HPT is used in [26] to obtain the periodic solution of the fractional sine-Gordon equation beside the Riemann-Liouville fractional derivative. The authors obtained a relationship between the frequency and amplitude, and the impact of the order of fractional derivative on the vibration property is investigated. In [27], the authors investigated the periodic solution of the nonlinear Duffing oscillator with fractional order utilizing a modification of HPT which is the insertion of an auxiliary parameter and using two homotopy parameters. A nonlinear packaging system has been solved analytically using HPT of Li-He’s in [28], in which the energy method is utilized to progress the frequency and the maximal displacement of the system.

In this paper, the solution of a nonlinear oscillating dynamical system is investigated. This system consists of a mass m₁ connected with a spring of linear stiffness and with other mass m₂ through a massless string of length l. HPT is utilized to obtain the solution of the equation of motion. This solution is graphically represented for different values of the system parameters and compared with the numerical solution of the governing equation of motion using the Runge-Kutta method [29] from fourth order. This comparison reveals high consistency between them which emphasizes the accuracy of the results obtained by HPT. The stability of the investigated model is presented and analyzed.

This paper is designated as follows. In Section 2, a description of the investigated problem and the derivation of the equation of motion are presented. Section 3 sheds light on the basic idea of HPT. Section 4 is devoted to reduce the equation of motion into appropriate equation and to obtain the solution of this equation analytically using HPT. In Section 5, we are going to represent the attained solution graphically and to obtain the numerical solution using the Runge-Kutta method. The stability of the obtained solution is discussed in Section 6. Finally, the manuscript is finished with some concluding remarks.

2. Description of the Problem

In this section, we are going to obtain the governing equation of motion of a nonlinear oscillation system using HPT of a dynamical model. This model consists of two masses; a first one m₁ moves horizontally in which it is attached to a spring of linear stiffness k and connected with the second mass m₂ with a massless string of length l see (Figure 1). Therefore let us consider that S₀, x and y are the natural length of the given spring, the horizontal coordinate of the centroid

Figure 1. The dynamical model.

of m₁ and the vertical coordinate of the centroid of m₂, respectively. Therefore, the potential and kinetic energies V and T of the system can be written in the forms

$V=-m_{2}g\sqrt{l^2-x^2}+\frac{k{x}^2}{2},T=\frac{1}{2}\left(m_1{\stackrel{\dot{}}{x}}^2+m_2{\stackrel{\dot{}}{y}}^2\right),$ (1)

where g is the gravitational acceleration, x is the extension of string after time t, dots denote to the differentiation with respect to time and $\left(\stackrel{\dot{}}{x}\mathrm{,}\stackrel{\dot{}}{y}\right)$ is the Cartesian velocity of the point $\left(x\mathrm{,}y\right)$.

According to (1), the Lagrangian $L=T-V$ has the form

$L=\frac{1}{2}\left(m_1+\frac{m_2{x}^2}{l^2-x^2}\right){\stackrel{\dot{}}{x}}^2+m_{2}g\sqrt{l^2-x^2}-\frac{k{x}^2}{2}.$ (2)

An inspection of the Lagrange’s function (2) shows that the investigated system has only one degree of freedom. Therefore, Lagrange’s equation for conservative system may be written as

$\frac{\text{d}}{\text{d}t}\left(\frac{\partial L}{\partial \stackrel{\dot{}}{x}}\right)-\frac{\partial L}{\partial x}=0.$ (3)

Here, x and $\stackrel{\dot{}}{x}$ are the generalized coordinate and velocity of the system respectively. Making use of (2) and (3) yields to the following form of the governing equation of motion

$\left(m_1+\frac{m_2{x}^2}{l^2-x^2}\right)\ddot{x}+\frac{m_2{l}^{2}x}{{\left(l^2-x^2\right)}^2}{\stackrel{\dot{}}{x}}^2+m_{2}g\frac{x}{{\left(l^2-x^2\right)}^{\frac{1}{2}}}+kx=0$ (4)

3. Homotopy Perturbation Technique

This section is devoted to illustrate HPT [7] through solving the following general nonlinear differential equation

$K\left(u\right)-f\left(r\right)=\mathrm{0,}r\in \Omega \mathrm{,}$ (5)

beside the following boundary condition

$B\left(u\mathrm{,}\frac{\partial u}{\partial n}\right)=\mathrm{0,}r\in \Gamma \mathrm{.}$ (6)

Here K and B represent the general differential operator and the boundary operator respectively, $f\left(r\right)$ denotes a known analytical function, $\Gamma$ is the boundary of a domain $\Omega$ and $\frac{\partial u}{\partial n}$ refers to differential along the normal drawn outwards from $\Omega$.

An inspection of Equation (5), broadly speaking, the operator K can be separated into two parts; which are a linear part L and a nonlinear one N. Therefore Equation (5) can be rewritten in the form

$L\left(u\right)+N\left(u\right)-f\left(r\right)=0$ (7)

It is worthwhile to notice that according to HPT, we can construct the homotopy $v\left(r\mathrm{,}\rho \right)\mathrm{<mo>\; :\; </mo>}\Omega \times \left[\mathrm{0,1}\right]\to R$, which satisfies

$H\left(v,\rho \right)=\left(1-\rho \right)\left(L\left(v\right)-L\left(U\right)\right)+\rho \left(K\left(v\right)-f\left(r\right)\right)=0,\rho \in \left[0,1\right]$ (8)

or in an equivalent form as

$H\left(v,\rho \right)=L\left(v\right)-L\left(U\right)+\rho L\left(U\right)+\rho \left(N\left(v\right)-f\left(r\right)\right)=\mathrm{0,}\rho \in \left[0,1\right]$ (9)

where $\rho \in \left[0,1\right]$ is a homotopy parameter and U (initial guess) is an initial approximation of Equation (5), in which it satisfies the boundary conditions.

In order to investigate the solution of (8) or (9), we express about this solution as a power series of $\rho$ as

$v=v_0+\rho v_1+{\rho }^2{v}_2+\cdots \mathrm{.}$ (10)

At $\rho \to 1$, Equations (8) or (9) corresponds to Equation (5) and the results in the approximation to the solution of Equation (5) can be expressed as

$u=\underset{\rho \to 1}{\lim }v=v_0+v_1+v_2+\cdots .$ (11)

It is important to note that, series (11) is convergent for more cases. Some criteria are suggested for convergence of this series, see [7].

4. Method of Solution

Dividing both sides of (4) by m₁ and consider that

${\omega }_0^2=\frac{k}{m_1}+\frac{Rg}{l},R=\frac{m_2}{m_1},u=\frac{x}{l},\left|u\right|\ll 1$ (12)

to reduce the equation of motion (4) to a more appropriate as

$\left[1+\frac{m_2{x}^2}{m_1\left(l^2-x^2\right)}\right]\ddot{x}+\left[\frac{m_2{l}^{2}x}{m_1{\left(l^2-x^2\right)}^2}\right]{\stackrel{\dot{}}{x}}^2+\frac{g{m}_{2}x}{m_1{\left(l^2-x^2\right)}^{\frac{1}{2}}}+\frac{kx}{m_1}=0.$

On the use of (12), the previous equation can be rewritten in the form

$\left[1+\frac{R{u}^2}{1-u^2}\right]\ddot{u}+\left[\frac{Ru}{{\left(1-u^2\right)}^2}\right]{\stackrel{\dot{}}{u}}^2+\frac{Rgu}{l{\left(1-u^2\right)}^{\frac{1}{2}}}+\frac{ku}{m_1}=0.$

Expanding the previous equation to obtain

$\begin{array}{l}\left[1+R{u}^2\left(1+u^2+\cdots \right)\right]\ddot{u}+Ru\left(1+2u^2+\cdots \right){\stackrel{\dot{}}{u}}^2\\ +\frac{Rg}{l}u\left(1+\frac{1}{2}{u}^2+\cdots \right)+\frac{ku}{m_1}=\mathrm{0\; <mo>\; ;\; </mo>}\left(u\ll 1\right)\mathrm{.}\end{array}$

Therefore, we obtain the following equation

$\left(1+R{u}^2\right)\ddot{u}+Ru{\stackrel{\dot{}}{u}}^2+{\omega }_0^{2}u+\frac{Rg}{2l}{u}^3+\cdots =0.$ (13)

A closer look of this equation reveals that it is a second order differential equation with high nonlinearity.

The aim of this section is to obtain the approximate solution of the governing equation of motion utilizing HPT in the presence of the following initial conditions

$u\left(0\right)=A,\stackrel{\dot{}}{u}\left(0\right)=0.$ (14)

By virtue of Equations (13) and (7), the linear part $L\left(u\right)$ and nonlinear one N(u) have the forms

$L\left(u\right)=\ddot{u}+{\omega }_0^{2}u\mathrm{,}$ (15)

$N\left(u\right)=R{u}^2\ddot{u}+Ru{\stackrel{\dot{}}{u}}^2+\frac{Rg}{2l}{u}^3\mathrm{,}$ (16)

where

$f\left(r\right)=0.$ (17)

Equation (8) can be rewritten in the form

$H\left(v,\rho \right)=\left(1-\rho \right)\left(L\left(v\right)-L\left(U\right)\right)+\rho \left(L\left(v\right)+N\left(v\right)-f\left(r\right)\right)=0,\rho \in \left[0,1\right]$ (18)

Substituting (15)-(17) into (18) to obtain

$\begin{array}{c}H\left(v\mathrm{,}\rho \right)=\left(1-\rho \right)\left(\ddot{v}+{\omega }_0^{2}v-\ddot{U}-{\omega }_0^{2}U\right)\\ +\rho \left(\ddot{v}+{\omega }_0^{2}v+R{v}^2\ddot{v}+Rv{\stackrel{\dot{}}{v}}^2+\frac{Rg}{2l}{v}^3\right)\\ =0.\end{array}$ (19)

Let $U=0$ ; the previous equation have the form

$H\left(v,\rho \right)=\left(1-\rho \right)\left(\ddot{v}+{\omega }_0^{2}v\right)+\rho \left(\ddot{v}+{\omega }_0^{2}v+R{v}^2\ddot{v}+Rv{\stackrel{\dot{}}{v}}^2+\frac{Rg}{2l}{v}^3\right)=0.$ (20)

Making use of (10) and (20), then equating the coefficients of similar powers of $\rho$ in both sides to obtain

Coefficient of ${\rho }^0$ :

${\ddot{v}}_0+{\omega }_0^2{v}_0=0,$ (21)

Coefficient of $\rho$ :

${\ddot{v}}_1+R{v}_0^2{\ddot{v}}_0+R{v}_0{\stackrel{\dot{}}{v}}_0{}^2+{\omega }_0^2{v}_1+\frac{Rg}{2l}{v}_0^3=0,$ (22)

Coefficient of ${\rho }^2$ :

${\ddot{v}}_2+R{v}_0^2{\ddot{v}}_1+2R{v}_0{\ddot{v}}_0{v}_1+2R{v}_0{\stackrel{\dot{}}{v}}_0{\stackrel{\dot{}}{v}}_1+R{v}_1{\stackrel{\dot{}}{v}}_0{}^2+{\omega }_0^2{v}_2+\frac{3Rg}{2l}{v}_1{v}_0^2=0.$ (23)

The previous Equations (21)-(23) can be solved subsequently with the aid if the following conditions

$\begin{array}{l}{v}_0\left(0\right)=A,{\stackrel{\dot{}}{v}}_0\left(0\right)=0,\\ v_1\left(0\right)=0,{\stackrel{\dot{}}{v}}_1\left(0\right)=0,\\ v_2\left(0\right)=0,{\stackrel{\dot{}}{v}}_2\left(0\right)=0,\end{array}$ (24)

to get

$v_0=Acos\left({\omega }_{0}t\right)\mathrm{,}$ (25)

$v_1=\frac{A^{3}Rsin\left({\omega }_{0}t\right)}{32l{\omega }_0^2}\left[2{\omega }_0\left(4l{\omega }_0^2-3g\right)t+\left(4l{\omega }_0^2-g\right)sin\left(2{\omega }_{0}t\right)\right]\mathrm{,}$ (26)

$\begin{array}{c}{v}_2=\frac{A^5{R}^2}{1024l^2{\omega }_0^4}\{g^2\left[22\cos \left({\omega }_{0}t\right)-\cos \left(3{\omega }_{0}t\right)\right]{\sin }^2\left({\omega }_{0}t\right)+3g^2{\omega }_{0}t[8\sin \left({\omega }_{0}t\right)\\ -3\sin \left(3{\omega }_{0}t\right)]-2g{\omega }_0^2\left[2l+9g{t}^2-6l\cos \left(2{\omega }_{0}t\right)+4l\cos \left(4{\omega }_{0}t\right)\right]\cos \left({\omega }_{0}t\right)\\ +24gl{\omega }_0^{3}t\left[\sin \left({\omega }_{0}t\right)+2\sin \left(3{\omega }_{0}t\right)\right]+24l{\omega }_0^4[-2l+2g{t}^2+l\cos \left(2{\omega }_{0}t\right)\\ +l\cos \left(4{\omega }_{0}t\right)]\cos \left({\omega }_{0}t\right)\}-16l^2{\omega }_0^{5}t\{3\left[2\sin \left({\omega }_{0}t\right)+\sin \left(3{\omega }_{0}t\right)\right]\\ -2{\omega }_{0}t\cos \left({\omega }_{0}t\right)\}.\end{array}$ (27)

Since $\rho \in \left[0,1\right]$, one gets directly the desired solution when $\rho \to 1$ in the form

$\begin{array}{c}u=\underset{\rho \to 1}{\lim }v=A\cos \left({\omega }_{0}t\right)+\frac{A^{3}R\sin \left({\omega }_{0}t\right)}{32l{\omega }_0^2}[2{\omega }_0\left(4l{\omega }_0^2-3g\right)t\\ +\left(4l{\omega }_0^2-g\right)\sin \left(2{\omega }_{0}t\right)]+\frac{A^5{R}^2}{1024l^2{\omega }_0^4}\{g^2[22\cos \left({\omega }_{0}t\right)\\ -\cos \left(3{\omega }_{0}t\right)]{\sin }^2\left({\omega }_{0}t\right)+3g^2{\omega }_{0}t[8\sin \left({\omega }_{0}t\right)-3\sin \left(3{\omega }_{0}t\right)]\\ -2g{\omega }_0^2\left[2l+9g{t}^2-6l\cos \left(2{\omega }_{0}t\right)+4l\cos \left(4{\omega }_{0}t\right)\right]\cos \left({\omega }_{0}t\right)\end{array}$

$\begin{array}{c}+24gl{\omega }_0^{3}t\left[\sin \left({\omega }_{0}t\right)+2\sin \left(3{\omega }_{0}t\right)\right]+24l{\omega }_0^4[-2l+2g{t}^2\\ +l\cos \left(2{\omega }_{0}t\right)+l\cos \left(4{\omega }_{0}t\right)]\cos \left({\omega }_{0}t\right)\}\\ -16l^2{\omega }_0^{5}t\left\{3\left[2\sin \left({\omega }_{0}t\right)+\sin \left(3{\omega }_{0}t\right)\right]-2{\omega }_{0}t\cos \left({\omega }_{0}t\right)\right\}.\end{array}$ (28)

5. Results and Discussion

In this section, we are going to shed light on the great accuracy of the results obtained by HPT when they are compared with the numerical results of the governing equation of motion (4) using the fourth order Runge-Kutta method [26].

Figures 2-4 are calculated at $m_1=3\text{kg}\mathrm{,5}\text{kg}$ and $m_1=7\text{kg}$ respectively for different values of $l=0.1\text{m},0.6\text{m}$ and 1 m, in which their parts (a), (b) and (c) are plotted when $m_2=0.2\text{kg}\mathrm{,0.4}\text{kg}$ and $m_2=0.6\text{kg}$ respectively. It is

Figure 2. Describes the variation of the solution v via time t for different values of l at $m_1=3\text{kg}$ when (a) $m_2=0.6\text{kg}$, (b) $m_2=0.4\text{kg}$, (c) $m_2=0.2\text{kg}$.

Figure 3. Illustrates the variation of the solution v via time t for different values of l at $m_1=5\text{kg}$ when (a) $m_2=0.6\text{kg}$, (b) $m_2=0.4\text{kg}$, (c) $m_2=0.2\text{kg}$.

Figure 4. Shows the time history of the solution v for different values of l at $m_1=7\text{kg}$ when (a) $m_2=0.6\text{kg}$, (b) $m_2=0.4\text{kg}$, (c) $m_2=0.2\text{kg}$.

worthwhile to notice that these drawings have periodic forms and therefore the attained solution has a stable manner.

An inspection of the corresponding parts of these figures reveals that when l increases from 0.1 m to 1 m passing the value 0.6 m; the number of oscillations decreases and the wavelength of the ripples increases while the amplitudes of these ripples remain unchanged.

When parts (a) of Figures 2-4 are generally compared to parts (b) and (c) of the same figures, we observe that when m₂ increases from 0.2 kg to 0.6 m through the value 0.4 m; the number of oscillations increases and the wavelength of waves decreases beside the constancy of their amplitudes.

Moreover, these results are plotted in some figures for the same considered parameters; see Figures 5-7 when $m_1=5\text{kg}$ and $m_2=\left(\mathrm{0.2,0.4,0.6}\right)\text{kg}$. Figures 5-7 are calculated at $l=0.1\text{m},0.6\text{m}$ and 1 m respectively. It is not difficult to notice from the parts of Figure 5 that, the difference between both results seems to be small as in Figure 5(a) in which this difference increases in Figure 5(b) and Figure 5(c) to some extent during the time interval $\left[\mathrm{5.5,10}\right]$.

On the other side, this difference becomes very slightly which can be neglected as in Figure 6 and Figure 7. Then we can conclude that the results become better when the length l equals 0.6 m and 1 m, than the value 0.1 m which gives more opportunity of the extension x.

Tables 1-8 reveal a comparison between the results obtained by HPT with the

Figure 5. Shows a comparison between the numerical solution (Num) and the approximate one obtained by HPT at $m_1=5\text{kg}$ and $l=0.1\text{m}$ when (a) $m_2=0.2\text{kg}$, (b) $m_2=0.4\text{kg}$, (c) $m_2=0.6\text{kg}$.

Figure 6. Describe the comparison between the numerical solution (Num) and the approximate one obtained by HPT at $m_1=5\text{kg}$ and $l=0.6\text{m}$ when (a) $m_2=0.2\text{kg}$, (b) $m_2=0.4\text{kg}$, (c) $m_2=0.6\text{kg}$.

Table 1. Error percentage of HPT for $m_1=5,m_2=0.2,l=0.1$.

Table 2. Error percentage of HPT for $m_1=5,m_2=0.4,l=0.1$.

Table 3. Error percentage of HPT for $m_1=5,m_2=0.6,l=0.1$.

Table 4. Error percentage of HPT for $m_1=5,m_2=0.2,l=0.6$.

Table 5. Error percentage of HPT for $m_1=5,m_2=0.4,l=0.6$.

Table 6. Error percentage of HPT for $m_1=5,m_2=0.6,l=0.6$.

Figure 7. Shows a comparison between the numerical solution (Num) and the approximate one obtained by HPT at $m_1=5\text{kg}$ and $l=1\text{m}$ when (a) $m_2=0.2\text{kg}$, (b) $m_2=0.4\text{kg}$, (c) $m_2=0.6\text{kg}$.

Table 7. Error percentage of HPT for $m_1=5,m_2=0.2,l=1$.

numerical ones of the governing equation of motion (4) that obtained using the fourth order Runge-Kutta method and the corresponding error percentage of HPT; for different values of parameters of the considered dynamical model. This compression shows high consistency between them which expresses the great accuracy of the obtained solutions using HPT.

Table 8. Error percentage of HPT for $m_1=5,m_2=0.6,l=1$.

6. Stability Analysis

In this section, we investigate the stability of the governing equation of motion (13). It is obvious from the preceding section that, this investigation will be unsuccessful in view of Equation (28). Therefore, we are going to obtain a periodic solution of (13).

It should be noticed that Equation (13) is transformed into linear and nonlinear parts as indicated in Equations (15) and (16) respectively in which ${\omega }_0$ denotes a natural frequency of Equation (15). It is clear that the linear part represents a simple harmonic equation. Therefore, the stability of this part depends upon the frequency ${\omega }_0$ which is always positive and consequently, the represented figures have periodic forms as expected. Therefore the system is always stable.

Now, let us focus attention on the stability of a nonlinear part in which we consider a nonlinear frequency analysis. Therefore, a nonlinear frequency ${\Omega }^2$ is assumed to be in the following form

${\Omega }^2={\omega }_0^2+\rho {\varpi }_1+{\rho }^2{\varpi }_2+\cdots \mathrm{,}$ (29)

where ${\varpi }_1\mathrm{,}{\varpi }_2\mathrm{,}\cdots$ are arbitrary parameters can be estimated.

According to the reported work [30] and HPT, we can write the approximate nonlinear frequency in the form

${\Omega }^2={\omega }_0^2+\underset{\rho \to 1}{\lim }{\displaystyle {\sum }_{i=1}^{\infty }{\rho }^i{\varpi }_i}$ (30)

Substitution of (29) into (20) yields

$\ddot{v}+{\Omega }^{2}v+\rho R\left[v^2\ddot{v}+v{\stackrel{\dot{}}{v}}^2+\frac{g}{2l}{v}^3\right]-{\displaystyle {\sum }_{i=1}^{\infty }{\rho }^i{\varpi }_{i}v}=0$ (31)

Making use of (10) and (31), then equating the coefficients of like powers of $\rho$ in both sides to obtain

Coefficient of ${\rho }^0$ :

${\ddot{v}}_0+{\Omega }^2{v}_0=\mathrm{0,}$ (32)

Coefficient of $\rho$ :

${\ddot{v}}_1+R\left[v_0^2{\ddot{v}}_0+v_0\left({\stackrel{\dot{}}{v}}_0^2+\frac{g}{2l}{v}_0^2\right)\right]+{\Omega }^2{v}_1-{\varpi }_1{v}_0=\mathrm{0,}$ (33)

Coefficient of ${\rho }^2$ :

${\ddot{v}}_2+R\left[v_0^2{\ddot{v}}_1+2v_0\left({\ddot{v}}_0{v}_1+{\stackrel{\dot{}}{v}}_0{\stackrel{\dot{}}{v}}_1\right)+v_1\left({\stackrel{\dot{}}{v}}_0^2+\frac{3g}{2l}{v}_0^2\right)\right]+{\Omega }^2{v}_2-{\varpi }_1{v}_1-{\varpi }_2{v}_0=\mathrm{0,}$ (34)

Taking into account conditions (24), one can solve Equations (32)-(34) subsequently to get

$v_0=Acos\left(\Omega t\right)\mathrm{.}$ (35)

It is worthy to mention that in order to get a uniform to expand solution, the terms that produce secular terms in Equations (33) and (34) must be deleted. Substituting (35) into (33) and (34), then expanding the trigonometric functions to obtain

$\begin{array}{l}{\ddot{v}}_1+{\Omega }^2{v}_1=A\left({\varpi }_1-R{A}^2{\Omega }^2\right)\cos \left(\Omega t\right)\\ +\frac{1}{8l}R{A}^3\left(4l{\Omega }^2-g\right)\left[3\cos \left(\Omega t\right)+\cos \left(3\Omega t\right)\right],\\ {\ddot{v}}_2+{\Omega }^2{v}_2=\frac{1}{8l}A\left[8l{\varpi }_2+QAR\left(16l{\Omega }^2-3g\right)\right]\cos \left(\Omega t\right)\\ +\frac{21}{4}Q{A}^{2}R{\Omega }^2\left[\cos \left(2\Omega t\right)+\cos \left(4\Omega t\right)\right]\\ +\frac{1}{4l}Q\left[4l{\varpi }_1-A^{2}R\left(2l{\Omega }^2+3g\right)\right]\cos \left(3\Omega t\right)\\ +\frac{1}{8l}Q{A}^{2}R\left[\left(16l{\Omega }^2-3g\right)\cos \left(5\Omega t\right)+2l{\Omega }^2\left(1+\cos \left(6\Omega t\right)\right)\right].\end{array}$ (36)

Omitting terms that lead to secular terms in (36) to get

${\varpi }_1=\frac{1}{8l}{A}^{2}R\left(3g-4l{\Omega }^2\right),{\varpi }_2=\frac{1}{8l}QAR\left(3g-16l{\Omega }^2\right).$ (37)

According to (37), one can write the solutions $v_1\left(t\right)$ and $v_2\left(t\right)$ of (36) in the form

$\begin{array}{l}{v}_1=Q\left[\cos \left(3\Omega t\right)-\cos \left(\Omega t\right)\right],\\ v_2=\frac{Q}{8}\left[\frac{{\varpi }_1}{{\Omega }^2}+\frac{A^{2}R}{2}\left(\frac{631}{21}-\frac{7g}{4l{\Omega }^2}\right)\right]\cos \left(\Omega t\right)-\frac{Q{A}^{2}R}{4\Omega }\sin \left(\Omega t\right)\\ -\frac{7}{20}Q{A}^{2}R\left[5\cos \left(2\Omega t\right)+\cos \left(4\Omega t\right)\right]\\ +\frac{Q}{32l{\Omega }^2}\left[3R{A}^2\left(g+7l{\Omega }^2\right)-4l{\varpi }_1\right]\cos \left(3\Omega t\right)\\ +\frac{RQ{A}^2}{16}\left[\frac{1}{12}\left(3g-16l{\Omega }^2\right)\cos \left(5\Omega t\right)-\frac{1}{35}\cos \left(6\Omega t\right)+1\right],\end{array}$ (38)

where

$Q=\frac{R{A}^3}{16l{\Omega }^2}\left(g-l{\Omega }^2\right)\mathrm{.}$

Making use of (10), (35) and (38), then considering $\rho \to 1$, to obtain the approximate periodic solution in the form

$\begin{array}{l}u=A\cos \left(\Omega t\right)+Q\left[\cos \left(3\Omega t\right)-\cos \left(\Omega t\right)\right]\\ +\frac{Q}{8}\left[\frac{{\varpi }_1}{{\Omega }^2}+\frac{A^{2}R}{2}\left(\frac{631}{21}-\frac{7g}{4l{\Omega }^2}\right)\right]\cos \left(\Omega t\right)-\frac{Q{A}^{2}R}{4\Omega }\sin \left(\Omega t\right)\\ -\frac{7}{20}Q{A}^{2}R\left[5\cos \left(2\Omega t\right)+\cos \left(4\Omega t\right)\right]\\ +\frac{Q}{32l{\Omega }^2}\left[3R{A}^2\left(g+7l{\Omega }^2\right)-4l{\varpi }_1\right]\cos \left(3\Omega t\right)\\ +\frac{RQ{A}^2}{16}\left[\frac{1}{12}\left(3g-16l{\Omega }^2\right)\cos \left(5\Omega t\right)-\frac{1}{35}\cos \left(6\Omega t\right)+1\right].\end{array}$ (39)

An inspection of the previous solution u is given as a function of time t and has a periodic form. Therefore, the arguments of the trigonometric functions must be real values. To achieve this aim substituting (37) into (29) and considering $\rho \to 1$, we obtain

$l\left[A^{2}R\left(A^{2}R-4\right)-8\right]{\Omega }^4+\left[A^{2}gR\left(3-\frac{7}{16}{A}^{2}R\right)+8l{\omega }_0^2\right]{\Omega }^2+\frac{3A^4{R}^2{g}^2}{64l}=0.$ (40)

Under the present circumstances, the stability conditions require that $\Omega$ must be taken a real and positive quantity. Therefore, the necessary and sufficient conditions for the stability have the forms

$A^{2}R\left(A^{2}R-4\right)-8>0,A^{2}gR\left(3-\frac{7}{16}{A}^{2}R\right)+8l{\omega }_0^2<0$ (41)

To gain more insight into the existence of real roots, the distinction of (40) must be positive or becomes worthless i.e.,

$64l^2{\left({\omega }_0^2\right)}^2+gR{A}^{2}l\left(48-7A^{2}R\right){\omega }_0^2+\frac{1}{8}{A}^4{R}^2{g}^2\left[A^{2}R\left(\frac{1}{32}{A}^{2}R-15\right)+84\right]\ge 0$ (42)

Therefore, one obtains the restrictions on the initial angular velocity $\overline{{\omega }_0}$ have the forms

${\omega }_0^2\ge \frac{1}{128l}{A}^{2}gR\left(7A^{2}R-48+4\sqrt{3A^{2}R\left(A^{2}R-4\right)-24}\right)\mathrm{,}$ (43)

or

${\omega }_0^2\le \frac{1}{128l}{A}^{2}gR\left(7A^{2}R-48-4\sqrt{3A^{2}R\left(A^{2}R-4\right)-24}\right)$ (44)

Beside the first condition in (41), the stability region requires that

$\begin{array}{l}{\omega }_0^2<\frac{1}{8l}{A}^{2}gR\left(\frac{7}{16}{A}^{2}R-3\right),\\ \text{with}\\ {\omega }_0^2\ge \frac{1}{128l}{A}^{2}gR\left(7A^{2}R-48+4\sqrt{3A^{2}R\left(A^{2}R-4\right)-24}\right),\end{array}$ (45)

or

$\begin{array}{l}{\omega }_0^2<\frac{1}{8l}{A}^{2}gR\left(\frac{7}{16}{A}^{2}R-3\right)\\ \text{with}\\ {\omega }_0^2\le \frac{1}{128l}{A}^{2}gR\left(7A^{2}R-48+4\sqrt{3A^{2}R\left(A^{2}R-4\right)-24}\right)\end{array}$ (46)

Based on the above inequalities we can obtain another condition of the stability between the parameters A and R as follows: From (45), one can deduce the following inequality easily

$\frac{1}{16}\left[7A^{2}R-48+4\sqrt{3A^{2}R\left(A^{2}R-4\right)-24}\right]<\left(\frac{7A^{2}R}{16}-3\right)$

Therefore, one gets the stability condition between A and R in the form

$A^{2}R\left(4-A^{2}R\right)>8$ (47)

7. Conclusion

The motion of a nonlinear oscillating dynamical system is studied. HPT is used to achieve the solution of the governing equation of motion. The graphical representations of the obtained solution are represented for some different values of the physical parameters of the studied system. The numerical results of the governing equation of motion are obtained utilizing the Runge-Kutta method from fourth order and compared with the obtained ones by HPT. The comparison between them reveals high consistency in both results which emphasize the accuracy of the obtained results by HPT. The stability criteria is investigated through a fourth order equation in terms of the initial frequency ${\omega }_0$.

Funding

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

Conflicts of Interest

The authors declare that they have no conflict of interest.

References

[1]	Nayfeh, A.H. and Mook, D.T. (1979) Nonlinear Oscillations. Wiley, New York.
[2]	Minorsky, N. and Krieger, R.E. (1974) Nonlinear Oscillations. Huntington, New York.
[3]	Nayfeh, A.H. (1973) Perturbation Methods. Wiley, New York.
[4]	Nayfeh, A.H. (1985) Problems in Perturbation. Wiley, New York.
[5]	Awrejcewicz, J., Andrianov, I.V. and Manevitch, L.I. (1998) Asymptotic Approaches in Nonlinear Dynamics (New Trends and Applications). Springer, Berlin. https://doi.org/10.1007/978-3-642-72079-6
[6]	El-Dib, Y.O. and Moatimid, G.M. (2018) On the Coupling of the Homotopy Perturbation and Frobenius Method for Exact Solutions of Singular Nonlinear Differential Equations. Nonlinear Science Letters A, 9, 220-230.
[7]	He, J.H. (1999) Homotopy Perturbation Technique. Computer Methods in Applied Mechanics and Engineering, 178, 257-262 https://doi.org/10.1016/S0045-7825(99)00018-3
[8]	He, J.H. (2002) Modified Lindstedt-Poincaré Methods for Some Strongly Nonlinear Oscillations, Part I—Expansion of a Constant. International Journal of Non-Linear Mechanics, 37, 309-314. https://doi.org/10.1016/S0020-7462(00)00116-5
[9]	He, J.H. (2008) Max-Min Approach to Nonlinear Oscillator. International Journal of Nonlinear Sciences and Numerical Simulation, 9, 207-210. https://doi.org/10.1515/IJNSNS.2008.9.2.207
[10]	Ganji, S.S., Ganji, D.D., Karimpour, S. and Babazadeh, H. (2009) Applications of He’s Homotopy Perturbation Method to Obtain Second-Order Approximations of the Coupled Two-Degree-of-Freedom Systems. International Journal of Nonlinear Sciences and Numerical Simulation, 10, 303-312. https://doi.org/10.1515/IJNSNS.2009.10.3.305
[11]	Yildirim, A. (2010) Determination of Periodic Solutions for Nonlinear Oscillators with Fractional Powers by He’s Modified Lindstedt-Poincaré Method. Meccanica, 45, 1-6. https://doi.org/10.1007/s11012-009-9212-4
[12]	El-Dib, Y.O. (2017) Homotopy Perturbation for Excited Nonlinear Equations. Science and Engineering Applications, 2, 96-108. https://doi.org/10.26705/SAEA.2017.2.1.96-108
[13]	El-Dib, Y.O. (2017) Multiple Scales Homotopy Perturbation Method for Nonlinear Oscillators. Nonlinear Science Letters A, 8, 352-364.
[14]	El-Dib, Y.O. (2018) Periodic Solution and Stability Behavior for Nonlinear Oscillator Having a Cubic Nonlinearity Time-Delayed. International Annals of Science, 5, 12-25. https://doi.org/10.21467/ias.5.1.12-25
[15]	Bayat, M., Pakar, I. and Bayat, M. (2015) Nonlinear Vibration of Mechanical Systems by Means of Homotopy Perturbation Method. Kuwait Journal of Science, 42, 64-85.
[16]	He, J.H. (2004) Comparison of Homotopy Perturbation Method and Homotopy Analysis Method. Applied Mathematics and Computation, 156, 527-539. https://doi.org/10.1016/j.amc.2003.08.008
[17]	Biazar, J. and Ghazvini, H. (2008) Numerical Solution for Special Non-Linear Fredholm Integral Equation by HPM. Applied Mathematics and Computation, 195, 681-687. https://doi.org/10.1016/j.amc.2007.05.015
[18]	Shaher, M., Erjaee, G.H. and Alnasr, M.H. (2009) The Modified Homotopy Perturbation Method for Solving Strongly Nonlinear Oscillators. Computers and Mathematics with Applications, 58, 2209-2220. https://doi.org/10.1016/j.camwa.2009.03.082
[19]	Öziş, T. and Akçı, C. (2011) Periodic Solutions for Certain Non-Smooth Oscillators by Iterated Homotopy Perturbation Method Combined with Modified Lindstedt-Poincaré Technique. Meccanica, 46, 341-347. https://doi.org/10.1007/s11012-010-9312-1
[20]	He, J.H. (2012) Homotopy Perturbation Method with an Auxiliary Term. Abstract and Applied Analysis, 2012, Article ID: 857612. https://doi.org/10.1155/2012/857612
[21]	Ghasemi, S.E., Zolfagharian, A. and Ganji, D.D. (2014) Study on Motion of Rigid Rod on a Circular Surface Using MHPM. Propulsion and Power Research, 3, 159-164. https://doi.org/10.1016/j.jppr.2014.07.003
[22]	Hosen, M.A. (2014) Approximate Solutions of the Equation of Motion’s of the Rigid Rod Which Rocks on a Circular Surface without Slipping. Ain Shams Engineering Journal, 5, 895-899. https://doi.org/10.1016/j.asej.2014.01.005
[23]	Belndez, A., Mndez, D., Belndez, T., Hernndez, A. and lvarez, M.L. (2008) Harmonic Balance Approaches to the Nonlinear Oscillators in Which the Restoring Force Is Inversely Proportional to the Dependent Variable. Journal of Sound and Vibration, 314, 775-782. https://doi.org/10.1016/j.jsv.2008.01.021
[24]	El-Dib, Y.O. and Moatimid, G.M. (2019) Stability Configuration of a Rocking Rigid Rod over a Circular Surface Using the Homotopy Perturbation Method and Laplace Transform. Arabian Journal for Science and Engineering, 44, 6581-6591. https://doi.org/10.1007/s13369-018-03705-6
[25]	Enayati, S.G., Azimi, M. and Jouya, M. (2017) Application of Modified Homotopy Perturbation Method and Amplitude Frequency Formulation to Strongly Nonlinear Oscillators. NTMSCI, 5, 66-82. https://doi.org/10.20852/ntmsci.2017.127
[26]	Shen, Y. and El-Dib, Y.O. (2020) A Periodic Solution of the Fractional Sine-Gordon Equation Arising in Architectural Engineering. Journal of Low Frequency Noise, Vibration and Active Control. https://doi.org/10.1177/1461348420917565
[27]	El-Dib, Y.O. and Elgazery, N.S. (2020) Effect of Dractional Derivative Properties on the Periodic Solution of the Nonlinear Oscillations. Fractals. https://doi.org/10.1142/S0218348X20500954
[28]	Ji, Q.-P., Wang, J., Lu, L.-X. and Ge, C.-F. (2020) Li-He’s Modified Homotopy Perturbation Method Coupled with the Energy Method for the Dropping Shock Response of a Tangent Nonlinear Packaging System. Journal of Low Frequency Noise, Vibration and Active Control. https://doi.org/10.1177/1461348420914457
[29]	Gilat, A. (2013) Numerical Methods for Engineers and Scientists. Wiley, Hoboken.
[30]	El-Dib, Y.O. (2018) Stability of a Strongly Displacement Time-Delayed Duffing Oscillator by the Multiple Scales-Homotopy Perturbation Method. Journal of Applied Mathematics and Computational Mechanics, 4, 260-274.