# Oscillator with a Sum of Noninteger-Order Nonlinearities

Authors: L. Cveticanin; T. Pogány
Journal: Journal of Applied Mathematics
 (2012)

---

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


---

## Body

## 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(1.1)ẍ+R(x)=0,
with initial conditions(1.2)x(0)=A,ẋ(0)=0,
where R(x) is an odd function of x, that is, R(x)=-R(-x) and A is the initial amplitude of vibration. Besides, the function R(x) is of polynomial type, whose terms are of integer and/or noninteger order(1.3)R(x)=∑αcα2x|x|α-1,
where α∈ℝ+ and cα2 are constants. Differential equation(1.4)ẍ+∑αcα2x|x|α-1=0
has a wide application as describes many oscillatory systems and processes which are considered in [1–6]. 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(1.5)ẍ+cα2x|x|α-1=0,
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 [8–16]. 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 [20–22]). 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.1)ẋ2+2cα2α+1|x|α+1=C.
Being both addends on the left nonnegative, the associated phase paths represent a generalized Lamé superellipse in the x—ẋ phase plane. There is a single equilibrium point x=ẋ=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)ẋ2+2cα2α+1|x|α+1=2cα2α+1Aα+1,A>0,
such that it has to be solved in the sequel. It is not hard to see that the previous equation can be rewritten into(2.3)ẋ=±K1-(|x|A)α+1,K=2|cα|A(α+1)/2α+1.
Now, we are ready to solve (1.5) analytically. Let us choose the positive right-hand-side expression in (1.5). It will be(2.4)L=∫dx1-(|x|/A)α+1=Kt+C,
where C denotes the integration constant. Expanding the integrand in L into a binomial series and then integrating termwise, we conclude(2.5)L=∫∑n=0∞(-1)n(-12n)(|x|A)(α+1)ndx=A∑n=0∞(-12n)1(α+1)n+1(|x|A)(α+1)n+1.
Employing the Pochhammer symbol (a)n=a(a+1)···(a+n-1),n∈ℕmutatis mutandis(2.6)(-12n)=(-1)n(1/2)nn!,1(α+1)n+1=1/(α+1)n+1/(α+1)=(1/(α+1))n((α+2)/(α+1))n;
hence(2.7)L=|x|∑n=0∞(1/2)n(1/(α+1))n((α+2)/(α+1))nn!(|x|A)(α+1)n=|x|2F1[12,1α+1α+2α+1(|x|A)α+1].
Now, by the well-known formula (see [23, 24]) that connects the Gaussian hypergeometric 2F1 and the incomplete Beta function Bz,(2.8)2F1[a,1-ba+1z]=azaBz(a,b),|z|<1.
Letting(2.9)a=1α+1,b=12,
we get(2.10)B(|x|/A)α+1(1α+1,12)=2(α+1)|cα|A(α-1)/2t+C.
Finally, the initial condition x(0)=A clearly gives(2.11)(α+1)|x|A2F1[1α+1,121+1α+1(|x|A)α+1]=B(|x|/A)α+1(1α+1,12)=B(1α+1,12)+2(α+1)|cα|A(α-1)/2t.
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 x(t) is(2.12)Tα=22A(1-α)/2cα2(α+1)B(1α+1,12),
and, in general, x(t)=x(t+T) is asynchronous. (In the linear case, T1=2πc1-1).


### 2.1. The Caseα=1

Equation (1.5) reduces the autonomous equation(2.13)ẍ+c12x=0.
By virtue of the initial condition x(0)=A, formula (2.11) reduces to(2.14)|x|2F1[12,1232(xA)2]=Aarcsin(xA)=A(π2+|c1|t);
therefore,(2.15)x(t)=±Asin⁡(π2+|c1|t)=Acos⁡(|c1|t).


### 2.2. The Caseα=3

Let us recall [27, 28] that(2.16)2F1[14,1254z]=1z4F(arcsinz4∣-1)=1z4sn-1(z∣-1),|z|<1,
where
(2.17)F(z∣m)=∫0z=am(z∣m)dt1-msin⁡2t,sn(z∣m)=sin⁡⁡am(z∣m),
and F,am,sn denote the incomplete elliptic integral of the first kind, the Jacobi amplitude, the Jacobi elliptic sn function and sn-1 stands for the inverse Jacobi elliptic sn functions, respectively.Thus, forα=3 we have(2.18)|x|2F1[14,1254(xA)4]=A⋅sn-1(|x|A|-1)=C+|c3|A22t.
Since the initial condition x(0)=A, we conclude(2.19)C=A⋅sn-1(1∣-1)=A⋅K(-1),
where K(m)=F(π/2∣m) denotes the celebrated complete elliptic integral of the first kind. The Jacobi elliptic sn(z∣m) has period ω=4K(m); for now(2.20)4K(-1)=B(14,12)=Γ2(1/4)2π≈5.244116
is the period of the function sn(z∣-1). Employing the quarter-period transformation formula for the Jacobi amplitude [29]:(2.21)am(K(m)-z∣m)=π2-am(1-mz|mm-1),m≤1
for m=-1, one deduces by (2.18) that(2.22)x(t)=A⋅sn(K(-1)+|c3|A2t|-1)=A⋅sin⁡am(K(-1)+|c3|A2t|-1)=A⋅sin⁡{π2-am(K(-1)+|c3|A2t|-1)}=A⋅cos⁡am(|c3|At|12).
Thus(2.23)x(t)=A⋅cn(|c3|At|12).
Here cn(z∣m)=cosam(z∣m) denotes the Jacobi ellipticcnfunction(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(0∣1/2)=1, formula (2.23) has the interpolatory property x(0)=A. 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 [z,m] 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-function(2.24)z↦12Bz(a,b)=12∫00≤z≤1ta-1(1-t)b-1dt.
Obviously, we are interested in case a=1/(α+1),b=1/2  (Senik [26] constructed inverses for the incomplete Beta-function Bz(a,b) for a=1/(n+1),  b=1/(m+1), n,m=(2μ+1)/(2ν+1),  μ,ν∈ℕ0). The Ateb (kind of inverse Beta) functions in the focus of our interest are solutions of the system of ODE(2.25)v̇-uα=0,u̇+2α+1v=0.
consult for example, Senik’s famous article [26, page 274, Equation (1.21)]. We write(2.26)v(z)=sa(1,α,z),u(z)=ca(α,1,z).
It can be easily verified that the inverse of (1/2)Bz(1/2,1/(α+1)) and v(z) coincide on [-(1/2)Πα,(1/2)Πα], where(2.27)Πα:=B(1α+1,12)
Having in mind the following set of properties:(2.28)sa(1,α,z)={-sa(1,α,-z),∓ca(α,1,12Πα±x),±sa(1,α,Πα±z),∓sa(1,α,2Πα∓z),
we see that sa(α,1,z) is an odd function of z∈ℝ; it is the so-called 2Πα-periodic sine Ateb, that is, sa-function. Also there holds [26, page 273, Equation (1.15)](2.29)sa2(1,α,z)+caα+1(α,1,z)=1,
and cosine Ateb, that is ca(1,α,z) function, is even and 2Πα-periodic having properties [26, page 273, Equations  (1.18-19), (1.24)]:(2.30)ca(α,1,z)={ca(α,1,-z),sa(1,α,12Πα±z),-ca(α,1,Πα±z),ca(α,1,2Πα±z).
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,z)=sinz,ca(1,1,z)=cosz.Now, inverting the half of the incomplete Beta-function in (2.11):(2.31)12B(|x|/A)α+1(1α+1,12)=Πα2+α+1|cα|2A(α-1)/2t,
we clearly deduce(2.32)x(t)=A⋅sa(1,α,Πα2+α+1|cα|2A(α-1)/2t).
Having in mind the quarter-period expansion formula (2.30), we arrive at(2.33)x(t)=A⋅ca(α,1,α+1|cα|2A(α-1)/2t),t∈R.
By ca(α,1,0)=1, we see that the initial condition x(0)=A 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].“Exponentsn=(α+1)/2and1/nbehave like odd integers,” while Senik’s restriction to some rational values of α constitutes the set of permitted α-values (see Table 1).Table 1
The values ofα=(2ν+1)/(2μ+1),μ,ν∈ℕ0, by Senik’s traces.
μ∖ν01234567…01357911131511/315/37/3311/313/3521/53/517/59/511/513/5331/73/75/719/711/713/715/741/91/35/97/9111/913/95/351/113/115/117/119/11113/1115/1161/133/135/137/139/1311/13115/1371/151/51/37/153/511/1513/151⋮⋱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 Atebca(n,m,z) 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(n,m,z) has been approximated by its sine Fourier series. Applying the there exposed method to ca(1,α,z), we conclude that(2.34)ca(α,1,z)=∑n=1∞ancos⁡πnzΠα
since obviously a0=0 and(2.35)an=2Πα∫0Παca(α,1,z)⋅cos⁡πnzΠαdz=2Πα∫0Παcos⁡πnzΠα{∫z1du¯(1-u¯2)α/(α+1)}dz.
We compute the an values numerically according to the prescribed accuracy. Of course, it is enough to compute ca(α,1,z) for z∈[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.36)ẍ-cα2x|x|α-1≡ẍ+(icα)2x|x|α-1=0
that is the coefficient cα in (1.5) becomes purely imaginary, see [25, page 47, Equation (43)], [26, Section 2].


## 2.1. The Caseα=1

Equation (1.5) reduces the autonomous equation(2.13)ẍ+c12x=0.
By virtue of the initial condition x(0)=A, formula (2.11) reduces to(2.14)|x|2F1[12,1232(xA)2]=Aarcsin(xA)=A(π2+|c1|t);
therefore,(2.15)x(t)=±Asin⁡(π2+|c1|t)=Acos⁡(|c1|t).


## 2.2. The Caseα=3

Let us recall [27, 28] that(2.16)2F1[14,1254z]=1z4F(arcsinz4∣-1)=1z4sn-1(z∣-1),|z|<1,
where
(2.17)F(z∣m)=∫0z=am(z∣m)dt1-msin⁡2t,sn(z∣m)=sin⁡⁡am(z∣m),
and F,am,sn denote the incomplete elliptic integral of the first kind, the Jacobi amplitude, the Jacobi elliptic sn function and sn-1 stands for the inverse Jacobi elliptic sn functions, respectively.Thus, forα=3 we have(2.18)|x|2F1[14,1254(xA)4]=A⋅sn-1(|x|A|-1)=C+|c3|A22t.
Since the initial condition x(0)=A, we conclude(2.19)C=A⋅sn-1(1∣-1)=A⋅K(-1),
where K(m)=F(π/2∣m) denotes the celebrated complete elliptic integral of the first kind. The Jacobi elliptic sn(z∣m) has period ω=4K(m); for now(2.20)4K(-1)=B(14,12)=Γ2(1/4)2π≈5.244116
is the period of the function sn(z∣-1). Employing the quarter-period transformation formula for the Jacobi amplitude [29]:(2.21)am(K(m)-z∣m)=π2-am(1-mz|mm-1),m≤1
for m=-1, one deduces by (2.18) that(2.22)x(t)=A⋅sn(K(-1)+|c3|A2t|-1)=A⋅sin⁡am(K(-1)+|c3|A2t|-1)=A⋅sin⁡{π2-am(K(-1)+|c3|A2t|-1)}=A⋅cos⁡am(|c3|At|12).
Thus(2.23)x(t)=A⋅cn(|c3|At|12).
Here cn(z∣m)=cosam(z∣m) denotes the Jacobi ellipticcnfunction(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(0∣1/2)=1, formula (2.23) has the interpolatory property x(0)=A. 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 [z,m] 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-function(2.24)z↦12Bz(a,b)=12∫00≤z≤1ta-1(1-t)b-1dt.
Obviously, we are interested in case a=1/(α+1),b=1/2  (Senik [26] constructed inverses for the incomplete Beta-function Bz(a,b) for a=1/(n+1),  b=1/(m+1), n,m=(2μ+1)/(2ν+1),  μ,ν∈ℕ0). The Ateb (kind of inverse Beta) functions in the focus of our interest are solutions of the system of ODE(2.25)v̇-uα=0,u̇+2α+1v=0.
consult for example, Senik’s famous article [26, page 274, Equation (1.21)]. We write(2.26)v(z)=sa(1,α,z),u(z)=ca(α,1,z).
It can be easily verified that the inverse of (1/2)Bz(1/2,1/(α+1)) and v(z) coincide on [-(1/2)Πα,(1/2)Πα], where(2.27)Πα:=B(1α+1,12)
Having in mind the following set of properties:(2.28)sa(1,α,z)={-sa(1,α,-z),∓ca(α,1,12Πα±x),±sa(1,α,Πα±z),∓sa(1,α,2Πα∓z),
we see that sa(α,1,z) is an odd function of z∈ℝ; it is the so-called 2Πα-periodic sine Ateb, that is, sa-function. Also there holds [26, page 273, Equation (1.15)](2.29)sa2(1,α,z)+caα+1(α,1,z)=1,
and cosine Ateb, that is ca(1,α,z) function, is even and 2Πα-periodic having properties [26, page 273, Equations  (1.18-19), (1.24)]:(2.30)ca(α,1,z)={ca(α,1,-z),sa(1,α,12Πα±z),-ca(α,1,Πα±z),ca(α,1,2Πα±z).
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,z)=sinz,ca(1,1,z)=cosz.Now, inverting the half of the incomplete Beta-function in (2.11):(2.31)12B(|x|/A)α+1(1α+1,12)=Πα2+α+1|cα|2A(α-1)/2t,
we clearly deduce(2.32)x(t)=A⋅sa(1,α,Πα2+α+1|cα|2A(α-1)/2t).
Having in mind the quarter-period expansion formula (2.30), we arrive at(2.33)x(t)=A⋅ca(α,1,α+1|cα|2A(α-1)/2t),t∈R.
By ca(α,1,0)=1, we see that the initial condition x(0)=A 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].“Exponentsn=(α+1)/2and1/nbehave like odd integers,” while Senik’s restriction to some rational values of α constitutes the set of permitted α-values (see Table 1).Table 1
The values ofα=(2ν+1)/(2μ+1),μ,ν∈ℕ0, by Senik’s traces.
μ∖ν01234567…01357911131511/315/37/3311/313/3521/53/517/59/511/513/5331/73/75/719/711/713/715/741/91/35/97/9111/913/95/351/113/115/117/119/11113/1115/1161/133/135/137/139/1311/13115/1371/151/51/37/153/511/1513/151⋮⋱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 Atebca(n,m,z) 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(n,m,z) has been approximated by its sine Fourier series. Applying the there exposed method to ca(1,α,z), we conclude that(2.34)ca(α,1,z)=∑n=1∞ancos⁡πnzΠα
since obviously a0=0 and(2.35)an=2Πα∫0Παca(α,1,z)⋅cos⁡πnzΠαdz=2Πα∫0Παcos⁡πnzΠα{∫z1du¯(1-u¯2)α/(α+1)}dz.
We compute the an values numerically according to the prescribed accuracy. Of course, it is enough to compute ca(α,1,z) for z∈[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.36)ẍ-cα2x|x|α-1≡ẍ+(icα)2x|x|α-1=0
that is the coefficient cα 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(3.1)ẍ+cα2x|x|α-1=-∑j∈Λcj2x|x|j-1,
under initial conditions (1.2), where cα2x|x|α-1 is the so-called “dominant term” in polynomial, (1.3) which satisfies the relation(3.2)cα2A|A|α-1≫cj2A|A|j-1,
for all j∈Λ and Λ⊂ℝ+ is some set of indices. For condition (3.2) the terms cj2x|x|j-1 are significantly smaller than the dominant one cα2x|x|α-1 and for all values of x we have(3.3)cα2x|x|α-1≫cj2x|x|j-1.
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 (cj2≪1) in comparison to cα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
(3.4)x=a(t)⋅ca(α,1,ψ(t)),ψ̇(t):=|cα|α+12[a(t)](α-1)/2+θ̇(t),
where the amplitude a and the phase angle ψ are time dependent. Let us recall that [26, page 274, Equation (1.21)]
(3.5)ddzca(α,1,z)=-2α+1sa(1,α,z),ddzsa(1,α,z)=caα(α,1,z).
So, the first time derivative of such x is
(3.6)ẋ=-2|cα|α+1a(α+1)/2sa(1,α,ψ(t)),
with
(3.7)ȧca(α,1,ψ)-2aθ̇α+1sa(1,α,ψ)=0.
Substituting the second derivative of ẍ and the function x (3.4) into ODE (3.1), we get
(3.8)ẍ+cα2x|x|α-1=-|cα|α+12a(α-1)/2sa(1,α,ψ)⋅ȧ-|cα|2α+1a(α+1)/2caα(α,1,ψ)⋅θ̇=-∑j∈Λcj2ajcaj(1,α,ψ)=:-Ξα,
such that it becomes a first-order differential equation
(3.9)sa(1,α,ψ)⋅ȧ+2aα+1caα(α,1,ψ)⋅θ̇=2α+1a(1-α)/2Ξα|cα|.
Coupled with (3.7), this system of ODEs in ȧ,θ̇ results in
(3.10)ȧ=2a(1-α)/2α+1|cα|sa(1,α,ψ)∑j∈Λcj2ajcaj(α,1,ψ)(3.11)θ̇=α+1a(1-α)/22|cα|∑j∈Λcj2aj-1caj+1(α,1,ψ).
Now, the solution of the averaged ODE (3.10) reads
(3.12)a(t)=A.
Indeed, bearing on mind that any T-periodic (f(t+T)=f(t)) integrable function f possesses the property
(3.13)∫γT+γf′(z)fβ(z)dz=0,β>0,γ∈R.
By letting f(z)=ca(α,1,z),T=2Πα, the averaging procedure immediately gives
(3.14)ȧ=a(1-α)/22(α+1)|cα|Πα∑j∈Λcj2aj∫02Παsa(1,α,ψ)caj(α,1,ψ)dψ=0.
Hence, by the initial condition x(0)=A, we deduce the claim (3.12). Now, it remains to solve (3.11). Because ca(α,1,z) is even and periodic, by averaging the right-hand-side expression in (3.11), we conclude
(3.15)θ̇=α+1a(1-α)/222|cα|π∑j∈Λcj2aj-1∫02Παcaj+1(α,1,ψ)dψ.
To find the value of the inner-most integral, we apply a result by Drogomirecka [36, page 108, Teorema], which reads as follows:
(3.16)∫02Παsap(n,m,z)caq(m,n,z)dz=12(1+(-1)p+(-1)q+(-1)p+q)B(p+1n+1,q+1m+1)
for all p,q∈{r/s:r∈ℤ,s=2k-1,k∈ℕ}, that is, the powers are rational, having integer numerator and odd positive denominator. Hence
(3.17)θ̇=α+1a(1-α)/22|cα|π∑j∈Λoddcj2aj-1B(j+2α+1,12),
so, according to a(t)≡A, we get
(3.18)θ(t)={α+12|cα|π∑j∈Λoddcj2Aj-(α+1)/2B(j+2α+1,12)}t.
Thus, the exact averaged perturbed solution of the Cauchy problem (3.1) with the initial conditions (1.2) reads
(3.19)x(t)=A⋅ca(α,1,α+1A(α-1)/22{|cα|+c*|cα|π}t),
with the correction constant
(3.20)c*=∑j∈Λoddcj2Aj-αB(j+2α+1,12).
Correction constant depends not only on the coefficients cj2 andcα2, but also on the initial amplitude and the order of the binoms in polynomial.Example 3.2.
Let us consider the Cauchy problem(3.21)ẍ+x|x|-2/3=0.01x,x(0)=0.1,ẋ(0)=0,
that is, c1/3=1,A=0.1,c1=-1. According to (2.33), the exact solution of the homogeneous ODE
(3.22)ẍ+x|x|-2/3=0
becomes
(3.23)x(t)=0.1⋅ca(13,1,23103t)≈0.1⋅ca(13,1,1.75909t),
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:
(3.24)x(t)=0.1⋅ca(13,1,23103{1-B(9/4,1/2)π108/3}t)≈0.1⋅ca(13,1,1.75758t).
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:(3.25)θ̇=α+1a(1-α)/222|cα|π∑j∈Λcj2aj-1∫02Παcaj+1(α,1,ψ)dψ=α+1a(1-α)/22|cα|π∑j∈Λcj2aj-1(1+(-1)j+1)∫0(1/2)Παcaj+1(α,1,ψ)dψ=2(α+1)a(1-α)/2|cα|π∑j∈Λoddcj2aj-1∫0(1/2)Παcaj+1(α,1,ψ)dψ.
Now, we point out the paper by Gendelman and Vakakis [18, page 1537, Equation (6b)], where the approximation formula(3.26)ca(α,1,z)≈1-2α+1ln⁡cosh⁡(α+12z),|z|≤12Πα,
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 1-lncoshz≈1-(1/2)z2,z∈ℝ. This procedure results in(3.27)I=∫0(1/2)Παcaj+1(α,1,ψ)dψ≈∫0(1/2)Πα{1-2α+1ln⁡cosh⁡(α+12ψ)}j+1dψ≈2α+1∫0(Πα/4)α+1(1-ψ2)j+1dψ=Πα22F1[12,-(j+1)32α+116Πα2],
such that it holds for all |α|≤α0, where α0 is the unique positive nil of the equation(3.28)B(1α+1,12)=4α+1⇔sa(1,α,2α+1)=1.
Let us remark here that α0=0.83073245123719… Therefore, we have(3.29)θ̇≈α+1Πα2|cα|π∑j∈Λoddcj2Aj-((α+1)/2)2F1[12,-(j+1)32α+116Πα2].
Also, in the case that j∈ℕ0, the hypergeometric term one reduces to the polynomial Pj+1(((α+1)/16)Πα2) of degree (j+1).Finally, we conclude that the approximate solutionx̃(t) of the averaged perturbed Cauchy problem (3.1) with the initial conditions (1.2) reads as follows:(3.30)x(t)≈x̃(t)=A⋅ca(α,1,α+1A(α-1)/22{|cα|+Παc̃|cα|π}t),
with the correction constant(3.31)c̃=∑j∈Λoddcj2Aj-α2F1[12,-(j+1)32α+116Πα2].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
(3.32)c̃=102/32F1[12,-4332112B2(34,12)]≈0.798431,
and the averaged, approximated system of ODEs is
(3.33)ȧ=0,θ̇=-0.012103c̃3πB(34,12).
Therefore, by (3.19) we have
(3.34)x(t)≈x̃(t)=0.1⋅ca(13,1,0.993910t).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(4.1)ẍ+c12x=-∑j∈Λcj2x|x|j-1,x(0)=A,ẋ(0)=0.
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(4.2)x=acos⁡(c1t+θ),
where a,θ stand for the related integration constants. Assuming that the perturbation term is a small one (cj2|A|j-1≪c12), 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,(4.3)x=a(t)cos⁡ψ(t)=:acos⁡ψ,
where ψ=ψ(t)=c1t+θ(t). First time derivative of x is(4.4)ẋ=-ac1sin⁡ψ
for(4.5)ȧcos⁡ψ-aθ̇sin⁡ψ=0.
Substituting (4.3) and the time derivative of (4.4) into initial ODE (4.1), we obtain(4.6)ȧsin⁡ψ+aθ̇cos⁡ψ=1c1∑j∈Λcj2ajcos⁡jψ=:Φ.
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 ȧ, we get(4.7)ψ̇=c1+Φa⋅cos⁡ψ,ȧ=Φ⋅sin⁡ψ.
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 T1=2π/c1. Assuming that ℜ{j}>-1, for the initial conditions (1.2) the averaged differential equation are(4.8)ȧ=0θ̇=12π∑j∈Λcj2aj-1∫02π/c1cos⁡j+1(c1t+θ)dt=12c1π∑j∈Λcj2aj-1(1+(-1)j+1)Γ((j+2)/2)Γ((j+3)/2).
So, for the initial conditions (1.2) the solutions of the averaged differential equation are(4.9)a(t)=A=const.ψ(t)=(c1+1c1π∑j∈Λoddcj2Aj-1Γ((j+2)/2)Γ((j+3)/2))t,
where Λodd denotes the set of indices which behave like odd integers and ψ(t) 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(4.10)x(t)=Acos⁡{(c1+c**c1)t},
with the correction term(4.11)c**=1π∑j∈Λcj2Aj-1Γ((j+2)/2)Γ((j+3)/2),
which depends on the coefficients of non-linearity, initial amplitude, and the order of non-linearity.Example 4.1.
Consider DE with initial conditions(4.12)ẍ+x=-c5/32x|x|2/3,x(0)=0.1,ẋ(0)=0.
For the case when c5/32=0.01, the related first approximation solution is
(4.13)xA=0.1cos⁡[(1+0.01⋅0.12/3Γ(11/6)πΓ(7/3))t]≈0.1cos⁡(1.00096t).
And for the case when c5/32=1, the related first approximation solution is
(4.14)xA1=0.1cos⁡[(1+1⋅0.12/3Γ(11/6)πΓ(7/3))t]≈0.1cos⁡(1.096t).
In Figure 1, the approximate analytic solutions xA (4.13) and xA1 (4.14) and the exact numerical ones xN and xN1 for (4.12) with c5/32=0.01 and c5/32=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 c5/32 is of the same order as of the linear term c12, 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/3≡0.021544.Figure 1
Time history diagrams obtained: (a) analyticallyxA-tand numerically xN-t for c5/32=0.01 and (b) analytically xA1-t and numerically xN1-t for c5/32=1.Example 4.2.
For the oscillator and initial conditions(4.15)ẍ+x=-0.5x|x|2/3,x(0)=3.0,ẋ(0)=0.
the analytical approximate solution is
(4.16)xA=3cos⁡[(1+0.5⋅32/3Γ(11/6)πΓ(7/3))t]≈3cos⁡(1.46343t).
In Figure 2, the analytical solution xA (4.16) and the numerically obtained solution xN 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.5·35/3≡3.12. The criterion that the coefficient c12 is higher than c5/32 is not valid for application of the approximate solving procedures in the strong non-linear oscillators.Figure 2
Time-history diagrams obtained analyticallyxA  and numerically xN.


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



---

*Source: 649050-2012-02-22.xml*