10464 FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 20, No 1, 2022, pp. 1 - 20 https://doi.org/10.22190/FUME220118004H © 2022 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper FORCED NONLINEAR OSCILLATOR IN A FRACTAL SPACE Ji-Huan He1,2,3, Galal M. Moatimid4, Marwa H. Zekry5 1School of Civil Engineering, Xi’an University of Architecture & Technology, Xi’an, China 2National Engineering Laboratory for Modern Silk, College of Textile and Clothing Engineering, Soochow University, Suzhou, China 3School of Science, Xi'an University of Architecture and Technology, Xi’an, China 4Department of Mathematics, Faculty of Education, Ain Shams University, Roxy, Cairo, Egypt 5Department of Mathematics and Computer Science, Faculty of Science, Beni-Suef University, Beni-Suef, Egypt Abstract. A critical hurdle of a nonlinear vibration system in a fractal space is the inefficiency in modelling the system. Specifically, the differential equation models cannot elucidate the effect of porosity size and distribution of the periodic property. This paper establishes a fractal-differential model for this purpose, and a fractal Duffing-Van der Pol oscillator (DVdP) with two-scale fractal derivatives and a forced term is considered as an example to reveal the basic properties of the fractal oscillator. Utilizing the two-scale transforms and He-Laplace method, an analytic approximate solution may be attained. Unfortunately, this solution is not physically preferred. It has to be modified along with the nonlinear frequency analysis, and the stability criterion for the equation under consideration is obtained. On the other hand, the linearized stability theory is employed in the autonomous arrangement. Consequently, the phase portraits around the equilibrium points are sketched. For the non-autonomous organization, the stability criteria are analyzed via the multiple time scales technique. Numerical estimations are designed to confirm graphically the analytical approximate solutions as well as the stability configuration. It is revealed that the exciting external force parameter plays a destabilizing role. Furthermore, both of the frequency of the excited force and the stiffness parameter, execute a dual role in the stability picture. Key words: Fractal vibration theory, Duffing-Van der Pol oscillator, Homotopy perturbation method, Expanded frequency analysis, Linearized stability theory, Multiple timescales technique Received January 18, 2022 / Accepted February 12, 2022 Corresponding author: ji-Huan He School of Civil Engineering, Xi’an University of Architecture & Technology, Xi’an 710055, China E-mail: hejihuan@suda.edu.cn 2 J.-H. HE, G. M. MOATIMID, M. H. ZEKRY 1. INTRODUCTION The fractal vibration theory has become a significant topic in both mathematics and mechanical engineering. The smart receptor system of a three-dimensional printing processing can be exactly controlled by the gecko-like fractal vibration system [1]. The fractal Toda oscillator exhibits some special properties, which the traditional vibration theory cannot account for [2]. The pull-in instability can be completely eliminated for an electrostatic nano/micro-electromechanical (N/MEMS) system by suitably controlling the fractal vibration system [3-5]. The nonlinear wave travelling in a fractal space can be converted into a fractal oscillator, and some fascinating properties are revealed [6]. The fractal vibration theory provides an effective tool for studying various unusual phenomena under the micro-gravity condition [7-9]. The fractal vibration can effectively model the vibration arising in porous media, and the low frequency property plays an important role in the anti-vibration design [10, 11]. Furthermore, the fractal Fangzhu oscillator explores the mechanism of an ancient technology for the collection of water from air [12]. Much literature focused on simple fractal oscillators [13-15], this paper considers a more generalized fractal vibration system in the following form: 2 2 3 2 ( ) cos , d z dz b cz az dz f d d       − − + + = (1) where dz/dτα is the two-scale fractal derivative [16 and 17]; z,τ,b,c,a,d,f and ω are displacement, time, the coefficient of the damping term, the coefficient of the nonlinear damping term, the root of natural frequencies, the cubic stiffness parameter, the amplitude of the excitation force, and frequency of this excitation, respectively. Eq. (1) displays a nonlinear vibration system in a fractal space as discussed by Tian et al. [3-5], where the air is considered as a porous medium and the fractal dimension α reflects the air distribution. When it vibrates in a vacuum, one has α=1. When f=b=c=0, Eq. (1) becomes the fractal Duffing oscillator [1 and 13-15], which can model the nonlinear vibration arising in a three-dimensional printing technology, concrete vibration and other applications. For more convenience, the graph x() = cos α will be plotted for various values of the parameter α in Fig. 1. It is obvious that the periodic solution becomes an asymptotic one. When α=1, Eq. (1) reflects the forced Duffing -Van der Pol driven oscillator. In view of the two-scale transforms [16], one gets z=x√b/c, τα=t/√a, b=µ√a and ω=Ω√a to convert Eq. (1) into the following form: 2 3 (1 ) cosx x x x x F t − − + + =  , (2) where β=bd/ac, F=(f/a)√c/b, , and the dot is the derivative with respect to t. The two-scale fractal theory is now widely used to model various phenomena arising in porous media or un-smooth boundaries. The fractal rheological model gives a quick physical insight into the rheological property of a non-Newtonian fluid [18 and 19]. The fractal approach provides a simple and good tool for mechanical and electrical properties of composite materials [20]. The fractal two-phase flow model can figure out the hidden properties in a polymer filling process, which cannot be revealed by the different models [21-23]. Now, the fractal variational principle has become a predominant topic in the Forced Nonlinear Oscillator in a Fractal Space 3 establishment of governing equations in a fractal space [24 and 25]. It is worthy to observe that when µ=0, Eq. (2) yields the Duffing equation. Simultaneously, the cancellation of the parameter β results in the VdP oscillator. In 1920, Van der Pol performed experimental studies on the dynamics of an oscillator electrical circuit. It was formulated by a second-order nonlinear differential equation, which was later known as the Van der Pol oscillator. Alghassab et al. [26] examined the Duffing oscillator together with the Van der Pol oscillator with a nonlinear damping term (DVdP). Apart from the linearized analysis, they exploited the Lyapunov theory in their analysis. Zhihong and Shaopu [27] presented a novel weak signal detection approach grounded on the DVdP oscillator. They found that the planned method was stronger than the Duffing oscillator. Under certain parametric conditions, Gao and Feng [28] deduced the first integrals without complicated computations of the DVdP. Their methodology was constructed of a series of variable transformations with the Preller-Singer method. Motsa and Sibanda [29] presented a novel application of the successive linearized method to the classical DVdP. They obtained accurate values of frequency and amplitude. Utilizing the homotopy analysis method, Cui et al. [30] explored a high accuracy curve and solutions of the DVdP. Their stability analysis adopted the Floquet theory. Additionally, their results were validated in light of the spectra analysis and bifurcation theory. Moatimid [31] examined the stability criteria of the parametric Duffing oscillator. The numerical calculations showed that the damper and the stiffness parameters have a destabilizing influence. The study of the parametric excited systems is actually imperative in light of external excited fields. Wei et al. [32] examined the major response of the DVdP to collective deterministic and arbitrary parametric excitations. They employed the method of multiple time scales to govern the equations of amplitude and phase. Dao et al. [33] considered the DVdP under the parametric exciting force. They evaluated the cases of a small parameter quasi-linear as well as the general case. Han et al. [34] scrutinized the effects of a gradually changing parametric excitation of the DVdP equation. They examined the periodic bifurcation delay behaviors when there are no spinning points. Their study showed that, at the first delay behavior, the delay time is decided by the initial time. The phenomenon of the entertaining of a generalized VdP was considered by Kovacic [35]. By regulating the method of averaging, the approximate amplitude-frequency response was obtained. The inspiration of the controls of the reinstating and damping force in the Fig. 1 Basic property of the fractal vibration, x() = cosα 4 J.-H. HE, G. M. MOATIMID, M. H. ZEKRY existence of this occurrence was examined. Shen et al. [36] investigated the primary and sub-harmonic simultaneous resonance of the Duffing oscillator with a fractional order derivative. The approximate analytical solution of the resonance was obtained. Furthermore, the rightness and acceptable exactness of the analytical solution were verified by a numerical simulation. Recently, Wu et al. [37] have examined periodic explanations of a harmonic forced Duffing oscillator with a time-delay state feedback by using the incremental harmonic balance method. Fractal differential equations are convenient in describing the numerous occurrences in porous media or on un-smooth boundaries. It is to be noted that the exact solutions are generally unreachable. Therefore, these nonlinear equations should be solved by expanding approximate methods. Many efforts were exerted to analytically treat these phenomena, ranging from the traditional perturbation methods to the multiple time scale technique. As known, the straightforwardness of the perturbation methods renders the nonlinear equations into linear ones. Certainly, all these perturbation methods are primarily contingent on a small parameter. Therefore, the presence of such a parameter is necessary to attain an approximate solution. Frequently, the solution of such a problem might be limited by the presence of this parameter. Fundamentally, the exact solution of this equation is difficult and a perturbation technique is required. All classical perturbation methods ranging from the traditional methods to the multiple-scales method depend mainly on a small parameter. Consequently, the problem cannot be solved without the occurrence of such a parameter. In recent years, many promising and powerful methods have been introduced to construct analytical solutions for such equations. One of the pioneer methods is the Homotopy perturbation method (HPM). This promising approach was first proposed by the Chinese mathematician He [38]. Homotopy is constructed in terms of a small embedded parameter. The HPM is an amalgamation of the traditional perturbation performance and the Homotopy method, which has removed the dominant restrictions of the perturbation methods. This HPM brings about an actual quick convergence of the solution series in most cases. Therefore, only after a few iterations, one reaches accurate solutions. HPM is a universal technique which can solve various kinds of nonlinear equations. This technique surpasses the traditional perturbation techniques. Gómez-Aguilar et al. [39] introduced the fractional equations of the mass-spring-damper system with Caputo and Caputo–Fabrizio derivatives. Their results indicated that the mechanical components display visco-elastic performances manufacturing temporal fractal diverse scales. Yang et al. [40] considered boundary-value problems in a fractal heat transfer. The exact solutions of non- differentiable type were attained by the local fractional differential transform method. Recently, Elías-Zúñiga et al. [41] introduced a new methodology for descending the careful, steady-state solutions of fractal damped and forced differential equations. Motivated by the aforementioned aspects, the aim of this examination is to investigate the DVdP since it is widely used in different fields. As previously shown, the dynamics of this equation have been widely investigated in recent years. Therefore, the task now is to investigate an analytical approximate solution together with the stability profile. To develop the existing problem, the rest of the paper is organized as follows. A periodic solution, based on an adaptation of the expanded natural frequency analysis, is given in Section 2. The stability analysis of the autonomous system, based on the linearized theory, is realized in Section 3. Additionally, the phase portraits are drawn in this Section. In Section 4, the stability analysis of the resonance as well as the non-resonance cases are accomplished by utilizing a combination of the Homotopy concept and the multiple time Forced Nonlinear Oscillator in a Fractal Space 5 scale analysis. Finally, concluding remarks grounded on the consequences are summarized in Section 5. 2. AN APPROXIMATE BOUNDED SOLUTION THROUGH RESOURCES OF THE NONLINEAR FREQUENCY ANALYSIS As known, the classical HPM fails to solve the damping nonlinear oscillator, especially for the Van der Pol oscillator because it contains linear and nonlinear damping forces. Therefore, the main purpose of this Section is to attain an analytical bounded approximate solution for the DVdP equation. The the classical HPM does not enable us to satisfy this end properly. Consequently, one necessarily seeks another new technique to realize a bounded solution of the considered equation. To this end, consider the Homotopy equation that is given by Eq. (3), where the natural frequency of the current model is the unity. 2 3 ( ( 1) cos ) 0; [0,1]x x x x x F t   + + − + −  =  (3) Along with this approach, the time dependent function may be determined as follows:   = = 0 )();( n n n txtx  (4) The following stability examination will be constructed on the nonlinear frequency analysis given by Moatimid [31]. In accordance with this approach, an expanded artificial frequency σ2 may be formulated as follows:   = += 1 2 1 j j j  (5) where the parameters σj (j=1,2,…) will be determined from the invalidation of the sources of the secular terms. By combining Eqs. (3) and (4), the Homotopy equation may be rewritten as: 2 2 2 3 1 2 3( ( 1) ( ) cos ) 0x x x x x x F t       + + − − + + + −  = (6) It is worth mentioning that if the standard initial conditions are considered as approached by Moatimid [36], the stiffness parameter µ will disappear. Therefore, to preserve the presence of this parameter in the following analysis, an adjustment of the initial conditions must take place. Consequently, a modification of the initial conditions may be written as follows: ,0)0( =x and   = = 0 )0( j j j x  (7) where the constants λi will be determined later as a function of the characteristics of the given equation. 6 J.-H. HE, G. M. MOATIMID, M. H. ZEKRY This will be done, together with the parameters σj as given in Eq. (5), in light of the cancellation of the secular terms. Taking the Laplace transforms to Eq. (6), while considering the initial conditions given in Eq. (7), the result becomes 2 30 31 2 2 2 2 2 2 2 2 2 2 2 3 1 2 32 2 { ( ; )} .... 1 { ( ( 1) ( ) cos )}. T T L x t s s s s L x x x x F t s                    = + + + + + + + + − − − + + + −  + (8) Employing the inverse transform of both sides of Eq. (8), one gets: 2 30 31 2 1 2 2 3 1 2 32 2 ( ; ) ...... sin 1 { ( ( 1) ( ) cos )}T T x t t L L x x x x F t s                     −   = + + + + −      − − + + + −   +  (9) By utilizing the expansion of the dependent parameter x(t;ρ) as given in Eq. (4), and by equating the coefficients of like powers ρ on both sides, one gets ,s in)(: 00 0 ttx     = (10) 1 2 31 1 0 0 1 0 02 2 1 : ( ) sin( ) { ( 1) cos }T Tx t t L L x x x x F t s         −   = − − − + −   +  (11) Characteristically, the uniform effective expansion emerges from the invalidation of the secular terms. For this purpose, one gets the following results: ▪ The cancellation of the secular terms, i.e. excluding the coefficients of the functions sin σt and cos σt gives  20 = (12) and  31 = (13) The second-order of Eq. (9) yields 2 1 2 22 2 1 0 0 0 1 1 1 2 0 0 12 2 1 : ( ) sin { ( ( 1) 2 ) 3 }T Tx t t L L x x x x x x x x x s          −   = − − + − − +  +  (14) ▪ Once more, the cancellation of the secular terms needs  4/51 −= (15) and 2222 2 8/)9(  +−= (16) The third-order of Eq. (9) results in 2 2 2 0 0 1 0 2 1 0 13 13 3 2 2 2 2 1 2 2 1 3 0 0 1 2 0 ( ( 1) ( 2 ) 2 )1 : ( ) sin 3 ( ) T T x x x x x x x x x x t t L L s x x x x x x x          −   − + + + −   = −    + − − + +     (17) Forced Nonlinear Oscillator in a Fractal Space 7 ▪ Once more, the cancellation of the secular terms necessitates 2 2 2 2 2 2 2 2 2 4 2 2 3 4 2 2 4 2 3 2 2 2 2 2 (81 5 )( ) (9 ) 72 (9 ) 24 (69 6 ) , 96 ( ) (9 ) F F              + − − − − + −  + = − − (18) and 4222 3 64/)9(9  += (19) It follows that the bounded solution x1(t) then becomes 2 2 2 2 21 1 2 2 2 2 2 2 2 1 ( ) sin ( ( 4 ( )) cos 3 ( ) sin 4 cos 4 ( ) ( ) cos 3 ( ) sin 3 ). x t t F t t F t t t                     = + − + − + − +  − − − − − (20) As mentioned before, for an easy follow-up of the analysis, the function x2(t) will be given in the Appendix. Therefore, the periodic analytic approximate solution of the governing equation of motion given in Eq. (2) may be written as follows: 2 0 1 2 1 ( ) lim ( ( ) ( ) ( ))x t x t x t x t    → = + + (21) By the inverse of two-scale transforms, t=√aτα, we can obtain the following solution of the fractal DVdP equation along with x0(√aτ α)+ x1(√aτ α)+ x2(√aτ α). In fact, the bounded approximate solution as given in Eq. (21) Necessitates that the arguments of the trigonometric functions are real. For this purpose, combining Eqs. (13), (16) and (19) into Eq. (5), it follows that the synthetic frequency satisfies a certain characteristic equation. The calculation showed that this equation represents a polynomial of sixth order. Once more, by means of HPM, one may obtain the approximate solution as follows: 2 4 6 ( ) , [0, 1]a b c    = + −  (22) where the constants a, b and c, may be listed as: )8/(9 22  −=a , ,)8(9/)819264( 22  −−+=b .)8(9/64 2  −=c (23) Following a similar procedure as given above, an approximate solution up to the second-order, of the simulated frequency, may be written as: 2 2 10 sss  ++ (24) where the constants si are given as: 2/)(, 2/3 10 acbasas +== (25) and )11187( 8 1 2222/5 2 cacbabas ++= (26) In accordance with the real value of the expanded frequency, the approximate bounded solution requires the following criterion: .08 2 − (27) 8 J.-H. HE, G. M. MOATIMID, M. H. ZEKRY Eq. (27) reads that the analytic approximate bounded solution cannot be realized where the stiffness parameter becomes negative. For more opportunities, the attained analytic approximate analytic solution must validate the numerical solution. For this objective, the fourth-order Runge–Kutta method (RK4) it utilized to confirm the approximate solution. Therefore, Fig. 2 is designed for demonstrating the analytic bounded solution given from the expanded frequency approach together with the numerical solution that is given by RK4. Through this figure, the usage of the sample chosen system showed that the approximate value of the expanded frequency equals σ≌ 3.29. Additionally, the calculations revealed that the initial condition becomes ẋ(0)=5.821. As shown from this figure, the amplitudes of the solution in the two cases are nearly coincident. This indicates a good agreement between the two methodologies. Fig. 2 Comparison of the approximate solution that is given in Eq. (21) with the numerical one for the case: µ=0.2, β=2, F=0.05 and Ω=0.4 Additionally, the bounded solution of the considered DVdP will be graphed. This solution is given by Eq. (21). Simultaneously, it depends mainly on the value of the expanded frequency (σ) that is given in Eq. (24). It is obvious that the validation of the value σ requires that 8β >µ2 as specified in the inequality (27). To prove the accuracy of the method, we depict in Fig. 3, the percent relative error E(t)=x2(t)/(x0(t)+x1(t) +x2(t)), i.e., the stepwise error, from the results depicted in this figure, we assure the accuracy as well as the stability of raising the order of the approximation. Fig. 3 The error function E(t) for a system having the particulars: µ=0.2, β=2, F=0.05 and Ω=0.4 Forced Nonlinear Oscillator in a Fractal Space 9 Along with these restrictions, in what follows two graphs are plotted for some chosen systems to indicate the bounded solutions, as shown in Figs. 4, and 5. Fig. 3 displays the analytical bounded approximate solution for a sample chosen system. The inspection of this figure indicates the behavior of the solution, which is a bounded one. Throughout this figure, the calculations showed that the approximate value of the expanded frequency equals σ≌ 8.716. Additionally, the calculations illustrated that the initial condition developed as ẋ(0)=13.017. Fig. 4 The approximate solution that is given in Eq. (21) for a system having the particulars: µ=0.2, β=2, F=50 and Ω=4 By contrast, Fig. 4 represents a system whose particulars are held fixed, except for the values of the parameter F. As shown from this figure, the amplitude of the wave solution increases as F also increases. This shows a destabilizing influence of this parameter. Physically, it is a convenient evidence because of the excitation of the external force. It should be noticed that the domain of the solution function is much shorter than the previous one. This occurs only to clearly indicate the influence. Fig. 5 The approximate solution that is given in Eq. (21), for a system having the particulars given in Fig. 2 but for the variation of the parameter F. 10 J.-H. HE, G. M. MOATIMID, M. H. ZEKRY 3. LINEARIZED STABILITY OF THE AUTONOMOUS SYSTEM The analysis here is concerned with the governing DVdP equation at which the external force is absent as (F→ 0). For this purpose, consider the following transformations: ,yx = and .yx  = (28) It follows that the governing equation that is given by Eq. (2) may be represented by the following first -order nonlinear ordinary differential equations: ),( yxfx = and ,),( yxhy = (29) where f(x,y)=y and h(x,y)=-x+µ(1−x2)y−β x3. The details of the current approach are given in our previous recent work [42]. To restrict the length of the paper, they will be omitted. For more opportunities, some numerical calculations are given here to indicate the nature of the characteristic eigenvalues and consequences. Therefore, the phase portraits are plotted. Table 1 Illustration of some sample values of the eigenvalues and the corresponding stability/instability Sample Chosen System Fixed Point Roots of the Eigenvalues Stability/Instability 1 0,1 ==  (0, 0) Pure imaginary .2,1 i= A stable center See Fig. 6 2 1,1 =−=  ( 1, 0) Real, equal, and of different signs .414.1,414.1 21 −== Saddle point (Unstable) See Fig. 7 3 5.0,2 =−=  1 , 0 2       Real, distinct, and of different signs .29473.1 ,54473.1 2 1 −= = Saddle point (Unstable) See Fig. 8 4 1,1 −=−=  (0, 0) complex conjugates with negative real part )31( 2 1 1 i−= Stable focus See Fig. 9 Unfortunately, the previous criterion ignores the contribution of the external excitation force. Therefore, the following analysis will be included to indicate this influence throughout the multiple scale technique. Forced Nonlinear Oscillator in a Fractal Space 11 Fig. 6 The phase portrait of a system having the particular: β =1 and µ=0. Fig. 7 The phase portrait of a system having the particular: β =−1 and µ=1. Fig. 8 The phase portrait of a system having the particular: β =−2 and µ=0.5. Fig. 9 The phase portrait of a system having the particular: β =−1 and µ=−1. 4. THE STABILITY ANALYSIS ALONG THE MULTIPLE-TIME SCALE TECHNIQUE As known the Poincare-Lindstedt method is an approach to postulate an asymptotic approximation of the periodic solution. On the other hand, it cannot be used to attain solutions that develop periodically on a slow timescale. The method of multiple scales is a more general approach that involves two key tricks. Therefore, the method of multiple scales encompasses a technique that is used to establish a uniform valid approximation to the general solution of perturbation problems in which the solution depends on extensively diverse scales. Nayfeh [43] divided the time into several time scales which will be considered to be independent variables. Based on this conception, Nayfeh formulated his theory. In agreement with the transformation of the derivatives, the original differential equation is divided into a set of partial differential equations. Subsequently, it seems that 12 J.-H. HE, G. M. MOATIMID, M. H. ZEKRY the problem has been complicated. The results obtained from this approach have shown that there are disadvantages of this complexity, which are far overbalanced by more compensations. This technique does not only offer a uniform expansion, but it is also still valid to be applied to numerous extremely nonlinear problems in extensive several areas of practical applications and engineering. Now, one seeks a uniform approximate solution Eq. (2) in the following form: )(),(),();( 2 101100  OTTxTTxtx ++= (30) It should be noted that the number of the independent time scales generates the appropriate order of expansion. Specifically, if the expansion is carried out up to the second-order, then only three scales T0, T1 and T2 are required. Replacing Eq. (30) as well as the transformation of the derivatives given below in the Homotopy equation Eq. (5), and then comparing the coefficient of like powers ρ, one finds the following equations: 0 2 0 0: ( 1) 0D x + = (31) and 0 0 2 2 3 0 1 0 1 0 0 0 0 0: ( 1) 2 ( 1) ( ) 2 i T i TF D x D D x x D x x e e    −  + = − − − − + + (32) where ....... 2 2 10 2 2 1 1 0 0 +++=+   +   +    DDD Tdt dT Tdt dT Tdt dT dt d  and 2 2 2 2 0 0 1 1 0 22 2 ( 2 ) ..., n n d D D D D D D D Tdt     + + + +   Therefore, it is convenient to write the solution of Eq. (32) in the following form: .,.)(),( 01100 cceTATTx iT += (33) where A(T1) is an unknown complex function of the time scale T1, and ..cc characterizes the complex conjugate of the foregoing terms. The governing equation of the complex function A(T1) is obtained in such a way that the time-dependent functions x1 must be of an oscillatory nature. Substituting Eq. (33) into Eq. (32), one finds 0 0 032 2 3 0 1 1 1 ( 1) ( 2 ( 3 ) ) ( ) . ., 2 iT iT i T D x iD A i A i A A e i A e Fe c c      + = − + − + − + + + (34) where the function Ā represents the complex conjugate of the function A. The uniform valid expansion of the function x1 requires a cancellation of the secular term. Accordingly, in order to attain a uniform valid expansion, the secular terms must be omitted. The invalidation of these secular terms arises from the coefficients of the exponentials Exp(±iT0). Therefore, the uniform valid expansion requires 0)3(2 2 1 =+−+− AAiAiAiD  (35) Forced Nonlinear Oscillator in a Fractal Space 13 Eq. (35) is a nonlinear first-order differential equation with complex coefficients. It is known as the amplitude equation. Occasionally, it is called the solvability condition. It follows that uniform solution of the particular solution of Eq. (34) may be formulated as .. )1(2 1 )( 8 1 ),( 00 2 33 101 ccFeeAiTTx TiiT + − ++=   (36) 4.1. Stability Analysis in the Non-Resonance Case To complete the stability profile, along with the non-resonance case, one goes back to the amplitude Eq. (35). Essentially, this equation might be utilized to regulate the indefinite function A in terms of the time-independent variable T1. Additionally, the stability presentation is contingent principally on the nature of this function. For this determination, one may apply the polar form solution as: )( 11 1)( 2 1 )( Ti eTTA  = (37) where the functions ζ (T1) and 𝛾(T1) are two real functions on the time scale 1T . Substituting Eq. (37) into Eq. (35), and then associating the real and imaginary parts on the two sides, one gets the following differential equations: 3 1 82     −= dT d (38) and 2 1 8 3   = dT d (39) Eqs. (38) and (39) are coupled first-order differential equations with real coefficients. As ρ→1, the parameter T1 will be transformed into t. Certainly, the stability conditions need that both of the functions ζ(t) and 𝛾(t) become bounded functions. For straightforwardness, these equations may be solved by using a perturbation approach as follows: Consider that the function ζ (t) may be perturbed around a steady-state response as: ζ (t)= ζ 0+ ζ 1(t). In this case, the zero-order solution yields 20 = (40) It follows that the first-order solution becomes t et   − =)(1 (41) In order to obtain a finite solution, the coefficient of the damped parameter (µ) must be positive. Combining equations (40), (41) and (39), the first-order solution of the function 𝛾1(t) may be written as: 1 3 ( ) (1 ), 2 t t e    − = − (42) where the initial condition: 𝛾1(0)=0 is used. 14 J.-H. HE, G. M. MOATIMID, M. H. ZEKRY It is worthy to notice that the stability occurs regardless of the sign of the parameter β, provided that 𝜇 >0. As a final result in the case of non-resonant, the approximate periodic solution of the DVdP given in Eq. (2) may be written as follows: 0 1 1 ( ) lim( ( ) ( ))x t x t x t   → = + (43) where x0(t) and x1(t) are the time-dependent functions obtained from Eqs. (33) and (36), respectively. Because of the presence of the non-homogeneous oscillatory term through the non- resonance case, the analysis will examine the resonance cases. Therefore, the following Subsection discusses the resonance case. 4.2. Stability Analysis of the Resonance Case As seen from the previous subsection, the non-resonance case fails to contain all parameters of the problem at hand during the stability criteria. Consequently, the following discussion is devoted to presenting resonance cases. Typically, only the resonance case is obtained. This case is definitely harmonic resonance. It is defined here as a sharp resonance. This benefit helps in introducing extra secular terms. Accordingly, one may write += 1 (44) where 𝛿 is some detuning small parameter. It follows that the governing equation x1 has an extra secular term. Inserting Eq. (44) into Eq. (34), one gets 0 01 32 2 3 0 1 1 1 ( 1) ( 2 ( 3 ) ) ( ) . ., 2 iT iTi T D x iD A i A i A A Fe e i A e c c      + = − + − + + − + + (45) The removal of the secular term from Eq. (45) requires 02/)3(2 1 2 1 =++−+− Ti FeAAiAiAiD   (46) Really, the particular integral of Eq. (45) is different from the previous corresponding equation as given in Eq. (34). Subsequently, the periodic approximate solution of the first- order will be transformed. The particular solution, in case of the harmonic resonance, becomes ..)( 8 1 ),( 0 33 101 cceAiTTx iT ++=  (47) Eq. (46) is a first-order nonlinear differential equation with complex coefficients. As previously mentioned, its solution may be considered as: 1)()( 11 Ti eTTA  = , (48) where 𝜓 (T1) is a real and time-dependent function of the independent parameter T1. Once more, as ρ→1, the parameter T1 will be converted to t. Substituting Eq. (48) into Eq. (45), and then separating the real and imaginary parts, one gets: , 22 3      −= dt d (49) .0 2 32 3 =+− F  (50) Forced Nonlinear Oscillator in a Fractal Space 15 The combination of Eqs. (49) and (50) yields (2 3 ) . 6 12 d F dt         + − = − (51) Eq. (51) is a linear differential equation in the function 𝜓. Its solution may be written as ,)32( 6 1 32 )(                 −+− − = tExp F t      (52) where the initial condition: 𝜓 (0)=0 is used. It follows that the stability occurs in such a way that (2 2 3 ) / 0,  − −  (53) where Eq. (44) is used. The instability criterion shown in the inequality (53) is independent of the parameter F. This happens because the parameter F lies away from the argument of the exponential function, as seen in Eq. (52). On the other hand, as Ω=1, the frequency of the exciting force has no implication in the stability criterion. This occurs because the natural frequency of the given DVdP is the unity. In the latter case, regardless of the sign of the parameter β, the stability criterion yields µ>0, which gives the same condition given in the non-resonance case. It is appropriate to corroborate this stability criterion numerically. Therefore, in what follows, two stability profiles are plotted. In the following figures, the white areas symbolize the unstable regions, whereas the dark areas stand for the stable regions. Fig. 10 represents the stability diagram, at which the parameter µ is plotted versus the parameter Ω, to show the influence of the parameter β on the stability picture. As seen, the plane is partitioned into white and dark regions. This figure is interpreted as follows: for positive values of the parameter µ, the increase of the parameter β matches larger stable Fig. 10 The stability criterion as given in the inequality (53), with variation of β 16 J.-H. HE, G. M. MOATIMID, M. H. ZEKRY regions. On the other hand, at the negative values of the parameter µ, the increase of the parameter β meets smaller stable regions. Therefore, the parameter β plays a dual role in the stability configuration. It should be noticed that this influence does not appear in the preceding analysis. Fig. 11 plots the stability profile, at which the parameter µ is plotted versus the parameter β to indicate the influence of the parameter Ω on the stability picture. As seen, the plane is partitioned into white and dark regions. This figure illustrates that for positive values of the parameter µ, the increase of the parameter β meets larger stable regions. On the other hand, at negative values of the parameter µ, one finds that the increase of the parameter β matches smaller stable regions. Once more, the parameter Ω plays a dual role in the stability diagram. This is the first mechanism for the parameter Ω in the entire analysis of the current paper. Fig. 11 The stability criterion as given in the inequality (53), with variation of Ω. 4. DISCUSSION AND CONCLUDING REMARKS Generally, a conservative vibration system must be free of the damping forces, and its amplitude will be fixed as t tends to infinity. As pointed out in Ref. [44], the inertia force in a fractal space can be approximately expressed as 2 2 2 2 = +(1 ) d z d z dz p p dd d    − (54) where p is a function of the fractal dimension. In our paper, a weak damping effect is considered. Its amplitude changes extremely slow, and we assume its value not to change for simplicity. Forced Nonlinear Oscillator in a Fractal Space 17 By the inverse two-scale transform [45-47], x=z√c/b, t=√a τα, the above results are converted to the fractal cases. Because of the great importance of the motivation of the DVdP in extensive applications in several areas in the sciences, many researchers have examined this problem. Really, the DVdP is one of the significant mathematical models for demonstrating a dynamical system having a single unstable fixed point and a single stable limit cycle. The forced DVdP does not have global analytical first integrals. This classical technique of He-Laplace does not have any authority to avoid the presence of the secular terms that have been assimilated during obtaining the approximate solution. Accordingly, the resulting approximate solution has an increased amplitude with the increase of time. On the other hand, the bounded solution is accomplished by means of nonlinear frequency analysis. Actually, this approach is easy, straightforward, powerful, and promising. It may be applied to analyze different types of highly nonlinear problems. Through this method, a dispersion relation that guarantees the periodic nature of the solution is reached. Once more, the HPM is utilized to realize an approximate solution of this characteristic equation. Simultaneously, stability restrictions are addressed. The periodic solution, in light of the second method, is graphed. By means of the linearized stability analysis, the stability criterion of the autonomous system is carried out. Finally, the multiple time scales together with the Homotopy concept are used to govern those stability criteria of the whole system. The latter analyses include the resonance as well as the non-resonance cases. It is concluded that the multiple scales method is one the most powerful mathematical tool in analyzing the nonlinear oscillation system that arises in physical applications and engineering. Overall, the concluding remarks of the present work may be summarized as follows: ▪ From the mathematical analysis of the nonlinear expanded frequency, one finds ▪ To maintain the presence of the cubic stiffness parameter µ in the stability criterion, a modification of the initial conditions has been made. ▪ The criterion of the analytical bounded approximate solution is given as 36β−µ2>0. ▪ In accordance with the multiple scales method, the following outcomes are reached: ▪ The stability criterion in the non-resonance case, regardless of the sign of the parameter β, requires that µ>0. ▪ The stability criterion, in the resonance case, needs µ(2Ω−2−3β) β <0. ▪ Because the natural frequency of the considered system is the unity, the parameter Ω has no implication on the stability configuration. ▪ The numerical calculations of the previous stability criterion show a dual role of both the parameters β and Ω. These mechanisms depend mainly on the sign of the parameter µ. ▪ Unfortunately, Through the multiple time scale method, the stability criteria are independent of the influence of the parameter F. 18 J.-H. HE, G. M. MOATIMID, M. H. ZEKRY REFERENCES 1. Zuo, Y.-T., 2021, A gecko-like fractal receptor of a three-dimensional printing technology: A fractal oscillator, Journal of Mathematical Chemistry, 59(4), pp. 1-10. 2. Elias-Zuniga, A., Palacios-Pineda, L.M., Jimenez-Cedeno, I.H., 2021, Equivalent power-form representation of the fractal Toda oscillator, Fractals, 29(1), 2150019. 3. Tian, D., Ain, Q.T., Anjum, N., He, C-H, Cheng, B., 2020, Fractal N/MEMS: From the pull-in instability to pull-in stability, Fractals, 29(2), 2150030. 4. Tian, D., He, C.H., 2021, A fractal micro-electromechanical system and its pull-in stability, Journal of Low Frequency Noise, Vibration and Active Control, 40(3), pp. 1380-1386. 5. He, J.H., Na, Q., Chun-Hui, H., Khaled, G., 2022, Fast identification of the pull-in voltage of a nano/micro-electromechanical system, Journal of Low Frequency Noise, Vibration and Active Control, doi:10.1177/14613484211068252 6. Liu, C. X., 2021, Periodic solution of Fractal phi-4 equation, Thermal Science, 25(2B), pp. 1345-1350. 7. Wang, K.L., 2021, He's frequency formulation for fractal nonlinear oscillator arising in a microgravity space, Numerical Methods for Partial Differential Equations, 37(2), pp. 1374-1384. 8. Wang, K.L., Wei C. F., 2020, A powerful and simple frequency formula to nonlinear fractal oscillators, Journal of Low Frequency Noise, Vibration and Active Control, 40(3), pp. 1373-1379. 9. Elias-Zuniga, A., Martinez-Romero, O., Trejo, D.O., Palacios-Pineda, L.M., 2021, An efficient ancient Chinese algorithm to investigate the dynamics response of a fractal microgravity forced oscillator, Accepted in Fractals. 10. He, C.H., Liu, C., Gepreel K. A., 2021, Low frequency property of a fractal vibration model for a concrete beam, Accepted in Fractals. 11. He, C.H., Liu. C., 2022, A modified frequency-amplitude formulation for fractal vibration systems, Fractals, doi: 10.1142/S0218348X22500463 12. He, C.H., Liu, C., He, J. H., Shirazi, A. H., Mohammad-Sedighi, H., 2021, Passive Atmospheric water harvesting utilizing an ancient Chinese ink slab and its possible applications in modern architecture, Facta Universitatis-Series Mechanical Engineering, 19(2), pp. 229-239. 13. Feng, G.Q., 2021, He’s frequency formula to fractal undamped Duffing equation, Journal of Low Frequency Noise, Vibration and Active Control, 40(4), pp. 1671-1676. 14. Elías-Zúñiga, A., Palacios-Pineda, L.M., Jiménez-Cedeño, I.H., Martínez-Romero, Olvera-Terjo, D., 2021, Analytical solution of the fractal cubi-quintic Duffing equation, Fractals, 29(4), 2150080. 15. Elias-Zuniga, A., Palacios-Pineda, L.M., Martinez-Romero, O., and Trejo, D.O., 2021, Dynamics response of the forced Fangzhu devices for water collection from air, Accepted in Fractals. 16. He, J.H., El-Dib, Y.O., 2021, A tutorial introduction to the two-scale fractal calculus and its application to the fractal Zhiber-Shabat Oscillator, Fractals, doi: 10.1142/S0218348X21502686. 17. He, C.H., Liu, S.H., Liu, C., Mohammad-Sedighi, H., 2021, A novel bond stress-slip model for 3-D printed concretes, Discrete and Continuous dynamical Systems S, doi: 10.3934/dcdss.2021161 18. Wang, K.L., 2022, Exact solitary wave solution for fractal shallow water wave model by He’s variational method, Modern Physics Letters B, 2150602, doi: 10.1142/S0217984921506028. 19. Wang, K.L., 2022, Solitary wave solution of nonlinear Bogoyavlenskii system by variational analysis method, International Journal of Modern Physics B, doi: 10.1142/S0217979222500151. 20. Zuo, Y.T., 2021, Effect of SiC particles on viscosity of 3D print paste: A fractal rheological model and experimental verification, Thermal Science, 25(3), pp. 2403-2407. 21. Zuo, Y.T., Liu, H.J., 2021, A fractal rheological model for sic paste using a fractal derivative, Journal of Applied and Computational Mechanics, 7, pp. 13-18. 22. Zuo, Y.-T., Liu, H.-J., 2021, Fractal approach to mechanical and electrical properties of graphene/sic composites, Facta Universitatis-Series Mechanical Engineering, 19(2), pp. 271-284. 23. Li, X.J., Liu, Z., He, J.H., 2020, A fractal two-phase flow model for the fiber motion in a polymer filling process, Fractals, 28(4), 2050093. 24. Kang-Le Wang, 2022, New variational theory for coupled nonlinear fractal Schrodinger system, International Journal of Numerical Methods for Heat and Fluid Flow, 32(2), pp. 589-597. 25. Liu, X.Y., Liu, Y.P., Wu, Z.W., 2021, Fractal calculus for modeling electrochemical capacitors under dynamical cycling, Thermal Science, 25 (2), pp. 1317-1320. 26. Alghassab, M., Mahmoud, A., Zohdy, M.Z., 2017, Nonlinear control of chaotic forced Duffing and Van der Pol oscillators, International Journal of Modern Nonlinear Theory and Application, 6(1), pp. 26-37. 27. Zhihong, Z., Shaopu, Y., 2015, Application of van der Pol–Duffing oscillator in weak signal detection, Computers & Electrical Engineering, 41, pp. 1-8. https://www.worldscientific.com/worldscinet/fractals https://www.worldscientific.com/toc/fractals/29/01 https://doi.org/10.1142/S0218348X22500463 https://www.worldscientific.com/toc/fractals/29/04 https://www.sciencedirect.com/science/journal/00457906 Forced Nonlinear Oscillator in a Fractal Space 19 28. Gao, G., Feng, Z., 2010, First integral for the Duffing Van der Pol type oscillator, UAB Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, pp. 123–133. 29. Motsa, S.S., Sibanda, P., 2012, A note on the solutions of the Van der Pol and Duffing equations using a linearization method, Mathematical Problems in Engineering, 2012, 693453. 30. Cui, J., Liang, J., Lin, Z., 2016, Stability analysis for periodic solutions of the Van der Pol-Duffing forced oscillator, Physica Scripta, 91, 015201. 31. Moatimid, G.M., 2020, Stability Analysis of a Parametric Duffing Oscillator, Journal of Engineering Mechanics, 146(5), 05020001. 32. Wei, X., Xiang-dong, W., Guang, M., Tong, F., 2002, Principal response of Van der Pol-Duffing oscillator under combined deterministic and random parametric excitation, Applied Mathematics and Mechanics, 23, pp. 299-310. 33. Dao, N.V., Dinh, N.V., Chi, T. K., 2007, Van der Pol oscillator under parametric and forced excitation, Ukrainian Mathematical Journal, 59(2), pp. 215-228. 34. Han, X., Bi, Q., Chun, Z., Yue, Y., 2014, Study of mixed-mode oscillations in parametrically excited Van der Pol, Nonlinear Dynamics, 77(4), pp. 1285-1296. 35. Kovacic, I., 2013, Harmonically excited generalized Van der Pol oscillators: Entertainment phenomenon, Meccanica, 48, pp. 2415-2425. 36. Shen, Y., Li, H., Yang, S., Peng, M., Han, Y., 2020, Primary and subharmonic simultaneous resonance of fractional-order Duffing oscillator, Nonlinear Dynamics, 102, pp. 1485-1497. 37. Wu, H., Zeng, X., Liu, Y., Lai, J., 2021, Analysis of Harmonically Forced Duffing Oscillator with Time Delay State Feedback by Incremental Harmonic Balance Method, Journal of Vibration Engineering & Technologies, 9, pp. 1239-1251. 38. He, J.‐H., El‐Dib, Y.O., Mady, A.A., 2021, Homotopy Perturbation Method for the Fractal Toda Oscillator, Fractal Fractional, 5(3), pp. 1-8. 39. Gómez-Aguilar, J.F., Yépez-Martínez, H., Calderón-Ramón, C., Cruz-Orduña, I., Escobar-Jiménez, R.F., Olivares-Peregrino, V. H., 2015, Modeling of a mass-spring-damper system by fractional derivatives with and without a singular kernel, Entropy, 17, pp. 6289-6303. 40. Yang, C-Y., Zhang, Y-D., Yang, X-J., 2016, Exact solutions for the differential equations in fractal heat transfer, Thermal Science, 20(3), pp. 747-750. 41. Elías-Zúñiga, A., Martínez-Romero, O., Trejo, D.O. Palacios-Pineda, L.M., 2021, Exact steady-state solution of fractals damped, and forced systems, Results in Physics, 28, 104580. 42. Ghaleb, A.F., Abou-Dina, M.S., Moatimid, G.M., Zekry, M.H., 2021, Analytic approximate solutions of the cubic–quintic Duffing–van der Pol equation with two-external periodic forcing terms: Stability analysis, Mathematics and Computers in Simulation, 180, pp. 129–151. 43. Nayfeh, A.H., 1973, Perturbation Methods, John Wiley-Interscience, New York, 425 p. 44. He, C.H., Liu, C., He, J.H., Mohammad-Sedighi, H., Shokri, A., Gepreel K.A., 2021, A fractal model for the internal temperature response of a porous concrete, Applied and Computational Mathematics, 20, pp. 1871-1875. 45. Liu, X.Y., Liu, Y.P., Wu, Z.W., 2021, Optimization of a fractal electrode-level charge transport model, Thermal Science, 25(3), pp. 2213-2220. 46. Liu, Y.P., Wang, C.C., Li, S.J., 2021, A fractal Langmuir kinetic equation and its solution structure, Thermal Science, 25(2), pp. 1351-1354. 47. Wang, K.J., 2021, Generalized variational principle and periodic wave solution to the modified equal width- Burgers equation in nonlinear dispersion media, Physics Letters A, 419(17), doi: 10.1016/j.physleta. 2021.127723. 20 J.-H. HE, G. M. MOATIMID, M. H. ZEKRY APPENDIX The function x2(t) that appears in Eq. (16) may be written as follows: 2 2 2 21 2 2 2 222 2 2 2 2 2 2 2 8 1536 ( ) sin 36 3 ( ) sin 3 (3 5 ) sin 5 sin 1 ( 1) 768 (2 ) 768 ( 2 ) sin( (2 )) sin( ( 2 )) ( ) cos ( 1)( 1) (3 ) ( 1)( 1) ( 3 )1 ( ) 3072 8 36 3 ( ) ( F F x t t t t t F F t t x t t x x F t               + + + + − −  −   −  −  +  − +  +  − − +  + +  −  − +   −  − − +  = − + + −  − 2 2 2 2 2 2 , 4608 2304 ) cos 3 cos Cos( (2 )) 1 ( 1) ( 1) ( 3 2 ) 2304 Cos( ( 2 )) ( 1) ( 3 2 ) F F t t t F t                    −  + +  +   −  + − +  +       − +    − − +  +   2 4 2 2 2 21 2 2 2 2 24 96 (3 78 11 ) ( ) (813 216 ) 48 (13 ) (73 72 ) 1 ( 9)( 1) F F x t t t t      −  +  = − − + + + +    −  −  −  , and 2 4 2 4 2 2 2 22 2 2 2 288(( 183 10 ) 4 (9 10 ) ) ( ) 44 864 48 (24 ) ( 9)( 1) t x t t t F       − −  +  + −  +  = − + − + +  −  − . The function x2(t) that appears in Eq. (33) may be written as follows: 2 2 2 21 22 23 24 4 2 2 2 2 2 2 2 2 1 cos sin cos 3 sin 3 (8 cos 5 (3 5 ) sin 5 ) 1 96 ( 3 cos sin ) 96 ( 3 cos( ( 2 )) ( 2 ) sin( ( 2 ))) ( ) sin 96 ( ) (3 )( ) 96 ( 3 cos( ( 2 )) ( 2 X t X t X t X t t t F t t F t t x t t F t                            + + + + + − + −  +   −  − +  −  − = + + + −  −  −  −  + +  + 2 2 2 , ) sin( ( 2 ))) (3 )( ) t                    +    +  −   2 2 2 2 2 4 2 2 4 21 4 2 2 2 2 2 ( (9 )( ) 36(111 18 ) ) (9 )( ) X F            = − − + +  −  − − ,                         −+− −−− +−−− −− = )18111(12 ))(9(29 ))(9(42 ))(9( 1 4224 22222 2222222 22222422     F X 2 2 23 4 2 2 2 9 (4 ( )) ( ) X F       = − − − , and 2 2 2 2 2 2 24 4 2 2 2 9 ( ( ) 2 ( 2 ( ))) ( ) X F        = − + − + − − .