^{1}

^{1}

^{*}

^{2}

^{2}

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.

Duffing oscillators have received remarkable attention in the recent decades due to the variety of their Engineering applications, for example magneto-elastic mechanical system [

Various researchers have used different methods in obtaining solutions, for instance, Ueda [

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 [

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 [

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.

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 [

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;

x ˙ = f ( x ) + ε g ( x ) . (1)

with x ∈ 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 ∈ R 2 . For periodically perturbed planar systems, we have the form;

x ˙ = f ( x ) + ε g ( x , t ) . (2)

where x ∈ R 2 and g is periodic with period t in T. We assume that f ∈ C ′ ( R 2 ) and g ∈ C ′ ( R 2 × R ) and we make the assumption;

1) For ε = 0 , the system (2) has a homoclinic orbit;

Γ 0 : X = γ 0 ( t ) , − ∞ < t < ∞ at a hyperbolic saddle point X 0 and

2) For ε = 0 , the system has a non-parametric family of periodic orbit. Then the Melnikov function M ( t 0 ) is defined as;

M ( t 0 ) = ∫ − ∞ ∞ f ( γ 0 ( t ) ) ^ g ( γ 0 ( t ) , t + t 0 ) 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 ε .

Lemma 2.1

Under assumption 1) and 2), for ε sufficiently small, the system (2) has a unique hyperbolic periodic orbit; γ ε ( t ) = X 0 + 0 ( ε ) of period T. Correspondingly, the Poincare map P ε has a unique hyperbolic fixed point of saddle type;

X ε = X 0 + 0 ( ε ) . (4)

Theorem 2.1

Under the assumption 1) and 2), if the Melnikov function M ( t 0 ) has a simple zero in [0,1], then for all sufficiently small ε ≠ 0 , the stable and unstable manifold of the Poincare map P ε intersect transversally, that is, P ε has a transverse homoclinic point.

This theorem was established by Melnikov [

Theorem 2.2 (Smale-Birkhoff Homoclinic Theorem) [

Let f be a diffeomorphism ( C ′ ) and suppose p is a hyperbolic fixed point. A homoclinic point is a point q ≠ 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 γ ( q ) = q n such that lim n → ∞ q n = lim n → ∞ q n = p . Since the stable and unstable manifolds are invariant, we have; q n ∈ W s ( p ) ∩ W u ( p ) for all n ∈ ℤ . Moreover, if q is transversal, so are all q n since f is diffeomorphism.

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 [

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 [

δ x ˙ i = ∂ u i ∂ x j δ 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 ( t ) = f ′ ( x 0 ) and x ( t ) + δ x ( t ) = f ′ ( x 0 + δ x 0 ) that starts over very close to each other separate exponentially with time. This sensitive initial condition can be quantified as;

‖ δ x ( t ) ‖ = e λ t ‖ δ x 0 ‖ (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;

λ = 1 t ln ‖ δ x ( t ) ‖ ‖ δ x 0 ‖ (6)

λ T ( X ( t ) , δ x 0 ) = 1 T ln ‖ δ X ( t + T ) ‖ ‖ δ x ( t ) ‖ (7)

λ l o c a l ( X ( t ) ) = lim 1 T ln ‖ δ X ( t + T ) ‖ ‖ δ X ( t ) ‖ (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 λ < 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 λ = 0 , the system is conservative and converges to a stable cycle limit. If λ > 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 λ = ∞ , 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 ( t ) and x 2 ( t ) which at some point t 0 are arbitrary close together and their difference in time tracked by the function;

δ x ( t ) = | x 1 ( t ) − x 2 ( t ) | . 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.

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

x ¨ + ε k x ˙ + α x + β x 3 = ε γ cos ( ω t ) (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;

x ˙ = ∂ H ∂ y + ε p ( x , y ) y ˙ = ∂ H ∂ x + ε q ( x , y ) } (10)

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

Let ( x , y ) = x ε ( t ) , y ε ( t ) be the solution of (3.1). Then the solution of the unperturbed system at ( ε = 0 ) is ( x , y ) = ( x 0 ( t ) , y 0 ( t ) ) . Further, we note that the unperturbed system at ε = 0 has one equilibrium point i.e. the center surrounded by a closed trajectories.

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

x ¨ + ε k x ˙ + α x + β x 3 = ε γ cos (ωt)

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

x ˙ = y + ε ( k y + γ cos ω t ) y ˙ = − α x − β x 3 } (11)

where α , β > 0 .

Let ( x , y ) = x ε ( t ) , y ε ( t ) be the solution of (11). The unperturbed system (11) has a Hamiltonian;

H = 1 2 y 2 + α 2 x 2 + β 4 x 4 . (12)

and one equilibrium point surrounded by a closed trajectories.

The solution of the unperturbed system is expressed as;

x 0 ( t ) = 2 a k 1 − k 2 b 1 − 2 k 2 s d ( a 1 − 2 k 2 t , k ) . (13)

where sd is a Jacobian function.

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

M ( t 0 ) = ∫ 0 T 0 ( − k x 0 ( t ) + γ cos ω t ) d t = ∫ 0 T 0 − k x 0 ( t ) d t + ∫ 0 T 0 γ cos ω t d t = − k L 1 + γ L 2 (14)

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

L 1 = ∫ 0 T 0 x 0 ( t ) d t = ∫ 0 T 0 2 a k 1 − k 2 b 1 − 2 k 2 s d ( a 1 − 2 k 2 t , k ) (15)

Then, using the following properties;

T 0 = 4 k ( k 2 ) 1 − 2 k 2 a (A3)

S d ( z , k ) = 0 (A2)

∫ 0 4 k s d 2 n ( z , k ) d z = 4 ∫ 0 k s d 2 n ( z , k ) d z (A3)

The Melnikov function becomes;

M ( t 0 ) = 4 1 − 2 k 2 a [ ∫ 0 k 2 a k 1 − k 2 b 1 − 2 k 2 s d 2 n ( z , k ) ] d z (16)

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

∫ 0 k s d 2 d z = 1 k 2 ( 1 − k 2 )

∫ 0 k s d 4 d z = 1 3 k 4 ( 1 − k 2 )

We obtain the final expression as;

M ( t 0 ) = 8 λ a 1 − 2 m ( 1 − 2 m ) 2

where λ = a / b .

Consider the Duffing equation below;

x ¨ + k x ˙ + α x + β x 3 = γ cos ( ω t ) (17)

where x ¨ and x ˙ are second-order and first-order derivative, α ∈ R n , β ∈ R n , δ is the damping, γ is the amplitude of the circle, ω 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 = x ˙ , g ( x , y ) = − k y − α x − β x 3

Then Equation (1) is equivalent to;

x ˙ = y y ˙ = g ( x , y ) + f ( t ) } (18)

which is written in matrix form as;

Y ˙ ( x ) = F ( Y ) (19)

According to the variational principle, its variational equations are;

Y ˙ ( t ) = J ( t ) Y ( t ) , Y ( 0 ) = I (20)

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

( ∂ f ∂ x ∂ f ∂ y ∂ g ∂ x ∂ g ∂ y ) = ( 0 1 − α − 3 β x 2 − k )

Then, QR factorization of Y ( t ) can be written as;

Y ( t ) = Q R (21)

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

Q ˙ R + Q R ˙ = J Q R . (22)

Q ( 0 ) R ( 0 ) = I

Left multiply Equation (22) by Q T and right multiply by R − 1 , we have;

Q T Q ˙ + R ˙ R − 1 = Q T J Q (23)

Q ( 0 ) = I , R ( 0 ) = 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 θ .

Q = [ cos θ sin θ − sin θ cos θ ]

The upper triangular matrix R can be expressed as;

R = [ e λ 1 ( t ) r 12 0 e λ 2 ( t ) ]

where θ is the angle variable, λ i ( t ) is the value associated with the Lyapunov exponent. Then,

Q T = [ cos θ − sin θ sin θ cos θ ]

R − 1 = [ e − λ 1 ( t ) − r 12 e λ 1 + λ 2 0 e − λ 2 ( t ) ]

Then putting Q T , R − 1 , Q and R into Equation (23), we have;

[ cos θ − sin θ sin θ cos θ ] [ − sin θ cos θ − cos θ − sin θ ] + [ de λ 1 ( t ) d t d r 12 d t 0 de λ 2 ( t ) d t ] [ e − λ 1 ( t ) − r 12 e λ 1 + λ 2 0 e − λ 2 ( t ) ] = [ cos θ − sin θ sin θ cos θ ] [ 0 1 − α − 3 β x 2 − k ] [ cos θ sin θ − sin θ cos θ ]

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

d λ 1 ( t ) d t = ∂ g ∂ y sin 2 θ − 1 2 [ 1 + ∂ g ∂ x ] sin ( 2 θ ) d λ 2 ( t ) d t = ∂ g ∂ y cos 2 θ + 1 2 [ 1 + ∂ g ∂ x ] sin ( 2 θ ) d θ ( t ) d t = − 1 2 ∂ g ∂ y sin 2 θ + sin 2 θ − ∂ g ∂ x cos 2 θ } (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;

d v 1 d t = ∂ g ∂ y d v 2 d t = ∂ g ∂ y cos 2 θ − 1 2 [ 1 + ∂ g ∂ x ] sin ( 2 θ ) d θ ( t ) d t = − 1 2 ∂ g ∂ y sin 2 θ + sin 2 θ − ∂ g ∂ x cos 2 θ } (25)

Then from;

d v 1 d t = d λ 1 d t + d λ 2 d t d v 2 d t = d λ 1 d t − d λ 2 d t } (26)

We obtain;

λ 1 ( t ) = [ v 1 ( t ) + v 2 ( t ) ] 2 λ 2 ( t ) = [ v 1 ( t ) − v 2 ( t ) ] 2 } (27)

The time evolution of the Lyapunov exponent is;

f 1 ( t ) = λ 1 ( t ) t f 2 ( t ) = λ 2 ( t ) t } (28)

Then, the Lyapunov exponent is;

λ 1 = lim t → ∞ λ 1 ( t ) t λ 2 = lim t → ∞ λ 2 ( t ) t } (29)

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

represent our chaotic solutions.

good agreement over the time interval shown. Also, at t = 0.067 , the maximum trajectory is x = 0.998 .

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,

x ¨ + ε k x ˙ + a x + β x 3 = ε y cos (ωt)

ε : = 0.01 , a : = 1 , β : = 0.5 , k : = 0.5 , ω : = 0.1

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

D ( t , X ) : = [ X 1 ε ⋅ y ⋅ cos ( ω ⋅ t ) − ε ⋅ k ⋅ X 1 − a ⋅ X 0 − b ⋅ ( X 0 ) 3 ]

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 : = ( 0 1 ) Vector of initial function values

N : = 1500 Numbers of solution values on [t_{0}, t_{1}]

S : = ( X 0 , t 0 , t 1 , N , D )

t : = S 〈 0 〉 Independent variables values

X 1 : = S 〈 1 〉 First solution function values

X 2 : = S 〈 2 〉 Second solution function values

0 | 1 | 2 | |
---|---|---|---|

0 | 0 | 0 | 0 |

1 | 0.067 | 0.067 | 0.998 |

2 | 0.133 | 0.133 | 0.991 |

3 | 0.2 | 0.199 | 0.98 |

4 | 0.267 | 0.264 | 0.964 |

5 | 0.333 | 0.327 | 0.944 |

6 | 0.4 | 0.389 | 0.918 |

7 | 0.467 | 0.449 | 0.888 |

8 | 0.533 | 0.508 | 0.853 |

9 | 0.6 | 0.563 | 0.812 |

10 | 0.667 | 0.616 | 0.766 |

12 | 0.8 | 0.711 | 0.658 |

13 | 0.867 | 0.753 | 0.596 |

14 | 0.933 | 0.79 | 0.53 |

15 | 1 | 0.823 | … |

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.

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

Eze, E.O., Goodluck, O.C., Ngozi, U.R. and Oko, T.O. (2019) Existence of Chaotic Phenomena in a Single-Well Duffing Oscillator under Parametric Excitations. World Journal of Mechanics, 9, 67-80. https://doi.org/10.4236/wjm.2019.94005