Abstract

Free and self-excited vibrations of conservative oscillators with polynomial nonlinearity are considered. Mathematical model of the system is a second-order differential equation with a nonlinearity of polynomial type, whose terms are of integer and/or noninteger order. For the case when only one nonlinear term exists, the exact analytical solution of the differential equation is determined as a cosine-Ateb function. Based on this solution, the asymptotic averaging procedure for solving the perturbed strong non-linear differential equation is developed. The method does not require the existence of the small parameter in the system. Special attention is given to the case when the dominant term is a linear one and to the case when it is of any non-linear order. Exact solutions of the averaged differential equations of motion are obtained. The obtained results are compared with “exact” numerical solutions and previously obtained analytical approximate ones. Advantages and disadvantages of the suggested procedure are discussed.

1. Introduction

We study free and self-excited vibrations of a conservative nonlinear oscillator having mathematical model of a so-called second-order autonomous differential equation̈𝑥+𝑅(𝑥)=0,(1.1) with initial conditions𝑥(0)=𝐴,̇𝑥(0)=0,(1.2) where 𝑅(𝑥) is an odd function of 𝑥, that is, 𝑅(𝑥)=𝑅(𝑥) and 𝐴 is the initial amplitude of vibration. Besides, the function 𝑅(𝑥) is of polynomial type, whose terms are of integer and/or noninteger order𝑅(𝑥)=𝛼𝑐2𝛼𝑥|𝑥|𝛼1,(1.3) where 𝛼+ and 𝑐2𝛼 are constants. Differential equation̈𝑥+𝛼𝑐2𝛼𝑥|𝑥|𝛼1=0(1.4) has a wide application as describes many oscillatory systems and processes which are considered in [16]. To find the exact analytical solution of (1.4) with initial conditions, (1.2) is usually impossible and various approximate analytical solving methods are developed. One of the possible solving methods is the variable scale procedure developed by Bondar (see [7, Section 1, pages 9–12] and the references therein). For the case when the nonlinearity is described with only one binomial̈𝑥+𝑐2𝛼𝑥|𝑥|𝛼1=0,(1.5) the period of vibration is exactly calculated [6] and the approximate solution of (1.5) is usually assumed in the form of a cosine trigonometric function [816]. In spite of the simplicity of the solution form, and the fact that the amplitude and the frequency of vibration are exactly obtained, difference between the approximate solution and the exact numerical one may be significant. It was the reason that the exact solution of (1.5) in the form of cosine Ateb periodic function being an inversion of the incomplete Beta function [17] is introduced [18, 19]. The Ateb functions are more complex than the trigonometric functions and, at the moment, are not familiar to engineers and technicians who have to use them. Usually, simple approximation of Ateb functions by elementary functions is proposed (see [2022]). In this paper, the exact Ateb periodic functions will be applied and widely discussed. In Section 2, the knowledge about Ateb function is extended and adopted for application in engineering. Based on the exact solution of (1.5), in Sections 3 and 4, the approximate averaging procedure for solving (1.4) with (1.2) is developed. Main property of the suggested method is that the obtained solution is the perturbed version of the solution of the differential equation with the “dominant” term which may be linear or non-linear of integer or non-integer order. The method does not require the existence of the small parameter. Accuracy of the procedure is proved on several examples.

2. Qualitative and Quantitative Analysis of a Pure Nonlinear Oscillator

First integral of (1.5) reads as follows:̇𝑥2+2𝑐2𝛼𝛼+1|𝑥|𝛼+1=𝐶.(2.1) Being both addends on the left nonnegative, the associated phase paths represent a generalized Lamé superellipse in the 𝑥---̇𝑥 phase plane. There is a single equilibrium point 𝑥=̇𝑥=0 such that is a centre; therefore, the solutions of (1.5) are periodic in time.

Simple application of the initial boundary conditions (1.2) leads to the equation [6, page 1067, Equation (19)]̇𝑥2+2𝑐2𝛼𝛼+1|𝑥|𝛼+1=2𝑐2𝛼𝐴𝛼+1𝛼+1,𝐴>0,(2.2) such that it has to be solved in the sequel. It is not hard to see that the previous equation can be rewritten intȯ𝑥=±𝐾1|𝑥|𝐴𝛼+1,𝐾=2||𝑐𝛼||𝐴(𝛼+1)/2.𝛼+1(2.3) Now, we are ready to solve (1.5) analytically. Let us choose the positive right-hand-side expression in (1.5). It will be𝐿=d𝑥1(|𝑥|/𝐴)𝛼+1=𝐾𝑡+𝐶,(2.4) where 𝐶 denotes the integration constant. Expanding the integrand in 𝐿 into a binomial series and then integrating termwise, we conclude𝐿=𝑛=0(1)𝑛12𝑛|𝑥|𝐴(𝛼+1)𝑛d𝑥=𝐴𝑛=012𝑛1(𝛼+1)𝑛+1|𝑥|𝐴(𝛼+1)𝑛+1.(2.5) Employing the Pochhammer symbol (𝑎)𝑛=𝑎(𝑎+1)(𝑎+𝑛1),𝑛 mutatis mutandis 12𝑛=(1)𝑛(1/2)𝑛,1𝑛!=(𝛼+1)𝑛+11/(𝛼+1)=𝑛+1/(𝛼+1)(1/(𝛼+1))𝑛((𝛼+2)/(𝛼+1))𝑛;(2.6) hence𝐿=|𝑥|𝑛=0(1/2)𝑛(1/(𝛼+1))𝑛((𝛼+2)/(𝛼+1))𝑛𝑛!|𝑥|𝐴(𝛼+1)𝑛=|𝑥|2𝐹112,1𝛼+1𝛼+2𝛼+1|𝑥|𝐴𝛼+1.(2.7) Now, by the well-known formula (see [23, 24]) that connects the Gaussian hypergeometric 2𝐹1 and the incomplete Beta function B𝑧,2𝐹1𝑎,1𝑏𝑎+1𝑧=𝑎𝑧𝑎B𝑧(𝑎,𝑏),|𝑧|<1.(2.8) Letting1𝑎=1𝛼+1,𝑏=2,(2.9) we getB(|𝑥|/𝐴)𝛼+11,1𝛼+12=||𝑐2(𝛼+1)𝛼||𝐴(𝛼1)/2𝑡+𝐶.(2.10) Finally, the initial condition 𝑥(0)=𝐴 clearly gives(𝛼+1)|𝑥|𝐴2𝐹11,1𝛼+1211+𝛼+1|𝑥|𝐴𝛼+1=B(|𝑥|/𝐴)𝛼+11,1𝛼+121=B,1𝛼+12+||𝑐2(𝛼+1)𝛼||𝐴(𝛼1)/2𝑡.(2.11) This relation will be our main tool in determining the explicit solution of Cauchy problem (1.5) with (1.2). Namely, depending on 𝛼 we conclude different types of solutions, expressed (in a closed form) via trigonometric function, Jacobi elliptic integral, and finally periodic Ateb and hyperbolic Ateb functions (see Rosenberg’s introductionary paper [25] and Senik’s exhaustive article [26], in which the author constructed the periodic Ateb and hyperbolic Ateb functions). Rosenberg reported [25, page 39, Equation (9)] that the period of vibration for 𝑥(𝑡) is𝑇𝛼=22𝐴(1𝛼)/2𝑐2𝛼B1(𝛼+1),1𝛼+12,(2.12) and, in general, 𝑥(𝑡)=𝑥(𝑡+𝑇) is asynchronous. (In the linear case, 𝑇1=2𝜋𝑐11).

2.1. The Case 𝛼=1

Equation (1.5) reduces the autonomous equation̈𝑥+𝑐21𝑥=0.(2.13) By virtue of the initial condition 𝑥(0)=𝐴, formula (2.11) reduces to|𝑥|2𝐹112,1232𝑥𝐴2𝑥=𝐴arcsin𝐴𝜋=𝐴2+||𝑐1||𝑡;(2.14) therefore,𝜋𝑥(𝑡)=±𝐴sin2+||𝑐1||𝑡||𝑐=𝐴cos1||𝑡.(2.15)

2.2. The Case 𝛼=3

Let us recall [27, 28] that2𝐹114,1254𝑧=14𝑧𝐹arcsin4=1𝑧14𝑧sn1(𝑧1),|𝑧|<1,(2.16) where 𝐹(𝑧𝑚)=𝑧=am0(𝑧𝑚)d𝑡1𝑚sin2𝑡,sn(𝑧𝑚)=sinam(𝑧𝑚),(2.17) and 𝐹,am,sn denote the incomplete elliptic integral of the first kind, the Jacobi amplitude, the Jacobi elliptic sn function and sn1 stands for the inverse Jacobi elliptic sn functions, respectively.

Thus, for 𝛼=3 we have|𝑥|2𝐹114,1254𝑥𝐴4=𝐴sn1|𝑥|𝐴|||||𝑐1=𝐶+3||𝐴22𝑡.(2.18) Since the initial condition 𝑥(0)=𝐴, we conclude𝐶=𝐴sn1(11)=𝐴𝐾(1),(2.19) where 𝐾(𝑚)=𝐹(𝜋/2𝑚) denotes the celebrated complete elliptic integral of the first kind. The Jacobi elliptic sn(𝑧𝑚) has period 𝜔=4𝐾(𝑚); for now14𝐾(1)=B4,12=Γ2(1/4)2𝜋5.244116(2.20) is the period of the function sn(z1). Employing the quarter-period transformation formula for the Jacobi amplitude [29]:𝜋am(𝐾(m)zm)=2am|||𝑚1𝑚𝑧𝑚1,𝑚1(2.21) for 𝑚=1, one deduces by (2.18) that||𝑐𝑥(𝑡)=𝐴sn𝐾(1)+3||𝐴2𝑡|||||||𝑐1=𝐴sinam𝐾(1)+3||𝐴2𝑡|||||𝜋1=𝐴sin2||𝑐am𝐾(1)+3||𝐴2𝑡|||||||𝑐1=𝐴cosam3||||1𝐴𝑡2.(2.22) Thus||𝑐𝑥(𝑡)=𝐴cn3||||1𝐴𝑡2.(2.23) Here cn(𝑧𝑚)=cosam(𝑧𝑚) denotes the Jacobi elliptic cn function (The elliptic functions sn,cn were defined by Jacobi under the names sin am, cos am; the notation sn, cn was introduced by Gudermann and Glashier, as stated by Du Val [30]). Being cn(01/2)=1, formula (2.23) has the interpolatory property 𝑥(0)=𝐴. It is worth to say that Lyapunov in his classical paper [31] introduced the Jacobi elliptic functions (cn and sn) which are the special case of Ateb cosinus and Ateb sinus functions for 𝛼=3. The same functions are used for solving the third order nonlinear differential equation of Duffing type by Bravo Yuste and Diaz Bejarano [32] but also Chen and Cheung [33]. The in-built subroutine JacobiCN [𝑧,𝑚] in Mathematica enables the exact curve plotting of the Jacobi elliptic functions (see [30]).

2.3. The General Case 𝛼>0 Solved by Periodic Ateb Functions

In his classical paper [25], Rosenberg introduced the so-called periodic Ateb-functions concerning the problem of inversion of the half of the incomplete Beta-function1𝑧2B𝑧1(𝑎,𝑏)=200𝑧1𝑡𝑎1(1𝑡)𝑏1d𝑡.(2.24) Obviously, we are interested in case 𝑎=1/(𝛼+1),𝑏=1/2 (Senik [26] constructed inverses for the incomplete Beta-function B𝑧(𝑎,𝑏) for 𝑎=1/(𝑛+1), 𝑏=1/(𝑚+1), 𝑛,𝑚=(2𝜇+1)/(2𝜈+1), 𝜇,𝜈0). The Ateb (kind of inverse Beta) functions in the focus of our interest are solutions of the system of ODĖ𝑣𝑢𝛼2=0,̇𝑢+𝛼+1𝑣=0.(2.25) consult for example, Senik’s famous article [26, page 274, Equation (1.21)]. We write𝑣(𝑧)=sa(1,𝛼,𝑧),𝑢(𝑧)=ca(𝛼,1,𝑧).(2.26) It can be easily verified that the inverse of (1/2)B𝑧(1/2,1/(𝛼+1)) and 𝑣(𝑧) coincide on [(1/2)Π𝛼,(1/2)Π𝛼], whereΠ𝛼1=B,1𝛼+12(2.27) Having in mind the following set of properties:1sa(1,𝛼,𝑧)=sa(1,𝛼,𝑧),ca𝛼,1,2Π𝛼,±𝑥±sa1,𝛼,Π𝛼,±𝑧sa1,𝛼,2Π𝛼,𝑧(2.28) we see that sa(𝛼,1,𝑧) is an odd function of 𝑧; it is the so-called 2Π𝛼-periodic sine Ateb, that is, sa-function. Also there holds [26, page 273, Equation (1.15)]sa2(1,𝛼,𝑧)+ca𝛼+1(𝛼,1,𝑧)=1,(2.29) and cosine Ateb, that is ca(1,𝛼,𝑧) function, is even and 2Π𝛼-periodic having properties [26, page 273, Equations (1.18-19), (1.24)]:1ca(𝛼,1,𝑧)=ca(𝛼,1,𝑧),sa1,𝛼,2Π𝛼,±𝑧ca𝛼,1,Π𝛼,±𝑧ca𝛼,1,2Π𝛼.±𝑧(2.30) By these two sets of relations, we see that functions sa,ca are defined on the whole range of . For 𝛼=1 sa, ca one reduce to the sine and cosine functions. That means sa(1,1,𝑧)=sin𝑧,ca(1,1,𝑧)=cos𝑧.

Now, inverting the half of the incomplete Beta-function in (2.11):12B(|𝑥|/𝐴)𝛼+11,1𝛼+12=Π𝛼2+||𝑐𝛼+1𝛼||2𝐴(𝛼1)/2𝑡,(2.31) we clearly deduceΠ𝑥(𝑡)=𝐴sa1,𝛼,𝛼2+||𝑐𝛼+1𝛼||2𝐴(𝛼1)/2𝑡.(2.32) Having in mind the quarter-period expansion formula (2.30), we arrive at𝑥(𝑡)=𝐴ca𝛼,1,||𝑐𝛼+1𝛼||2𝐴(𝛼1)/2𝑡,𝑡.(2.33) By ca(𝛼,1,0)=1, we see that the initial condition 𝑥(0)=𝐴 is satisfied as well.

Moreover, we have to point out a restricting characteristics of Rosenberg’s and Senik’s inversion. Rosenberg pointed out the rule [25, page 40].

Exponents 𝑛=(𝛼+1)/2 and 1/𝑛 behave like odd integers,” while Senik’s restriction to some rational values of 𝛼 constitutes the set of permitted 𝛼-values (see Table 1).

Approximate methods by numerically obtained evaluations of Ateb-functions have been performed by Droniuk, Nazarkevich, and coworkers for information protection, see [34, 35] and the references therein. Further study on Ateb-function integral was realized by Drogomirecka in [36].

Being cosine Ateb ca(𝑛,𝑚,𝑧) even 2Π𝛼-periodic function, it is a perfect candidate for a cosine Fourier series expansion. Let us mention that finite Fourier series approximation has been discussed in [34, Section 4], while in [35, pages 122, 124] the sine Ateb sa(𝑛,𝑚,𝑧) has been approximated by its sine Fourier series. Applying the there exposed method to ca(1,𝛼,𝑧), we conclude thatca(𝛼,1,𝑧)=𝑛=1𝑎𝑛cos𝜋𝑛𝑧Π𝛼(2.34) since obviously 𝑎0=0 and𝑎𝑛=2Π𝛼Π𝛼0ca(𝛼,1,𝑧)cos𝜋𝑛𝑧Π𝛼2d𝑧=Π𝛼Π𝛼0cos𝜋𝑛𝑧Π𝛼1𝑧d𝑢1𝑢2𝛼/(𝛼+1)d𝑧.(2.35) We compute the 𝑎𝑛 values numerically according to the prescribed accuracy. Of course, it is enough to compute ca(𝛼,1,𝑧) for 𝑧[0,Π𝛼/2], another values we calculate by means of formula (2.30), see also [37].

Another model in approximating Ateb-functions is the Taylor series expansion, which corresponds to the investigations by Gricik and Nazarkevich [38].

2.4. The Case 𝛼>0 Solved by Hyperbolic Ateb Function

This kind functions arise equation (1.5) one transform into ̈𝑥𝑐2𝛼𝑥|𝑥|𝛼1̈𝑥+i𝑐𝛼2𝑥|𝑥|𝛼1=0(2.36) that is the coefficient 𝑐𝛼 in (1.5) becomes purely imaginary, see [25, page 47, Equation (43)], [26, Section 2].

3. Differential Equation with a Nonlinear Dominant Term Application of Gendelman and Vakakis [18] Approximation

Let us consider the mathematical model of the oscillator in the form̈𝑥+𝑐2𝛼𝑥|𝑥|𝛼1=𝑗Λ𝑐2𝑗𝑥|𝑥|𝑗1,(3.1) under initial conditions (1.2), where 𝑐2𝛼𝑥|𝑥|𝛼1 is the so-called “dominant term” in polynomial, (1.3) which satisfies the relation𝑐2𝛼𝐴||𝐴||𝛼1𝑐2𝑗𝐴||𝐴||𝑗1,(3.2) for all 𝑗Λ and Λ+ is some set of indices. For condition (3.2) the terms 𝑐2𝑗𝑥|𝑥|𝑗1 are significantly smaller than the dominant one 𝑐2𝛼𝑥|𝑥|𝛼1 and for all values of 𝑥 we have𝑐2𝛼𝑥|𝑥|𝛼1𝑐2𝑗𝑥|𝑥|𝑗1.(3.3) The right-hand-side expression in (3.1) is the so-called perturbation term which represents a small value. Due to (3.3), we assume the solution of (3.1) as the perturbed version of the solution of the differential equation (1.5).

Remark 3.1. Usually, in the literature it is stated that if the right-hand-side coefficients are small (𝑐2𝑗1) in comparison to 𝑐2𝛼 the solution of the differential equation (3.1) differs from the solution of the differential equation (1.5) just for a bit. Unfortunately, this conclusion is not valid for strong non-linear differential equations where the initial conditions have a significant influence. Then, the limitation (3.2) has to be introduced and satisfied.
Having in mind the solution (2.33) of the differential equation (1.5) with only one nonlinear term, the perturbed solution of (3.1) is assumed as ||𝑐𝑥=𝑎(𝑡)ca(𝛼,1,𝜓(𝑡)),̇𝜓(𝑡)=𝛼||𝛼+12[]𝑎(𝑡)(𝛼1)/2+̇𝜃(𝑡),(3.4) where the amplitude 𝑎 and the phase angle 𝜓 are time dependent. Let us recall that [26, page 274, Equation (1.21)] d2d𝑧ca(𝛼,1,𝑧)=d𝛼+1sa(1,𝛼,𝑧),d𝑧sa(1,𝛼,𝑧)=ca𝛼(𝛼,1,𝑧).(3.5) So, the first time derivative of such 𝑥 is ̇𝑥=2||𝑐𝛼||𝑎𝛼+1(𝛼+1)/2sa(1,𝛼,𝜓(𝑡)),(3.6) with ̇𝜃̇𝑎ca(𝛼,1,𝜓)2𝑎𝛼+1sa(1,𝛼,𝜓)=0.(3.7) Substituting the second derivative of ̈𝑥 and the function 𝑥 (3.4) into ODE (3.1), we get ̈𝑥+𝑐2𝛼𝑥|𝑥|𝛼1||𝑐=𝛼||𝛼+12𝑎(𝛼1)/2||𝑐sa(1,𝛼,𝜓)̇𝑎𝛼||2𝑎𝛼+1(𝛼+1)/2ca𝛼(̇𝜃𝛼,1,𝜓)=𝑗Λ𝑐2𝑗𝑎𝑗ca𝑗(1,𝛼,𝜓)=Ξ𝛼,(3.8) such that it becomes a first-order differential equation sa(1,𝛼,𝜓)̇𝑎+2𝑎𝛼+1ca𝛼̇(𝛼,1,𝜓)𝜃=2𝑎𝛼+1(1𝛼)/2Ξ𝛼||𝑐𝛼||.(3.9) Coupled with (3.7), this system of ODEs in ̇𝜃̇𝑎, results in ̇𝑎=2𝑎(1𝛼)/2||𝑐𝛼+1𝛼||sa(1,𝛼,𝜓)𝑗Λ𝑐2𝑗𝑎𝑗ca𝑗(𝛼,1,𝜓)(3.10)̇𝜃=𝛼+1𝑎(1𝛼)/22||𝑐𝛼||𝑗Λ𝑐2𝑗𝑎𝑗1ca𝑗+1(𝛼,1,𝜓).(3.11) Now, the solution of the averaged ODE (3.10) reads 𝑎(𝑡)=𝐴.(3.12) Indeed, bearing on mind that any 𝑇-periodic (𝑓(𝑡+𝑇)=𝑓(𝑡)) integrable function 𝑓 possesses the property 𝛾𝑇+𝛾𝑓(𝑧)𝑓𝛽(𝑧)d𝑧=0,𝛽>0,𝛾.(3.13) By letting 𝑓(𝑧)=ca(𝛼,1,𝑧),𝑇=2Π𝛼, the averaging procedure immediately gives 𝑎̇𝑎=(1𝛼)/22||𝑐(𝛼+1)𝛼||Π𝛼𝑗Λ𝑐2𝑗𝑎𝑗2Π𝛼0sa(1,𝛼,𝜓)ca𝑗(𝛼,1,𝜓)d𝜓=0.(3.14) Hence, by the initial condition 𝑥(0)=𝐴, we deduce the claim (3.12). Now, it remains to solve (3.11). Because ca(𝛼,1,𝑧) is even and periodic, by averaging the right-hand-side expression in (3.11), we conclude ̇𝜃=𝛼+1𝑎(1𝛼)/222||𝑐𝛼||𝜋𝑗Λ𝑐2𝑗𝑎𝑗12Π𝛼0ca𝑗+1(𝛼,1,𝜓)d𝜓.(3.15) To find the value of the inner-most integral, we apply a result by Drogomirecka [36, page 108, Teorema], which reads as follows: 2Π𝛼0sa𝑝(𝑛,𝑚,𝑧)ca𝑞1(𝑚,𝑛,𝑧)d𝑧=21+(1)𝑝+(1)𝑞+(1)𝑝+𝑞B𝑝+1,𝑛+1𝑞+1𝑚+1(3.16) for all 𝑝,𝑞{𝑟/𝑠𝑟,𝑠=2𝑘1,𝑘}, that is, the powers are rational, having integer numerator and odd positive denominator. Hence ̇𝜃=𝛼+1𝑎(1𝛼)/22||𝑐𝛼||𝜋𝑗Λodd𝑐2𝑗𝑎𝑗1B𝑗+2,1𝛼+12,(3.17) so, according to 𝑎(𝑡)𝐴, we get 𝜃(𝑡)=𝛼+12||𝑐𝛼||𝜋𝑗Λodd𝑐2𝑗𝐴𝑗(𝛼+1)/2B𝑗+2,1𝛼+12𝑡.(3.18) Thus, the exact averaged perturbed solution of the Cauchy problem (3.1) with the initial conditions (1.2) reads 𝑥(𝑡)=𝐴ca𝛼,1,𝛼+1𝐴(𝛼1)/22||𝑐𝛼||+𝑐||𝑐𝛼||𝜋𝑡,(3.19) with the correction constant 𝑐=𝑗Λodd𝑐2𝑗𝐴𝑗𝛼B𝑗+2,1𝛼+12.(3.20) Correction constant depends not only on the coefficients 𝑐2𝑗 and𝑐2𝛼, but also on the initial amplitude and the order of the binoms in polynomial.

Example 3.2. Let us consider the Cauchy problem ̈𝑥+𝑥|𝑥|2/3=0.01𝑥,𝑥(0)=0.1,̇𝑥(0)=0,(3.21) that is, 𝑐1/3=1,𝐴=0.1,𝑐1=1. According to (2.33), the exact solution of the homogeneous ODE ̈𝑥+𝑥|𝑥|2/3=0(3.22) becomes 1𝑥(𝑡)=0.1ca3,1,233110𝑡0.1ca3,1,1.75909𝑡,(3.23) such that is 2Π1/3=2B(3/4,1/2)4.79256-periodic function.
The exact solution, achieved by the averaged perturbation method, according to (3.19) is: 1𝑥(𝑡)=0.1ca3,1,233101B(9/4,1/2)𝜋108/3𝑡10.1ca3,1,1.75758𝑡.(3.24)
Amplitude of vibration remains constant and is equal to initial one, but the period of vibration has to be corrected. Frequency of perturbed vibration depends on the initial amplitude, coefficients of non-linearity, and order of non-linear terms.

Further discussion in evaluating the integral in (3.15) by an approximative method leads to the following considerations:̇𝜃=𝛼+1𝑎(1𝛼)/222||𝑐𝛼||𝜋𝑗Λ𝑐2𝑗𝑎𝑗12Π𝛼0ca𝑗+1=(𝛼,1,𝜓)d𝜓𝛼+1𝑎(1𝛼)/22||𝑐𝛼||𝜋𝑗Λ𝑐2𝑗𝑎𝑗11+(1)𝑗+1(1/2)Π𝛼0ca𝑗+1=(𝛼,1,𝜓)d𝜓2(𝛼+1)𝑎(1𝛼)/2||𝑐𝛼||𝜋𝑗Λodd𝑐2𝑗𝑎𝑗1(1/2)Π𝛼0ca𝑗+1(𝛼,1,𝜓)d𝜓.(3.25) Now, we point out the paper by Gendelman and Vakakis [18, page 1537, Equation (6b)], where the approximation formula2ca(𝛼,1,𝑧)1𝛼+1lncosh𝛼+12𝑧1,|𝑧|2Π𝛼,(3.26) has been used in approximative treatment of the two degree-of-freedom damped nonlinear system governed by the equation of motion. The approximation (3.26) is periodically continued to the whole . Approximation (3.26) and the corresponding period are accurate up to 𝒪((𝛼+1)2) which will be enough sharp for our considerations (see [18]). So, we apply first the Gendelman-Vakakis approximation in the integrand, then one takes the first two terms in the Maclaurin series for 1lncosh𝑧1(1/2)𝑧2,𝑧. This procedure results in=(1/2)Π𝛼0ca𝑗+1(𝛼,1,𝜓)d𝜓(1/2)Π𝛼021𝛼+1lncosh𝛼+12𝜓𝑗+12d𝜓𝛼+1(Π𝛼/4)0𝛼+11𝜓2𝑗+1=Πd𝜓𝛼22𝐹1123,(𝑗+1)2𝛼+1Π162𝛼,(3.27) such that it holds for all |𝛼|𝛼0, where 𝛼0 is the unique positive nil of the equationB1,1𝛼+12=42𝛼+1sa1,𝛼,𝛼+1=1.(3.28) Let us remark here that 𝛼0=0.83073245123719 Therefore, we havė𝜃𝛼+1Π𝛼2||𝑐𝛼||𝜋𝑗Λodd𝑐2𝑗𝐴2𝑗((𝛼+1)/2)𝐹1123,(𝑗+1)2𝛼+1Π162𝛼.(3.29) Also, in the case that 𝑗0, the hypergeometric term one reduces to the polynomial 𝑃𝑗+1(((𝛼+1)/16)Π2𝛼) of degree (𝑗+1).

Finally, we conclude that the approximate solution ̃𝑥(𝑡) of the averaged perturbed Cauchy problem (3.1) with the initial conditions (1.2) reads as follows:𝑥(𝑡)̃𝑥(𝑡)=𝐴ca𝛼,1,𝛼+1𝐴(𝛼1)/22||𝑐𝛼||+Π𝛼̃𝑐||𝑐𝛼||𝜋𝑡,(3.30) with the correction constant̃𝑐=𝑗Λodd𝑐2𝑗𝐴2𝑗𝛼𝐹1123,(𝑗+1)2𝛼+1Π162𝛼.(3.31)

Example 3.3. Let us study the previous example. Because 𝛼=1/3<𝛼0, it is in the permitted range for 𝛼, so the correction constant will be ̃𝑐=1022/3𝐹1124,3321B12234,120.798431,(3.32) and the averaged, approximated system of ODEs is ̇̇𝑎=0,𝜃=0.012310̃𝑐B33𝜋4,12.(3.33) Therefore, by (3.19) we have 1𝑥(𝑡)̃𝑥(𝑡)=0.1ca3.,1,0.993910𝑡(3.34)

The discrepancy (Δ=43.5%) between the coefficients of the arguments in exact (3.24) and the approximate (3.34) solutions is explained by the rough triple approximation for the cosine-Ateb function appearing in (3.26). To make precise the Gendelman-Vakakis approximation method, they propose to take further terms in the estimate (3.26), compare [18]. However, this significantly complicated our exposition.

4. Differential Equation with a Linear Dominant Term

Let us consider the case when the dominant term is linear and the second order nonlinear Cauchy problem is̈𝑥+𝑐21𝑥=𝑗Λ𝑐2𝑗𝑥|𝑥|𝑗1,𝑥(0)=𝐴,̇𝑥(0)=0.(4.1) Let us mention that this system one reduces to the system (3.1) associated with the linear case 𝛼=1.

As (2.15) suggests, the general solution of (4.1) can be taken in the form𝑐𝑥=𝑎cos1𝑡+𝜃,(4.2) where 𝑎,𝜃 stand for the related integration constants. Assuming that the perturbation term is a small one (𝑐2𝑗|𝐴|𝑗1𝑐21), it is supposed that the vibrations of the perturbed system are close to that of the linear system. So, we are looking for an approximate solution of (4.1) in the form of (2.15), but with time variable amplitude and phase,𝑥=𝑎(𝑡)cos𝜓(𝑡)=𝑎cos𝜓,(4.3) where 𝜓=𝜓(𝑡)=𝑐1𝑡+𝜃(𝑡). First time derivative of 𝑥 iṡ𝑥=𝑎𝑐1sin𝜓(4.4) foṙ̇𝑎cos𝜓𝑎𝜃sin𝜓=0.(4.5) Substituting (4.3) and the time derivative of (4.4) into initial ODE (4.1), we obtaiṅ1̇𝑎sin𝜓+𝑎𝜃cos𝜓=𝑐1𝑗Λ𝑐2𝑗𝑎𝑗cos𝑗𝜓=Φ.(4.6) Hence, second-order nonlinear ODE (4.1) is rewritten into a system of two coupled first order ODEs (4.5) and (4.6). Solving this system with respect to ̇𝜃 and ̇𝑎, we geṫ𝜓=𝑐1+Φ𝑎cos𝜓,̇𝑎=Φsin𝜓.(4.7) To solve this system it is not an easy task. At this point we introduce the averaging procedure over the period of vibration equal to 𝑇1=2𝜋/𝑐1. Assuming that {𝑗}>1, for the initial conditions (1.2) the averaged differential equation arė1̇𝑎=0𝜃=2𝜋𝑗Λ𝑐2𝑗𝑎𝑗12𝜋/𝑐10cos𝑗+1𝑐11𝑡+𝜃d𝑡=2𝑐1𝜋𝑗Λ𝑐2𝑗𝑎𝑗11+(1)𝑗+1Γ((𝑗+2)/2).Γ((𝑗+3)/2)(4.8) So, for the initial conditions (1.2) the solutions of the averaged differential equation are𝜓𝑐𝑎(𝑡)=𝐴=const.(𝑡)=1+1𝑐1𝜋𝑗Λodd𝑐2𝑗𝐴𝑗1Γ((𝑗+2)/2)Γ((𝑗+3)/2)𝑡,(4.9) where Λodd denotes the set of indices which behave like odd integers and 𝜓(𝑡) takes its final form by virtue of reasons similar to Rosenberg-Senik’s rules.

Having in mind the initial conditions (1.2), the solution of (4.1) in the first approximation is obtained as𝑐𝑥(𝑡)=𝐴cos1+𝑐𝑐1𝑡,(4.10) with the correction term𝑐=1𝜋𝑗Λ𝑐2𝑗𝐴𝑗1Γ((𝑗+2)/2),Γ((𝑗+3)/2)(4.11) which depends on the coefficients of non-linearity, initial amplitude, and the order of non-linearity.

Example 4.1. Consider DE with initial conditions ̈𝑥+𝑥=𝑐25/3𝑥|𝑥|2/3,𝑥(0)=0.1,̇𝑥(0)=0.(4.12) For the case when 𝑐25/3=0.01, the related first approximation solution is 𝑥𝐴=0.1cos1+0.010.12/3Γ(11/6)𝑡𝜋Γ(7/3)0.1cos(1.00096𝑡).(4.13) And for the case when 𝑐25/3=1, the related first approximation solution is 𝑥𝐴1=0.1cos1+10.12/3Γ(11/6)𝑡𝜋Γ(7/3)0.1cos(1.096𝑡).(4.14) In Figure 1, the approximate analytic solutions 𝑥𝐴 (4.13) and 𝑥𝐴1 (4.14) and the exact numerical ones 𝑥𝑁 and 𝑥𝑁1 for (4.12) with 𝑐25/3=0.01 and 𝑐25/3=1, respectively, are plotted. Numerical solutions are obtained by applying of the Runge-Kutta procedure. It can be seen that the difference is negligible. In spite of the fact that the coefficient of the non-linear term 𝑐25/3 is of the same order as of the linear term 𝑐21, the averaging solving procedure gives very accurate results due to the fact that the relation (3.2) is satisfied, that is,0.1>0.15/30.021544.

Example 4.2. For the oscillator and initial conditions ̈𝑥+𝑥=0.5𝑥|𝑥|2/3,𝑥(0)=3.0,̇𝑥(0)=0.(4.15) the analytical approximate solution is 𝑥𝐴=3cos1+0.532/3Γ(11/6)𝑡𝜋Γ(7/3)3cos(1.46343𝑡).(4.16) In Figure 2, the analytical solution 𝑥𝐴 (4.16) and the numerically obtained solution 𝑥𝑁 of (4.15) are plotted. It can be seen that the difference is significant. Namely, the suggested criteria (3.2) are not satisfied: 3̸0.535/33.12. The criterion that the coefficient 𝑐21 is higher than 𝑐25/3 is not valid for application of the approximate solving procedures in the strong non-linear oscillators.

5. Conclusion

Due to previous consideration the following is concluded.(1)The periodical cosinus-Ateb function is the exact analytical solution of the second-order pure non-linear differential equation, which describes the vibration of the conservative oscillator. Using the Fourier series expansion of the cosinus-Ateb function and taking the first term, the approximate description of the motion with a cosinus periodical function is obtained. Such description of the motion is usually in application.(2)Dominant binomial in the polynomial which describes the non-linear property of the conservative oscillator is obtained using not only the coefficient of the term but also the initial amplitude of the system. For the initial case, the dominant term has the highest value in comparison to other terms in polynomial. The higher the difference between the dominant and other terms, the higher the accuracy of vibration calculation based on the oscillator with dominant term.(3)The approximate solving procedures have to be based on the solution of the differential equation with the dominant term. The additional terms in the differential equation give the perturbation of that solution. Motion of the oscillator is the perturbed version of vibration of the pure non-linear system. The smaller the perturbation terms, the higher the accuracy of the approximate solution.(4)The approximate averaging method is suitable for application if the perturbation of the oscillatory motion is small and the criteria suggested in the paper (3.2) are satisfied. If the coefficient of perturbation is small, it does not mean that the method is applicable for the strong non-linear systems. The suggested solving method does not require the existence of the small parameter.(5)For the non-linear conservative oscillator, the averaging procedure gives the amplitude of vibration that corresponds to the initial one, and the frequency of vibration depends not only on the initial amplitude but also on the coefficients and orders of the nonlinearity.

Acknowledgments

The investigation is supported by the Ministry of Science and Technological Development of Serbia under Grants ON 174028 and IT 41007 and Provincial Secretariat for Science and Technological Development of Vojvodina under Grant 114-451-2094/2011.