Existence of Chaotic Phenomena in a Single-Well Duffing Oscillator under Parametric Excitations

Everestus Obinwanne Eze¹, Okorafor Chinenye Goodluck¹, Ujumadu Rosary Ngozi², Toochukwu Ogbonnia Oko¹

¹Department of Mathematics, Michael Okpara University of Agriculture, Umudike, Nigeria.
²Department of Mathematics, Chukwuemeka Odimegwu Ojukwu University, Uli, Nigeria.

DOI: 10.4236/wjm.2019.94005                

Abstract

In this paper, the existence of chaotic behavior in the single-well Duffing Oscillator was examined under parametric excitations using Melnikov method and Lyapunov exponents. The minimum and maximum values were obtained and the dynamical behaviors showed the intersections of manifold which was illustrated using the MATCAD software. This extends some results in the literature. Simulation results indicate that the single-well oscillator is sensitive to sinusoidal signals in high frequency cases and with high damping factor, the amplitude of the oscillator was reduced.

Keywords

Single-Well Duffing Oscillator, Melnikov Method, Positive Lyapunov Exponent, Threshold Values, MATCAD Software

1. Introduction

Duffing oscillators have received remarkable attention in the recent decades due to the variety of their Engineering applications, for example magneto-elastic mechanical system [1] , large amplitude oscillator of centrifugal governing system [2] , nonlinear vibration beans and plates [3] [4] and fluid flow induced vibration [5] . It is famous for the existence of chaos behavior in recent decades [6] . In 1979, the chaotic phenomena in Duffing equation had been investigated by Ueda [7] . Chaos occurs when the behavior of the dynamical system is extremely sensitive to initial conditions. In mechanical system, it means a motion which trajectories starting from slightly different initial conditions diverge exponentially.

Various researchers have used different methods in obtaining solutions, for instance, Ueda [7] used the numerical simulation where changes in attractors were obtained under various parameters.

The problem in Duffing-type system still remains a puzzle to so many scientists for instance, suppressing and inducing of chaos, influence of time delay, fractional dynamics [8] - [13] .

Melnikov method and Lyapunov exponents are very significant analytical techniques for determining chaos. The main idea of this method is to measure the distance between the stable and unstable manifolds and if the stable and unstable manifolds intensively intersect once, they will intersect infinite times [6] . Thus, according to Smale-Birkhoff theorem in [1] , it implies the existence of the chaotic behavior in the Smale-horseshoe sense. The Melnikov theory was firstly used to study chaos in Duffing system by Holmes [2] and generalized Melnikov function was developed by Wiggins [13] [14] . This criterion is just the necessary condition for chaos but not sufficient for chaos, therefore, it must be sufficient conditions for the suppression of chaos [15] . The Lyapunov exponent is an important indicator in determining the sensitivity of chaotic behavior which characterizes the average rate of the system in phase space between adjacent tracks of convergence and divergence [16] . Whether the Lyapunov exponent is greater than zero or not is one of the most straight forward criterions to distinguish the chaotic systems [17] . In other to calculate the Lyapunov exponent, some methods of solutions includes [16] , nonlinear adaptive filtering method, QR matrix factorization method and its improvement methods. This paper makes use of two methods, the Melnikov method and improved QR matrix factorization.

The objectives of this paper therefore are to investigate the existence of chaotic behavior in a single-well Duffing oscillator forced by parametric excitations.

The rest of the paper is organized as follows: Section 2, explained the preliminaries to the results, Section 3 gives the main results using the Melnikov method and the calculation of the Lyapunov exponent and Section 4 presents the numerical simulations and finally some conclusions are given in Section 5.

2. Preliminaries

2.1. Melnikov Method

One of the main tools for determining the existence or non-existence of chaos in a perturbed Hamiltonian system is Melnikov. In his theory, the distance between stable and unstable manifolds of the perturbed system were calculated up to the first order term.

Melnikov method is a procedure which gives a bound on the parameters of a system such that chaos is predicted not to occur. The Melnikov method investigate the homoclinic bifurcation in the forced Duffing oscillator system with linear and non-damping. It measures the distance between stable and unstable manifolds in the Poincare section [6] and to preserve the homoclinic loops under a perturbation requires that at $t_0$ , if $M\left(t_0\right)$ , that is the Melnikov function has a simple zero, then a homoclinic bifurcation occurs, implying that the chaotic motion occurs.

2.2. Melnikov Method for Predicting Chaos

Melnikov method gives an analytic tool for establishing the existence of transverse homoclinic points of the Poincare map for a periodic orbit of a perturbed dynamical system of the form;

$\stackrel{\dot{}}{x}=f\left(x\right)+\epsilon g\left(x\right)$ . (1)

with $x\in R^n$ . It can also be used to establish the existence of sub-harmonic periodic orbits of perturbed system of the form in (1). Furthermore, it can be used to show the existence of limit cycles and separatix cycles of perturbed planar system with $x\in R^2$ . For periodically perturbed planar systems, we have the form;

$\stackrel{\dot{}}{x}=f\left(x\right)+\epsilon g\left(x,t\right)$ . (2)

where $x\in R^2$ and g is periodic with period t in T. We assume that $f\in C^{\prime }\left(R^2\right)$ and $g\in C^{\prime }\left(R^2\times R\right)$ and we make the assumption;

1) For $\epsilon =0$ , the system (2) has a homoclinic orbit;

${\Gamma }_0:X={\gamma }_0\left(t\right),-\infty <t<\infty$ at a hyperbolic saddle point $X_0$ and

2) For $\epsilon =0$ , the system has a non-parametric family of periodic orbit. Then the Melnikov function $M\left(t_0\right)$ is defined as;

$M\left(t_0\right)={\displaystyle {\int }_{-\infty }^{\infty }f\left({\gamma }_0\left(t\right)\right)^g\left({\gamma }_0\left(t\right),t+t_0\right)\text{d}t}$ . (3)

The Melnikov method can be interpreted as a derivation in energy from the value on the perturbed separatix. Before stating main result established by Melnikov concerning the existence of transverse homoclinic point of the Poincare section, we need the following lemma and theory which establish the existence of a periodic orbit and hence the existence of the Poincare map with sufficient $\epsilon$ .

Lemma 2.1

Under assumption 1) and 2), for $\epsilon$ sufficiently small, the system (2) has a unique hyperbolic periodic orbit; ${\gamma }_{\epsilon }\left(t\right)=X_0+0_{\left(\epsilon \right)}$ of period T. Correspondingly, the Poincare map $P_{\epsilon }$ has a unique hyperbolic fixed point of saddle type;

$X_{\epsilon }=X_0+0_{\left(\epsilon \right)}$ . (4)

Theorem 2.1

Under the assumption 1) and 2), if the Melnikov function $M\left(t_0\right)$ has a simple zero in [0,1], then for all sufficiently small $\epsilon \ne 0$ , the stable and unstable manifold of the Poincare map $P_{\epsilon }$ intersect transversally, that is, $P_{\epsilon }$ has a transverse homoclinic point.

This theorem was established by Melnikov [1] . The idea of the proof is that $M\left(t_0\right)$ is a measure of the separation of the stable and unstable manifold of the Poincare map. The theory is an important result because it establishes the existence of transverse homoclinic point for $P_{\epsilon }$ . It implies the existence of strange invariant set for some iterate of $P_{\epsilon }$ and the same type of chaotic dynamics for system (2) as for the Smale horseshoe map. Generally, the Melnikov method is very useful for detecting the presence of transverse homoclinic orbits and the occurance of homoclinic bifurcations.

Theorem 2.2 (Smale-Birkhoff Homoclinic Theorem) [18]

Let f be a diffeomorphism ( $C^{\prime }$ ) and suppose p is a hyperbolic fixed point. A homoclinic point is a point $q\ne p$ which is in the stable and unstable manifolds. If the stable and unstable manifolds intersect transversally at q, then q is called transverse. This implies that there is a homoclinic orbit $\gamma \left(q\right)=q_n$ such that ${\lim }_{n\to \infty }{q}_n={\lim }_{n\to \infty }{q}_n=p$ . Since the stable and unstable manifolds are invariant, we have; $q_n\in W^s\left(p\right)\cap W^u\left(p\right)$ for all $n\in \mathbb{Z}$ . Moreover, if q is transversal, so are all $q_n$ since f is diffeomorphism.

2.3. Method of Lyapunov Exponent

The method of Lyapunov exponent serves as a useful tool to qualify chaos. Specifically, Lyapunov exponent measures the rate of convergence or divergence of nearby trajectories [16] [18] . Negative Lyapunov exponents indicate convergence while positive Lyapunov exponents demonstrate divergence and chaos. The magnitude of the Lyapunov exponents is an indicator of the time scale on which chaotic behavior can be predicted or transients for the positive and negative cases respectively [19] .

Physically, the Lyapunov exponent measures average exponential divergence or convergence between trajectories that differ only in having an infinitesimally small difference in their initial condition. The system is said to be chaotic if the trajectories remain within a bounded set of the dynamics. If one considers a ball of points in N-dimensional phase space in which each point follows its own trajectory based upon the system equations of motion over time, the ball of points will collapse to a simple point, will stay a ball or will become ellipsoid in shape [20] . The measure of the rate at which this infinitesimal ball collapse or expands is the Lyapunov exponent. For a system written in the state-space form $\stackrel{\dot{}}{x}=u\left(x\right)$ , small derivation from trajectory can be expressed by the equation

$\delta {\stackrel{\dot{}}{x}}_i=\frac{\partial u_i}{\partial x_j}\delta x_j$ . The maximal Lyapunov exponent is then defined by this equation.

Other useful quantities are the short time Lyapunov exponent and the local Lyapunov exponent. A short time Lyapunov exponent is simply a Lyapunov exponent defined over some finite time interval. The local Lyapunov exponent is a short time Lyapunov exponent in the limit where the time interval approaches zero. Both are dependent on starting points and the short time Lyapunov exponent is also independent on the magnitude of the time interval. If all points in the neighborhood of a trajectory converge towards the same orbit, the attractor is a fixed point or a limit cycle. However, if the attractor is strange, any two trajectories $x\left(t\right)=f^{\prime }\left(x_0\right)$ and $x\left(t\right)+\delta x\left(t\right)=f^{\prime }\left(x_0+\delta x_0\right)$ that starts over very close to each other separate exponentially with time. This sensitive initial condition can be quantified as;

$\Vert \delta x\left(t\right)\Vert ={\text{e}}^{\lambda t}\Vert \delta x_0\Vert$ (5)

where λ, the mean rate of separation of trajectories of the system is called the Lyapunov exponent, which can be estimated for long time t as;

$\lambda =\frac{1}{t}\ln \frac{\Vert \delta x\left(t\right)\Vert }{\Vert \delta x_0\Vert }$ (6)

${\lambda }_T\left(X\left(t\right),\delta x_0\right)=\frac{1}{T}\ln \frac{\Vert \delta X\left(t+T\right)\Vert }{\Vert \delta x\left(t\right)\Vert }$ (7)

${\lambda }_{local}\left(X\left(t\right)\right)=\lim \frac{1}{T}\ln \frac{\Vert \delta X\left(t+T\right)\Vert }{\Vert \delta X\left(t\right)\Vert }$ (8)

Equations (7) and (8) are for short and local Lyapunov exponent. The exponent can be positive or negative but at least one must be positive for an attractor to be classified as chaotic. In particular, if $\lambda <0$ , the system converges to a stable fixed point or periodic orbits. A negative value of the Lyapunov exponent is characteristic of dissipative or non-conservative systems. If $\lambda =0$ , the system is conservative and converges to a stable cycle limit. If $\lambda >0$ , the system is unstable and chaotic. Hence, if the system is chaotic, it will have at least one positive Lyapunov exponent. Thus, the definition of chaotic system is based on a positive Lyapunov exponent. Finally, If $\lambda =\infty$ , the system is random.

Generally, the most used measure of sensitive initial condition is a system characterization by the Lyapunov exponent, which quantifies the rate of separation of infinitesimal close trajectories. For example, consider a one-dimensional system with two trajectories $x_1\left(t\right)$ and $x_2\left(t\right)$ which at some point $t_0$ are arbitrary close together and their difference in time tracked by the function;

$\delta x\left(t\right)=\left|x_1\left(t\right)-x_2\left(t\right)\right|$ . The sign of the lyapunov exponent characterizes whether or not the system is exhibiting chaotic behavior. If the exponent is negative, the system, at least in that set of initial conditions is said to be stable (like trajectories go to like trajectories). A Lyapunov exponent of zero implies an unstable system which is essentially on the edge stable and chaotic. And of course a positive exponent implies the system is chaotic where trajectories exhibit exponential divergence.

3. Results

3.1. The Single-Well Duffing Oscillator

The single-well Duffing equation under parametrical excitation is shown below;

$\ddot{x}+\epsilon k\stackrel{\dot{}}{x}+\alpha x+\beta x^3=\epsilon \gamma \cos \left(\omega t\right)$ (9)

The System (9) has a unique hyperbolic limit cycle. Using the Melnikov theory, an analysis has been performed of the limit circles in oscillator systems described by single-well Duffing equation under perturbation.

Briefly, we describe Melnikov function and the bifurcations in perturbed Hamiltonian system as;

$\begin{array}{l}\stackrel{\dot{}}{x}=\frac{\partial H}{\partial y}+\epsilon p\left(x,y\right)\\ \stackrel{\dot{}}{y}=\frac{\partial H}{\partial x}+\epsilon q\left(x,y\right)\end{array}\}$ (10)

where H the Hamiltonian $H=H\left(x,y\right)$ is the analytic function. Also the perturbation functions $p\left(x,y\right)$ and $q\left(x,y\right)$ are analytic, ε is a small parameter.

Let $\left(x,y\right)=x_{\epsilon }\left(t\right),y_{\epsilon }\left(t\right)$ be the solution of (3.1). Then the solution of the unperturbed system at ( $\epsilon =0$ ) is $\left(x,y\right)=\left(x_0\left(t\right),y_0\left(t\right)\right)$ . Further, we note that the unperturbed system at $\epsilon =0$ has one equilibrium point i.e. the center surrounded by a closed trajectories.

3.2. Melnikov Function for the Perturbed Single-Well Duffing Equation

In this work, the single-well Duffing equation is represented by;

$\ddot{x}+\epsilon k\stackrel{\dot{}}{x}+\alpha x+\beta x^3=\epsilon \gamma \cos (\omega t)$

This equation can be rewritten in the following perturbed Hamiltonian system;

$\begin{array}{l}\stackrel{\dot{}}{x}=y+\epsilon \left(ky+\gamma \cos \omega t\right)\\ \stackrel{\dot{}}{y}=-\alpha x-\beta x^3\end{array}\}$ (11)

where $\alpha ,\beta >0$ .

Let $\left(x,y\right)=x_{\epsilon }\left(t\right),y_{\epsilon }\left(t\right)$ be the solution of (11). The unperturbed system (11) has a Hamiltonian;

$H=\frac{1}{2}{y}^2+\frac{\alpha }{2}{x}^2+\frac{\beta }{4}{x}^4$ . (12)

and one equilibrium point surrounded by a closed trajectories.

The solution of the unperturbed system is expressed as;

$x_0\left(t\right)=\frac{\sqrt{2}\sqrt{ak}\sqrt{1-k^2}}{\sqrt{b}\sqrt{1-2k^2}}sd\left(\frac{\sqrt{a}}{\sqrt{1-2k^2}}t,k\right)$ . (13)

where sd is a Jacobian function.

Then, the Melnikov function for the System (3) is given as;

$\begin{array}{c}M\left(t_0\right)={\displaystyle {\int }_0^{T_0}\left(-k{x}_0\left(t\right)+\gamma \cos \omega t\right)\text{d}t}\\ ={\displaystyle {\int }_0^{T_0}-k{x}_0\left(t\right)\text{d}t}+{\displaystyle {\int }_0^{T_0}\gamma \cos \omega t\text{d}t}\\ =-k{L}_1+\gamma L_2\end{array}$ (14)

now taking (13) into (14), we get;

$\begin{array}{c}{L}_1={\displaystyle {\int }_0^{T_0}{x}_0\left(t\right)\text{d}t}\\ ={\displaystyle {\int }_0^{T_0}\frac{\sqrt{2}\sqrt{ak}\sqrt{1-k^2}}{\sqrt{b}\sqrt{1-2k^2}}sd\left(\frac{\sqrt{a}}{\sqrt{1-2k^2}}t,k\right)}\end{array}$ (15)

Then, using the following properties;

$T_0=4k\left(k^2\right)\frac{\sqrt{1-2k^2}}{\sqrt{a}}$ (A3)

$Sd\left(z,k\right)=0$ (A2)

${\displaystyle {\int }_0^{4k}s{d}^{2n}\left(z,k\right)\text{d}z}=4{\displaystyle {\int }_0^{k}s{d}^{2n}\left(z,k\right)\text{d}z}$ (A3)

The Melnikov function becomes;

$M\left(t_0\right)=\frac{4\sqrt{1-2k^2}}{\sqrt{a}}\left[{\displaystyle {\int }_0^k\frac{\sqrt{2}\sqrt{ak}\sqrt{1-k^2}}{\sqrt{b}\sqrt{1-2k^2}}s{d}^{2n}\left(z,k\right)}\right]\text{d}z$ (16)

after a long calculation and introducing the notation $m=k^2$ and the following identities;

${\displaystyle {\int }_0^{k}s{d}^2\text{d}z}=\frac{1}{k^2\left(1-k^2\right)}$

${\displaystyle {\int }_0^{k}s{d}^4\text{d}z}=\frac{1}{3k^4\left(1-k^2\right)}$

We obtain the final expression as;

$M\left(t_0\right)=\frac{8\lambda \sqrt{a}\sqrt{1-2m}}{{\left(1-2m\right)}^2}$

where $\lambda =a/b$ .

3.3. The Lyapunov Exponent of a Single-Well Duffing Oscillator

Consider the Duffing equation below;

$\ddot{x}+k\stackrel{\dot{}}{x}+\alpha x+\beta x^3=\gamma \cos \left(\omega t\right)$ (17)

where $\ddot{x}$ and $\stackrel{\dot{}}{x}$ are second-order and first-order derivative, $\alpha \in {\mathcal{R}}^n$ , $\beta \in {\mathcal{R}}^n$ , $\delta$ is the damping, $\gamma$ is the amplitude of the circle, $\omega$ is the angular frequency of the driven circle. In other to Type equation here, determine whether the system is in chaotic state, we need to calculate the Lyapunov exponent using the QR factorization method.

Let $y=\stackrel{\dot{}}{x}$ , $g\left(x,y\right)=-ky-\alpha x-\beta x^3$

Then Equation (1) is equivalent to;

$\begin{array}{l}\stackrel{\dot{}}{x}=y\\ \stackrel{\dot{}}{y}=g\left(x,y\right)+f\left(t\right)\end{array}\}$ (18)

which is written in matrix form as;

$\stackrel{\dot{}}{Y}\left(x\right)=F\left(Y\right)$ (19)

According to the variational principle, its variational equations are;

$\stackrel{\dot{}}{Y}\left(t\right)=J\left(t\right)Y\left(t\right),Y\left(0\right)=I$ (20)

where $Y\left(t\right)$ is a 2 by 2 matrix, I is a 2 by 2 unit matrix, $J\left(t\right)$ is the Jacobian matrix of the system and its expression is;

$\left(\begin{array}{cc}\frac{\partial f}{\partial x}& \frac{\partial f}{\partial y}\\ \frac{\partial g}{\partial x}& \frac{\partial g}{\partial y}\end{array}\right)=\left(\begin{array}{cc}0& 1\\ -\alpha -3\beta x^2& -k\end{array}\right)$

Then, QR factorization of $Y\left(t\right)$ can be written as;

$Y\left(t\right)=QR$ (21)

where Q is orthogonal matrix, R is upper triangular matrix. Substituting (21) in (20), we obtain the variational equation;

$\stackrel{\dot{}}{Q}R+Q\stackrel{\dot{}}{R}=JQR$ . (22)

$Q\left(0\right)R\left(0\right)=I$

Left multiply Equation (22) by $Q^{\text{T}}$ and right multiply by $R^{-1}$ , we have;

$Q^{\text{T}}\stackrel{\dot{}}{Q}+\stackrel{\dot{}}{R}{R}^{-1}=Q^{\text{T}}JQ$ (23)

$Q\left(0\right)=I,R\left(0\right)=I$

The orthogonal matrix Q is written as a function of angle variables. To the Duffing equation, its orthogonal matrix Q can be expressed by one angle $\theta$ .

$Q=\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]$

The upper triangular matrix R can be expressed as;

$R=\left[\begin{array}{cc}{\text{e}}^{{\lambda }_1\left(t\right)}& r_{12}\\ 0& {\text{e}}^{{\lambda }_2\left(t\right)}\end{array}\right]$

where $\theta$ is the angle variable, ${\lambda }_i\left(t\right)$ is the value associated with the Lyapunov exponent. Then,

$Q^{\text{T}}=\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right]$

$R^{-1}=\left[\begin{array}{cc}{\text{e}}^{-{\lambda }_1\left(t\right)}& \frac{-r_{12}}{{\text{e}}^{{\lambda }_1+{\lambda }_2}}\\ 0& {\text{e}}^{-{\lambda }_2\left(t\right)}\end{array}\right]$

Then putting $Q^{\text{T}},R^{-1},Q$ and R into Equation (23), we have;

$\begin{array}{l}\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right]\left[\begin{array}{cc}-\sin \theta & \cos \theta \\ -\cos \theta & -\sin \theta \end{array}\right]+\left[\begin{array}{cc}\frac{{\text{de}}^{{\lambda }_1\left(t\right)}}{\text{d}t}& \frac{\text{d}{r}_{12}}{\text{d}t}\\ 0& \frac{{\text{de}}^{{\lambda }_2\left(t\right)}}{\text{d}t}\end{array}\right]\left[\begin{array}{cc}{\text{e}}^{-{\lambda }_1\left(t\right)}& \frac{-r_{12}}{{\text{e}}^{{\lambda }_1+{\lambda }_2}}\\ 0& {\text{e}}^{-{\lambda }_2\left(t\right)}\end{array}\right]\\ =\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right]\left[\begin{array}{cc}0& 1\\ -\alpha -3\beta x^2& -k\end{array}\right]\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]\end{array}$

The correspondent matrix elements on both sides of (23) are equal, so we get;

$\begin{array}{l}\frac{\text{d}{\lambda }_1\left(t\right)}{\text{d}t}=\frac{\partial g}{\partial y}{\sin }^2\theta -\frac{1}{2}\left[1+\frac{\partial g}{\partial x}\right]\sin \left(2\theta \right)\\ \frac{\text{d}{\lambda }_2\left(t\right)}{\text{d}t}=\frac{\partial g}{\partial y}{\cos }^2\theta +\frac{1}{2}\left[1+\frac{\partial g}{\partial x}\right]\sin \left(2\theta \right)\\ \frac{\text{d}\theta \left(t\right)}{\text{d}t}=-\frac{1}{2}\frac{\partial g}{\partial y}\sin 2\theta +{\sin }^2\theta -\frac{\partial g}{\partial x}{\cos }^2\theta \end{array}\}$ (24)

We add and subtract the first two differential equations and get a new differential equation. Together with the third differential equation, we obtain three new equations;

$\begin{array}{l}\frac{\text{d}{v}_1}{\text{d}t}=\frac{\partial g}{\partial y}\\ \frac{\text{d}{v}_2}{\text{d}t}=\frac{\partial g}{\partial y}\cos 2\theta -\frac{1}{2}\left[1+\frac{\partial g}{\partial x}\right]\sin \left(2\theta \right)\\ \frac{\text{d}\theta \left(t\right)}{\text{d}t}=-\frac{1}{2}\frac{\partial g}{\partial y}\sin 2\theta +{\sin }^2\theta -\frac{\partial g}{\partial x}{\cos }^2\theta \end{array}\}$ (25)

Then from;

$\begin{array}{l}\frac{\text{d}{v}_1}{\text{d}t}=\frac{\text{d}{\lambda }_1}{\text{d}t}+\frac{\text{d}{\lambda }_2}{\text{d}t}\\ \frac{\text{d}{v}_2}{\text{d}t}=\frac{\text{d}{\lambda }_1}{\text{d}t}-\frac{\text{d}{\lambda }_2}{\text{d}t}\end{array}\}$ (26)

We obtain;

$\begin{array}{l}{\lambda }_1\left(t\right)=\frac{\left[v_1\left(t\right)+v_2\left(t\right)\right]}{2}\\ {\lambda }_2\left(t\right)=\frac{\left[v_1\left(t\right)-v_2(t)\right]}{2}\end{array}\}$ (27)

The time evolution of the Lyapunov exponent is;

$\begin{array}{l}{f}_1\left(t\right)=\frac{{\lambda }_1\left(t\right)}{t}\\ f_2\left(t\right)=\frac{{\lambda }_2\left(t\right)}{t}\end{array}\}$ (28)

Then, the Lyapunov exponent is;

$\begin{array}{l}{\lambda }_1={\lim }_{t\to \infty }\frac{{\lambda }_1\left(t\right)}{t}\\ {\lambda }_2={\lim }_{t\to \infty }\frac{{\lambda }_2\left(t\right)}{t}\end{array}\}$ (29)

4. Numerical Simulation of Single-Well Duffing Oscillator

In this section, we compare the numerical solution of Equation (9) using MATCAD simulation. In Figures 1-6, the trajectory versus time response curves are plotted for different sets of parameter values noted in the figure captions. In all figures, the solid lines represent the numerical solution and the dashed lines

Figure 1. Trajectory-time response curves and parameter values of $\epsilon =0.01$ , $\alpha =1$ , $\beta =0.5$ , $k=0.5$ , $\omega =0.1$ . The solid lines represent the numerical solution and the dash lines represent the chaotic behavior.

Figure 2. Velocity-time response curves and parameter values of $\epsilon =0.01$ , $\alpha =1$ , $\beta =0.5$ , $k=0.5$ , $\omega =0.1$ . The solid lines represent the numerical solution and the dash lines represent the chaotic behavior.

Figure 3. The phase portrait orbits in the chaotic state at $k=0.5$ . $\epsilon =0.01$ , $\alpha =1$ , $\beta =0.5$ , $k=0.5$ , $\omega =0.1$ .

represent our chaotic solutions.

Figure 1 and Figure 2 compare solutions by considering a strong nonlinearity value of $\epsilon =0.01$ . The periodic solution of the Duffing’s equation were shown by the relationship between the first solution function values and the independent variables values as shown in Table 1. The values were generated using the vector initial function values and the constant. However, the solutions are in

Figure 4. Trajectory-time response curves and parameter values of $\epsilon =0.01$ , $\alpha =1$ , $\beta =0.5$ , $k=0.5$ , $\omega =0.1$ . The solid lines represent the numerical solution and the dash lines represent the chaotic behavior.

Figure 5. Velocity-time response curves and parameter values of $\epsilon =0.01$ , $\alpha =1$ , $\beta =0.5$ , $k=0.5$ , $\omega =0.1$ . The solid lines represent the numerical solution and the dash lines represent the chaotic behavior.

Figure 6. The phase portrait orbits in the chaotic state at $k=2$ .

good agreement over the time interval shown. Also, at $t=0.067$ , the maximum trajectory is $x=0.998$ . Figure 3 is the phase diagram of the chaotic system at $\epsilon =0.01$ and damping factor $k=1$ .

The same conclusion can be drawn from Figures 4-6 but with damping factor at $k=2$ . The solutions are in excellent agreement over the time interval shown. However, Figure 6 is the phase portrait of the chaotic system.

$\ddot{x}+\epsilon k\stackrel{\dot{}}{x}+ax+\beta x^3=\epsilon y\cos (\omega t)$

$\epsilon :=0.01,a:=1,\beta :=0.5,k:=0.5,\omega :=0.1$

Define a function that determines a vector of derivatives values at any solution point $\left(t,Y\right)$ :

$D\left(t,X\right):=\left[\begin{array}{c}{X}_1\\ \epsilon \cdot y\cdot \cos \left(\omega \cdot t\right)-\epsilon \cdot k\cdot X_1-a\cdot X_0-b\cdot {\left(X_0\right)}^3\end{array}\right]$

Define an additional argument for the ODE solver:

$t_0:=0$ Initial value of independent variable

$t_1:=100$ Final value of independent variable

$X_0:=\left(\begin{array}{c}0\\ 1\end{array}\right)$ Vector of initial function values

$N:=1500$ Numbers of solution values on [t₀, t₁]

$S:=\left(X_0,t_0,t_1,N,D\right)$

$t:=S^{\langle 0\rangle }$ Independent variables values

$X_1:=S^{\langle 1\rangle }$ First solution function values

$X_2:=S^{\langle 2\rangle }$ Second solution function values

Table 1. Solution matrix table for solution functions.

5. Conclusions

In the present study, the chaotic behavior in single-well Duffing oscillator is investigated using Melnikov approach and Lyapunov exponent. The distance between the stable and unstable manifold of the nonlinear system is calculated by Melnikov approach. The Lyapunov exponent of the nonlinear system is evaluated by QR factorization to determine whether the chaotic phenomenon of the nonlinear system actually occurs.

As a result, threshold values were obtained and the dynamical behaviors showing the intersections of manifold were illustrated. To detect the chaotic phenomena of the nonlinear system, the Melnikov approach, Lyapunov exponent, the time history, phase portrait of the nonlinear system were presented for various cases.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Guckenheimer, J. and Holmes, P. (1983) Global Bifurcations. In: Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, New York, 289-352. https://doi.org/10.1007/978-1-4612-1140-2_6
[2]	Holmes, P.J. (1979) A Nonlinear Oscillator with a Strange Attractor. Philosophical Transaction of the Royal Society of London, Series A: Mathematical and Physical Science, 292, 419-448. https://doi.org/10.1098/rsta.1979.0068
[3]	Kim, S. and Kim, Y. (2000) Dynamic Stabilization in the Double-Well Duffing Oscillator. Physical Review, 61, 6517-6520. https://doi.org/10.1103/physreve.61.6517
[4]	Venkatesan, A. and Lakshmanan, M. (1997) Bifurcation and Chaos in the Double-Well Duffing-van der Pol Oscillator: Numerical and Analytical Studies. Physical Review, 56, 6321-6330. https://doi.org/10.1103/PhysRevE.56.6321
[5]	McCue, L. and Troesch, A.W. (2011) Use of Lyapunov Exponents to Predict Chaotic Vessel Motions. In: Almeida Santos Neves, M., Belenky, V., de Kat, J., Spyrou, K. and Umeda, N., Eds., Contemporary Ideas on Ship Stability and Capsizing in Waves, Springer, Dordrecht, 415-432. https://doi.org/10.1007/978-94-007-1482-3_23
[6]	Zhang, Y., Yue, X., Du, L., Wang, L. and Fang, T. (2016) Generation and Evolution of Chaos in Double-Well Duffing Oscillator under Parametrical Excitation. Journal of Shock and Vibration, 2016, Article ID: 6109062. https://doi.org/10.1155/2016/6109062
[7]	Ueda, Y. (1979) Randomly Transitional Phenomena in the System Governed by Duffing’s Equation. Journal of Statistical Physics, 20, 181-196. https://doi.org/10.1007/BF01011512
[8]	Sun, Z.-K., Xu, W., Yang, X.-L. and Fang, T. (2006) Inducing or Suppressing Chaos in a Double-Well Duffing Oscillator by Time Delay Feedback. Chaos, Solitons and Fractals, 27, 705-714. https://doi.org/10.1016/j.chaos.2005.04.041
[9]	Liu, D., Li, J. and Xu, Y. (2014) Principle Resonance Responses of SDOF Systems with Small Fractional Derivative Damping under Narrow-Band Random Parametric Excitation. Communication in Nonlinear Science and Numerical Simulation, 19, 3642-3652. https://doi.org/10.1016/j.cnsns.2014.03.018
[10]	Syta, A., Litak, G., Lenci, S. and Scheffler, M. (2014) Chaotic Vibrations of the Duffing System with Fractional Damping. Chaos, 24, Article ID: 013107. https://doi.org/10.1063/1.4861942
[11]	Tehrani, M.G., Wilmshurst, L. and Elliot, S.J. (2013) Reacceptance Method for Active Vibration Control of a Nonlinear System. Journal of Sound and Vibration, 332, 4440-4449. https://doi.org/10.1016/j.jsv.2013.04.002
[12]	Zhao, X.R., Xu, W., Yang, Y.G. and Wang, X.Y. (2016) Stochastic Responses of a Viscoelastic-Impact System under Additive and Multiplicative Random Excitations. Communications in Nonlinear Science and Numerical Simulation, 35, 166-176. https://doi.org/10.1016/j.cnsns.2015.11.008
[13]	Wiggins, S. (1988) Global Bifurcation and Chaos: Analytical Methods. Springer, New York. https://doi.org/10.1007/978-1-4612-1042-9
[14]	Wiggins, S. (1990) Introduction to Applied Nonlinear Dynamical Systems and Chaos. Springer, New York. https://doi.org/10.1007/978-1-4757-4067-7
[15]	Garcia-Margallo, J. and Bejarano, J.D. (1998) Melnikov’s Method for Non-Linear Oscillators with Non-Linear Excitations. Journal of Sound and Vibration, 212, 311-319. https://doi.org/10.1006/jsvi.1997.1443
[16]	Wolf, A. (1986) Quantifying Chaos with Lyapunov Exponents. Chaos. Chapter, 13, Princeton University Press, Princeton, NJ. https://doi.org/10.1515/9781400858156.273
[17]	Zeni, A.Z. and Gallas, A.C. (1995) Lyapunov Exponent for a Duffing Oscillator. Physica D: Nonlinear Phenomena, 89, 71-82. https://doi.org/10.1016/0167-2789(95)00215-4
[18]	Teschl, G. (2012) Ordinary Differential Equations and Dynamical Systems. American Mathematical Society, Providence, Rhode Island.
[19]	Haken, H. (1981) Chaos and Order in Nature. Proceedings of the International Symposium on Synergetics at Schloβ Elmau, Bavaria, Germany, April 27-May 2, 1981, 2-11. https://doi.org/10.1007/978-3-642-68304-6
[20]	Glass, L. and Mackey, M.C. (1998) From Clock to Chaos: The Rhythms of Life. Princeton University Press, Princeton, NJ.