BIFURCATIONS OF PERIODIC SOLUTIONS AND CHAOS IN DUFFING-VAN DER POL EQUATION WITH ONE EXTERNAL FORCING∗

The Duffing-Van der Pol equation with fifth nonlinear-restoring force and one external forcing term is investigated in detail: the existence and bifurcations of harmonic and second-order subharmonic, and third-order subharmonic, third-order superharmonic and m-order subharmonic under small perturbations are obtained by using second-order averaging method and subharmonic Melnikov function; the threshold values of existence of chaotic motion are obtained by using Melnikov method. The numerical simulation results including the influences of periodic and quasi-periodic and all parameters exhibit more new complex dynamical behaviors. We show that the reverse period-doubling bifurcation to chaos, period-doubling bifurcation to chaos, quasi-periodic orbits route to chaos, onset of chaos, and chaos suddenly disappearing, and chaos suddenly converting to period orbits, different chaotic regions with a great abundance of periodic windows (periods:1,2,3,4,5,7,9,10,13,15,17,19,21,25,29,31,37,41, and so on), and more wide period-one window, and varied chaotic attractors including small size and maximum Lyapunov exponent approximate to zero but positive, and the symmetry-breaking of periodic orbits. In particular, the system can leave chaotic region to periodic motion by adjusting the parameters p, β, γ, f and ω, which can be considered as a control strategy.


Introduction
In this paper, we consider the following Duffing-Van der Pol (DVP) equation with fifth nonlinear force and one external forcing x + p(x 2 − 1)ẋ + ω 2 0 x + βx 3 + γx 5 = f cos ωt, (1.1) where p, ω 0 , β, γ, f, ω are real parameters. Physically, p can be regarded as dissipation or damping, β and γ are the strength of nonlinearity, f and ω as the amplitude and frequency of the external force.
The Duffing-Van der Pol equation is a combination of Duffing and Van der Pol equations and has application in the simulation of nonlinear oscillation systems [9,20]. The dynamical behaviors of the Duffing-Van der Pol equation as the parameters varying are considered in [8,14,15,19]. In the previously work [7], we investigate the Duffing-Van der Pol equation with two external forcings, and give the threshold values of existence of chaos under the periodic perturbation, and the criterion of existence of chaos in averaged system under quasi-periodic perturbation for ω 2 = nω 1 + ϵσ, n = 1, 3, 5, and can't prove the criterion of existence of chaos in second-order averaged system under quasi-periodic perturbation for ω 2 = nω 1 + ϵσ, n = 2, 4, 6,7,8,9,10, where σ is not rational to ω 1 , and obtain rich dynamical behaviors by numerical simulations. In this paper, the bifurcations of periodic orbits, and chaos for (1.1) are investigated by using second-order averaging method and Melnikov methods in [1, 3-5, 11-13, 16-18, 21, 23, 24]. We give the existence and bifurcations of harmonics and second-order subharmonic, and third-order subharmonic, third-order superharmonic and m-order subharmonic under small perturbations, and the criterion of existence of chaos. We also give numerical simulations including bifurcation diagrams of fixed points and system, computation of maximum Lyapunov exponents, phase portraits, and Poincaré map, and consider the influences of periodic and quasi-periodic and all parameters. We show that the reverse period-doubling bifurcation to chaos, period-doubling bifurcation to chaos, quasi-period route to chaos, onset of chaos, and chaos suddenly disappearing and chaos suddenly converting to period orbits, different chaotic regions with a great abundance of periodic windows (periods: 1,2,3,4,5,7,9,10,13,15,17,19,21,25,29,31,37,41, and so on), and more wide period-one window, and varied chaotic attractors including small size and maximum Lyapunov exponent approximate to zero but positive, and the symmetry-breaking of periodic orbits. In particular, the system can leave chaotic region to periodic motion by adjusting the parameters p, β, γ, f and ω, which can be considered as a control strategy.
The paper is organized as follows. Analytical results for the conditions of existence and bifurcations of harmonic and second-order subharmonics are given in section 2 and 3, respectively. In section 4, we provide the conditions of existence and bifurcation for the third-order subharmonic resonance. In section 5, the criterion for the existence of m-order subharmonics is proved by using subharmonic Melnikov function. In section 6, we present the conditions of existence and bifurcation for the third-order superharmonic resonance. In section 7, Melnikov's method is used to prove the existences of homoclinic bifurcation and heteroclinic bifurcation. The numerical simulations including bifurcation diagrams in (x − f ), (x − ω 0 ), (x − γ), (x − β), (x − p), and (x − ω) planes, the computation of maximum Lyapunov exponent corresponding to bifurcation diagram, the phase portrait and Poincaré map at neighborhood of critical values, are given in section 8. Finally, we give remark in section 9.

Primary Resonance and Bifurcation
In this section, we consider primary resonance using the second-order averaging method. Introduce a small parameter ϵ, such that 0 < ϵ ≪ 1 and replace p and f by ϵp and ϵ 3 2 f respectively, then (1.1) can be rewritten as Let (x 0 , 0) be a center of (2.1) for ϵ = 0 and the frequency of periodic orbit near the center is approximately given by If the ratio of ω and a 1 is a rational number, then the resonance behavior may occur in (2.1). We now consider the case of primary resonance ω ≈ a 1 for (2.1). Assume that Then (2.1) can be rewritten as We use the Van der Pol transformation and carry out averaging up to second-order for (2.5). We obtain the following averaged equationu In polar coordinates r = √ u 2 + v 2 and θ = arctan (v/u), (2.7) becomeṡ The fixed points of (2.8) satisfy the following equation where y = r 2 . Let and y = x − 2Ω 3b0 , then (2.9) becomes The discriminant of (2.11) is given by 13) and ∆ ′ denotes the discriminant of (2.13), then ∆ ′ = 16 We have the following conclusion as b 0 Ω < 0.
(i) If ∆ ′ > 0, then there exist two real roots of (2.13) as following and (ii) If ∆ ′ < 0, then there exists not any real root of (2.13).
(ii) If (Ω, f 2 ) lies in the curve f 2 1 or f 2 2 , there are three real roots of (2.9), but two of them coincide; on f 2 1 the coincidence occurs for the root y =
|), there are three coincident real roots, which correspond to one nontrivial fixed point of (2.8).
1 , there is one real root of (2.9), which corresponds to one nontrivial fixed point of (2.8). The stability of fixed point (r s , θ s ) of (2.8) is determined by the characteristic values where y s = r 2 s and F ′ 0 (y s ) = 3b 2 0 y 2 s + 4b 0 Ωy s + Ω 2 + p 2 α 2 a 2 1 . Furthermore, we know that the averaged system (2.8) (or (2.7)) has no closed orbit by the Dulac's criterion.
Thus, by considering the above stability conditions and Lemma 2.1, we can get the following conclusion for α > 0: (ii) For Ω 2 − 3p 2 α 2 a 2 1 > 0, b 0 Ω < 0 and 0 < ϵ ≪ 1, there exist two stable foci-node and one saddle in region (II); there exists a stable foci-node in region (I) and (III). On the curves f 2 1 and f 2 2 there is one stable foci-node and a fixed point which has one zero eigenvalue.
|) at which there is a fixed point with one zero eigenvalue.
(iv) When f 2 increases, the fixed point changes from one to three, passing through f 2 2 for b 0 < 0 or f 2 1 for b 0 > 0; on f 2 2 or f 2 1 , there are two fixed points and one of them has a zero eigenvalue, so f 2 2 or f 2 1 is supercritical saddle-node bifurcation. When f 2 decreases, the fixed point changes from one to three, passing through f 2 1 for b 0 < 0 or f 2 2 for b 0 > 0; on f 2 1 or f 2 2 , there are two fixed points and one of them has a zero eigenvalue, so f 2 1 or f 2 2 is subcritical saddle-node bifurcation. For α < 0, the stable foci-node becomes unstable and the stability of other fixed points does not change in Lemma 2.2.
(iv) The harmonic solution of (2.1) is approximately given by where (r s , θ s ) is given by the equilibrium solution of averaged (2.8). The other solutions in (2.8) correspond to the almost periodic solutions or chaotic motions in (2.1).

Second-order Subharmonic Resonance and Bifurcation
In this section we consider the second-order subharmonic resonance ω ≈ 2a 1 and set ϵΩ = (ω 2 − 4a 2 1 )/4. Replace p and f by ϵp and ϵf (0 < ϵ ≪ 1), respectively, and then (1.1) can be written as (3.1) Using regular perturbation methods, one obtains harmonic of (1.1) as To investigate stability of the harmonicx(t), one can set in (3.1) and carrying out averaging up to second order, one haṡ In polar coordinates, (3.6) becomeṡ By analyzing the fixed points of (3.7), we have the following conclusion.
(a) (b) Figure 3. Super-critical and sub-critical saddle-node bifurcation curves: By the averaging theorem we get the following conclusion. . Each pair of nontrivial fixed point (r ± , θ ± ) and (r ± , θ ± + π) corresponds to a single second-order subharmonic of (3.1), which is approximately given by (ii) There are two bifurcation curves: saddle-node bifurcations of subharmonics occur near the curve and period doubling bifurcations of harmonics occur near the curve (iii) In region (I), there is a stable fixed points (0, 0) of (3.6), which corresponds to a stable (resp. unstable) non-resonant harmonic of (3.1) for α > 0 and an unstable fixed points (0, 0) of (3.6), which corresponds to an unstable nonresonant harmonic of (3.1) for α < 0; in region (II), there are three fixed points and one of them corresponds to a unstable non-resonant harmonic and a single stable resonant subharmonic of period-two of (3.1) for α > 0 and a single unstable resonant subharmonic of period-two of (3.1) for α < 0; in region (III), there are five fixed points which correspond to a stable nonresonant harmonic and a stable resonant subharmonic for α > 0 or an unstable non-resonant harmonic and an unstable resonant subharmonic for α < 0 and an unstable resonant subharmonic of (3.1), where the regions(I)-(III) are in Fig.3.

The m-Order Subharmonics and Bifurcation
In this section we investigate the existence of m-order subharmonics of (1.1) by using Melnikov's method for subharmonic which is defined in [23]. Consider the perturbed system {ẋ = y, Let q µ (t) = (x µ (t), y µ (t)) (µ ∈ (µ 1 , µ 2 )) denote a one-parameter family of periodic orbits with period 2πm nω of (5.1) for ϵ = 0, where µ 1 and µ 2 are constants, and m and n are relatively prime. In [24], it has been proved that M m (t 0 ) can have simple zero only if n = 1, so the Melnikov function for q µ (t) of (5.1) is given by then M m (t 0 ) has a simple zero and a necessary condition for the occurrence of subharmonics of period 2πm ω of (5.1) is given by (5.4). The bifurcation curve of subharmonic is created and occurs at
, there exists one resonant third-order superharmonic.

Chaos in Equation (1.1)
In this section, we discuss the chaotic behaviors of equation (1.1) in which f and p are assumed to be small parameters with order ε. Rewriting (1.1) as an autonomous system gives   ẋ = y, where εf = f, εp 1 = p.
The unperturbed system of system (7.1) has two homoclinic orbits Γ ± hom and two heteroclinic orbits Γ ± het for γ > 0, β < 0 and β 2 − 4γω 2 0 > 0. The phase portrait of the unperturbed system is shown in Fig.5(a) for ω 0 = 1, β = −3, and γ = 2. When the perturbation is added, the homoclinic or heteroclinic orbits break, and may have transverse homoclinic or heteroclinic orbits. By the Smale-Birkhoff homoclinic theorem [3,23], the existence of such orbits results in chaotic dynamics. We therefore apply the Melnikov method to equation (7.1) for finding the criteria of the existence of homoclinic or heteroclinic bifurcation and chaos.
Suppose that the homoclinic or heteroclinic orbits of the unperturbed systems are written as (x 0 (t), y 0 (t)), then the Melnikov function for system (7.1) can be given where t 0 is the cross-section time of the Poincaré map and t 0 can be interpreted as the initial time of the forcing term. Because it is difficult to give analytical expression of (x 0 (t), y 0 (t)), we will compute x 0 (t) and y 0 (t) numerically. We note that y 0 (t) is a function of time from −∞ to +∞. We therefore choose a starting point P 1 which is an intersecting point of homoclinic orbit Γ ± hom1 with x-axis, and a starting point P 2 which is an intersecting point of heteroclinic orbit Γ ± het with y-axis, and y 0 (t) would be an odd function of time for the homoclinic orbit and an even function of time for the heteroclinic orbit (see Fig.5(a)). For the homoclinic orbits Γ ± hom1 , the Melnikov function can be simplified as where A = ∫ +∞ 0 (x 2 0 (t) − 1)y 2 0 (t)dt is a constant once x 0 (t) and y 0 (t) are given, I hom1 (ω) = ∫ +∞ 0 y 0 (t) sin(ωt)dt is a function of the frequency ω when the homoclinic orbits (x 0 (t), y 0 (t)) are given. Thus, if (7.4) then there is at 0 such that M 1 (t 0 ) = 0 and ∂M1 ∂t0 | t=t0 ̸ = 0,t 0 ∈ [0, 2π ω ] and the following theorem can be obtained.
Theorem 7.1. The homoclinic bifurcation will occur at f = R 1 (ω), (7.5) this implies that if ε > 0 is sufficiently small, the transverse homoclinic orbits exist and system (7.1) may be chaotic. For the heteroclinic orbits, y 0 (t) is even, then the Melnikov function on Γ ± het can be written as y 0 (t) cos(ωt)dt is a function of the frequency ω when the heteroclinic orbits (x 0 (t), y 0 (t)) are given. Thus, if then there is at ′ 0 such that M 2 (t ′ 0 ) = 0 and ∂M1 2π ω ] and the following theorem can be obtained.
Theorem 7.2. The heteroclinic bifurcation will occur at f = R 2 (ω), (7.8) this implies that if ε > 0 is sufficiently small, the transverse heteroclinic orbits exist and system (1.1) may be chaotic.
Remark. In order to compare the dynamical behavior of system (1.1) with that in [7], we take up with the same values for each case. And we show the differences of bifurcation diagrams and dynamical behaviors of system (1.1) and the system in [7] from the following results.
For case (5). The bifurcation diagram of system (1.1) in (p, x) plane and the corresponding maximum Lyapunov exponents for ω 0 = 1 are shown in Fig.15(a) and (b), the local amplified bifurcation diagrams of Fig.15(a) for p ∈ (2.95, 3.35) and p ∈ (3.95, 4.2) are given in Fig.15(c)and (d), respectively. We show that the period-doubling bifurcation leading to a small chaotic region, and the inverse perioddoubling bifurcation leading to another small chaotic region, and two more wide period-one regions.  The bifurcation diagram of system (1.1) in (p, x) plane and corresponding maximum Lyapunov exponents for ω 0 = √ 2 are shown in Fig.16(a) and (b), the local amplifications of Fig.16(a) for p ∈ (0, 0.4) is given in Fig.16(c). There are the periodone orbit suddenly converting to quasi-periodic orbits and the chaotic region with small positive Lyapunov exponents, and complex periodic windows (periods 3,11,23,11,23,21,11,19,29,5,21,39,21 and so on) as p increasing, and they appear alternatively, and the final chaotic region, ending in a periodic-one orbit. The quasiperiodic orbit at f = 0.5(L = 0) and chaotic attractor at f = 4.75(L = 0.110068) are given in Fig.17(a) and (b), respectively.
For case (6). (i) The bifurcation diagram of system (1.1) in (ω, x) plane and maximum Lyapunov exponents for ω 0 = √ 2 are shown in Fig.18(a) and (b). The local amplified bifurcation diagram of Fig.18(a) for ω ∈ (0, 0.4) and ω ∈ (0.136, 0.15) are given in Fig.18(c) and (d). The system (1.1) exhibits the quasi-periodic route to chaos in Fig.18(c), and another two routes to chaos: period-doubling bifurcation to chaos, and inverse period-doubling bifurcation to chaos in Fig.18(a) and (d). We also show that the interleaving occurrence of chaotic behaviors and quasi-periodic, and the chaotic regions with period-windows, and the jumping behaviors of periodic orbit and the symmetry-breaking of one pair period-orbit, and the chaos ending in a periodic orbit in Fig.18(a)-(d). The phase portrait of quasi-periodic orbits is given in Fig.18(e).
Remark. The bifurcation diagram in (ω, x) plane and maximum Lyapunov exponents for ω = 1 (in case (6)) are similar to the Fig.18  (ii) The bifurcation diagram of system (1.1) in (ω, x) plane and the corresponding maximum Lyapunov exponents for ω 0 = 1, and f = 2 are shown in Fig.19(a) and (b), the local amplified bifurcation diagram of (a) for ω ∈ (1.3, 1.5) is given in (c). There are the chaotic regions with period windows (period 3, 7 and so on), and the chaos suddenly converting to periodic one orbit.
(iii) The bifurcation diagram of system (1.1) in (ω, x) plane for β = 0.2, γ = 1, ω 0 = √ 2, p = 3 and f = 0.5 and maximum Lyapunov exponents are shown in Fig.20(a) and(b), where we show that the maximum Lyapunov exponents are equal or approximate to zero except for a small interval of ω. The system (1.1) therefore presents fundamentally the quasiperiodic behaviors and chaotic behaviors except for a small period-window. The phase portraits of the chaotic behavior and quasiperiodic behavior for ω = 4.693 (L = 0.00039129) and ω = 5.95 (L = −0.0020617) are given in Fig.20(c) and (d), respectively.

Remark
In this paper and previously paper [7], we have investigated the dynamical behaviors of Duffing-Van der Pol equation with an external forcing and two external forcings, respectively. We give the conditions of existences and bifurcations for harmonic, subharmonics, superharmonic and chaos for the system with one forcing in the paper. And we give the threshold values of existence of chaotic motion under the periodic perturbation for original system and under quasiperiodic perturbation for averaged system for the system with two external forcings, but haven't to give the condition of various subharmonics in [7]. Moreover, comparing Fig.10, Fig.15, Fig.18, Fig.21, Fig.24, Fig.25 in [7] with Fig.6, Fig.10, Fig.11, Fig.13, Fig.15, Fig.18 in the paper, respectively, we shown the dynamical behaviors of one external forcing are different from that of two external forcings. In particular, there are the more small chaotic regions (Lyapunov exponent approximating to zero and small size) with a great abundance of periodic windows, and the system can be leave chaotic region to periodic motion by adjusting the parameters p, β, γ, f and ω, which can be considered as a control strategy for the system with one external forcing. But there are the more large chaotic regions (larger Lyapunov exponent and wide region) with quasi-periodic, periodic widows, and the system can leave chaotic region to periodic motion by adjusting the parameter f 2 and there are always the chaotic behavior or quasi-periodic behavior as adjusting other parameters for two external forcings. These results are useful for understanding dynamical behaviors for the Duffing-Van der Pol system.