# Heteroclinic Bifurcation Behaviors of a Duffing Oscillator with Delayed Feedback

Authors: Shao-Fang Wen; Ju-Feng Chen; Shu-Qi Guo
Journal: Shock and Vibration
 (2018)

---

## Abstract

The heteroclinic bifurcation and chaos of a Duffing oscillator with forcing excitation under both delayed displacement feedback and delayed velocity feedback are studied by Melnikov method. The Melnikov function is analytically established to detect the necessary conditions for generating chaos. Through the analysis of the analytical necessary conditions, we find that the influences of the delayed displacement feedback and delayed velocity feedback are separable. Then the influences of the displacement and velocity feedback parameters on heteroclinic bifurcation and threshold value of chaotic motion are investigated individually. In order to verify the correctness of the analytical conditions, the Duffing oscillator is also investigated by numerical iterative method. The bifurcation curves and the largest Lyapunov exponents are provided and compared. From the analysis of the numerical simulation results, it could be found that two types of period-doubling bifurcations occur in the Duffing oscillator, so that there are two paths leading to the chaos in this oscillator. The typical dynamical responses, including time histories, phase portraits, and Poincare maps, are all carried out to verify the conclusions. The results reveal some new phenomena, which is useful to design or control this kind of system.


---

## Body

## 1. Introduction

It is well known that many kinds of nonlinearities exist in engineering systems, such as parametric excitation, nonsmoothness, time delay, discontinuity, and large deformation [1–6]. These nonlinearities will make the system responses more complicated and the system performance deteriorated. In recent years, the complicated dynamics of typical nonlinear systems have been more and more concerned. Duffing oscillator is one of the most common and typical models in nonlinear dynamical systems. Among these nonlinear systems the well-known Duffing oscillator is quite suitable to model the large deformation structure in many physics and engineering fields [7–10]. Moreover, time delay is an unavoidable and very common problem when those systems were controlled. It can influence the system dynamical characteristic and even destroy the system stability [11, 12]. The complicated Duffing nonlinear system with time delay is possible to generate more complex bifurcation and chaotic dynamic phenomena [11–14]. For example, Luo [15] studied the bifurcation trees of time-delay Duffing oscillator by a semianalytical method. Amer et al. [16] investigated the Duffing oscillator with parametric excitation under time-delay feedback based on the multiple scales perturbation method and analyzed the influences of the system parameters. Ji and Leung [17] studied the bifurcation control of a Duffing oscillator with parametrical excitation by a linear time-delayed feedback and discovered that the stable region could be broadened through choosing an appropriate feedback control. Zhang et al. [18] investigated the multipulse global bifurcations and chaotic behaviors of a cantilever beam by extended Melnikov method. Theodossiades and Natsiavas [19] and Van Dooren [20] studied the period-doubling bifurcations for gear-pair systems with periodic stiffness and backlash. Rusinek et al. [21, 22] investigated the dynamics of the time-delay Duffing oscillator. Besides, they also studied the chaos and feedback control of the time-delay Duffing system and found the suitable time delay and feedback gain would destroy the chaotic attractor for Duffing system. Nana Nbendjo et al. [23] studied the effects of the control parameters in a double-well Duffing oscillator with the time delay by Melnikov theory.The heteroclinic bifurcation and chaos behaviors are two of the most important characteristics in nonlinear dynamical systems. Many scholars paid attention to the heteroclinic bifurcation and chaos phenomena of nonlinear Duffing oscillator [24–26]. For example, Yang and Sun et al. [27, 28] investigated the necessary condition for the generating chaos of a double-well Duffing oscillator with bounded noise excitations and time-delay feedback by Melnikov theory. Siewe et al. [29] investigated the necessary condition for chaotic behavior of Duffing-Rayleigh system based on Melnikov theory. Cao et al. [30] presented a novel construction of homoclinic and heteroclinic orbits for nonlinear systems with a perturbation-incremental method. Chen and Yan [31] studied the heteroclinic bifurcation behavior in the Duffing-VDP oscillator by the hyperbolic Lindstedt-Poincare method and obtained the analytical heteroclinic solution. Chacon [32] and Maki et al. [33] investigated heteroclinic bifurcation phenomenon for different types of nonlinear systems based on Melnikov method, respectively. Lei and Zhang [34] investigated the chaotic motion of the Duffing system and came to the conclusion that the threshold for generating chaos could be changed by choosing the internal parameters of trichotomous noise.All the above analyses mainly focused on qualitative, numerical, or simplified analytical solutions about the necessary condition for generating chaos. In this paper, the first-order exact analytical solution of the necessary condition for generating chaos in sense of Smale horseshoes in a Duffing oscillator with both delayed displacement feedback and delayed velocity feedback is obtained based on Melnikov theory. Besides, the numerical heteroclinic bifurcation results were presented and some new phenomena were found in the Duffing oscillator with time-delay feedback. The basic structure of this paper is arranged as follows. The Melnikov function is obtained based on Melnikov method, and the analytical necessary condition for generating chaos is also obtained in Section2. In Section 3, the bifurcation curves and the largest Lyapunov exponents by numerical method are investigated. It can be found that there are two paths leading to chaos via period-doubling bifurcation in this system. Then the time histories, phase portraits, and Poincare maps with the typical system parameters are all presented to verify the new phenomenon. The influences of the delayed displacement and velocity feedback parameters on generating chaos are analyzed. The comparisons of the numerical results with the analytical results are also fulfilled in this section. Finally, the conclusions of this paper are summarized.


## 2. Analytical Necessary Condition for Chaos in Sense of Smale Horseshoes

Duffing oscillator is one of the most familiar systems in nonlinear dynamics. Under the function of both delayed displacement feedback and delayed velocity feedback, a Duffing oscillator with forcing excitation would be investigated in this section. The dynamic equation is(1)x¨+kx-αx3+cx˙=fcos⁡ωt+uxt-τ1+vx˙t-τ2,where k is the linear stiffness coefficient, α is nonlinear stiffness coefficient, c is linear damping coefficient, f and ω are excitation amplitude and frequency, respectively, u and v are displacement and velocity feedback coefficients, respectively, and τ1 and τ2 are time delays of displacement and velocity feedback, respectively. Here all the system parameters are positive and dimensionless.Introducing the transformation(2)c=εc1,f=εf1,u=εu1,v=εv1,where ε is a small real parameter, (1) turns into(3)x¨+kx-αx3=εf1cos⁡ωt+u1xt-τ1+v1x˙t-τ2-c1x˙.Supposing ε=0, the unperturbed system is (4)x¨+kx-αx3=0.There are three equilibrium points, where (0,0) is a center, and ±k/α,0 are two saddle points. Generally speaking, if the unperturbed system is conservative system and the number of saddle points is 1 larger than that of the centers, there may be heteroclinic orbit in this system. Here the heteroclinic orbit connecting the two saddle points satisfies the formula(5)12x˙2+k2x2-α4x4=k24α.Supposing x=0 at t=0, the calculating result is(6)x˙0=±k2α.Integrating (5), one could obtain(7)∫0xdx±k2/2α-kx2+α/2x4=t.Calculating (7), it yields(8)x±t=±kαtanh⁡k2t.So (3) can be rewritten as(9a)x⃑˙=f⃑x⃑+εg⃑x⃑,t,where(9b)x⃑=x1x2=xx˙,(9c)f⃑x⃑=x2-kx1+αx13,(9d)g⃑x⃑,t=ε0-c1x2+f1cos⁡ωt+u1x1t-τ1+v1x2t-τ2.Then one can establish the heteroclinic orbit as follows:(10a)x1±t=±kαtanh⁡k2t,(10b)x2±t=±k2αsech2⁡k2t.Melnikov theory [35], proposed by Melnikov, is a perturbation method originally for proving the existence of transverse homoclinic or heteroclinic orbits. Melnikov theory has been applied in many researches. The Melnikov technique was firstly applied to study a periodically driven Duffing oscillator by Holmes [36]. The Melnikov technique is also valid in time-delay systems, which was proved in [37]. According to Melnikov theory, the perturbed function must be a periodicity of function. In (9a), the perturbed function g⃑(x⃑,t) with time delay is also a periodic function which is satisfied with the basic condition of Melnikov theory. Therefore, many scholars have studied the dynamical characteristics of the time-delay Duffing system by Melnikov theory such as [23, 27, 28, 37]. All the existing analyses were mainly focused on qualitative, numerical, or simplified analytical solutions about the necessary condition for generating chaos. Here, Melnikov method is also applied to yield the Melnikov function as follows:(11a)Mt0=f⃑x⃑∧g⃑x⃑,t=ε∫-∞∞x2±t-c1x2±t+f1cos⁡ωt+t0+u1x1±t-τ1+v1x2±t-τ2dt=M1t0+M2t0+M3t0+M4t0,where(11b)M1t0=-εc1∫-∞∞x2±t2dt=-22εc1kk3α,11c1M2t0=∓εkf12α∫-∞∞sech2⁡k2tcos⁡ωt+t0dt.Making use of the odevity of integrand, 11c1 becomes11c2M2t0=∓εkf12αcos⁡ωt0∫-∞∞sech2⁡k2tcos⁡ωtdt=∓2εf1πω2αcos⁡ωt0csch⁡πω2k.And the other parts in (11a) are (11d)M3t0=εu1∫-∞∞x2±tx1±t-τ1dt=εkku12α∫-∞∞sech2⁡k2ttanh⁡k2t-τ1dt=2εku12αcsch2⁡k2τ12kτ1-sinh⁡2kτ1,(11e)M4t0=εv1∫-∞∞x2±tx2±t-τ2dt=εk2v12α∫-∞∞sech2⁡k2tsech2⁡k2t-τ2dt=2εkkv1αcsch3⁡k2τ22kτ2cosh⁡k2τ2-2sinh⁡k2τ2.Through the analysis of the above results, the necessary condition for generating chaos in sense of Smale horseshoes can be obtained as follows:(12)2εf1πω2αcsch⁡πω2k>-22εc1kk3α+2εku12αcsch2⁡k2τ12kτ1-sinh⁡2kτ1+2εkkv1αcsch3⁡k2τ22kτ2cosh⁡k2τ2-2sinh⁡k2τ2.Replacing the parameters in (12) with the original system parameters, one can get(13)2fπω2αcsch⁡πω2k>-22ckk3α+2ku2αcsch2⁡k2τ12kτ1-sinh⁡2kτ1+2kkvαcsch3⁡k2τ22kτ2cosh⁡k2τ2-2sinh⁡k2τ2.


## 3. Numerical Simulation and Influences of Delayed Feedback on Bifurcation and Chaotic Behaviors

If the delayed displacement feedback coefficientu and velocity feedback coefficient v are all small enough, then (13) becomes(14)2fπω2αcsch⁡πω2k>22ckk3α-2ku2αcsch2⁡k2τ12kτ1-sinh⁡2kτ1-2kkvαcsch3⁡k2τ22kτ2cosh⁡k2τ2-2sinh⁡k2τ2.Through the analysis of (14), we could find that the analytical necessary condition for generating chaos is influenced by the delayed displacement feedback and delayed velocity feedback, respectively. That is to say, the coupling relationship does not exist between delayed displacement feedback and delayed velocity feedback. According to the above analysis of the relationship between the two kinds of feedback, the influences of the two kinds of delayed feedback on bifurcation and chaos behaviors will be investigated separately in the following sections.


### 3.1. The Influences of Delayed Displacement Feedback Parameters

If there is only the delayed displacement feedback in (1), (14) becomes(15)2fπω2αcsch⁡πω2k>22ckk3α-2ku2αcsch2⁡k2τ12kτ1-sinh⁡2kτ1.


#### 3.1.1. The Influence of Feedback Coefficientu

By the analysis of (15), we can find(16)csch2⁡k2τ1>0,limτ1→0⁡2kτ1-sinh⁡2kτ1=0,2kτ1-sinh⁡2kτ1<0,so that(17)-2ku2αcsch2⁡k2τ12kτ1-sinh⁡2kτ1>0.It could be found that the right-hand parts of (15) would become larger with the increasing of the delayed displacement feedback coefficient u. Accordingly, the critical excitation amplitude fmin for generating chaos will become larger with the increasing of u. In other words, the threshold value for chaotic motion in (1) will become larger, so that it is harder to generate chaos with larger u in (1).In order to verify the validity of the analytical necessary condition for generating chaos, a set of illustrative system parameters is chosen ask=1, c=0.4, α=1, ω=0.8, and τ1=0.5. The numerical solutions for fmin of (1) are also investigated by numerical iterative method. When u=0.1, the bifurcation diagram and the corresponding largest Lyapunov exponents are shown in Figures 1(a) and 1(b), respectively. The numerical simulation method for the largest Lyapunov exponents about delayed differential equation can be found in [38]. From the analysis of Figures 1(a) and 1(b), it could be found that the fmin values of generating chaos in two figures are consistent, which implies that the results by numerical iterative method are correct. Through the numerical simulation, the authors also observe some special dynamical phenomenon in this complex system. We find that there exist two paths leading to chaos via different period-doubling bifurcation in (1) which are shown in Figures 2(a) and 3(a), respectively. If the two paths are drawn on the same diagram which is shown in Figure 1, it could be found that the intersection point of the two paths is f=0.3774. When f<0.3774, there is a single periodic solution. If f>0.3774, there will be two periodic-1 solutions, which depend on the initial values. The path shown in Figure 2(a) is simply named as “Type 1” and the other path shown in Figure 3(a) is “Type 2.” In order to verify the conclusion, the time history, phase portrait, and Poincare maps at different typical points are presented and analyzed in the following part.Figure 1
The system responses for different excitation amplitudef.
(a)
Bifurcation diagram
(b)
The largest Lyapunov exponentsFigure 2
(a)
Bifurcation diagram for Type 1
(b)
The system response atf=0.39 for Type 1. (b1) Time history; (b2) phase portrait; (b3) Poincare maps
(c)
The system response atf=0.402 for Type 1. (c1) Time history; (c2) phase portrait; (c3) Poincare maps
(d)
The system response atf=0.4044 for Type 1. (d1) Time history; (d2) phase portrait; (d3) Poincare mapsFigure 3
(a)
Bifurcation diagram for Type 2
(b)
The system response atf=0.39 for Type 2. (b1) Time history; (b2) phase portrait; (b3) Poincare maps
(c)
The system response atf=0.402 for Type 2. (b1) Time history; (b2) phase portrait; (b3) Poincare maps
(d)
The system response atf=0.404 for Type 2. (b1) Time history; (b2) phase portrait; (b3) Poincare mapsThe transitions observed in Type 1 occur atf≈0.3998 (1P→2P) and 0.4040 (2P→4P), respectively. Some typical samples in different ranges are taken to analyze the dynamical phenomenon.(1)
Whenf=0.39 and the initial value [x0,x˙0]=[-0.3956,0.5858], the time history, phase portrait, and Poincare maps are shown in Figures 2(b1), 2(b2), and 2(b3), respectively. It could be found that it is period-1 motion in this case.(2)
Whenf=0.402, the time history, phase portrait, and Poincare maps are shown in Figures 2(c1), 2(c2), and 2(c3), respectively. It is shown the dynamical motion is period-2 in this case.(3)
Whenf=0.4044, the time history, phase portrait, and Poincare maps are shown in Figures 2(d1), 2(d2), and 2(d3), respectively. It could be found that it is period-4 motion in this case.The transitions observed in Type 2 occur atf≈0.4 (1P→2P) and 0.4038 (2P→4P), respectively.(1)
Whenf=0.39 and the initial value [x0,x˙0]=[-0.6,0.5], the time history, phase portrait, and Poincare maps are shown in Figures 3(b1), 3(b2), and 3(b3), respectively. It could be found that the dynamical motion is period-1 in this case.(2)
Whenf=0.402, the time history, phase portrait, and Poincare maps are shown in Figures 3(c1), 3(c2), and 3(c3), respectively. It could be found that it is period-2 motion in this case.(3)
Whenf=0.404, the time history, phase portrait, and Poincare maps are shown in Figures 3(d1), 3(d2), and 3(d3), respectively. It could be found that the dynamical motion is period-4 in this case.At last,f is selected as 0.41 and the initial values are selected randomly. It could be found that it is a chaotic motion in this case. The time history, phase portrait, and Poincare maps are shown in Figures 4(a), 4(b), and 4(c), respectively. Through the analysis of these dynamical characteristics at different typical points, the conclusion that there exist two paths leading to chaos in (1) can be proved.Figure 4
The system response atf=0.41. (a) Time history; (b) phase portrait; (c) Poincare maps.
(a)
(b)
(c)If the same system parameters are substituted into (15), the relation curve of displacement feedback coefficient u and the critical excitation amplitude fmin could be got and shown in Figure 5 with solid line. According to the above numerical analysis process, the critical excitation amplitude fmin with different displacement feedback coefficient u could be found, and the relationship of u and fmin is also shown in Figure 5 with dots. From Figure 5, it could be found that the tendency of analytical results is similar to the numerical simulation. fmin will get smaller with the decreasing of u. That means the increase of u would make the threshold value for generating chaos of (1) larger. The larger u is able to prevent the generation of chaos.Figure 5
The relation curves ofu and fmin.


#### 3.1.2. The Influence of the Time Delayτ1

From the analysis of (15), it could be found that it is an increasing function when τ1∈[0~1.5] and -2ku/2αcsch2⁡k/2τ12kτ1-sinh(2kτ1)>0. System typical parameters are chosen as k=1, c=0.2, α=0.5, ω=1, and u=0.01, so that the relation curve of kτ1 and fmin based on (15) is given in Figure 6 with solid line. The numerical results based on the above-mentioned numerical process are shown in Figure 6 with dots. From the observation of Figure 6, it could be obtained that increase of kτ1 would make the threshold value for generating chaos larger, and it is against the generation of chaos.Figure 6
The relation curves ofkτ1 and fmin.


### 3.2. The Influences of Delayed Velocity Feedback Parameters

If there is only the delayed velocity feedback in (1), (14) turns into(18)2αfπωsech⁡πω2k>4ckk3α-2kkvαcsch3⁡kτ22sinh⁡2kτ2-kτ23+cosh⁡2kτ2.


#### 3.2.1. The Influence of Feedback Coefficientv

Here system typical parameters are chosen ask=1, c=0.4, α=1, ω=0.8, τ2=0.5 and v is chosen from 0 to 0.16. According to (18), the relation curve of v and fmin is obtained and shown in Figure 7 with solid line. The numerical simulation results based on the above-mentioned numerical method are shown in Figure 7 with dots. From the observation of Figure 7, it could be found that the increase of v would make the threshold value for generating chaos smaller, and it is beneficial to generate chaos.Figure 7
The relation curves ofv and fmin.


#### 3.2.2. The Influence of the Time Delayτ2

The system parameters are chosen ask=1, c=0.4, α=1, ω=0.8, v=0.1 and τ2 is chosen from 0 to 0.16. According to (18), the relation curve of τ2 and fmin is drawn in Figure 8 with solid line. The numerical simulation results are shown in Figure 8 with dots. From the observation of Figure 8, it could be found that the increase of τ2 would make the threshold for generating chaos larger, and it is against generating chaos.Figure 8
The relation curves ofτ2 and fmin.Through all the above analysis, it could be found that tendencies of the analytical solutions for the influences of all the delayed feedback parameters are consistent with the numerical iterative simulations. That is to say, a qualitative agreement between the numerical and analytical solutions is obtained. It is generally well known that the analytical necessary condition for the chaos in sense of Smale horseshoes by Melnikov method is the first-order approximate result, so that the quantitative differences between the numerical results and analytical solutions are acceptable. Although there are some differences of analytical solutions with numerical results, the conclusions of the analytical necessary condition for generating chaos are helpful to design this kind of delayed system or control the chaos by choosing appropriate system parameters.


## 3.1. The Influences of Delayed Displacement Feedback Parameters

If there is only the delayed displacement feedback in (1), (14) becomes(15)2fπω2αcsch⁡πω2k>22ckk3α-2ku2αcsch2⁡k2τ12kτ1-sinh⁡2kτ1.


### 3.1.1. The Influence of Feedback Coefficientu

By the analysis of (15), we can find(16)csch2⁡k2τ1>0,limτ1→0⁡2kτ1-sinh⁡2kτ1=0,2kτ1-sinh⁡2kτ1<0,so that(17)-2ku2αcsch2⁡k2τ12kτ1-sinh⁡2kτ1>0.It could be found that the right-hand parts of (15) would become larger with the increasing of the delayed displacement feedback coefficient u. Accordingly, the critical excitation amplitude fmin for generating chaos will become larger with the increasing of u. In other words, the threshold value for chaotic motion in (1) will become larger, so that it is harder to generate chaos with larger u in (1).In order to verify the validity of the analytical necessary condition for generating chaos, a set of illustrative system parameters is chosen ask=1, c=0.4, α=1, ω=0.8, and τ1=0.5. The numerical solutions for fmin of (1) are also investigated by numerical iterative method. When u=0.1, the bifurcation diagram and the corresponding largest Lyapunov exponents are shown in Figures 1(a) and 1(b), respectively. The numerical simulation method for the largest Lyapunov exponents about delayed differential equation can be found in [38]. From the analysis of Figures 1(a) and 1(b), it could be found that the fmin values of generating chaos in two figures are consistent, which implies that the results by numerical iterative method are correct. Through the numerical simulation, the authors also observe some special dynamical phenomenon in this complex system. We find that there exist two paths leading to chaos via different period-doubling bifurcation in (1) which are shown in Figures 2(a) and 3(a), respectively. If the two paths are drawn on the same diagram which is shown in Figure 1, it could be found that the intersection point of the two paths is f=0.3774. When f<0.3774, there is a single periodic solution. If f>0.3774, there will be two periodic-1 solutions, which depend on the initial values. The path shown in Figure 2(a) is simply named as “Type 1” and the other path shown in Figure 3(a) is “Type 2.” In order to verify the conclusion, the time history, phase portrait, and Poincare maps at different typical points are presented and analyzed in the following part.Figure 1
The system responses for different excitation amplitudef.
(a)
Bifurcation diagram
(b)
The largest Lyapunov exponentsFigure 2
(a)
Bifurcation diagram for Type 1
(b)
The system response atf=0.39 for Type 1. (b1) Time history; (b2) phase portrait; (b3) Poincare maps
(c)
The system response atf=0.402 for Type 1. (c1) Time history; (c2) phase portrait; (c3) Poincare maps
(d)
The system response atf=0.4044 for Type 1. (d1) Time history; (d2) phase portrait; (d3) Poincare mapsFigure 3
(a)
Bifurcation diagram for Type 2
(b)
The system response atf=0.39 for Type 2. (b1) Time history; (b2) phase portrait; (b3) Poincare maps
(c)
The system response atf=0.402 for Type 2. (b1) Time history; (b2) phase portrait; (b3) Poincare maps
(d)
The system response atf=0.404 for Type 2. (b1) Time history; (b2) phase portrait; (b3) Poincare mapsThe transitions observed in Type 1 occur atf≈0.3998 (1P→2P) and 0.4040 (2P→4P), respectively. Some typical samples in different ranges are taken to analyze the dynamical phenomenon.(1)
Whenf=0.39 and the initial value [x0,x˙0]=[-0.3956,0.5858], the time history, phase portrait, and Poincare maps are shown in Figures 2(b1), 2(b2), and 2(b3), respectively. It could be found that it is period-1 motion in this case.(2)
Whenf=0.402, the time history, phase portrait, and Poincare maps are shown in Figures 2(c1), 2(c2), and 2(c3), respectively. It is shown the dynamical motion is period-2 in this case.(3)
Whenf=0.4044, the time history, phase portrait, and Poincare maps are shown in Figures 2(d1), 2(d2), and 2(d3), respectively. It could be found that it is period-4 motion in this case.The transitions observed in Type 2 occur atf≈0.4 (1P→2P) and 0.4038 (2P→4P), respectively.(1)
Whenf=0.39 and the initial value [x0,x˙0]=[-0.6,0.5], the time history, phase portrait, and Poincare maps are shown in Figures 3(b1), 3(b2), and 3(b3), respectively. It could be found that the dynamical motion is period-1 in this case.(2)
Whenf=0.402, the time history, phase portrait, and Poincare maps are shown in Figures 3(c1), 3(c2), and 3(c3), respectively. It could be found that it is period-2 motion in this case.(3)
Whenf=0.404, the time history, phase portrait, and Poincare maps are shown in Figures 3(d1), 3(d2), and 3(d3), respectively. It could be found that the dynamical motion is period-4 in this case.At last,f is selected as 0.41 and the initial values are selected randomly. It could be found that it is a chaotic motion in this case. The time history, phase portrait, and Poincare maps are shown in Figures 4(a), 4(b), and 4(c), respectively. Through the analysis of these dynamical characteristics at different typical points, the conclusion that there exist two paths leading to chaos in (1) can be proved.Figure 4
The system response atf=0.41. (a) Time history; (b) phase portrait; (c) Poincare maps.
(a)
(b)
(c)If the same system parameters are substituted into (15), the relation curve of displacement feedback coefficient u and the critical excitation amplitude fmin could be got and shown in Figure 5 with solid line. According to the above numerical analysis process, the critical excitation amplitude fmin with different displacement feedback coefficient u could be found, and the relationship of u and fmin is also shown in Figure 5 with dots. From Figure 5, it could be found that the tendency of analytical results is similar to the numerical simulation. fmin will get smaller with the decreasing of u. That means the increase of u would make the threshold value for generating chaos of (1) larger. The larger u is able to prevent the generation of chaos.Figure 5
The relation curves ofu and fmin.


### 3.1.2. The Influence of the Time Delayτ1

From the analysis of (15), it could be found that it is an increasing function when τ1∈[0~1.5] and -2ku/2αcsch2⁡k/2τ12kτ1-sinh(2kτ1)>0. System typical parameters are chosen as k=1, c=0.2, α=0.5, ω=1, and u=0.01, so that the relation curve of kτ1 and fmin based on (15) is given in Figure 6 with solid line. The numerical results based on the above-mentioned numerical process are shown in Figure 6 with dots. From the observation of Figure 6, it could be obtained that increase of kτ1 would make the threshold value for generating chaos larger, and it is against the generation of chaos.Figure 6
The relation curves ofkτ1 and fmin.


## 3.1.1. The Influence of Feedback Coefficientu

By the analysis of (15), we can find(16)csch2⁡k2τ1>0,limτ1→0⁡2kτ1-sinh⁡2kτ1=0,2kτ1-sinh⁡2kτ1<0,so that(17)-2ku2αcsch2⁡k2τ12kτ1-sinh⁡2kτ1>0.It could be found that the right-hand parts of (15) would become larger with the increasing of the delayed displacement feedback coefficient u. Accordingly, the critical excitation amplitude fmin for generating chaos will become larger with the increasing of u. In other words, the threshold value for chaotic motion in (1) will become larger, so that it is harder to generate chaos with larger u in (1).In order to verify the validity of the analytical necessary condition for generating chaos, a set of illustrative system parameters is chosen ask=1, c=0.4, α=1, ω=0.8, and τ1=0.5. The numerical solutions for fmin of (1) are also investigated by numerical iterative method. When u=0.1, the bifurcation diagram and the corresponding largest Lyapunov exponents are shown in Figures 1(a) and 1(b), respectively. The numerical simulation method for the largest Lyapunov exponents about delayed differential equation can be found in [38]. From the analysis of Figures 1(a) and 1(b), it could be found that the fmin values of generating chaos in two figures are consistent, which implies that the results by numerical iterative method are correct. Through the numerical simulation, the authors also observe some special dynamical phenomenon in this complex system. We find that there exist two paths leading to chaos via different period-doubling bifurcation in (1) which are shown in Figures 2(a) and 3(a), respectively. If the two paths are drawn on the same diagram which is shown in Figure 1, it could be found that the intersection point of the two paths is f=0.3774. When f<0.3774, there is a single periodic solution. If f>0.3774, there will be two periodic-1 solutions, which depend on the initial values. The path shown in Figure 2(a) is simply named as “Type 1” and the other path shown in Figure 3(a) is “Type 2.” In order to verify the conclusion, the time history, phase portrait, and Poincare maps at different typical points are presented and analyzed in the following part.Figure 1
The system responses for different excitation amplitudef.
(a)
Bifurcation diagram
(b)
The largest Lyapunov exponentsFigure 2
(a)
Bifurcation diagram for Type 1
(b)
The system response atf=0.39 for Type 1. (b1) Time history; (b2) phase portrait; (b3) Poincare maps
(c)
The system response atf=0.402 for Type 1. (c1) Time history; (c2) phase portrait; (c3) Poincare maps
(d)
The system response atf=0.4044 for Type 1. (d1) Time history; (d2) phase portrait; (d3) Poincare mapsFigure 3
(a)
Bifurcation diagram for Type 2
(b)
The system response atf=0.39 for Type 2. (b1) Time history; (b2) phase portrait; (b3) Poincare maps
(c)
The system response atf=0.402 for Type 2. (b1) Time history; (b2) phase portrait; (b3) Poincare maps
(d)
The system response atf=0.404 for Type 2. (b1) Time history; (b2) phase portrait; (b3) Poincare mapsThe transitions observed in Type 1 occur atf≈0.3998 (1P→2P) and 0.4040 (2P→4P), respectively. Some typical samples in different ranges are taken to analyze the dynamical phenomenon.(1)
Whenf=0.39 and the initial value [x0,x˙0]=[-0.3956,0.5858], the time history, phase portrait, and Poincare maps are shown in Figures 2(b1), 2(b2), and 2(b3), respectively. It could be found that it is period-1 motion in this case.(2)
Whenf=0.402, the time history, phase portrait, and Poincare maps are shown in Figures 2(c1), 2(c2), and 2(c3), respectively. It is shown the dynamical motion is period-2 in this case.(3)
Whenf=0.4044, the time history, phase portrait, and Poincare maps are shown in Figures 2(d1), 2(d2), and 2(d3), respectively. It could be found that it is period-4 motion in this case.The transitions observed in Type 2 occur atf≈0.4 (1P→2P) and 0.4038 (2P→4P), respectively.(1)
Whenf=0.39 and the initial value [x0,x˙0]=[-0.6,0.5], the time history, phase portrait, and Poincare maps are shown in Figures 3(b1), 3(b2), and 3(b3), respectively. It could be found that the dynamical motion is period-1 in this case.(2)
Whenf=0.402, the time history, phase portrait, and Poincare maps are shown in Figures 3(c1), 3(c2), and 3(c3), respectively. It could be found that it is period-2 motion in this case.(3)
Whenf=0.404, the time history, phase portrait, and Poincare maps are shown in Figures 3(d1), 3(d2), and 3(d3), respectively. It could be found that the dynamical motion is period-4 in this case.At last,f is selected as 0.41 and the initial values are selected randomly. It could be found that it is a chaotic motion in this case. The time history, phase portrait, and Poincare maps are shown in Figures 4(a), 4(b), and 4(c), respectively. Through the analysis of these dynamical characteristics at different typical points, the conclusion that there exist two paths leading to chaos in (1) can be proved.Figure 4
The system response atf=0.41. (a) Time history; (b) phase portrait; (c) Poincare maps.
(a)
(b)
(c)If the same system parameters are substituted into (15), the relation curve of displacement feedback coefficient u and the critical excitation amplitude fmin could be got and shown in Figure 5 with solid line. According to the above numerical analysis process, the critical excitation amplitude fmin with different displacement feedback coefficient u could be found, and the relationship of u and fmin is also shown in Figure 5 with dots. From Figure 5, it could be found that the tendency of analytical results is similar to the numerical simulation. fmin will get smaller with the decreasing of u. That means the increase of u would make the threshold value for generating chaos of (1) larger. The larger u is able to prevent the generation of chaos.Figure 5
The relation curves ofu and fmin.


## 3.1.2. The Influence of the Time Delayτ1

From the analysis of (15), it could be found that it is an increasing function when τ1∈[0~1.5] and -2ku/2αcsch2⁡k/2τ12kτ1-sinh(2kτ1)>0. System typical parameters are chosen as k=1, c=0.2, α=0.5, ω=1, and u=0.01, so that the relation curve of kτ1 and fmin based on (15) is given in Figure 6 with solid line. The numerical results based on the above-mentioned numerical process are shown in Figure 6 with dots. From the observation of Figure 6, it could be obtained that increase of kτ1 would make the threshold value for generating chaos larger, and it is against the generation of chaos.Figure 6
The relation curves ofkτ1 and fmin.


## 3.2. The Influences of Delayed Velocity Feedback Parameters

If there is only the delayed velocity feedback in (1), (14) turns into(18)2αfπωsech⁡πω2k>4ckk3α-2kkvαcsch3⁡kτ22sinh⁡2kτ2-kτ23+cosh⁡2kτ2.


### 3.2.1. The Influence of Feedback Coefficientv

Here system typical parameters are chosen ask=1, c=0.4, α=1, ω=0.8, τ2=0.5 and v is chosen from 0 to 0.16. According to (18), the relation curve of v and fmin is obtained and shown in Figure 7 with solid line. The numerical simulation results based on the above-mentioned numerical method are shown in Figure 7 with dots. From the observation of Figure 7, it could be found that the increase of v would make the threshold value for generating chaos smaller, and it is beneficial to generate chaos.Figure 7
The relation curves ofv and fmin.


### 3.2.2. The Influence of the Time Delayτ2

The system parameters are chosen ask=1, c=0.4, α=1, ω=0.8, v=0.1 and τ2 is chosen from 0 to 0.16. According to (18), the relation curve of τ2 and fmin is drawn in Figure 8 with solid line. The numerical simulation results are shown in Figure 8 with dots. From the observation of Figure 8, it could be found that the increase of τ2 would make the threshold for generating chaos larger, and it is against generating chaos.Figure 8
The relation curves ofτ2 and fmin.Through all the above analysis, it could be found that tendencies of the analytical solutions for the influences of all the delayed feedback parameters are consistent with the numerical iterative simulations. That is to say, a qualitative agreement between the numerical and analytical solutions is obtained. It is generally well known that the analytical necessary condition for the chaos in sense of Smale horseshoes by Melnikov method is the first-order approximate result, so that the quantitative differences between the numerical results and analytical solutions are acceptable. Although there are some differences of analytical solutions with numerical results, the conclusions of the analytical necessary condition for generating chaos are helpful to design this kind of delayed system or control the chaos by choosing appropriate system parameters.


## 3.2.1. The Influence of Feedback Coefficientv

Here system typical parameters are chosen ask=1, c=0.4, α=1, ω=0.8, τ2=0.5 and v is chosen from 0 to 0.16. According to (18), the relation curve of v and fmin is obtained and shown in Figure 7 with solid line. The numerical simulation results based on the above-mentioned numerical method are shown in Figure 7 with dots. From the observation of Figure 7, it could be found that the increase of v would make the threshold value for generating chaos smaller, and it is beneficial to generate chaos.Figure 7
The relation curves ofv and fmin.


## 3.2.2. The Influence of the Time Delayτ2

The system parameters are chosen ask=1, c=0.4, α=1, ω=0.8, v=0.1 and τ2 is chosen from 0 to 0.16. According to (18), the relation curve of τ2 and fmin is drawn in Figure 8 with solid line. The numerical simulation results are shown in Figure 8 with dots. From the observation of Figure 8, it could be found that the increase of τ2 would make the threshold for generating chaos larger, and it is against generating chaos.Figure 8
The relation curves ofτ2 and fmin.Through all the above analysis, it could be found that tendencies of the analytical solutions for the influences of all the delayed feedback parameters are consistent with the numerical iterative simulations. That is to say, a qualitative agreement between the numerical and analytical solutions is obtained. It is generally well known that the analytical necessary condition for the chaos in sense of Smale horseshoes by Melnikov method is the first-order approximate result, so that the quantitative differences between the numerical results and analytical solutions are acceptable. Although there are some differences of analytical solutions with numerical results, the conclusions of the analytical necessary condition for generating chaos are helpful to design this kind of delayed system or control the chaos by choosing appropriate system parameters.


## 4. Conclusion

In this paper, the heteroclinic bifurcation behavior of a Duffing oscillator under forcing excitation with both delayed displacement feedback and delayed velocity feedback is investigated. The Melnikov function is established and the analytical necessary condition for generating chaos in (1) is obtained. The numerical simulation based on numerical iterative method is used to verify the correctness of the analytical necessary condition for generating chaos. The bifurcation and the large Lyapunov exponent figures are presented. From the analysis of numerical simulations, it could be found that two paths leading to chaos via period-doubling bifurcation exist in (1). The influences of the delayed feedback parameters are studied, respectively. From the analysis, we find that the tendencies of the analytical solutions for the influences of all the delayed feedback parameters are consistent with the numerical iterative simulations, which verifies the correctness of the analytical necessary condition. The systematic and comprehensive results are obtained. Those results will be helpful to design or control this kind of delayed feedback system.



---

*Source: 7213606-2018-01-21.xml*