#  Chaos analysis in the un-perturbed system for Quintic Duffing equation

Zeraoulia Rafik111Zeraoulia  University Batna2
53, Route de Constantine. Fésdis, Batna 05078, Algeria
<EMAIL_ADDRESS>,khenchela, Lycee Mourri toufana Alvaro
humberto Salas  Group Fizmako Research
universidade national de colombia,Bogota
<EMAIL_ADDRESS>
((to be inserted by publisher))

###### Abstract

In this paper we have investigated the driven cubic-quintic Duffing equation.
Using the Melnikov analysis technique, we are able to predict the number of
limit cycles around the equilibrium and to develop a theoretical approach to
chaos suppression in damped driven systems. The analytical solution and the
bifurcation of limit cycles of the investigated system have been studied., We
discussed our nonlinear dynamics, namely, the cubic quintic Duffing equation
in the absence of noise and we may discusse it in the presence of noise in the
futur research, investigating for the Hamiltonian of the system (not
perturbed) in both cases is presented.

###### keywords:

Homoclinic orbits- chaos theory- duffing equation-Hamiltonian.

##  1 Introduction

The harmonically driven damped pendulum is often used as a simple example of a
chaotic system, the equation is just

$\ddot{\phi}+\frac{1}{q}\dot{\phi}+\sin\phi=A\cos(\omega t)$ (1)

As long as $A$ and $\omega$ are small it behaves like a driven harmonic
oscillator, and asymptotically settles into regular oscillations with a fixed
period. However, as $A$ (or $\omega$) are increased, with the rest of
parameters fixed, the system undergoes a cascade of period doubling
bifurcations leading to chaotic behavior, which then gives way to regular
oscillations again when it is increased further. For example, when $q=2$ and
$\omega=2/3$ the first period doubling (”symmetry breaking”) occurs at
$A\approx 1.07$ and the first chaos at $A\approx 1.08$. These rigorous results
seem to be obtained by numerical simulations. One can be actually interested
in situations where chaos does not occur H.W. Haslach [1982]. Are there known
rigorous conditions on $A,\omega$ and $q$ that put the system below the first
period doubling? However this question does not belong to the aim of this
paper but it would be very interesting to conclude somethings about chaotics
behaviors of some dynamics and to discover new ways to supress chaos in the
cubic-Quintic Duffing Equation which it is the aim of our research in this
paper. The use of Melnikov analysis (MA) techniques Melnikov [1963] has
allowed the development of a theoretical approach to chaos suppression in
damped driven systems, and involves adding periodic chaos-suppressing (CS)
excitations B.Palmero& F. Chacon, R. [2022]. This MA-based approach has been
shown to be reliable in suppressing chaos in a Duffing oscillator by a fine
choice of the shape of the external periodic excitation Melnikov [1963], a
generalized Duffing oscillator with fractional-order deflection Gilbert Lewis
& Frank Monasa [1982], coupled arrays of damped Alvaro Humberto Salas
[2022],Alvaro Humberto [2022], periodically forced, nonlinear oscillators , as
well as in starlike networks of dissipative nonlinear oscillators Cveticanin,
L [1993]. The Duffing equation (or Duffing oscillator) named after George
Duffing is a nonlinear second order differential equation used to model
certain damped and driven oscillators with a more complicated potential than
in simple harmonic motion (P. Holmes [1979]) The Duffing equation is an
example of a dynamical system that exhibits chaotic behaviour B.Palmero& F.
Chacon, R. [2022].The equation is given by :

$\ddot{x}+\delta\dot{x}+\rho x+\mu x^{3}=\lambda\cos(\omega t),\text{ }$ (2)

where the (unknown) function $x=x(t)$ is the displacement at time $t$. The
damping factor $\delta$ controls the size of the damping, the $\rho$ controls
the size of the stiffness and the $\mu$ controls the amount of nonlinearity in
the restoring force. If $\mu=0$, the Duffing equation describes a damped and
driven simple harmonic oscillator. The quantity $\lambda$ controls the
amplitude of the periodic driving force. If $\lambda=0$, we have a system
without driving force. The quantity $\omega$ controls the frequency of the
periodic driving force.G. Prathap & T., Varadan [1976]

In this paper the special case as the modified formula of (2) which is called
cubic-quintic Duffing equation A Elías & Zúñiga [2013] is considered :

$\ddot{x}-ax+bx^{3}+cx^{5}=\varepsilon(\gamma\cos\omega t\
-\delta\dot{x}).\text{ }$ (3)

The method of chaos control by delayed self-controlling feedback developed by
Pyragas (Pyragas [1992] ,Pyragas [1996], Pyragas [2001]) is applied for (3).
The cubic-quintic Duffing oscillator which is defined in (3) has been
investigated and make comparison with different theoretical approach to get
analytical solution in the absence of drivingEl-Dib, Yusry O., Elgazery,
Nasser S., Mady, Amal A. and Alyousef & Haifa A [2022],The cubic Duffing
equation (3) can as well be used to model the nonlinear spring-mass system
(Hassan Nayfa [1973], J. C. Amazigo [2011]) as well as the motion of a
classical particle in a double well potentialLo, C.C & Gupta, S.D. [1978] .
(3) with initial condition $x(0)=A,\dot{x}=0$ with
$\omega=(2k+1)\frac{\pi}{2},k\in\mathbb{Z}$ were proposed as a system by
Correig in A. M. Correig & M. Urquizu [2002] as a model of microseism time
series and have been used in M. O. Oyesanya [2008]to model the prediction of
earthquake occurrence. It was also used to model the transverse oscillation of
nonlinear beams in H. M. Sedighi, K. H. Shirazi & J. Zare [2012].

##  2 Analytical Solution

Let us consider the i.v.p :

$\ddot{x}-ax+bx^{3}+cx^{5}=0,x(0)=x_{0}\text{ and }x^{\prime}(0)=0.$ (4)

The exact solution may be written in the form

$x=x(t)=x_{0}\frac{\sqrt{1+\lambda+\mu}\cdot\text{cn}(\sqrt{\omega}t,m)}{\sqrt{1+\lambda\cdot\text{cn}^{2}(\sqrt{\omega}t,m)+\mu\cdot\text{cn}^{4}(\sqrt{\omega}t,m)}}.$
(5)

$\ddot{x}-ax+bx^{3}+cx^{5}=$
---
$\frac{x_{0}\sqrt{1+\lambda+\mu}\text{cn}}{\left(1+\lambda\text{cn}^{2}+\mu\text{cn}^{4}\right)^{5/2}}\left[\begin{array}[]{c}(-a-\omega+2m\omega-3\lambda\omega+3m\lambda\omega)+\\\
\\\
\left(-2a\lambda-2m\omega+2\lambda\omega-4m\lambda\omega-10\mu\omega+10m\mu\omega+bx_{0}^{2}+b\lambda
x_{0}^{2}+b\mu x_{0}^{2}\right)\text{cn}^{2}+\\\ \\\
\left(\begin{array}[]{c}-a\lambda^{2}-2a\mu+m\lambda\omega+10\mu\omega-20m\mu\omega-\lambda\mu\omega+m\lambda\mu\omega+b\lambda
x_{0}^{2}+\\\ \\\ b\lambda^{2}x_{0}^{2}+b\lambda\mu
x_{0}^{2}+cx_{0}^{4}+2c\lambda x_{0}^{4}+c\lambda^{2}x_{0}^{4}+2c\mu
x_{0}^{4}+2c\lambda\mu
x_{0}^{4}+c\mu^{2}x_{0}^{4}\end{array}\right)\text{cn}^{4}+\\\ \\\
-\mu\left(2a\lambda-10m\omega-2\lambda\omega+4m\lambda\omega-2\mu\omega+2m\mu\omega-
bx_{0}^{2}-b\lambda x_{0}^{2}-b\mu x_{0}^{2}\right)\text{cn}^{6}+\\\ \\\
-\mu(a\mu-3m\lambda\omega+\mu\omega-2m\mu\omega)\text{cn}^{8}.\end{array}\right]$
(6)

Equating to zero the coefficients of cnj ($j=0,2,4,6,8$) in the last
expression gives an algebraic system for determining the unknown constants
$\lambda,\mu$, $\omega$ and $m$. The solutions are :

$\begin{array}[]{ccc}\begin{array}[]{c}\omega=\frac{1}{12}\left(-12a+9bx_{0}^{2}+6cx_{0}^{4}+\sqrt{3}\sqrt{\Delta}\right)\\\
\end{array}&\begin{array}[]{c}m=\frac{x_{0}^{2}\left(3b+2cx_{0}^{2}\right)\left(-bx_{0}^{2}+2cx_{0}^{4}+\sqrt{3}\sqrt{\Delta}\right)-4a\left(4cx_{0}^{4}+\sqrt{3}\sqrt{\Delta}\right)}{4\left(6a-3bx_{0}^{2}-2cx_{0}^{4}\right)\left(a-x_{0}^{2}\left(b+cx_{0}^{2}\right)\right)}\\\
\end{array}&\begin{array}[]{c}\lambda=\frac{-3bx_{0}^{2}-6cx_{0}^{4}+\sqrt{3}\sqrt{\Delta}}{12\left(-a+bx_{0}^{2}+cx_{0}^{4}\right)}\\\
\end{array}\\\
\begin{array}[]{c}\omega=\frac{1}{12}\left(-12a+9bx_{0}^{2}+6cx_{0}^{4}-\sqrt{3}\sqrt{\Delta}\right)\\\
\end{array}&\begin{array}[]{c}m=\frac{4a\left(\sqrt{3}\sqrt{\Delta}-4cx_{0}^{4}\right)-x_{0}^{2}\left(3b+2cx_{0}^{2}\right)\left(bx_{0}^{2}-2cx_{0}^{4}+\sqrt{3}\sqrt{\Delta}\right)}{4\left(6a-3bx_{0}^{2}-2cx_{0}^{4}\right)\left(a-bx_{0}^{2}-cx_{0}^{4}\right)}\\\
\end{array}&\begin{array}[]{c}\lambda=\frac{3bx_{0}^{2}+6cx_{0}^{4}+\sqrt{3}\sqrt{\Delta}}{12\left(a-bx_{0}^{2}-cx_{0}^{4}\right)}\\\
\end{array}\\\
\Delta=x_{0}^{4}\left(16ac+3b^{2}-4bcx_{0}^{2}-4c^{2}x_{0}^{4}\right)>0.&\left(6a-3bx_{0}^{2}-2cx_{0}^{4}\right)\left(a-bx_{0}^{2}-cx_{0}^{4}\right)\neq
0.&\end{array}$ (7)
$\begin{array}[]{cc}\begin{array}[]{c}\omega=\frac{b(3\lambda+2)x_{0}^{2}(\lambda+\mu+1)-2a(\lambda(3\lambda+4)-5\mu+1)}{6\lambda(\lambda+1)+10\mu+2}.\\\
\end{array}&\begin{array}[]{c}m=\frac{2a(\lambda(3\lambda+2)-5\mu)-b(3\lambda+1)x_{0}^{2}(\lambda+\mu+1)}{2a(\lambda(3\lambda+4)-5\mu+1)-b(3\lambda+2)x_{0}^{2}(\lambda+\mu+1)}.\\\
\end{array}\\\
\begin{array}[]{c}\lambda=\frac{2\left(3b+2cx_{0}^{2}\right)\left(-12a+9bx_{0}^{2}+6cx_{0}^{4}\pm
2\sqrt{6}\sqrt{\delta}\right)}{3x_{0}^{2}\left(-16ac-3b^{2}+4bcx_{0}^{2}+4c^{2}x_{0}^{4}\right)}\\\
\\\
\delta=\left(6a-3bx_{0}^{2}-2cx_{0}^{4}\right)\left(a-bx_{0}^{2}-cx_{0}^{4}\right)>0.\end{array}&\begin{array}[]{c}\mu=\frac{96a^{2}+x_{0}^{4}\left(51b^{2}-112ac\right)-144abx_{0}^{2}+76bcx_{0}^{6}+28c^{2}x_{0}^{8}\pm
4\sqrt{6}\sqrt{\delta}\left(4a-3bx_{0}^{2}-2cx_{0}^{4}\right)}{x_{0}^{4}\left(16ac+3b^{2}-4bcx_{0}^{2}-4c^{2}x_{0}^{4}\right)}\\\
\\\ -16ac-3b^{2}+4bcx_{0}^{2}+4c^{2}x_{0}^{4}\neq 0.\end{array}\end{array}$
(8)

We have the following homoclinic orbits :

$x(t)=x_{0}\frac{\sqrt{1+\lambda}\text{tanh}(\sqrt{k}t)}{\sqrt{1+\lambda\text{tanh}^{2}(\sqrt{k}t)}}=\frac{A\text{tanh}(\sqrt{k}t)}{\sqrt{1+\lambda\text{tanh}^{2}(\sqrt{k}t)}}$,
$A=x_{0}\sqrt{1+\lambda}.$
---
for the choices
$k=\frac{1}{2}\left(-bx_{0}^{2}-2cx_{0}^{4}\right),\lambda=-\frac{2cx_{0}^{2}}{3\left(b+2cx_{0}^{2}\right)},x_{0}=\sqrt{\frac{-b\pm\sqrt{b^{2}+4ac}}{2c}}.$
$b^{2}+4ac>0.$
(9)
$x(t)=x_{0}\frac{\sqrt{1+\lambda}\text{sech}(\sqrt{k}t)}{\sqrt{1+\lambda\text{sech}^{2}(\sqrt{k}t)}}=\frac{A\text{sech}(\sqrt{k}t)}{\sqrt{1+\lambda\text{sech}^{2}(\sqrt{k}t)}}$,
$A=x_{0}\sqrt{1+\lambda}.$
---
for the choices
| $k=a,$ | $\lambda=\frac{bx_{0}^{2}-2a}{4a-bx_{0}^{2}},$ | $x_{0}=\frac{1}{2}\sqrt{\frac{-3b\pm\sqrt{48ac+9b^{2}}}{c}}$
---|---|---
$48ac+9b^{2}>0.$
(10)

We also have the following homoclinic orbits (Melnikov [1963]) :

$x(t)=x_{0}\frac{\sqrt{1+\lambda}\text{tanh}(\sqrt{k}t)}{\sqrt{1+\lambda\text{tanh}^{2}(\sqrt{k}t)}}=\frac{A\text{tanh}(\sqrt{k}t)}{\sqrt{1+\lambda\text{tanh}^{2}(\sqrt{k}t)}}$,
$A=x_{0}\sqrt{1+\lambda}.$
---
for the choices
$k=\frac{1}{2}\left(-bx_{0}^{2}-2cx_{0}^{4}\right),\lambda=-\frac{2cx_{0}^{2}}{3\left(b+2cx_{0}^{2}\right)},x_{0}=\sqrt{\frac{-b\pm\sqrt{b^{2}+4ac}}{2c}}.$
$b^{2}+4ac>0.$
(11)
$x(t)=x_{0}\frac{\sqrt{1+\lambda}\text{sech}(\sqrt{k}t)}{\sqrt{1+\lambda\text{sech}^{2}(\sqrt{k}t)}}=\frac{A\text{sech}(\sqrt{k}t)}{\sqrt{1+\lambda\text{sech}^{2}(\sqrt{k}t)}}$,
$A=x_{0}\sqrt{1+\lambda}.$
---
for the choices
| $k=a,$ | $\lambda=\frac{bx_{0}^{2}-2a}{4a-bx_{0}^{2}},$ | $x_{0}=\frac{1}{2}\sqrt{\frac{-3b\pm\sqrt{48ac+9b^{2}}}{c}}$
---|---|---
$48ac+9b^{2}>0.$
(12)

Example 1. Let us consider the Duffing equation

$\ddot{x}-x+x^{3}+x^{5}=0,x(0)=1\text{ and }x^{\prime}(0)=0.$ (13)

Exact solution :

$x(t)=\frac{\sqrt{\frac{1}{6}(3+\sqrt{3})}\text{cn}\left(\sqrt{1+\frac{1}{\sqrt{3}}}t,\sqrt{3}-1\right)}{\sqrt{1+\frac{1}{6}(3+\sqrt{3})\text{cn}^{2}\left(\sqrt{1+\frac{1}{\sqrt{3}}}t,\sqrt{3}-1\right)}}.$
(14)

See Figure 1

Example 2. Let

$\ddot{x}+x+2x^{3}+3x^{5}=0,x(0)=1\text{ and }x^{\prime}(0)=0.$ (15)

Exact solution :

$x(t)=x_{0}\frac{\sqrt{1+\lambda+\mu}\cdot\text{cn}(\sqrt{\omega}t,m)}{\sqrt{1+\lambda\cdot\text{cn}^{2}(\sqrt{\omega}t,m)+\mu\cdot\text{cn}^{4}(\sqrt{\omega}t,m)}},$
---
where
$\lambda=4-3\sqrt{2},\mu=12\sqrt{2}-17,\omega=3\sqrt{2},m=\frac{1}{6}\left(3-2\sqrt{2}\right).$
(16)

See Figure 2.

###  2.1 Approximate Analytical Solution for the General Case.

Suppose we are given that :

$\ddot{x}+\omega_{0}^{2}x+F(t,x,\dot{x})=0\text{, }x(0)=x_{0}\text{ and
}x^{\prime}(0)=\dot{x}_{0}$ (17)

Let us consider the following $p$-problem :

$\ddot{x}+\omega_{0}^{2}x+pF(t,x,\dot{x})=0\text{, }x(0)=x_{0}\text{ and
}x^{\prime}(0)=\dot{x}_{0}$ (18)

Let $x_{p}=x_{p}(t)$ be the solution to the $p$-problem. We seek approximate
analytical solution in the ansatz form

$x_{p}=\mathbf{a}\cos\left(\psi\right)+\sum_{n=1}^{N}p^{n}u_{n}(\mathbf{a},\psi)+\text{o}(p^{N+1}),$
(19)

where each $u_{n}$ is a periodic function of $\psi,$ and a and $\psi$ are
assumed to vary with time according to

$\frac{d\mathbf{a}}{dt}\equiv\mathbf{a}=\sum_{n=1}^{N}p^{n}\mathbf{A}_{n}(\mathbf{a})+\text{o}(p^{N+1}).$
(20)
$\frac{d\psi}{dt}\equiv\dot{\psi}=\omega_{0}+\sum_{n=1}^{N}p^{n}\psi_{n}(\mathbf{a})+\text{o}(p^{N+1})\text{
.}$
---
or
$\frac{d\psi}{dt}\equiv\dot{\psi}=\sqrt{\omega_{0}^{2}+\sum_{n=1}^{N}p^{n}\psi_{n}(\mathbf{a})+\text{o}(p^{N+1})}.$
(21)

Define

$H(x,t)=\ddot{x}+\omega_{0}^{2}x+pF(t,x,\dot{x}).$ (22)

The next step is to write the residual $H_{p}(x_{p},t)$ as a power series in
$p$

$H(x_{p},t)=p\Upsilon_{1}+p^{2}\Upsilon_{2}+p^{3}\Upsilon_{3}+\cdots.$ (23)

For the determination of the unknown functions $u_{n}$, $\psi_{n}$,
$\mathbf{A}_{n},$ and $\mathbf{a},$ we equate to zero the coefficients
$\Upsilon_{n}$ in Eq. (23) and then we can get a system of odes. To avoid the
so-called secularity, we choose only the solutions that do not contain
$\cos\psi$ nor $\sin\psi$.

In the case, for $N=2$ (the second order approximation), we may use the
following formulas (we neglected all terms containing $p^{j}$ for $j\geq 3$) :

$\ddot{x}_{p}+\omega_{0}^{2}x_{p}=$
---
$p\omega_{0}\left(-2\sin(\psi)\mathbf{A}_{1}(\mathbf{a})-2\mathbf{a}\cos(\psi)\varphi_{1}(\mathbf{a})+\omega_{0}\left(u_{1}(\mathbf{a},\psi)+u_{1}{}^{(0,2)}(\mathbf{a},\psi)\right)\right)+$
$p^{2}\left(\begin{array}[]{c}\left(-\mathbf{a}\left(\varphi_{1}(\mathbf{a}){}^{2}+2\omega_{0}\varphi_{2}(\mathbf{a})\right)+\mathbf{A}_{1}(\mathbf{a})\mathbf{A}_{1}(\mathbf{a})\right)\cos(\psi)-\\\
\\\
\left(2\omega_{0}\mathbf{A}_{2}(\mathbf{a})+\mathbf{A}_{1}(\mathbf{a})\left(2\varphi_{1}(\mathbf{a})+\mathbf{a}\varphi_{1}^{\prime}(\mathbf{a})\right)\right)\sin(\psi)+\\\
\\\
\omega_{0}\left(\omega_{0}u_{2}(\mathbf{a},\psi)+2\varphi_{1}(\mathbf{a})u_{1}{}^{(0,2)}(\mathbf{a},\psi)+\omega_{0}u_{2}{}^{(0,2)}(\mathbf{a},\psi)+2\mathbf{A}_{1}(\mathbf{a})u_{1}{}^{(1,1)}(\mathbf{a},\psi)\right)\end{array}\right).$
(24)

Let $x_{p}=x$ for sake of simplicity. Then

$\begin{array}[]{l}\begin{array}[]{l}p\dot{x}=p^{2}\left(\text{$\omega_{0}$}u_{1}{}^{(0,1)}(\mathbf{a},\psi)+\mathbf{A}_{1}(\mathbf{a})\cos(\psi)-\mathbf{a}\varphi_{1}(\mathbf{a})\sin(\psi)\right)-\mathbf{a}p\text{$\omega_{0}$}\sin(\psi).\\\
\\\
px(t)^{2}=2\mathbf{a}p^{2}\cos(\psi)u_{1}(\mathbf{a},\psi)+\mathbf{a}^{2}p\cos^{2}(\psi).\\\
\\\
px(t)^{3}=\frac{3}{2}\mathbf{a}^{2}p^{2}(\cos(2\psi)+1)u_{1}(\mathbf{a},\psi)+\frac{1}{4}\mathbf{a}^{3}p(3\cos(\psi)+\cos(3\psi)).\\\
\\\
px^{4}=\mathbf{a}^{3}p^{2}(3\cos(\psi)+\cos(3\psi))u_{1}(\mathbf{a},\psi)+\frac{1}{8}\mathbf{a}^{4}p(4\cos(2\psi)+\cos(4\psi)+3).\\\
\\\
px^{5}=\frac{5}{8}\mathbf{a}^{4}p^{2}(4\cos(2\psi)+\cos(4\psi)+3)u_{1}(\mathbf{a},\psi)+\frac{1}{16}p\left(10\mathbf{a}^{5}\cos(\psi)+5\mathbf{a}^{5}\cos(3\psi)+\mathbf{a}^{5}\cos(5\psi)\right).\\\
\\\
px^{6}=\frac{3}{8}\mathbf{a}^{5}p^{2}(10\cos(\psi)+5\cos(3\psi)+\cos(5\psi))u_{1}(\mathbf{a},\psi)+\frac{1}{32}\mathbf{a}^{6}p(15\cos(2\psi)+6\cos(4\psi)+\cos(6\psi)+10).\\\
\\\
px^{7}=\frac{7}{32}\mathbf{a}^{6}p^{2}(15\cos(2\psi)+6\cos(4\psi)+\cos(6\psi)+10)u_{1}(\mathbf{a},\psi)+\frac{1}{64}\mathbf{a}^{7}p(35\cos(\psi)+21\cos(3\psi)+7\cos(5\psi)+\cos(7\psi)).\\\
\\\
p\dot{x}x=-\frac{1}{2}\mathbf{a}p^{2}\left(2\text{$\omega_{0}$}\sin(\psi)u_{1}(\mathbf{a},\psi)-2\text{$\omega_{0}$}\cos(\psi)u_{1}{}^{(0,1)}(\mathbf{a},\psi)+\mathbf{A}_{1}(\mathbf{a})(-\cos(2\psi))-\mathbf{A}_{1}(\mathbf{a})+\mathbf{a}\varphi_{1}(\mathbf{a})\sin(2\psi)\right)-\frac{1}{2}\mathbf{a}^{2}p\text{$\omega_{0}$}\sin(2\psi).\\\
\\\
p\dot{x}x^{2}=-\frac{1}{4}\mathbf{a}^{2}p^{2}\left(\begin{array}[]{c}-2\text{$\omega_{0}$}u_{1}{}^{(0,1)}(\mathbf{a},\psi)+4\text{$\omega_{0}$}\sin(2\psi)u_{1}(\mathbf{a},\psi)-2\text{$\omega_{0}$}\cos(2\psi)u_{1}{}^{(0,1)}(\mathbf{a},\psi)-3\mathbf{A}_{1}(\mathbf{a})\cos(\psi)-\\\
\mathbf{A}_{1}(\mathbf{a})\cos(3\psi)+\mathbf{a}\varphi_{1}(\mathbf{a})\sin(\psi)+\mathbf{a}\varphi_{1}(\mathbf{a})\sin(3\psi)\end{array}\right)-\\\
\frac{1}{4}\mathbf{a}^{3}p\text{$\omega_{0}$}(\sin(\psi)+\sin(3\psi)).\\\ \\\
p\dot{x}^{2}=-\mathbf{a}p^{2}\text{$\omega_{0}$}\left(2\text{$\omega_{0}$}\sin(\psi)u_{1}{}^{(0,1)}(\mathbf{a},\psi)+\mathbf{A}_{1}(\mathbf{a})\sin(2\psi)+\mathbf{a}\varphi_{1}(\mathbf{a})\cos(2\psi)-\mathbf{a}\varphi_{1}(\mathbf{a})\right)-\frac{1}{2}\mathbf{a}^{2}p\text{$\omega_{0}$}^{2}(\cos(2\psi)-1).\\\
\\\
p\dot{x}^{3}=-\frac{3}{4}\mathbf{a}^{2}p^{2}\text{$\omega_{0}$}^{2}\left(\begin{array}[]{c}-2\text{$\omega_{0}$}u_{1}{}^{(0,1)}(\mathbf{a},\psi)+2\text{$\omega_{0}$}\cos(2\psi)u_{1}{}^{(0,1)}(\mathbf{a},\psi)+\\\
\mathbf{A}_{1}(\mathbf{a})(-\cos(\psi))+\mathbf{A}_{1}(\mathbf{a})\cos(3\psi)+3\mathbf{a}\varphi_{1}(\mathbf{a})\sin(\psi)-\mathbf{a}\varphi_{1}(\mathbf{a})\sin(3\psi)\end{array}\right)-\\\
\frac{1}{4}\mathbf{a}^{3}p\text{$\omega_{0}$}^{3}(3\sin(\psi)-\sin(3\psi)).\end{array}\end{array}$
(25)

Then, the solution to the original problem is obtained by letting $p=1.$ Let
us consider the i.v.p.

$\ddot{x}-ax+bx^{3}+cx^{5}=\varepsilon(\gamma\cos\omega t\
-\delta\dot{x})\text{, }x(0)=x_{0}\text{ and }x^{\prime}(0)=\dot{x}_{0}.\text{
}$ (26)

####  2.1.1 First Case. $a<0$

Let $\omega_{0}=\sqrt{-a}.$The associated $p-$problem reads

$\ddot{x}+\omega_{0}^{2}x+p\left[\epsilon\dot{x}+~{}bx^{3}+cx^{5}-\phi(t)\right]=0,~{}x(0)=x_{0}\text{
and }x^{\prime}(0)=\dot{x}_{0}.\text{ }$ (27)

where

$\epsilon=\varepsilon\delta\text{ and }\phi(t)=\varepsilon\gamma\cos\omega t.$
(28)

Using the above formulas (25) gives

$\ddot{x}+\omega_{0}^{2}x+p\left[\epsilon\dot{x}+~{}bx^{3}+cx^{5}-\phi(t)\right]=$
---
$\frac{1}{16}p\left(\begin{array}[]{c}16\omega_{0}^{2}u_{1}(\mathbf{a},\psi)+16\omega_{0}^{2}u_{1}{}^{(0,2)}(\mathbf{a},\psi)+10\mathbf{a}^{5}c\cos(\psi)+\\\
5\mathbf{a}^{5}c\cos(3\psi)+\mathbf{a}^{5}c\cos(5\psi)+12\mathbf{a}^{3}b\cos(\psi)+\\\
4\mathbf{a}^{3}b\cos(3\psi)-32\omega_{0}\mathbf{A}_{1}(\mathbf{a})\sin(\psi)-\\\
32\mathbf{a}\omega_{0}\varphi_{1}(\mathbf{a})\cos(\psi)-16\mathbf{a}\omega_{0}\epsilon\sin(\psi)-16\phi(t)\end{array}\right)+$
$\frac{1}{8}p^{2}\left(\begin{array}[]{c}15\mathbf{a}^{4}cu_{1}(\mathbf{a},\psi)+20\mathbf{a}^{4}c\cos(2\psi)u_{1}(\mathbf{a},\psi)+\\\
5\mathbf{a}^{4}c\cos(4\psi)u_{1}(\mathbf{a},\psi)+12\mathbf{a}^{2}bu_{1}(\mathbf{a},\psi)+12\mathbf{a}^{2}b\cos(2\psi)u_{1}(\mathbf{a},\psi)+\\\
16\omega_{0}\mathbf{A}_{1}(\mathbf{a})u_{1}{}^{(1,1)}(\mathbf{a},\psi)+16\omega_{0}\varphi_{1}(\mathbf{a})u_{1}{}^{(0,2)}(\mathbf{a},\psi)+8\omega_{0}^{2}u_{2}(\mathbf{a},\psi)+\\\
8\omega_{0}^{2}u_{2}{}^{(0,2)}(\mathbf{a},\psi)+8\omega_{0}\epsilon
u_{1}{}^{(0,1)}(\mathbf{a},\psi)+\\\
\sin(\psi)\left(-8\mathbf{aA}_{1}(\mathbf{a})\varphi_{1}^{\prime}(\mathbf{a})-16\mathbf{A}_{1}(\mathbf{a})\varphi_{1}(\mathbf{a})-16\omega_{0}\mathbf{A}_{2}(\mathbf{a})-8\mathbf{a}\epsilon\varphi_{1}(\mathbf{a})\right)+\\\
\cos(\psi)\left(8\epsilon\mathbf{A}_{1}(\mathbf{a})+8\mathbf{A}_{1}(\mathbf{a})\mathbf{A}_{1}(\mathbf{a})-16\mathbf{a}\omega_{0}\varphi_{2}(\mathbf{a})-8\mathbf{a}\varphi_{1}(\mathbf{a}){}^{2}\right)\end{array}\right)$
(29)

The required solutions reads

$\begin{array}[]{l}u_{1}(\mathbf{a},\psi)=\frac{1}{384\omega_{0}^{2}}\left[\mathbf{a}^{5}c\cos(5\psi)+3\mathbf{a}^{3}\cos(3\psi)\left(5\mathbf{a}^{2}c+4b\right)+384\phi(t)\right].\\\
\\\
u_{2}(\mathbf{a},\psi)=\frac{\mathbf{a}^{2}}{294912\omega_{0}^{4}}\left(\begin{array}[]{c}\mathbf{a}^{7}c^{2}(-5280\cos(3\psi)+160\cos(5\psi)+95\cos(7\psi)+3\cos(9\psi))+\\\
\\\ 72\mathbf{a}^{5}bc(-164\cos(3\psi)+4\cos(5\psi)+\cos(7\psi))+\\\ \\\
32\mathbf{a}^{3}\left(9b^{2}(\cos(5\psi)-21\cos(3\psi))+20c\omega_{0}\epsilon(27\sin(3\psi)+\sin(5\psi))\right)+\\\
\\\ 12288\mathbf{a}^{2}c\phi(t)(20\cos(2\psi)+\cos(4\psi)-45)+\\\ \\\
6912\mathbf{a}b\omega_{0}\epsilon\sin(3\psi)+147456b\phi(t)(\cos(2\psi)-3)\end{array}\right).\\\
\\\ \mathbf{A}_{1}(\mathbf{a})=-\frac{\mathbf{a}\epsilon}{2}\text{,
}\mathbf{A}_{2}(\mathbf{a})=\frac{5\mathbf{a}^{5}c\epsilon+3\mathbf{a}^{3}b\epsilon}{16\omega_{0}^{2}}.\\\
\\\
\varphi_{1}(\mathbf{a})=\frac{5\mathbf{a}^{4}c+6\mathbf{a}^{2}b}{16\omega_{0}}\text{,
}\varphi_{2}(\mathbf{a})=\frac{-55\mathbf{a}^{8}c^{2}-240\mathbf{a}^{6}bc-180\mathbf{a}^{4}b^{2}-384\omega_{0}^{2}\epsilon^{2}}{3072\omega_{0}^{3}}.\end{array}$
(30)

Then

$\begin{array}[]{l}x_{p}(t)=\mathbf{a}\cos(\psi)+p\frac{\mathbf{a}^{5}c(15\cos(3\psi)+\cos(5\psi))+12\mathbf{a}^{3}b\cos(3\psi)+384\phi(t)}{384\omega_{0}^{2}}+\\\
\\\
\frac{\mathbf{a}^{2}p^{2}}{294912\omega_{0}^{4}}\left(\begin{array}[]{c}\mathbf{a}^{7}c^{2}(-5280\cos(3\psi)+160\cos(5\psi)+95\cos(7\psi)+3\cos(9\psi))+\\\
\\\ 72\mathbf{a}^{5}bc(-164\cos(3\psi)+4\cos(5\psi)+\cos(7\psi))+\\\ \\\
32\mathbf{a}^{3}\left(-189b^{2}\cos(3\psi)+9b^{2}\cos(5\psi)+20c\omega_{0}\epsilon(27\sin(3\psi)+\sin(5\psi))\right)+\\\
\\\ 12288\mathbf{a}^{2}c\phi(t)(20\cos(2\psi)+\cos(4\psi)-45)+\\\ \\\
6912\mathbf{a}b\omega_{0}\epsilon\sin(3\psi)+147456b\phi(t)(\cos(2\psi)-3)\end{array}\right).\\\
\\\ \text{The odes for determining }\mathbf{a}=\mathbf{a}(t)\text{ and
}\psi=\psi(t)\text{ are : }\\\ \\\
\mathbf{a}=-p\frac{\mathbf{a}\epsilon}{2}+p^{2}\frac{1}{16\omega_{0}^{2}}\left[3b\epsilon\mathbf{a}^{3}+5c\epsilon\mathbf{a}^{5}\right]\\\
\\\
\dot{\psi}=\omega_{0}+\frac{p}{16\omega_{0}}\left(6\mathbf{a}^{2}b+5\mathbf{a}^{4}c\right)+\frac{p^{2}}{3072\omega_{0}^{3}}\left(-55\mathbf{a}^{8}c^{2}-240\mathbf{a}^{6}bc-180\mathbf{a}^{4}b^{2}-384\omega_{0}^{2}\epsilon^{2}\right).\end{array}$
(31)

Example 3. see Figure 3 Let

$\ddot{x}+x+0.025~{}\dot{x}+2x^{3}+x^{5}u=0.01\cos(0.1t)\wedge x(0)=0.25\wedge
x^{\prime}(0)=0.\text{ }$ (32)

The approximate analytical solution is given by

$x(t)=0.01\cos(0.1t)+0.23913e^{-0.0125t}\cos\left(1.t-0.0204369e^{-0.05t}-1.71549e^{-0.025t}+1.72421\right)+$
---
$\left(0.0000305442e^{-0.0625t}+0.000854635e^{-0.0375t}\right)\cos\left(3.t-0.0613106e^{-0.05t}-5.14647e^{-0.025t}+5.17263\right)+$
$2.03$e$-6e^{-0.0625t}\cos\left(5.t-0.102184e^{-0.05t}-8.57746e^{-0.025t}+8.62105\right).$

####  2.1.2 Second Case . $a>0.$

Let us consider the I.V.P

$x(t)=\eta+u(t)\text{, where }-a\eta+b\eta^{3}+c\eta^{5}=0\text{, }\eta\neq
0.$ (33)

Then

$\ddot{x}-ax+bx^{3}+cx^{5}-\varepsilon(\gamma\cos\omega t\ -\delta\dot{x})=$
---
$u^{\prime\prime}(t)+\omega_{0}^{2}u(t)+B~{}u(t)^{2}+~{}Cu(t)^{3}+Du(t)^{4}+Eu(t)^{5}-\varepsilon\left(\gamma\cos(\omega
t)-\delta u^{\prime}(t)\right),$
where
$\omega_{0}^{2}=-a+3b\eta^{2}+5c\eta^{4}$, $\
B=3b\eta+10c\eta^{3},~{}C=b+10c\eta^{2}$, $D=5c\eta$ and $E=c.$
(34)

The associated $p-$problem reads

$u^{\prime\prime}(t)+\omega_{0}^{2}u(t)+p\left[\epsilon\dot{x}+B~{}u(t)^{2}+~{}Cu(t)^{3}+Du(t)^{4}+Eu(t)^{5}-\phi(t)\right]=0,~{}u(0)=x_{0}-\eta\text{
and }u^{\prime}(0)=\dot{x}_{0}.\text{ }$ (35)

where

$\epsilon=\varepsilon\delta\text{ and }\phi(t)=\varepsilon\gamma\cos\omega t.$
(36)

The second order approximation is given by

$\begin{array}[]{l}u(t)=\mathbf{a}\cos(\psi)+\\\
\frac{p}{1920\text{$\mathbf{\omega}_{0}$}^{2}}\left(\begin{array}[]{c}5\mathbf{a}^{5}E(15\cos(3\psi)+\cos(5\psi))+16\mathbf{a}^{4}D(20\cos(2\psi)+\cos(4\psi)-45)+\\\
60\mathbf{a}^{3}C\cos(3\psi)+320\mathbf{a}^{2}B(\cos(2\psi)-3)+1920\phi(t)\end{array}\right)+\\\
\\\
\frac{\mathbf{a}^{2}p^{2}\left(\begin{array}[]{c}175\mathbf{a}^{7}E^{2}(-5280\cos(3\psi)+160\cos(5\psi)+95\cos(7\psi)+3\cos(9\psi))+\\\
640\mathbf{a}^{6}DE(-24710\cos(2\psi)-168\cos(4\psi)+198\cos(6\psi)+5\cos(8\psi)+38115)-\\\
280\mathbf{a}^{5}\left(36\cos(3\psi)\left(205CE+72D^{2}\right)-4\cos(5\psi)\left(45CE+184D^{2}\right)-\cos(7\psi)\left(45CE+16D^{2}\right)\right)+\\\
1792\mathbf{a}^{4}(-5\cos(2\psi)(2425BE+1458CD)+6\cos(4\psi)(27CD-20BE)+\\\
9\cos(6\psi)(5BE+2CD)+150(140BE+81CD))+\\\
1120\mathbf{a}^{3}\left(\begin{array}[]{c}-27\cos(3\psi)\left(16BD+35C^{2}\right)+\\\
\cos(5\psi)\left(176BD+45C^{2}\right)+100E\text{$\mathbf{\omega}_{0}$}\epsilon(27\sin(3\psi)+\sin(5\psi))\end{array}\right)+\\\
21504\mathbf{a}^{2}(-775BC\cos(2\psi)+25BC\cos(4\psi)+1500BC+800D\text{$\mathbf{\omega}_{0}$}\epsilon\sin(2\psi)+\\\
16D\text{$\mathbf{\omega}_{0}$}\epsilon\sin(4\psi)+100E\phi(t)(20\cos(2\psi)+\cos(4\psi)-45))+\\\
134400\mathbf{a}\left(8B^{2}\cos(3\psi)+9C\text{$\mathbf{\omega}_{0}$}\epsilon\sin(3\psi)+48D\phi(t)\cos(3\psi)\right)+\\\
2867200(2B\text{$\mathbf{\omega}_{0}$}\epsilon\sin(2\psi)+9C\phi(t)(\cos(2\psi)-3))\end{array}\right)}{51609600\text{$\mathbf{\omega}_{0}$}^{4}}.\\\
\\\ \text{The odes for determining }\mathbf{a=a(t)}\text{ and
}\psi=\psi(t)\text{ are : }\\\ \\\
\mathbf{\dot{a}}=-\frac{\mathbf{a}p\epsilon}{2}+p^{2}\frac{5\mathbf{a}^{5}E\epsilon+3\mathbf{a}^{3}C\epsilon}{16\text{$\mathbf{\omega}_{0}$}^{2}}.\\\
\\\
\dot{\psi}=\text{$\mathbf{\omega}_{0}+$}\frac{p\left(5\mathbf{a}^{4}E+6\mathbf{a}^{2}C\right)}{16\text{$\mathbf{\omega}_{0}$}}+\frac{p^{2}\left(-275\mathbf{a}^{8}E^{2}-1200\mathbf{a}^{6}CE-6048\mathbf{a}^{6}D^{2}-13440\mathbf{a}^{4}BD-900\mathbf{a}^{4}C^{2}-6400\mathbf{a}^{2}B^{2}+23040\mathbf{a}^{2}D\phi(t)+15360B\phi(t)-1920\text{$\mathbf{\omega}_{0}$}^{2}\epsilon^{2}\right)}{15360\text{$\mathbf{\omega}_{0}$}^{3}}\\\
\lx@intercol\hfil\hfil\lx@intercol\end{array}$ (37)

for the error of approximate analytical on can refer to figure 4

##  3 Homoclinic orbits in the unperturbed system

For the unperturbed system with fractional order displacement, when
$\varepsilon=0$, the differential equation (3) simplifies to

$\ddot{x}-ax+bx^{3}+cx^{5}=0.~{}$ (38)

Let

$\Delta:=b^{2}+4ac.$ (39)

Equilibrium points for $\Delta>0$ are :

$\left(x,\dot{x}\right)=\left(\pm\sqrt{\frac{-b+\sqrt{b^{2}+4ac}}{2c}},0\right)\text{
: centers}$ (40)

Define

$x_{e}^{+}=\sqrt{\frac{-b+\sqrt{b^{2}+4ac}}{2c}}\text{ and
}x_{e}^{-}=-\sqrt{\frac{-b+\sqrt{b^{2}+4ac}}{2c}}$ (41)

The energy function for (38) is

$\frac{1}{2}\dot{x}(t)^{2}-\frac{1}{2}ax(t)^{2}+\frac{1}{4}bx(t)^{4}+\frac{1}{6}cx(t)^{6}=K$
(42)

where $K$ is the energy constant dependent on the initial amplitude
$x(0)=x_{0}$ and initial velocity $x^{\prime}(0)=\dot{x}_{0}$ :

$K=\frac{1}{2}\dot{x}_{0}^{2}-\frac{1}{2}ax_{0}^{2}+\frac{1}{4}bx_{0}^{4}+\frac{1}{6}cx_{0}^{6}.$
(43)

Dependently on $K$, the level sets are different. For all of them it is common
that they form closed periodic orbits which surround the fixed points
$(\mathrm{x},\dot{\mathrm{x}})=(x_{e}^{+},\ 0)$ or
$(\mathrm{x},\dot{\mathrm{x}})=(x_{e}^{-},\ 0)$ or all the three fixed points
$(x_{e}^{\pm},\ 0)$ and $(0,\ 0)$. The boundary between these two groups of
orbits corresponds to $K=0$, when

$\dot{x}_{0}=\pm
x_{0}\sqrt{\frac{1}{6}\left(6a-3bx_{0}^{2}-2cx_{0}^{4}\right)}.$ (44)

The level set

$\frac{\dot{x}^{2}}{2}-\frac{1}{2}ax^{2}+\frac{1}{4}bx^{4}+\frac{1}{6}cx^{6}=0\text{
}\quad$ (45)

is composed of two homoclinic orbits

$\Gamma_{+}^{0}(\mathrm{t})\equiv(x_{+}^{0}(\mathrm{t}),\dot{x}_{+}^{0}(\mathrm{t})),\text{
}\quad$ (46)
$\Gamma_{-}^{0}(\mathrm{t})\equiv(x_{-}^{0}(\mathrm{t}),\dot{x}_{-}^{0}(\mathrm{t})),\text{
}$ (47)

which connect the fixed hyperbolic saddle point $(0,\ 0)$ to itself and
contain the stable and unstable manifolds. The functions
$x_{\pm}^{0}(\mathrm{t})$ may be evaluated using formulas (9)- (12). See
Figure 5

The homoclinic orbit P. Holmes& J. Marsden [1981] separates the phase plane
into two areas. Inside the separatrix curve the orbits are around one of the
centers, and outside the separatrix curve the orbits surround both the centers
and the saddle point. Physically it means that for certain initial conditions
the oscillations are around one steady-state position, and for others around
all the steady- state solutions (two stable and an unstable).

##  4  Melnikov’s criteria for chaos

Let us form Melnikov’s function for (3) and $\Gamma_{+}^{0}(\mathrm{t})$, i.e.
$\Gamma_{-}^{0}(\mathrm{x})$ given by (46) and (47)

$M(t_{0})=\int_{-\infty}^{+\infty}\dot{x}^{0}(t)[\gamma\cos\omega(t+t_{0})-\delta\dot{x}^{0}(t)]dt,\text{
}\quad$ (48)

Let

$x^{0}(t)=\frac{A\text{sech}\left(\sqrt{k}t\right)}{\sqrt{1+\lambda\cdot\text{sech}^{2}\left(\sqrt{k}t\right)}}.$
(49)

Then

$\begin{array}[]{l}\dot{x}^{0}(t)=-\frac{A\sqrt{k}\tanh\left(\sqrt{k}t\right)\text{sech}\left(\sqrt{k}t\right)}{\left(\lambda\text{sech}^{2}\left(\sqrt{k}t\right)+1\right)^{3/2}}\\\
\text{and}\\\
\dot{x}^{0}(t)^{2}=\frac{2A^{2}k\sinh^{2}\left(2\sqrt{k}t\right)}{\left(\cosh\left(2\sqrt{k}t\right)+2\lambda+1\right)^{3}}.\end{array}$
(50)

We have :

$M(t_{0})=\gamma\int_{-\infty}^{+\infty}\dot{x}^{0}(t)\cos\omega(t+t_{0})dt-\delta\int_{-\infty}^{+\infty}\dot{x}^{0}(t)^{2}dt$
---
$=\gamma I_{1}-\delta I_{2}$,
where
$I_{1}=-A\sqrt{k}\int_{-\infty}^{+\infty}\frac{\tanh\left(\sqrt{k}t\right)\text{sech}\left(\sqrt{k}t\right)}{\left(1+\lambda\text{sech}^{2}\left(\sqrt{k}t\right)\right)^{3/2}}\cos\omega(t+t_{0})dt$
and
$I_{2}=2A^{2}k\int_{-\infty}^{+\infty}\frac{\sinh^{2}\left(2\sqrt{k}t\right)}{\left(\cosh\left(2\sqrt{k}t\right)+2\lambda+1\right)^{3}}dt.$
(51)

The value $I_{2}$ is evaluated as

$I_{2}=\frac{A^{2}\sqrt{k}\left(2\sqrt{\lambda+1}\lambda^{3/2}+\sqrt{\lambda(\lambda+1)}-\tanh^{-1}\left(\sqrt{\frac{\lambda}{\lambda+1}}\right)\right)}{4(\lambda(\lambda+1))^{3/2}}.$
(52)

The first integral is hard to evaluate in closed form. For this reason, we
will approximate it by taking into account the folowing Chebyshev
approximation:

$\frac{x}{(1+\lambda x^{2})^{3/2}}\approx rx+sx^{3}$ for $-1\leq x\leq 1,$
---
where
$r=\frac{64\sqrt{2}\left(-\sin^{2}\left(\frac{\pi}{8}\right)\sqrt{\lambda\sin^{2}\left(\frac{\pi}{8}\right)+1}-\lambda\sin^{4}\left(\frac{\pi}{8}\right)\sqrt{\lambda\sin^{2}\left(\frac{\pi}{8}\right)+1}+\cos^{2}\left(\frac{\pi}{8}\right)\sqrt{\lambda\cos^{2}\left(\frac{\pi}{8}\right)+1}+\lambda\cos^{4}\left(\frac{\pi}{8}\right)\sqrt{\lambda\cos^{2}\left(\frac{\pi}{8}\right)+1}\right)}{\left(4-\left(\sqrt{2}-2\right)\lambda\right)^{3/2}\left(\left(2+\sqrt{2}\right)\lambda+4\right)^{3/2}}$
and
$s=\frac{64\sqrt{2}\left(\lambda\sin^{2}\left(\frac{\pi}{8}\right)\sqrt{\lambda\sin^{2}\left(\frac{\pi}{8}\right)+1}+\sqrt{\lambda\sin^{2}\left(\frac{\pi}{8}\right)+1}+\lambda\left(-\cos^{2}\left(\frac{\pi}{8}\right)\right)\sqrt{\lambda\cos^{2}\left(\frac{\pi}{8}\right)+1}-\sqrt{\lambda\cos^{2}\left(\frac{\pi}{8}\right)+1}\right)}{\left(4-\left(\sqrt{2}-2\right)\lambda\right)^{3/2}\left(\left(2+\sqrt{2}\right)\lambda+4\right)^{3/2}}..$
(53)

Let us evaluate $I_{1}$ in (3) using (53). We have :

$I_{1}=-A\sqrt{k}\int_{-\infty}^{+\infty}\frac{\tanh\left(\sqrt{k}t\right)\text{sech}\left(\sqrt{k}t\right)}{\left(1+\lambda\text{sech}^{2}\left(\sqrt{k}t\right)\right)^{3/2}}\cos\omega(t+t_{0})dt\approx$
---
$A\sqrt{k}\int_{-\infty}^{+\infty}\tanh\left(\sqrt{k}t\right)\left[r\text{sech}\left(\sqrt{k}t\right)+s\text{sech}^{3}\left(\sqrt{k}t\right)\right]\cos\omega(t+t_{0})dt$
$A\sqrt{k}\left[r\int_{-\infty}^{+\infty}\tanh\left(\sqrt{k}t\right)\text{sech}\left(\sqrt{k}t\right)\cos\omega(t+t_{0})dt+s\int_{-\infty}^{+\infty}\tanh\left(\sqrt{k}t\right)\text{sech}^{3}\left(\sqrt{k}t\right)\cos\omega(t+t_{0})dt\right]$
$A\sqrt{k}\left[-r\frac{\omega\pi}{k}\text{sech}\left(\frac{\omega\pi}{2\sqrt{k}}\right)\sin\left(\omega
t_{0}\right)+s\frac{\omega\pi\left(k+\omega^{2}\right)}{6k^{2}}\text{sech}\left(\frac{\omega\pi}{2\sqrt{k}}\right)\sin\left(\omega
t_{0}\right)\right].$
(54)

Thus the Melnikov function reads

$M(x^{0}(t),t_{0})=M(t_{0})=$
---
$\gamma
A\sqrt{k}\left[-r\frac{\omega\pi}{k}+s\frac{\omega\pi\left(k+\omega^{2}\right)}{6k^{2}}\right]\text{sech}\left(\frac{\pi}{2\sqrt{k}}\omega\right)\sin\left(\omega
t_{0}\right)-\delta\frac{A^{2}\sqrt{k}\left(2\sqrt{\lambda+1}\lambda^{3/2}+\sqrt{\lambda(\lambda+1)}-\tanh^{-1}\left(\sqrt{\frac{\lambda}{\lambda+1}}\right)\right)}{4(\lambda(\lambda+1))^{3/2}}$
for the orbit
$x^{0}(t)=\frac{A\text{sech}\left(\sqrt{k}t\right)}{\sqrt{1+\lambda\cdot\text{sech}^{2}\left(\sqrt{k}t\right)}}.$
The numbers $r$ and $s$ are found from (53).
(55)

Now, assume an orbit of the form

$x^{0}(t)=\frac{A\text{tanh}\left(\sqrt{k}t\right)}{\sqrt{1+\lambda\cdot\text{tanh}^{2}\left(\sqrt{k}t\right)}}.$
(56)

We have :

$M(t_{0})=\gamma\int_{-\infty}^{+\infty}\dot{x}^{0}(t)\cos\omega(t+t_{0})dt-\delta\int_{-\infty}^{+\infty}\dot{x}^{0}(t)^{2}dt$
---
$=\gamma J_{1}-\delta J_{2}$,
where
$J_{1}=A\sqrt{k}\int_{-\infty}^{+\infty}\frac{\text{sech}^{2}\left(\sqrt{k}t\right)\cos\left(\omega\left(t+t_{0}\right)\right)}{\left(1+\lambda\tanh^{2}\left(\sqrt{k}t\right)\right)^{3/2}}dt$
and
$J_{2}=A^{2}k\int_{-\infty}^{+\infty}\frac{\text{sech}^{4}\left(\sqrt{k}t\right)}{\left(1+\lambda\tanh^{2}\left(\sqrt{k}t\right)\right)^{3}}dt.$
(57)

The value $J_{2}$ is evaluated as

$J_{2}=\frac{A^{2}\sqrt{k}\left(\sqrt{\lambda}(3\lambda+1)+(\lambda+1)(3\lambda-1)\tan^{-1}\left(\sqrt{\lambda}\right)\right)}{4\lambda^{3/2}(\lambda+1)}.$
(58)

In order to evaluate the value of $J_{1}$ we will use the following Chebyshev
approximation :

$\begin{array}[]{l}\frac{1-x}{(1+\lambda x)^{3/2}}\approx
1+\bar{r}x+\bar{s}x^{2}\text{ for }-1\leq x\leq 1,\\\ \text{where}\\\
\begin{array}[]{c}\bar{r}=\frac{1}{3}\sqrt{2}\left(\frac{2\sqrt{3}}{\left(\sqrt{3}\lambda+2\right)^{3/2}}-\frac{3}{\left(\sqrt{3}\lambda+2\right)^{3/2}}-\frac{2\sqrt{3}}{\left(2-\sqrt{3}\lambda\right)^{3/2}}-\frac{3}{\left(2-\sqrt{3}\lambda\right)^{3/2}}\right)\\\
\bar{s}=\frac{2}{3}\left(-\frac{\sqrt{6}}{\left(\sqrt{3}\lambda+2\right)^{3/2}}+\frac{2\sqrt{2}}{\left(\sqrt{3}\lambda+2\right)^{3/2}}+\frac{\sqrt{6}}{\left(2-\sqrt{3}\lambda\right)^{3/2}}+\frac{2\sqrt{2}}{\left(2-\sqrt{3}\lambda\right)^{3/2}}-2\right)\end{array}\end{array}$
(59)

We have :

$J_{1}=A\sqrt{k}\int_{-\infty}^{+\infty}\frac{\text{sech}^{2}\left(\sqrt{k}t\right)\cos\left(\omega\left(t+t_{0}\right)\right)}{\left(1+\lambda\tanh^{2}\left(\sqrt{k}t\right)\right)^{3/2}}dt=A\sqrt{k}\int_{-\infty}^{+\infty}\frac{\left(1-\text{tanh}^{2}\left(\sqrt{k}t\right)\right)\cos\left(\omega\left(t+t_{0}\right)\right)}{\left(1+\lambda\tanh^{2}\left(\sqrt{k}t\right)\right)^{3/2}}dt\approx$
---
$A\sqrt{k}\int_{-\infty}^{+\infty}\left[1+~{}\bar{r}\cdot\text{tanh}^{2}\left(\sqrt{k}t\right)+\bar{s}\cdot\text{tanh}^{4}\left(\sqrt{k}t\right)\right]\cos\omega(t+t_{0})dt=$
$A\sqrt{k}\int_{-\infty}^{+\infty}\left[~{}\bar{r}\cdot\text{tanh}^{2}\left(\sqrt{k}t\right)+\bar{s}\cdot\text{tanh}^{4}\left(\sqrt{k}t\right)\right]\cos\omega(t+t_{0})dt=$
$A\sqrt{k}\left[-\bar{r}\frac{\omega\pi\text{csch}\left(\frac{\omega\pi}{2\sqrt{k}}\right)\cos\left(\omega
t_{0}\right)}{k}+\bar{s}\frac{\omega\pi\left(\omega^{2}-8k\right)\text{csch}\left(\frac{\omega\pi}{2\sqrt{k}}\right)\cos\left(\omega
t_{0}\right)}{6k^{2}}\right].$
(60)

Thus the Melnikov function for the orbit (53) is given by

$M(t_{0})=\gamma
A\sqrt{k}\left[-\bar{r}\frac{\omega\pi}{k}+\bar{s}\frac{\omega\pi\left(\omega^{2}-8k\right)}{6k^{2}}\right]\text{csch}\left(\frac{\omega\pi}{2\sqrt{k}}\right)\cos\left(\omega
t_{0}\right)-\delta\frac{A^{2}\sqrt{k}\left(\sqrt{\lambda}(3\lambda+1)+(\lambda+1)(3\lambda-1)\tan^{-1}\left(\sqrt{\lambda}\right)\right)}{4\lambda^{3/2}(\lambda+1)}.$
(61)

The numbres $\bar{r}~{}$ and $\bar{s}~{}$are found from (59). Other formulas
that may be useful for computing are given in the Appendix.

##  5 Chaos Control.

###  5.1 Poincaré Map.

Let $z_{0}(\omega,t)$ be the solution to the i.v.p.

$\ddot{x}+\delta\dot{x}-ax+bx^{3}+cx^{5}=\gamma\cos\omega t\text{,
}x(0)=x_{0}\text{ and }x^{\prime}(0)=\dot{x}_{0}\text{ on~{}}0\leq
t\leq\frac{2\pi}{\omega}.\text{ }$ (62)

Next, let $z_{1}(\omega,t)$ be the solution to the i.v.p.

$\ddot{x}+\delta\dot{x}-ax+bx^{3}+cx^{5}=\gamma\cos\omega t\text{,
}x\left(0\right)=z_{0}\left(\omega,\frac{2\pi}{\omega}\right)\text{ and
}x^{\prime}\left(0\right)=z_{0}^{\prime}\left(\omega,\frac{2\pi}{\omega}\right)\text{
on~{}}0\leq t\leq\frac{2\pi}{\omega}\text{ .}$ (63)

Suppose we already found $z_{j}(\omega,t)$ for $j=0,1,2,...,n-1$. Then, the
function $z_{n}(\omega,t)$ is defined to be the solution to the i.v.p.

$\ddot{x}+\delta\dot{x}-ax+bx^{3}+cx^{5}=\gamma\cos\omega t\text{,
}x\left(0\right)=z_{n-1}\left(\omega,\frac{2\pi}{\omega}\right)\text{ and
}x^{\prime}\left(0\right)=z_{n-1}^{\prime}\left(\omega,\frac{2\pi}{\omega}\right)\text{
on~{}}0\leq t\leq\frac{2\pi}{\omega}.\text{ }$ (64)

We obtain the sequences

$P_{n}=z_{n-1}\left(\omega,\frac{2\pi}{\omega}\right)\text{ and
}Q_{n}=z_{n-1}^{\prime}\left(\omega,\frac{2\pi}{\omega}\right)\text{,
}n=1,2,3,....$ (65)

Thus, for a given $\omega$ we find the respective $\gamma_{\omega}$ value such
that the oscillator is chaotic for $\gamma=\gamma_{\omega}$ and non chaotic
for $0<\gamma<\gamma_{\omega}$. The $\gamma_{\omega}$ value is determined
experimentally. Let us consider the particular values

$a=b=1\text{,~{}}c=0,\text{ }\delta=0.1\text{ }\omega=1.4$

The first $\gamma$ chaotic $\gamma~{}$value was estimated as
$\gamma=\gamma_{1.4}=0.34$. The transition to chaos appears to occur between
$\gamma=0.34$. and$~{}\gamma=0.35$. See Figure 6

For the values $0<\gamma<0.34$ we have several bifurcation values (S. Wiggins
[1988]) . See Figures a, Figure c and figure 7

The curve for chaos in Figure 7. The Chebyshev approximation is depicted in
figure 8

The results of computations are shown in Table 1 for different values of
$\omega.$

$\begin{array}[]{ccc}\omega&\gamma_{\omega}&\lambda_{\omega}\\\ 0.050&0.387&0.044\\\ 0.100&0.402&0.023\\\ 0.125&0.402&0.065\\\ 0.150&0.397&0.063\\\ 0.200&0.380&0.051\\\ 0.225&0.389&0.064\\\ 0.250&0.381&0.057\\\ 0.300&0.382&0.091\\\ 0.350&0.360&0.020\\\ 0.400&0.342&0.047\\\ 0.450&0.478&0.009\\\ 0.500&0.640&0.098\\\ 0.525&0.633&0.069\\\ 0.600&0.381&0.020\\\ 0.625&0.376&0.015\\\ 0.650&0.375&0.027\\\ 0.700&0.429&0.067\\\ 0.725&0.450&0.008\\\ 0.750&0.522&0.011\\\ 0.800&0.721&0.151\\\ 0.825&0.810&0.190\\\ 0.925&1.336&0.134\\\ 0.950&1.261&0.067\\\ 1.000&0.939&0.119\\\ 1.100&0.173&0.115\\\ 1.125&0.147&0.058\\\ 1.150&0.136&0.068\\\ 1.200&0.199&0.048\\\ 1.225&0.214&0.074\end{array}$ | $\begin{array}[]{ccc}\omega&\gamma_{\omega}&\lambda_{\omega}\\\ 1.250&0.233&0.164\\\ 1.300&0.268&0.088\\\ 1.325&0.285&0.133\\\ 1.350&0.304&0.170\\\ 1.400&0.340&0.140\\\ 1.425&0.358&0.003\\\ 1.450&0.381&0.002\\\ 1.500&0.423&0.015\\\ 1.525&0.447&0.131\\\ 1.550&0.471&0.023\\\ 1.600&0.510&0.091\\\ 1.625&0.518&0.167\\\ 1.650&0.529&0.163\\\ 1.700&0.548&0.044\\\ 1.725&0.556&0.065\\\ 1.750&0.573&0.137\\\ 1.800&0.605&0.094\\\ 1.825&0.609&0.071\\\ 1.850&0.618&0.062\\\ 1.900&0.636&0.303\\\ 1.925&0.643&0.306\\\ 1.950&0.655&0.121\\\ 2.000&0.684&0.186\\\ 2.025&0.688&0.054\\\ 2.050&0.695&0.255\\\ 2.100&0.705&0.106\\\ 2.125&0.706&0.110\\\ 2.150&0.707&0.213\\\ 2.200&0.703&0.111\end{array}$ | $\begin{array}[]{ccc}\omega&\gamma_{\omega}&\lambda_{\omega}\\\ 2.300&0.699&0.027\\\ 2.325&0.724&0.233\\\ 2.350&0.763&0.006\\\ 2.400&0.840&0.024\\\ 2.425&0.858&0.266\\\ 2.450&0.862&0.006\\\ 2.500&0.839&0.272\\\ 2.525&0.841&0.203\\\ 2.550&0.826&0.084\\\ 2.600&0.783&0.125\\\ 2.700&0.786&0.002\\\ 2.800&1.263&0.012\\\ 2.825&1.351&0.016\\\ 2.850&1.525&0.039\\\ 2.900&1.871&0.026\\\ 2.925&1.937&0.022\\\ 3.000&2.155&0.004\\\ 3.100&2.168&0.032\\\ 3.200&2.530&0.031\\\ 3.225&2.590&0.051\\\ 3.500&3.870&0.025\\\ 3.525&3.955&0.013\\\ 3.650&4.736&0.058\\\ 3.700&4.999&0.315\\\ 3.750&4.987&0.365\\\ 3.800&5.066&0.027\\\ 3.900&6.137&0.014\\\ 3.925&6.288&0.038\\\ 4.000&6.787&0.049\end{array}$
---|---|---
$\text{Table 1. Suspected for chaos }\gamma_{\omega}\text{values and their
Lyapunov exponents. }$

###  5.2 Chaos supression. Delayed feedback controller using Pyragas method .

Assume that the following oscillator is chaotic Md Abdur Razzak, [2016]:

$\ddot{x}+\delta\dot{x}-ax+bx^{3}+cx^{5}=\gamma\cos\omega t.$ (66)

In order to supress the chaos, we introduce two constants $\mu$ and $\tau$ as
follows :

$x^{\prime\prime}(t)+\delta
x^{\prime}(t)-ax(t)+bx^{3}(t)+cx^{5}(t)=\gamma\cos\omega
t+\mu[x^{\prime}(t-\tau)-x^{\prime}(t)].$

The function

$x_{\mu,\tau}(t)=\mu[x^{\prime}(t-\tau)-x^{\prime}(t)]$

is called a delayed feed-back controller. The constants $\mu$ and $\tau$ are
chosen so that the solution to the i.v.p.

$x^{\prime\prime}(t)+\delta
x^{\prime}(t)-ax(t)+bx^{3}(t)+cx^{5}(t)=\gamma\cos\omega
t+\mu[x^{\prime}(t-\tau)-x^{\prime}(t)],x(0)=x_{0}\text{ and
}x^{\prime}(0)=\dot{x}_{0}.$

is periodic with period $\tau$. The values of the constants $\mu$ and $\tau$
are determined experimentally.

Example 5. Let us consider the following chaotic oscillator
($\delta=0.1,a=b=1,c=0.2,\gamma=0.35$ and $\omega=1.4$) :

$x^{\prime\prime}(t)+0.1x^{\prime}(t)-x(t)+x(t)^{3}+0.2x(t)^{5}=0.35\cos(1.4t)$

The Pincarpñe map (P. Holmes [1979]) is displayed in Figure 9

Now, we introduce the controller :

$x^{\prime\prime}(t)+0.1x^{\prime}(t)-x(t)+x(t)^{3}+0.2x(t)^{5}=0.35\cos(1.4t)+\mu[x^{\prime}(t-\tau)-x^{\prime}(t)].$

See Figures 10 and figure 11 for different values of the parameters $\mu$ and
$\tau$. The optimal values are

$\mu=2.25311\text{ and }\tau=3.73093.$

Figure 1: *

The Pincarpñe map for different values of the parameters $\mu$ and $\tau$

The solution to ther i.v.p.

$x^{\prime\prime}(t)+0.1x^{\prime}(t)-x(t)+x(t)^{3}+0.2x(t)^{5}=0.35\cos(1.4t)+2.25311[x^{\prime}(t-3.73093)-x^{\prime}(t)]\text{,
}x(0)=0\text{ and }x^{\prime}(0)=0$

is periodic with period $T=3.73093$. The Chebyshev approximation for the
periodic solution 0n $0\leq t\leq 3.73093~{}$is given by

$x_{\text{Chebyshev}}(t)=-\frac{3t^{5}}{2038}+\frac{8t^{4}}{337}-\frac{517t^{3}}{4498}+\frac{514t^{2}}{2927}-\frac{14t}{4985}.$

See Figure 12

Figure 2: *

The Chebyshev approximation for the periodic solution 0n $0\leq t\leq 3.73093$

##  6 The corresponding Hamiltonian(un-perturbed)system of Quintic-Quibic
duffing equation

For investigation of our hamiltonian in the presdence of noise Lin, H., & Yim,
S. C. S. [1996] we may consider in this section the following noise
equation(stochaostic differentail equation):

$dX_{t}=a(X_{t},t)dt+b(X_{t},t)dW_{t},\quad X_{0}=x_{0}$ (67)

The Euler–Maruyama states Hakima Bessaih [1999] that SDE (stochaostic
differentail equation) defined in (67) can be approximated recursively by

$X_{t_{i+1}}=X_{t_{i}}+a(X_{t_{i}},t_{i})\Delta t+b(X_{t_{i}},t_{i})\Delta
W_{t_{i}},\quad X_{t_{0}}=x_{0}\mbox{ for }0\leq i\leq N-1\quad(1)$

where $0=t_{0}<t_{1}<\cdots<t_{N}=T$, $\Delta t=T/N$ and $\Delta
W_{t_{i}}=W_{t_{i+1}}-W_{t_{i}}$, The random variables $W_{t_{n}}$ are
independent and identically distributed normal random variables (Zeraoulia
Rafik & Alvaro Salas and D.ocampo [2018]) with expected value zero and
variance $\Delta t$. cubic-quintic stochastic differential equation (Duffing
equation) let

$\ddot{x}+\epsilon\gamma\delta\dot{x}-ax+bx^{3}+cx^{5}=\epsilon\gamma\delta\cos(\omega
t)+\eta[t]$ (68)

Here $\eta[t]$ is Gaussian white noise. The equation (68) describes the
stochastic motion of a particle in a harmonic potential. $\gamma$ is a
modulation factor determining the relative strength of deterministic and
stochastic forcings such that $0<\gamma<1$,The solutions of that problem in
the case of absence of white noise term for some values of $a=1,b=-1$ and
$\omega$ with $c=0$ , Equation(68) can be written as :

$\dot{x}=v\\\ \dot{v}=ax-bx^{3}-\gamma\epsilon
v+A\cos(wt)+\eta(t),A=\gamma\epsilon$ (69)

We may give the correspending Hamiltonian (un-perturbed)system of Quintic-
Quibic duffing equation in the absence of noise and in the presence of it .We
may start with absence of white Gaussian noise ($\eta(t)=0$),The Quintic-
Quibic duffing equation Duffing equation can be integrated upon multiplication
by the velocity:

$\dot{x}\left(\ddot{x}+\omega_{0}^{2}x+\beta\,x^{3}+c\,x^{4}\right)=\frac{\text{d}}{{\text{d}}t}\left(\frac{1}{2}\,\dot{x}^{2}+\frac{1}{2}\,\omega_{0}^{2}x^{2}+\frac{1}{4}\,\beta\,x^{4}+\frac{1}{6}\,c\,x^{6}\right)=0.$
(70)

Integration yields

$E(t)=\frac{1}{2}\,\dot{x}^{2}+\frac{1}{2}\,\omega_{0}^{2}x^{2}+\frac{1}{4}\,\beta\,x^{4}+\frac{1}{6}\,c\,x^{6}=\mbox{constant}.$
(71)

The function in parenthesis
$H=\frac{1}{2}\,\dot{x}^{2}+\frac{1}{2}\,\omega_{0}^{2}x^{2}+\frac{1}{4}\,\beta\,x^{4}+\frac{1}{6}\,c\,x^{6}$
s called the Hamiltonian for Quintic Quibic Duffing equation. Then

$\dot{x}=\frac{\partial H}{\partial y},\qquad\dot{y}=-\frac{\partial
H}{\partial x}.$

For positive coefficients $\omega_{0}^{2}$ and $\beta$ , the solution is
bounded:$|x|\leq\sqrt{2H/\omega_{0}^{2}}$ and $|\dot{x}|\leq\sqrt{2H}$ . When
$\gamma\geq 0$ the function $E(t)$ satisfies:

$\frac{{\text{d}}E(t)}{{\text{d}}t}=-\gamma\,\dot{x}^{2}\leq 0;$

therefore, $E(t)$ is a Lyapunov function, and every trajectory moves on the
surface of $E(t)$ toward the equilibrium position the origin. When the Duffing
equation has nonzero coefficients, there exist stationary solutions that are
obtained upon solving the Quintic equation:

$\omega_{0}^{2}x+\beta\,x^{3}+cx^{5}=0\qquad\mbox{or}\qquad
x\left(\omega_{0}^{2}+\beta\,x^{2}+cx^{4}\right)=0.$ (72)

So we get two other multiple equilibrium solutions
$\left\\{\left\\{x\to-\frac{\sqrt{-\text{cx}^{4}-w^{2}}}{\sqrt{B}}\right\\},\left\\{x\to\frac{\sqrt{-\text{cx}^{4}-w^{2}}}{\sqrt{B}}\right\\}\right\\}$,
To analyze their stability, we apply the linearization procedure, so we
calculate the Jacobian matrix, we may let this for readers. Now , we may
investigate for the correspending Hamiltonian in the presence of white
Gaussian noise.here we may use the correspending coupled system of first order
differential equation which is defined in (69),let
$\dot{x}=\dot{p},\dot{x}=\dot{q}$, such that $p$ and $q$ two states variable ,
The correspending Hamiltonian in the presence of white Gaussian noise can be
obtained using (69) with $a=1,b=c=-1$ by:

$\displaystyle\begin{bmatrix}\dot{p}\\\
\dot{q}\end{bmatrix}=f(q,p)=\begin{bmatrix}\frac{\partial H(q,p)}{\partial
p}\\\ \frac{\partial H(q,p)}{\partial q}\end{bmatrix}=\begin{bmatrix}p\\\
q-q^{3}-q^{5}\end{bmatrix}$ (73)

where $[q,p]^{T}=[x,\dot{x}]^{T}$ and $H(q,p)$ represent the Hamiltonian P.
Holmes [1979], equation (68) for $a=1,b=c=-1$ can be expressed as :

$\displaystyle\begin{bmatrix}\dot{p}\\\
\dot{q}\end{bmatrix}=f(q,p)+h(q,p,t),h(q,p,t)=\begin{bmatrix}0\\\
-cp+A\cos(\omega t)+\eta(t)\end{bmatrix}$ (74)

where $h(q,p,t)$ is the perturbation to the Hamiltonian system

##  7 Conclusion:

Control of chaos remains an area of intensive research. Reliable forecasting
of the dynamics of nonlinear systems with chaotic behaviour Zhu & Leung [1999]
is a challenging task. It can be solved in several ways. For example, by
localizing a chaotic attractor, while obtaining a rough forecasting, or by
introducing control of unstable periodic orbits embedded in a chaotic
attractor, thereby making the behavior of the system predictable for given
values of its parameters .We have proved using Melinikov Method that we can
predict the number of limit cycle around the equilibrium .We were able to
discuss Quibic quantic duffing equation in the stochastics sense (only
investigation for hamiltonian) such that noisly model probably it can be
considered as omnipresent in biomedical systems and signals. Conventional
views assume that its presence is detrimental to systems performance and
accuracy. Hence, various analytic approaches and instrumentation have been
designed to remove noise L. Cohen [2005]. On the contrary, recent
contributions have shown that noise can play a beneficial role in biomedical
systems. The results of this literature review indicate that noise is an
essential part of biomedical systems and often plays a fundamental role in the
performance of these systems. Furthermore R. Benzi, A. Sutera, & A. Vulpiani
[1981], in preliminary work, noise has demonstrated therapeutic potential to
alleviate the effects of various diseases. Further research into the role of
noise and its applications in medicine is likely to lead to novel approaches
to the treatment of diseases and prevention of disability.

Note Added Appendix of this paper were achieved in independent file.

Acknowledgments I would like to thank my friend Alvaro and Lorenzo for their
contribution and their helpful to improve the quality of this paper.And I also
thank refree for their comments .

## References

  * A Elías & Zúñiga [2013] , A Elías & Zúñiga,[2013],“Analytical Solution for the Cubic-Quintic Duffing Oscillator Equation with Physics Applications,”Applied Math. Modelling 37, 2574–2579 (2013)
  * Alvaro Humberto Salas [2022] , Alvaro,[2022],“Exact solution of the cubic-quintic Duffing oscillator,”Complexity Vol. 2022, 9269957, 14 pages (2022).
  * Alvaro Humberto [2022] , Alvaro,[2022],“Analytical Approximant to a Quadratically Damped Duffing Oscillator,”Complexity Vol. 2022, 9269957, 14 pages (2022).Research Article — Open Access Volume 2022 — Article ID 3131253 — https://doi.org/10.1155/2022/3131253
  * A.Lasota [1994] A. Lasota & M. C. Mackey. [1994] “Chaos, Fractals, and Noise.Stochastic Aspects of Dynamics,” IEEE J. Solid-State Circuits , Springer, Berlin p. 1994
  * A. Beléndez. M. L. Alvarez. J. Francés. [2012] A. Beléndez. M. L. Alvarez. J. Francés. S. Bleda. T. Beléndez. & A. Nájera. E. [2012] ,“Analytical Approximate Solutions for the Cubic-Quintic Duffing Oscillator in Terms of Elementary Functions,” J. Appl. Math. 2012 (SI02) 1 J.  , https://doi.org/10.1155/2012/286290
  * Hassan Nayfa [1973] A. Hassan Nayfeh,[1973], ”Perturbation Methods, John Wiley and Sons”, Inc., Hoboken, 1973.
  * B.Palmero& F. Chacon, R. [2022] B.Palmero & F. Chacon, R. [2022] “Suppressing chaos in damped driven systems by non-harmonic excitations, experimental robustness against potential mismatches,” Nonlinear Dynamics  108, 2643–2654,https://doi.org/10.1007/s11071-022-07329-2
  * Cveticanin, L [1993] Cveticanin ,L,[1993] ,“Extension of Melnikov criterion for the differential equation with complex function,” Nonlinear Dynamics 4  4, 139–152,https://doi.org/10.1007/BF00045251
  * E.N. Dudnik, Yu.I. Kuznetsov, I.I. Minakova, Yu.M. & Romanovskiii [1983] E.N. Dudnik, Yu.I. Kuznetsov, I.I. Minakova, Yu.M. & Romanovskiii(1983),‘Synchronization in systems with strange attractors,”Moscow University Physics Bulletin Series 3  ,24 ,84-87.
  * El-Dib, Yusry O., Elgazery, Nasser S., Mady, Amal A. and Alyousef & Haifa A [2022] El-Dib, Yusry O., Elgazery, Nasser S., Mady, Amal A. and Alyousef, Haifa A.[2022], ”On the modeling of a parametric cubic–quintic nonconservative Duffing oscillator via the modified homotopy perturbation method”,Zeitschrift für Naturforschung A, vol. 77, no. 5, pp. 475-486. https://doi.org/10.1515/zna-2021-0354
  * G. Prathap & T., Varadan [1976] G. Prathap, T., Varadan,[1976] “The inelastic large deformation of beams,”Journal of Applied Mechanics  ,43 ,689-690.
  * Gilbert Lewis & Frank Monasa [1982] Gilbert Lewis & Frank Monasa,[1982] “Large deflections of cantilever beams of non-linear materials of the Ludwick type subjected to an end moment,” International Journal of Non-Linear Mechanics, Volume 17, Issue 1,  1, ISSN 0020-7462,https://doi.org/10.1016/0020-7462(82)90032-4
  * H.W. Haslach [1982] H.W. Haslach [1982] “Post-buckling behavior of columns with non-linear constitutive equations,” In ternational Journal of Non-Linear Mechanics 20
  * H. M. Sedighi, K. H. Shirazi & J. Zare [2012] H. M. Sedighi, K. H. Shirazi & J. Zare,[2012] “An Analytic Solution of Transversal Oscillation of Quintic Nonlinear Beam with Homotopy Analysis Method,” JInternational Journal of Nonlinear Mechanics, ,Vol. 47, No. 10, 2012, pp. 777-784.
  * Hakima Bessaih [1999] , Hakima Bessaihi[1999],“Martingale solutions for stochastic Euler equations, Stochastic Analysis and Applications,” Stochastic Analysis and Applications, 17:5, 713-725, DOI: 10.1080/07362999908809631
  * Lin, H., & Yim, S. C. S. [1996] . ”Analysis of a Nonlinear System Lin, H., & Yim, S. C. S. [1996],”Exhibiting Chaotic, Noisy Chaotic, and Random Behaviors.” ,ASME. J. Appl. Mech. June 1996; 63(2): 509–516. https://doi.org/10.1115/1.2788897
  * L. Cohen [2005] L. Cohen[2005], “The history of noise,” IEEE Signal Processing Magazine, vol. 22, no. 6, pp. 20–45, Nov. 2005.
  * J. C. Amazigo [2011] J. C. Amazigo,[2001], ”Perturbation Methods, John Wiley and Sons”, Inc., Hoboken, 1973.
  * Lo, C.C & Gupta, S.D. [1978] Lo, C.C & Gupta, S.D,[1978]“Bending of a Nonlinear Rectangular Beam in Large Deflection,”Journal of Applied Mechanics  ,45 ,213-215
  * Md Abdur Razzak, [2016] Md Abdur Razzak,[2016], “An analytical approximate technique for solving cubic–quintic Duffing oscillator, Alexandria Engineering Journal ”, Volume 55, Issue 3, 2016, ,Pages 2959-2965, ISSN 1110-0168, ,https://doi.org/10.1016/j.aej.2016.04.036.
  * Melnikov [1963] Melnikov,[1963], ”On the stability of the center for time periodic perturbations”, Transac tions of the Moscow Mathematical Society”, 12, (1963) 1-57.
  * P. Holmes [1979] P. Holmes,[1979]“A nonlinear oscillator with a strange attractor,”Philosophical Transac tions of the Royal Society of London Series A  ,292 ,419-448.
  * P. Holmes& J. Marsden [1981] P. Holmes& J. Marsden,[1981]“A partial differential equation with infinitely many periodic orbits: chaotic oscillations of a forced beam,”Archives for Rational,Holmes, Philip and Marsden, Jerrold E. ,292 ,419-448.
  * Pyragas [1992] , Pyragas,[1992], ”Continuous control of chaos by self controlling feedback”, Physics Le tters A 170 (1992) 421-428.
  * Pyragas [1996] , Pyragas,[1996], ”Continuous control of chaos by self controlling feedback”, Academic Press, San Diego, , pp. 118-123.
  * Pyragas [2001] , Pyragas,[2001], ”Control of chaos via an unstable delayed feedback controller”, Physical Review Letters ,86, (2001) 2265-2268.
  * S. Wiggins [1988] S. Wiggins,[1988],“Global Bifurcations and Chaos: Analytical Methods,”Springer, New York, 1988. , Volume 73 , ISBN : 978-1-4612-1041-2
  * A. M. Correig & M. Urquizu [2002] A. M. Correig & M. Urquizu,[2002] “Some Dynamical Aspects of Microseism Time Series,” Geophysical Journal International, Vol. 149, No. 3, 2002, pp. 589-598.
  * M. O. Oyesanya [2008] M. O. Oyesanya,[2008] “Duffing Oscillator as a Model for Predicting Earthquake Occurrence 1,” Journal of Nigerian Association of Mathematical Physics, Vol. 12, 2008, pp. 133-142.
  * R. Benzi, A. Sutera, & A. Vulpiani [1981] , R. Benzi, A. Sutera, & A. Vulpiani[1981],“The mechanism of stochastic resonance,” Journal of Physics A: Mathematical and General, vol. 14, no. 11, pp. L453–L457, Nov. 1981.
  * Zhu & Leung [1999] Zhu, Z. & Leung, H. [1999] “Optimal synchronization of chaotic systems in noise,” IEEE Trans. Circ. Syst.-I$:$ Fund. Th. Appl. 46, 1320–1329.
  * Zeraoulia Rafik & Alvaro Salas and D.ocampo [2018] , Zeraoulia Rafik & Alvaro Salas and D.ocampo [2018], ”A New Special Function and Its Application in Probability”, https://www.hindawi.com/journals/ijmms/2018/5146794/,International Journal of Mathematics and Mathematical Sciences / 2018 / Article