Canards, heteroclinic and homoclinic orbits for a slow-fast predator-prey model of generalized Holling type III

Available online at www.sciencedirect.com

ScienceDirect J. Differential Equations 267 (2019) 3397–3441 www.elsevier.com/locate/jde

Canards, heteroclinic and homoclinic orbits for a slow-fast predator-prey model of generalized Holling type III Cheng Wang a , Xiang Zhang b,∗ a School of Mathematical Sciences, Shanghai Jiao Tong University, Shanghai 200240, PR China b School of Mathematical Sciences, MOE–LSC, Shanghai Jiao Tong University, Shanghai, 200240, PR China

Received 12 October 2018; revised 26 February 2019 Available online 6 May 2019

Abstract For a classical ratio-dependent predator-prey model with the generalized Holling type III functional response, it was previously investigated in [20] by Hsu and Huang for global stability of an equilibrium, and in [21] by Huang, Ruan and Song for subcritical Hopf and Bogdanov-Takens bifurcations. Here in this model when prey reproduces much faster than predator, by using geometric singular perturbation theory, we achieve much richer new dynamical phenomena than the existing ones, such as the existence of canard cycles, canard explosion and relaxation oscillations, heteroclinic and homoclinic orbits, cyclicity of slowfast cycles, and the coexistence of the Hopf cycle and the relaxation oscillation. On global stability of the equilibrium we also provide less restricted conditions than the existing ones. © 2019 Elsevier Inc. All rights reserved. MSC: 34C37; 34C26; 34C07; 37G15 Keywords: Predator-prey model; Slow-fast system; Geometric singular perturbation theory; Heteroclinic and homoclinic orbits; Canard cycle; Relaxation oscillation

* Corresponding author.

E-mail addresses: [email protected] (C. Wang), [email protected] (X. Zhang). https://doi.org/10.1016/j.jde.2019.04.008 0022-0396/© 2019 Elsevier Inc. All rights reserved.

3398

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Table 1 Biological meaning of the parameters of system (1.1).

1. Introduction In this paper, employing geometric singular perturbation theory (GSPT) laid by Fenichel [14–17], we present some new insights into the rich dynamics that can arise in a classical ratiodependent predator-prey model, suggested by Freedman and Mathsen in [18] and Collings in [4]. The model we are studying takes the following form  dx x = x˙ = rx 1 − − yp(x), dτ¯ K  y  dy = y˙ = sy 1 − , dτ¯ hx

(1.1)

with a generalized Holling type III functional response p(x) =

mx 2 , ax 2 + bx + 1

where x(τ¯ ) and y(τ¯ ) indicate the prey and predator densities respectively, at time τ¯ . The param√ eters r, K, s, h, m and a are positive constants, and b > −2 a (so that ax 2 + bx + 1 > 0 for all x ≥ 0 and hence p(x) > 0 for all x > 0). The predator’s numerical response is Leslie form originated by Leslie in [30]. The parameters a, b and m are fitting parameters of the functional response. The biological meaning of other parameters are shown in Table 1. An important feature of biological systems is that they often depend on multiple scales. Multiple scales problems of biological systems are usually modeled by slow-fast systems [19]. Next we provide necessary background on slow-fast systems and GSPT limited to two dimensional differential systems for our applications. Consider a planar system of differential equations ε

dx =ε x˙ = f (x, y, λ, ε), dτ dy =y˙ = g(x, y, λ, ε), dτ

(1.2)

where (x, y) ∈ R2 are the coordinates of the phase space, λ ∈ Rk are parameters, and 0 < ε  1 is a sufficiently small perturbation parameter. The functions f and g are assumed to be sufficiently smooth. System (1.2) is called a slow-fast system with a fast variable x and a slow variable y. By switching from the slow time scale τ to the fast time scale t = τ/ε, system (1.2) can be rescaled to

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3399

dx =x  = f (x, y, λ, ε), dt dy =y  = εg(x, y, λ, ε). dt

(1.3)

Taking the singular limit ε → 0 in system (1.2) gives 0 = f (x, y, λ, 0), y˙ = g(x, y, λ, 0),

(1.4)

which is a differential-algebraic equation called the slow subsystem. The flow generated by system (1.4) is called the slow flow. Similarly, letting ε → 0 in system (1.3) yields x  = f (x, y, λ, 0), y  = 0,

(1.5)

which can be seen as a one dimensional differential equation in the fast variable x via treating the slow variable y as a parameter. We call system (1.5) the fast subsystem, whose flow is called the fast flow. The slow subsystem is defined on the critical set C0 := {(x, y) ∈ R2 : f (x, y, λ, 0) = 0} If C0 is a submanifold, we refer to C0 as a critical manifold. Note that C0 does not have to be a manifold [26] in general. Here in this paper, the critical set will be a manifold. Observe also that C0 can also be interpreted as a manifold of equilibria for the fast subsystem (1.5). A subset S ⊂ C0 is called normally hyperbolic if ∂x f (p, λ, 0) = 0 for all p ∈ S. A point p on S is normally attracting (resp. normally repelling) in case that ∂x f (p, λ, 0) has a strictly negative (resp. strictly positive) eigenvalue. In the other case we call the points p contact points where normal hyperbolicity of the fast subsystem (1.5) is lost. If the ∂x f (p, λ, 0) = 0 at p ∈ S then the implicit function theorem describes S locally as a graph h0 : R → R,

f (h0 (y), y, λ, 0) = 0

near p. This allows us to write the slow subsystem (1.4) more concisely as y˙ = g(h0 (y), y, λ, 0). One goal of GSPT is to use these lower-dimensional flows of the slow subsystem (1.4) and of the fast subsystem (1.5) to understand the dynamics of the full system (1.2) or (1.3) for sufficiently small ε > 0. GSPT was pioneered by Fenichel [17,14–16] and popularized by Jones [22]. Today, GSPT encompasses not only the results of Fenichel but also a much wider range of geometric techniques, such as blow-up method [11,12,28,25,27,26], exchange lemma [22,32,24,23,35,36], and slow-fast normal form theory [1] and so on. Canards were discovered and first analyzed by French mathematicians Benoît et al. in [2]. By definition a trajectory segment of the slow-fast system (1.2) is a canard if it follows an attracting

3400

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Fig. 1.1. The critical manifold C0 := {(x, y)|y = F (x)} (dashed black) consists of three submanifolds: C0l and C0r attracting, and C0m repelling, together with two normally non-hyperbolic points (xm , F (xm )) and (xM , F (xM )). Double arrows indicate fast flow, and single arrows indicate slow flow. (a) The canard-type slow-fast cycle without head. (b) The transitory slow-fast cycle. (c) The canard-type slow-fast cycle with head. (d) The common slow-fast cycle.

slow manifold, passes close to a contact point of the critical manifold, and then follows a repelling slow manifold for O(1) time. An orbit lying in the intersection of the attracting slow manifold and the repelling slow manifold is called a maximal canard. The curves consisting of equilibria for the fast subsystem (1.5) are also called slow curves. The regular orbits of the fast subsystem (1.5) are called fast orbits. Given an ordered sequence of fast orbits (compact parts including α and ω limits) and compact pieces of slow arcs, oriented by the respective fast and slow subsystems, assuming that their union is a topological circle, we call the sequence a limit periodic set or a slow-fast cycle [9]. A common slow-fast cycle is a slow-fast cycle where the slow parts are either all attracting or all repelling. Fig. 1.1 illustrates partly the geometrical descriptions of the slow-fast cycles. Denote by Xε,λ the vector field associated to system (1.2). For a fixed λ = λ0 and a slow-fast cycle γ0 containing both attracting and repelling parts of the critical curve, if there are δ > 0 and ε0 > 0, such that for each ε ∈ (0, ε0 ) and λ = λ(ε) ∼ λ0 , the vector field Xε,λ has a limit cycle γελ in the δ-neighborhood of the slow-fast cycle γ0 , corresponding to (ε, λ) = (0, λ0 ), and γελ → γ0 (in Hausdorff sense) as ε → 0, then γελ is called a canard cycle, bifurcating from γ0 . The maximal number of such canard cycles, taking into account their multiplicities, is called the cyclicity of γ0 for Xε,λ at (ε, λ) = (0, λ0 ) and is denoted by Cycl(Xε,λ , γ0 , (0, λ0 )) [31]. Relaxation oscillation is a periodic orbit of the slow-fast system (1.2) which converges (with respect to Hausdorff distance) to a slow-fast cycle consisting of alternating fast and slow segments forming a closed loop, where the slow ones have the same stability (e.g. Fig. 1.1 (d)) as ε → 0. The very fast transition upon variation of a control parameter from a small amplitude limit cycle (Hopf cycle) via canard cycles to a large amplitude relaxation oscillation within an exponentially small range O(exp(−1/ε)) of the control parameter is called canard explosion. System (1.1), with the generalized Holling type III functional response, has been partially analyzed in [20], in which the authors studied the global stability of the positive locally asymptotically stable equilibrium for b > 0. In [21], the authors focused on the subcritical Hopf bifurcation and Bogdanov-Takens bifurcation of system (1.1). But except [20,21], there are few other interesting previous works on dynamics of systems (1.1). Relative works are those by Li and Zhu in [31] applying the GSPT to study the canard cycles of classical Gause prey-predator model with generalized Holling type III and IV functional response. In this paper, we assume that the prey reproduces much faster than the predator in system (1.1), i.e., s  r. We note that this is a natural assumption if predators live very long and encounter many different generations of prey. For example, hares and squirrels reproduce much faster than their predators, such as lynx and coyotes. Another example satisfying this assump-

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3401

tion is the spruce-budworm system [33]. Applying the GSPT we can obtain several new and interesting phenomena on dynamics of system (1.1), for instance the existence of canards, of heteroclinic orbits and of homoclinic orbits, which have never been observed previously in the existing references. Since system (1.1) is not defined on the y-axis and the origin represents extinction of both populations, from the biological point of view it is interesting to consider the dynamics near the origin. Using Poincaré compactification and successive blow-up technique, we show that the origin is a non-hyperbolic saddle, and so it implies that when the prey is at a low density, the predator density decreases and approaches extinction. In the interior of the first quadrant, the dynamics of system (1.1) depends on the next three different cases. (I) System (1.1) has only one positive equilibrium. We prove that the unique positive equilibrium is a global attractor when it lies on the normally hyperbolic attracting critical submanifold, see Theorem 5.1. When it lies on the contact points, system (1.1) exhibits canard explosion phenomenon, see Theorem 5.2. When it lies on the normally hyperbolic repelling critical submanifold, there exists a unique relaxation oscillation in a tubular neighborhood of the common slow-fast cycle, which is the only limit cycle in the first quadrant, see Theorem 5.3. (II) System (1.1) has two positive equilibria. We show that there are infinitely many heteroclinic orbits, see Theorem 5.5; there exists a unique homoclinic orbit and a limit cycle, see Theorem 5.6; or there exists a unique relaxation oscillation, see Theorem 5.7. (III) System (1.1) has three positive equilibria E1 (x1 , y1 ), E2 (x2 , y2 ) and E3 (x3 , y3 ) with x1 < x2 < x3 . There exist heteroclinic orbits connecting E1 and E2 , and E2 and E3 , see Theorem 5.8; there exist two heteroclinic orbits connecting equilibria E1 and E2 , and infinitely many heteroclinic orbits connecting E1 and E3 , see Theorem 5.9; there exists a unique relaxation oscillation, see Theorem 5.10; or there exists a unique homoclinic orbit which could bifurcate a unique canard cycle, see Theorem 5.11, 5.12 and 5.13. Finally, by computing slow divergence integrals, we prove with an additional condition that there is at most one limit cycle bifurcating from the slow-fast cycles, see Theorem 6.9. Moreover, we show via canard explosion that two limit cycles, which are the Hopf cycle and the relaxation oscillation or the canard cycle without head and the canard cycle with head, can coexist, see Proposition 6.6. The main tools in this paper are the slow-fast normal form theory [1] and the geometric singular perturbation theory including Fenichel’s theory [14–17] and the results near normally nonhyperbolic points developed by Krupa and Szmolyan [25,27]. Our studies are analytic together with numerical realization. The rest of this paper is organized as follows. Section 2 presents a normal form of the slow-fast predator-prey system, the conditions for existence of the S-shaped critical manifold and equilibria, and the invariant region. Section 3 presents dynamics near the origin and the boundary equilibrium (K, 0). Section 4 develops a slow-fast normal form for system (1.1) near a contact point. Section 5 describes dynamics of system (1.1) in the interior of the first quadrant. Section 6 investigates the cyclicity of slow-fast cycles that can occur in system (1.1). Section 7 describes the biological relevance of our results on relaxation oscillation.

3402

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

2. The slow-fast predator-prey model In this section, a normalized slow-fast predator-prey system is derived from system (1.1), and by which one discusses the shape of its critical manifold. Also, the results on the existence of equilibria and the positive invariant region of the slow-fast system are presented. Note that system (1.1) is defined on the set := {(x, y) ∈ R2 : x > 0, y ≥ 0}, taking into account the biological meaning of the model. In order to simplify calculations in studying dynamics of system (1.1) in the interior of the first quadrant, we first normalize system (1.1). Lemma 2.1. System (1.1) is topologically equivalent to the system    dx¯ x¯ ¯ = x¯ 2 x¯ 2 + b¯ x¯ + 1 1 − − x¯ 3 y, dt K¯   dy¯ = ε y( ¯ x¯ 2 + b¯ x¯ + 1) x¯ − h¯ y¯ , dt

(2.1)

√ ¯ h) ¯ ∈ R3+ and b¯ > −2, where b¯ = √b , K¯ = K a, h¯ = ra and ε = s . with (ε, K, mh r a ¯ × R := {(x, Proof. The invertible map  :

¯ y) ¯ : x¯ > 0, y¯ ≥ 0} × R → × R by ⎛

√

1 r a (x, y, τ¯ ) = (x, ¯ y, ¯ t) := ⎝ √ x, y, ¯ ¯ m a

(x¯ 2 +

√b x¯ a

+ 1)x¯

r

⎞ t⎠ ,

preserves orientation of the time, and sends system (2.1) to     dx¯ b x¯ 2 2 − x¯ 3 y, ¯ = x¯ + √ x¯ + 1 x¯ 1 − √ dt a aK    dy¯ b ra y¯ s = y¯ x¯ 2 + √ x¯ + 1 x¯ − , dt r mh a √ b ra s which by setting b¯ = √ , K¯ = K a, h¯ = and ε = gives system (2.1). mh r a

2

2.1. Slow-fast structure Dropping the bars for notational convenience in system (2.1), we obtain the following standard slow-fast predator-prey system:   dx x = x  = x 2 x 2 + bx + 1 1 − − x 3 y ≡ f (x, y, η), dt K dy = y  = εy(x 2 + bx + 1) (x − hy) ≡ εg(x, y, η), dt

(2.2)

where x > 0, y ≥ 0, η = (b, K, h) and 0 < ε  1. Setting ε = 0 in (2.2) yields the fast subsystem:

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3403

Fig. 2.1. (a) The curve C1 is tangent to the straight line L : {b = −2} at the point A(−2, 1). (b) y = G(x) with (b, K) ∈ U . (c) The S-shaped critical manifold C0 (blue curve) for system (2.2) consists of three submanifolds: C0l and C0r attracting parts, C0m repelling part. Two normally nonhyperbolic points (xm , F (xm )) and (xM , F (xM )) (black dot) are also shown. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

  x x  = x 2 x 2 + bx + 1 1 − − x 3 y, K y  = 0,

(2.3)

which is a parametrized differential equation, where the slow variable y acts as a parameter. Another way to understand dynamics of system (2.2) is on the slow time scale τ := εt, instead of the fast time scale t, i.e., the next system ε

  dx x = x˙ = x 2 x 2 + bx + 1 1 − − x 3 y, dτ K dy = y˙ = y(x 2 + bx + 1) (x − hy) . dτ

(2.4)

Taking the singular limit in system (2.4) leads to the slow subsystem:   x 0 = x 2 x 2 + bx + 1 1 − − x 3 y, K y˙ = y(x 2 + bx + 1) (x − hy) ,

(2.5)

which is a differential-algebraic equation defined on the critical manifold

1 x C0 := (x, y)|y = F (x) := (x 2 + bx + 1)(1 − ) . x K

(2.6)

2.2. Critical manifold and slow-fast cycles This subsection discusses shape of the critical manifold C0 as the parameters are varied. Some 1 3 2 easy calculations show that F  (x) = − Kx 2 G(x) with G(x) = 2x + (b − K)x + K. At the  root x = (K − b)/3 of G (x) = 0, G(x) takes the value −C1 (K, b)/27, where C1 (K, b) := (K − b)3 − 27K. Fig. 2.1 (a) demonstrates the curve C1 := {C1 (K, b) = 0}. Clearly G(x) has always a negative root. The next result summarizes the properties of G(x) and F (x).

3404

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Lemma 2.2. For the functions F (x) and G(x), the following statements hold. (a) If the parameters (b, K) ∈ U1 := {(b, K)| C1 (K, b) ≤ 0, b +2 > 0 and K > 0}, the function G(x) has a unique real zero, which is negative. Consequently, F (x) is monotone for x > 0. (b) If the parameters (b, K) ∈ U2 := {(b, K)| C1 (K, b) > 0, b + 2 > 0 and 0 < K < 1}, the function G(x) = 0 has two positive zeros, which are both greater than K, and consequently F (x) is negative at these two zeros. (c) If the parameters (b, K) ∈ U := {(b, K)| C1 (K, b) > 0, b + 2 > 0 and K > 1}, the function G(x) has two positive roots xm and xM satisfying 0 < xm < xM < K. We remark that for K = 1 and b > −2, it forces that C1 (K, b) < 0, which belongs to the case (a). Proof. Statement (a) follows easily. The reasons that statements (b) and (c) hold are the following: • F (x) has a unique zero, which is at x = K, where G(K) = K(K 2 + bK + 1) > 0 for K > 0. • The local minimum of G(x) is negative taking at x = (K − b)/3 (> 0). • In the region U2 of the case (b), which is located under the line K = −b/2, one has K < (K − b)/3. This implies that K < xm < (K − b)/3 < xM , statement (b) follows. • In the region U of the case (c), which is located above the line K = −b/2, one has K > (K − b)/3. This implies that 0 < xm < (K − b)/3 < xM < K, statement (c) follows. 2 In case (a) of Lemma 2.2, since the function F (x) is monotone for x > 0, system (2.2) cannot exhibit canard cycles. In case (b) of Lemma 2.2, the fact that the function F (x) is negative at the extreme points conflicts with the biological hypothesis. Therefore, in the following we need only to investigate the dynamics of system (2.2) restricted to (b, K) ∈ U , where the critical manifold C0 is of S-shaped curve with the y-axis as its vertical asymptote. The critical manifold C0 consists of three submanifolds: C0l , C0r and C0m , which are all normally hyperbolic, and the first two are attracting, and the last one is repelling. The normally hyperbolicity is lost at both non-degenerate contact points (xm , F (xm )) and (xM , F (xM ) (note: the non-degeneracy of a contact point (x0 , y0 ) means F  (x0 ) = 0 and F  (x0 ) = 0). Obviously,  K−b K > b for (b, K) ∈ U . Furthermore, it follows from G(0) = G( K−b 2 ) = K > 0 and G ( 3 ) = 0 that 0 < xm < (K − b)/3 < xM < (K − b)/2,

for (b, K) ∈ U.

See Fig. 2.1 (b) and (c) for a geometrical description of the previous analysis. It follows from system (2.5) that for x > 0 the dynamics of the slow subsystem are governed by F  (x)

dx = g(x, F (x), η), dτ

(2.7)

with g given in (2.2) and η = (b, K, h). We first study the local dynamics of system (2.2) near the contact point (xm , F (xm )). The similar arguments work well for the contact point (xM , F (xM )). For (xm , F (xm )) we distinguish two cases:

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3405

Fig. 2.2. The U -shaped critical manifold C0 (blue curve) for system (2.2) consists of two submanifolds: C0l is attracting, while C0m is repelling. The normally nonhyperbolic point (xm , F (xm )) (black dot) is shown. Double arrows (blue) indicate fast flow, and single arrows (blue) indicate slow flow. (a) The slow manifolds Cεl and Cεm (red curve) of system (2.2) near the jump point (xm , F (xm )). (b) The slow manifolds Cεl and Cεm (red curve) of system (2.2) near the canard point (xm , F (xm )).

(a) g(xm , F (xm ), η) < 0. This is equivalent to h > F x(xmm ) . Since F  (x) = (x − xm )F  (xm ) + O((x − xm )2 ), system (2.7) is singular at xm . In this case the contact point is called a jump point [25]. (b) g(xm , F (xm ), η0 ) = 0 for some η = η0 . It means that h = F x(xmm ) . Along the critical manifold one has dx g(x, F (x), η0 ) gx (xm , F (xm ), η0 ) + O((x − xm )) = = . dτ F  (x) F  (xm ) + O((x − xm )) 2 + bx + 1) > 0 and F  (x ) > 0. One can check easily that gx (xm , F (xm ), η0 ) = F (xm )(xm m m In this case the contact point is called a canard point [25].

The local dynamics of system (2.2) near the jump point and canard point are shown in Fig. 2.2, see also [25]. 2.3. The existence of equilibria The equilibria of system (2.2) are O = (0, 0), Eb = (K, 0), and the possible positive equilibria E ∗ = (x ∗ , y ∗ ), whose abscissa x ∗ satisfy P (x) := x 3 + (K/ h + b − K)x 2 + (1 − Kb)x − K = 0. Define A := (K/ h + b − K)2 + 3(Kb − 1)

(2.8)

and B := 2(K/ h + b − K)3 − 9(1 − Kb)(K/ h + b − K) − 27K. Some computations show that the extreme value of P (x) is equal to zero if and only if

:= B 2 − 4A3 = 0.

(2.9)

3406

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Based on the above analysis, we obtain the following results. Lemma 2.3 (Existence of equilibria). (a) If > 0, then system (2.2) has a unique positive equilibrium, which is an elementary antisaddle. (b) If = 0 and √ (b1 ) A > 0 and (K/ h + b − K) + A > 0, then system (2.2) has a unique positive equilibria, which is an elementary anti-saddle; (b2 ) A = 0, then system (2.2) has a unique positive equilibria (− 13 (K/ h + b − K), 1 − 3h (K/ h + b − K)), which a degenerate equilibrium. √ (c) If = 0, A > 0 and (K/ h + b − K) + A < 0, then system (2.2) has two different positive equilibria. More precisely, √ (c1 ) If B = −2A3/2 , A > 0 and (K/ h + b − K) + A < 0, then system (2.2) has two different positive equilibria: a degenerate equilibrium (− 13 (K/ h + b − K + √ √ 1 A), − 3h (K/ h + b − K + A)) and an elementary anti-saddle equilibrium √

√

1 (− 13 (K/ h + b − K) + 2 3 A , − 3h (K/ h + b − K) + 23hA ). √ 3/2 (c2 ) If B = 2A , A > 0 and (K/ h + b − K) + A < 0, then system (2.2) has two different positive equilibria: an elementary anti-saddle equilibrium (− 13 (K/ h + b − K) − √ √ 2 A 1 (K/ h + b − K) − 23hA ) and a degenerate 3 ,− √3h 1 √ K) + A), 3h (−(K/ h + b − K) + A)). √

equilibrium ( 13 (−(K/ h + b −

(d) If < 0, K/ h < K − b − 3(1 − Kb) and −2 < b < K1 , then system (2.2) has three different positive equilibria E1 (x1 , y1 ), E2 (x2 , y2 ) and E3 (x3 , y3 ) satisfying x1 < x2 < x3 , which are all elementary equilibria, and especially E2 is a saddle. Proof. It follows from roots of cubic algebraic equations, see also [21].

2

Remark 2.1. Note that in system (2.2), f and g are independent of ε. Therefore, the positive equilibria (x ∗ , h−1 x ∗ ) of system (2.2) coincide with the positive equilibria of the slow subsystem (2.5). 2.4. Positive invariance This subsection provides some positive invariant set, which will be used later on. Lemma 2.4. The set  = {(x, y) ∈ R × R : 0 < x ≤ K; 0 ≤ y ≤ K/ h} is positively invariant for system (2.2). If (b, K) ∈ U , the set is global attracting in the first quadrant of the (x, y) plane for 0 < ε  1. Proof. Note that the x-axis and y-axis are both invariant under the flow of system (2.2). Re dx  stricting the first equation of system (2.2) on x = K gives dt  = −K 3 y < 0 for y > 0. And x=K  K 2 restricting the second equation of system (2.2) on y = Kh yields dy dt  K = ε h (x + bx + 1)(x − y= h

K) ≤ 0 for 0 < x ≤ K. This proves the positive invariance of the set  under the flow of system (2.2).

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3407

Fig. 3.1. Directional blow-up at the origin O(0, 0) of system (3.1).

When (b, K) ∈ U , the global attraction of  for system (2.2) with 0 < ε  1 follows from Fig. 2.1 (c) and the dynamics of its associated fast subsystem. 2 3. Dynamics near O(0, 0) and Eb (K, 0) 3.1. Dynamics near O(0, 0) Note that system (2.2) can be written as      dx = x 2 + b − K1 x 3 + 1 − Kb x 4 − x 3 y − K1 x 5 , dt   dy = ε (xy − hy 2 ) + b(x 2 y − hxy 2 ) + (x 3 y − hx 2 y 2 ) . dt

(3.1)

The equilibrium O(0, 0) of system (3.1) is isolated and degenerate. The next result presents its property. Proposition 3.1. The origin O(0, 0) of system (3.1) is a non-hyperbolic saddle for all values of the parameters and 0 < ε  1. Proof. Applying the horizontal blow-up transformation (x, y, t) → (r, ry, τ1 /r), and the vertical blow-up transformation (x, y, t) → (rx, r, τ1 /r), we get one saddle for each of the blow–up systems as shown in Fig. 3.1. Then the proof follows from blow–down. 2 3.2. Dynamics near Eb (K, 0) Eb (K, 0) is a boundary equilibrium of system (2.2), and it is always a saddle. For details, see [21]. Set ym = F (xm ) and yM = F (xM ). 4. Slow-fast normal form Based on the results obtained in [25] and [27], we now study the slow-fast normal form of system (2.2) near the contact point (xm , ym ) (the similar arguments work also for the contact point (xM , yM )). Rewrite system (2.2) via time rescaling dt1 = x 3 dt , and moving its contact point (xm , ym ) to the origin gives

3408

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

dx¯ 3 3 −4 −5 2 = −y¯ + x¯ 2 ((K − xm )/(Kxm ) − xm x¯ + xm x¯ + O(x¯ 3 )), dt1 dy¯ −2 = εxm (p0 + p1 x¯ + p2 x¯ 2 + y(q0 + q1 x¯ + q2 y) ¯ + O(|(x, ¯ y)| ¯ 3 )), dt1

(4.1)

where   −1 2 p0 = ym 1 − hym xm + bxm + 1), (xm −2 2 2 ym (2bhxm ym + ym hxm − bxm + 3hym − 2xm ), p1 = x m −3 2 2 ym (3bhxm ym + ym hxm − bxm + 6hym − 3xm ), p2 = −xm −1 2 q0 = x m (xm − 2hym )(xm + bxm + 1),

(4.2)

−2 2 2 (4bhxm ym + 2ym hxm − bxm + 6hym − 2xm ), q 1 = xm −1 2 (xm + bxm + 1). q2 = −hxm

Note that p0 and p1 have the relationship −1 2 2 p1 = ym (b + xm ) + ym xm − 3p0 (xm (xm + bxm + 1))−1 − p0 (xm + 2b)(xm + bxm + 1)−1 .

When (b, K) ∈ U , we can choose suitable h such that |p0 | is small enough and p1 > 0. In order to obtain the standard slow-fast normal form (3.6) in [25], we set λ := −

3 )p (K − xm 0 3/2

2p Kxm 1

.

(4.3)

3 − K = x 3 − (−2x 3 + Then it will be sufficiently small if p0 is sufficiently small. Note that xm m m 2 ) = x 2 (3x − (K − b)) < 0. In order to use the theory developed in [27] and [25], we (K − b)xm m m write system (4.1) as

du = −vh1 (u, v, λ) + u2 h2 (u, v, λ) + εh3 (u, v, λ), dt2 dv = ε (uh4 (u, v, λ) − λh5 (u, v, λ) + vh6 (u, v, λ)) , dt2 where (u, v) take values near (0, 0) and h1 (u, v, λ) = 1; h2 (u, v, λ) = 1 −

√ 2p K 2 x m p1 K 3 xm 1 u + u2 + O(|u|3 ); 3 )2 3 )3 (K − xm (K − xm

2 p2 Kxm u + O(|u|2 ); 3) p1 (K − xm √ K p1 Kxm q0 u + q2 v + O(|(u, v)|2 ). h5 (u, v, λ) = 1; h6 (u, v, λ) = √ + q1 3 3 x m p1 K − xm K − xm

h3 (u, v, λ) = 0; h4 (u, v, λ) = 1 + √

System (4.4) is a slow-fast normal form of system (2.2) near the contact point (xm , ym ).

(4.4)

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3409

−1 = 0 ⇔ h = x y −1 . Hence, It follows from (4.2) and (4.3) that λ = 0 ⇔ p0 = 0 ⇔ 1 −hym xm m m by the formulae (3.12) and (3.13) of [27] we have

∂h3 ∂h1 (0, 0, 0) = 0; a2 = (0, 0, 0) = 0; ∂u ∂u  2 + bx + 1) ∂h2 K 2 xm ym (xm m ; (0, 0, 0) = − a3 = 3 )2 ∂u (K − xm  √ 2 + 2bx + 3) −K xm ym (xm x 2 + bxm + 1 ∂h4 m ; a5 = h6 (0, 0, 0) = − m a4 = , (0, 0, 0) =  2 + bx + 1(K − x 3 ) ∂u xm ym xm m m a1 =

and 2 3 −2 ˆ Aˆ = −a2 + 3a3 − 2a4 − 2a5 = (xm ym (xm + bxm + 1))− 2 (K − xm ) A1 , 1

where    2    3 2 2 4 5 3 Aˆ 1 =2 K − xm + bxm + 1 + Kxm ym bxm − xm + 3 K − 4bxm − 2xm − 6xm xm . (4.5) We can see from (4.2) and (4.3) that λ=−

3 )(x − hy )(x 2 + bx + 1) (K − xm m m m m . 1/2 2 2 Kym (2bhxm ym + ym hxm − bxm + 3hym − 2xm )3/2

(4.6)

According to formulae (3.15) and (3.16) in [27],√ we can calculate the singular Hopf bifur√ cation curve λH ( ε) and maximal canard curve λc ( ε) for the slow-fast normal form (4.4) as follows: 1/2 √ a 1 + a5 1 2 −1 −1 λH ( ε) = − ym ε + O(ε 3/2 ) ε + O(ε 3/2 ) = (xm + bxm + 1)xm 2 2 and √ a 1 + a5 ˆ λc ( ε) = −( + A/8)ε + O(ε3/2 ) 2    1/2 1 2 −1 −1 ˆ + bxm + 1)xm ym − A/8 ε + O(ε 3/2 ). (xm = 2 Some further calculations show that the asymptotic expansion of singular Hopf bifurcation √ √ curve hH ( ε) and maximal canard curve hc ( ε) for system (2.2) can be written as ⎞ ⎛ 3/2 ym 3/2 (y (b + x ) + ) Kx √ m m m xm ⎝ xm hH ( ε) = ε + O(ε 3/2 )⎠ 1 + 3/2 3 2 1/2 ym 2ym (K − xm )(xm + bxm + 1) and

(4.7)

3410

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Fig. 5.1. Five cases.

⎛ x 2 +bx +1 2 ˆ ( 12 ( m xm ymm )1/2 − A/8)Kx √ xm ⎝ m (ym (b + xm ) + 1+ hc ( ε) = 3 )y (x 2 + bx + 1) ym (K − xm m m m

ym 3/2 xm )

⎞ ε + O(ε 3/2 )⎠ , (4.8)

where

1 2



2 +bx +1 xm m xm ym

−

Aˆ 8

can be simplified as

2 + bx + 1)(K − x 3 )2 − x y K((bx − x 2 + 3)K − 4bx 4 − 2x 5 − 6x 3 ) 2(xm m m m m m m m m m  . 2 3 2 8 xm ym (xm + bxm + 1)(K − xm )

Hereafter, a singular Hopf bifurcation occurs when the eigenvalues become singular as ε → 0. That is, the√linearized system has a pair of ‘singular eigenvalues’ at the Hopf bifurcation point λ = λH ( ε), i.e. μ(λ, ε) = α(λ, ε) + iβ(λ, ε), so that α(λH , ε) = 0, β(λH , ε) = 0, with limε→0 |β(λH , ε)| = 0 on the fast time t2 and limε→0 |β(λH , ε)| = ∞ on the slow time t3 := εt2 . See [3] for more details. 5. Dynamics in the interior of the first quadrant According to Lemma 2.2 the first quadrant of the phase space of system (2.2) is separated in three domains defined in the cases (a), (b) and (c) of the lemma. In the cases (a) and (b), the dynamics of system (2.2) for 0 < ε  1 is simple. So, we restrict our attention to the case (c) of Lemma 2.2, i.e. (b, K) ∈ U = {(b, K)| (K − b)3 > 27K, b > −2, K > 1}. Then the critical manifold C0 of system (2.2) has the S-shaped structure with two contact points (xm , ym ) and (xM , yM ) satisfying 0 < xm < xM < K, which are located in the first quadrant. According to the number of positive equilibria of system (2.2), we separate our next discussion in four subsections. 5.1. The system has a unique positive equilibrium In this case the parameters (K, b, h) belong to the region S1 := {(K, b, h)|(K, √b) ∈ U and > 0}, or S2 := {(K, b, h)|(K, b) ∈ U, = 0, A > 0 and K/ h + b − K + A > 0} ∪ {(K, b, h)|(K, b) ∈ U, = 0 and A = 0}, with A and defined in (2.8) and (2.9) respectively, and the unique equilibrium is (x0 , h−1 x0 ), which by Remark 2.1 is also the unique positive equilibrium of the slow subsystem (2.2). Taking different value of x0 ∈ (0, K), there are five distinguished cases to be considered as shown in Fig. 5.1. Since the cases (d) and (e) in Fig. 5.1 could be investigated in the same way as the cases (b) and (a) respectively, we analysis only the first three cases. For cases (d) and (e), we just

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3411

Fig. 5.2. The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5). The critical manifold C0 consists of three submanifolds: C0l and C0r are attracting parts, while C0m is repelling. The equilibrium (x0 , h−1 x0 ) (black dot) of system (2.2) and the two contact points (xm , F (xm )) and (xM , F (xM )) are shown as well. Double arrows indicate fast flow, and single arrows indicate slow flow on the critical manifold.

study dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) and omit the similar description for system (2.2) or (2.4). 5.1.1. Case (a): x0 < xm We can check that the equilibrium (x0 , h−1 x0 ) is a stable hyperbolic node of the slow subsystem (2.5). Then Fenichel’s theory verifies that (x0 , h−1 x0 ) is a stable hyperbolic node of system (2.2). Furthermore, the fast subsystem (2.3) and the slow subsystem (2.5) have the dynamics as shown in Fig. 5.2 (a). The next result verifies the global attraction of (x0 , y0 ) with y0 = h−1 x0 . Theorem 5.1. If the parameters (K, b, h) ∈ S1 ∪ S2 and the unique positive equilibrium (x0 , y0 ) lies on the normally hyperbolic attracting critical manifold C0l , then (x0 , y0 ) is a global attractor for system (2.2) in the interior of the first quadrant. We remark that Hsu and Huang in [20] have obtained that the unique positive equilibrium is globally asymptotically stable when b > 0. Here we further prove that the statement holds for b > −2 and ε sufficiently small. Proof. The line x = xm divides the interior of the first quadrant into two disjoint regions which we label (from left to right) R1 and R2 . First, it follows from Lemma 2.4 that all orbits in the first quadrant of system (2.2) are positively bounded. Next we prove that all these orbits are attracted to (x0 , y0 ). For doing so, we define a Dulac function D(x, y) = (x 2 + bx + 1)x −2 y −2 ,

(x, y) ∈ R1 .

An easy computation shows that for (x, y) ∈ R1 and 0 ≤ ε  1 ∂f (x,y)D(x,y) ∂x

+

∂g(x,y)D(x,y) ∂y

= −x

2 +bx+1

x2y2

 x2 + ε x 2 +bx+1

1 (2x 3 Kx 2

 + (b − K)x 2 + K) < 0.

By Dulac’s criterion, system (2.2) has no closed orbits lying entirely in R1 .

3412

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Fig. 5.3. (a) Orbits of system (2.2) starting at points (0.1, 2), (0.1, 0.5), (0.1, 0.1), (3, 1) and (1.5, 2) for ε = 0.01. (b) Behavior of the prey and predator with the initial densities x(0) = 0.1 and y(0) = 0.1.

Since the contact point (xM , yM ) is a jump point (see Fig. 2.2 (a)), according to Fenichel’s theory and Theorem 2.1 of [25] on the analysis of the jump point, we get together with Fig. 5.2 (a) that the trajectory starting at a point in R2 must eventually cross the line x = xm , and it has (x0 , y0 ) as its ω limit. It completes the proof of the theorem. 2 Next example illustrates our Theorem 5.1. Example 5.1. Set ε = 0.01, K = 3, b = −9/5 and h = 1/3. Then system (2.2) has the unique positive equilibrium (x0 , y0 ) = (0.3706, 1.1119) and the two contact points (xm , ym ) = (1.0567, 0.1315) and (xM , yM ) = (2.0393, 0.2337) of the critical manifold C0 . Fig. 5.3 (a) shows some typical orbits of system (2.2), and the behavior of the prey and the predator that system (2.2) is modeling in Fig. 5.3 (b). 5.1.2. Case (b): x0 = xm Now the unique positive equilibrium (x0 , h−1 x0 ) of system (2.2) coincides with the contact point (xm , ym ), and it becomes a canard point [25] (see Fig. 2.2 (b)). Note that h = xymm at the moment. The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) are shown in Fig. 5.2 (b). Let s0 = yM − ym . For s ∈ (0, s0 ), let xl (s) < xm (s) < xr (s) be the three distinct roots of F (x) = h−1 x0 + s. Set xl (0) = xm (0) = xm and xm (s0 ) = xr (s0 ) = xM , respectively. Define the continuous family of slow-fast cycles γ (s) for s ∈ [0, 2s0 ] as follows: • γ (0) is the canard point (xm , ym ); γ (s0 ) is the transitory slow-fast cycle (see Fig. 1.1 (b)); and γ (2s0 ) is the common slow-fast cycle (see Fig. 1.1 (d)). • For s ∈ (0, s0 ), γ (s) = {(x, F (x)) : x ∈ [xl (s), xm (s)]} ∪ {(x, s) : x ∈ [xl (s), xm (s)]}, see Fig. 1.1 (a). • For s ∈ (s0 , 2s0 ), γ (s) = {(x, F (x)) : x ∈ [xl (s0 ), xm (2s0 − s)]} ∪ {(x, 2s0 − s) : x ∈ [xm (2s0 − s), xr (2s0 − s)]} ∪ {(x, F (x)) : x ∈ [xM , xr (2s0 − s)]} ∪ {(x, yM ) : x ∈ [xl (s0 ), xM ]}, see Fig. 1.1 (c).

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3413

The following results show the existence of canard cycles and of canard explosion of system (2.2). Theorem 5.2. Assume that (K, b, h) ∈ S1 ∪ S2 , and that the unique positive equilibrium (x0 , y0 ) of the slow subsystem of system (2.2) coincides with the contact point (xm , ym ). Let h and λ satisfy the relationship (4.6). For 0 < ε  1, the following statements hold for system (2.2). (a) There exists an h0 > 0 such that for |h − xymm | < h0 , in a suitable neighborhood of (x0 , y0 ), system (2.2) has precisely one equilibrium point E → (x0 , y0 ) as (h, ε) → ( xymm , 0). Further√ more, there exists a√curve h = hH ( ε)√depending on ε of Hopf bifurcations such that E is stable for h < hH ( ε), where h = hH ( ε) is given by (4.7). Moreover, the Hopf bifurcation is nondegenerate when Aˆ 1 = 0, and it is supercritical if Aˆ 1 < 0 and subcritical if Aˆ 1 > 0, where Aˆ 1 is given by (4.5). √ √ (b) There exists a parameterized smooth family s → (h(s, ε), γ (s, ε)) of periodic√orbits bifurcated from the slow-fast cycles γ (s) for each s ∈ (0, 2(yM −√ym )) and γ (s, ε) → γ (s, 0) = γ (s) as ε → 0. Furthermore, there exists a curve h = hc ( ε) depending on ε, for β ∈ (0, 1), the canard explosion occurs when s ∈ [ε β , 2(yM − ym ) − ε β ], where √ √ β |h(s, ε) − hc ( ε)| ≤ e−1/ε , √ and h = hc ( ε) is given by (4.8). Proof. Denote by η0 the quantity η = (b, K, h)|h= xm . We can see from (2.2) that f (x0 , y0 , η0 ) = ym

0, fx (x0 , y0 , η0 ) = 0, g(x0 , y0 , η0 ) = 0, fxx (x0 , y0 , η0 ) = −xK −1 (6x 2 + 2(b − K)x) > 0, fy (x0 , y0 , η0 ) = −x03 < 0, gx (x0 , y0 , η0 ) = y0 (x02 + bx0 + 1) > 0, gh (x0 , y0 , η0 ) = −y02 (x02 + bx0 + 1) < 0. Then the unique positive equilibrium (x0 , y0 ) is a generic canard point at h = xymm . As shown in Section 4, system (2.2) is equivalent to the slow-fast normal form (4.4) near the generic canard point (x0 , y0 ) = (xm , ym ). For system (4.4), (0, 0) is a generic canard point for λ = 0. (a) By Theorem 3.1 of [27] on the analysis of Hopf bifurcations, there exist ε0 > 0 and λ0 > 0 such that for 0 < ε < ε0 and |λ| < λ0 , in a suitable neighborhood of the origin, system (4.4) has there exists precisely one equilibrium point p∗ with p∗ → (0, 0) as (λ, ε) → (0, 0). Furthermore, √ √ a curve λ = λH ( ε) of Hopf bifurcations such that p∗ is stable for λ < λH ( ε) and 1/2 √ 1 2 −1 −1 λH ( ε) = ym ε + O(ε 3/2 ). (xm + bxm + 1)xm 2 Furthermore, the Hopf bifurcation is nondegenerate when Aˆ 1 = 0, and it is supercritical if Aˆ 1 < 0 and subcritical if Aˆ 1 > 0. Moreover, from (4.6) it follows that 3 )(x 2 + bx + 1) − (K − x 3 )(x − hy )(x 2 + bx + 1)p ∗ dλ ym (K − xm m m m m m m m 1 = 2 − bx 2 + 3hy − 2x )3 dh K 2 ym (2bhxm ym + ym hxm m m m 2 − bx 2 + 3hy − 2x )1/2 (2bx y + y x 2 + 3y ) > 0. with p1∗ = 32 Kym (2bhxm ym + ym hxm m m m m m m m m dλ > 0. Hence the algebraic equation (4.6) via λ = Then it is feasible to choose h such that dh 1/2

3414

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

√ √ √ λH (√ε) has a unique solution hH ( ε), which satisfies that λ < λH ( ε) if and only if h < hH ( ε). The statement follows. (b) By Theorems 3.3 and 3.5 of [27] on the canard explosion, system (4.4) has a family of periodic orbits √ √ s → (h(s, ε), γ (s, ε)),

s ∈ (0, 2(yM − ym ))

√ which is C k -smooth in (s, ε). For ν ∈ (0, 1), the canard explosion always occurs when s ∈ [ε ν , 2(yM − ym ) − ε ν ], √ √ ν | h(s, ε) − hc ( ε) |≤ O(e−1/ε ), √ where hc ( ε) is given by (4.8). Then the statement (b) holds for system (2.2). It completes the proof of the theorem. 2 Remark 5.1. In Theorem 5.2 the properties of the Hopf bifurcation depend on the sign of Aˆ 1 . Because of the complexity of the expression of Aˆ 1 we cannot give a general discrimination on the sign of Aˆ 1 . The next examples show that all three cases Aˆ 1 < 0, Aˆ 1 > 0 and Aˆ 1 = 0 could ˆ Set K = 225/32, b = 29/32 happen. Since Aˆ and Aˆ 1 have the same sign, we state them in A. and h = 0.6205, one has Aˆ < 0. Set K = 10, b = −11/10 and h = 1.3058, one has Aˆ > 0. Set K = 1.5, b = −1.9473 and h = 65.5777, one has Aˆ = 0. Next numerical example exhibits clearly the phenomena stated in Theorem 5.2. Example 5.2. Set ε = 0.0001, K = 225/32, b = 29/32 and h = 0.6205. Then, system (2.2) has a unique positive equilibrium (x0 , y0 ) = (1.5, 2.4174). The two contact points are (xm , ym ) = (x0 , y0 ) and (xM , yM ) = (2.5, 2.4529). When h increases, system (2.2) exhibits the supercritical singular Hopf bifurcation, the canard explosion and the relaxation oscillation as shown in Fig. 5.4. √ In Fig. 5.5, we √draw a sketch of the singular Hopf bifurcation curve hH ( ε) and maximal canard curve hc ( ε) in (ε, h) plane for Example 5.2. 5.1.3. Case (c): xm < x0 < xM The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) are shown in Fig. 5.6 (a). For 0 < ε  1, we have the next result. Theorem 5.3. Assume that (K, b, h) ∈ S1 ∪ S2 and x0 lies on the normally hyperbolic repelling critical submanifold C0m . Then for each fixed ε > 0 sufficiently small, system (2.2) has a unique limit cycle, say γε , which is located in a small tubular neighborhood of the common slow-fast cycle γ (2(yM − ym )). Furthermore, γε converges to γ (2(yM − ym )) in the Hausdorff distance as ε → 0. Proof. Now the unique positive equilibrium (x0 , y0 ) of the slow subsystem is an unstable singular node, and the critical submanifold C0m is normally hyperbolic and unstable, see Fig. 5.6 (a). Moreover, the two contact points (xm , ym ) and (xM , yM ) are both jump points. Fenichel’s theory shows that the equilibrium (x0 , y0 ) is still an unstable node of system (2.2) for 0 < ε  1, and

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3415

Fig. 5.4. Canards for system (2.2) when ε = 0.0001. (a) The unique positive equilibrium is stable. (b) A periodic orbit generated in a singular Hopf bifurcation. (c) Canard cycle without head. (d) Canard cycle with head. (d) Relaxation oscillation. (f) The amplitude (vertical axis) of the limit cycle when h is varied (denote the amplitude by the maximum value of x).

√ √ Fig. 5.5. The singular Hopf bifurcation curve hH ( ε) (solid blue) and maximal canard curve hc ( ε) (solid red) for system (2.2) when K = 225/32, b = 29/32 (Example 5.2). Consider ε > 0 fixed. By increasing h, a singular Hopf bifurcation takes place, giving rise to a small stable limit cycle called Hopf cycle with amplitude O(ε). The amplitude √ of the Hopf cycle is growing as h increases. When h reaches the dashed line hs ( ε), the Hopf cycle becomes a canard √ cycle without head. Along the curve hr ( ε) the family of canard cycles ends at a relaxation oscillation. The black dashed √ √ √ √ lines hs ( ε) and hr ( ε) mark the beginning and ending of the canard explosion. Furthermore, |hi ( ε) − hc ( ε)| = O(e−K0 /ε ), i = s, r, for some K0 > 0 as ε → 0.

3416

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Fig. 5.6. (a) The dynamics of fast subsystem (2.3) and of slow subsystem (2.5) for xm < x0 < xM . (b) The sketch of the relaxation oscillation for system (2.2) with 0 < ε  1. (c) The relaxation oscillation for system (2.2) when ε = 0.0001, K = 225/32, b = 29/32 and h = 14400/17549. (d) The behavior of x and y in Fig. 5.6 (c) depending on time t . (e) and (f) show the dynamics of the fast and of the slow subsystems (2.3) and (2.5) for x0 = xM and x0 > xM , respectively.

that the critical attracting submanifolds C0l and C0r perturb to nearby attracting slow manifolds Cεl and Cεr , respectively. By Theorem 2.1 in [25], the attracting slow manifold Cεl (resp. Cεr ) can be continued past the jump point (xm , ym ) (resp. (xM , yM )) and then follow approximately a layer of the fast subsystem (2.5) until it arrives at the attracting slow manifold Cεr (resp. Cεl ). Take a small horizontal section transversal to C0l , see Fig. 5.6 (b), and choose any two trajectories 1ε and 2ε of system (2.2) starting on . For ε > 0 sufficiently small, it follows from Fenichel’s theory that 1ε and 2ε are attracted to Cεl with exponential rate O(e−1/ε ). Theorem 2.1 in [25] shows that 1ε and 2ε pass by the jump point (xm , ym ) contracting exponentially toward each other until they reach a neighborhood of C0r . Then 1ε and 2ε are attracted toward Cεr with exponentially small rate, and pass by the jump point (xM , yM ) once more with exponentially small attracting rate, and finally they return to . Following the positive orbits of system (2.2) starting on , we can get an attracting map  :

→ with an exponential small contracting rate O(e−1/ε ). The contraction mapping theorem verifies that  has a unique fixed point on . Consequently, system (2.2) for 0 < ε  1 has a unique limit cycle, say γε , passing through , which is located in a small neighborhood of the common slow-fast cycle γ (2(yM − ym )). Applying again Fenichel’s theory we get that γε → γ (2(yM − ym )) as ε → 0. We claim that the limit cycle γε is a unique periodic orbit of system (2.2) in the first quadrant of the phase plane. Indeed, any orbit of system (2.2) starting above the line y = yM (resp. under the line y = ym )) will rapidly move to the vicinity of C0l (resp. of C0r ), then follows the slow manifold Cεl (resp. Cεr ). These facts imply that system (2.2) cannot have a limit cycle located

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3417

Fig. 5.7. Six cases: (a) x1 < xm < x2 < xM , (b) x1 = xm < x2 < xM , (c) xm < x1 < x2 < xM , (d) xm < x1 < x2 < xM , (e) xm < x1 < x2 = xM , (f) xm < x1 < xM < x2 . In cases (a), (b) and (c), the curve y = F (x) is tangent to the straight line y = h−1 x at the point (x2 , F (x2 )). In cases (d), (e) and (f), the curve y = F (x) is tangent to the straight line y = h−1 x at the point (x1 , F (x1 )).

outside a small neighborhood of the common slow-fast cycle γ (2(yM − ym )). Hence the limit cycle γε is the unique one for system (2.2) in the first quadrant. The claim follows. It completes the proof of the theorem. 2 By numerical simulation we exhibit the limit cycle obtained in Theorem 5.3 as shown in Fig. 5.6 (c) and (d). 5.1.4. Case (d) and (e): x0 = xM and x0 > xM In the cases x0 = xM and x0 > xM , the fast and slow subsystems (2.3) and (2.5) have the dynamics shown in Fig. 5.6 (e) and (f), respectively. Following similar arguments as those treated in the cases x0 = xm and x0 < xm , we get the same conclusions as stated in Theorem 5.2 and Theorem 5.1, respectively. The details are omitted. 5.2. The system has two positive equilibria Now the parameters (K, √ b, h) belong to the set S3 := {(K, b, h)|(K, b) ∈ U, = 0, A > 0 and Kh−1 + b − K + A < 0}, where A and are defined in (2.8) and (2.9) respectively, and system (2.2) has exactly two positive equilibria E1 (x1 , h−1 x1 ) and E2 (x2 , h−1 x2 ), with one degenerate, which are both located on the critical manifold C0 . Here we assume that x1 < x2 . Then the locations of these two equilibria can be distinguished into six cases, as shown in Fig. 5.7, where in the first three cases E2 (x2 , y2 ) is degenerate, whereas in the last three cases E1 (x1 , y1 ) is degenerate. Since the cases (d), (e) and (f ) in Fig. 5.7 can be investigated in the completely same way as the cases (c), (b) and (a), respectively, we consider only the first three cases, where E2 is a degenerate equilibrium. Lemma 5.4. Set H1 := −10hx23 + 6((K − b)h − K)x22 + 3(Kb − 1)hx2 + Kh,

(5.1)

√ where x2 = − 13 (K/ h + b − K) + 13 A in cases (a), (b) and (c) of Fig. 5.7 and x2 = √ − 13 (K/ h + b − K) − 13 A in cases (d), (e) and (f) of Fig. 5.7. For system (2.2), if H1 = 0, then E2 is a saddle-node, and its local structure is shown in Fig. 5.8 (b), where Ll and Lr are the two branches of the unstable manifold associated to the positive eigenvalue. Proof. It follows from normal form near (x2 , y2 ) via the techniques and Theorem 2.19 in [10]. 2

3418

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Fig. 5.8. (a) The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for 0 < x1 < xm < x2 < xM . (b) The phase portrait in the neighborhood of saddle-node E2 . (c) Two typical heteroclinic orbits connecting the equilibria E1 and E2 for the numerical example.

Remark 5.2. Next example illustrates that H1 = 0 can happen. For instance, K = 10, b = −9/5 and h = 2.7061 verifies H1 = −108.2759892 = 0. 5.2.1. Case (a): x1 < xm < x2 < xM In this case, the two positive equilibria E1 (x1 , h−1 x1 ) and E2 (x2 , h−1 x2 ) lie on the normally hyperbolic attracting submanifold C0l and repelling submanifold C0m , respectively. The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) are shown in Fig. 5.8, where E1 and E2 are respectively a stable node and a saddle-node of the slow subsystems (2.5). By Fenichel’s theory, E1 perturbs to a stable node of system (2.2) for 0 < ε  1. Recall from the above analysis that E2 is a saddle-node of system (2.2) if H1 = 0. Next result provides the existence of heteroclinic orbits connecting the equilibria E1 and E2 . Theorem 5.5. Assume (K, b, h) ∈ S3 and 0 < x1 < xm < x2 < xM . For 0 < ε  1, the following statements hold. (a) System (2.2) has neither canard cycle nor relaxation oscillation. (b) If H1 given by (5.1) is not equal to zero, system (2.2) has infinitely many heteroclinic orbits connecting E1 and E2 . Proof. (a) It follows from the facts that system (2.2) has neither canard points nor slow-fast cycles. (b) By Fenichel’s theory, C0r and C0l perturb to the nearby slow submanifolds Cεr and Cεl . By Theorem 2.1 in [25] the slow submanifold Cεr can be continued and passes the fold point (xM , F (xM )). Then it follows approximately a layer of the fast subsystem (2.3) until it arrives in a neighborhood of the slow submanifold Cεl , and it finally is attracted to the stable node E1 . The branch Lr of the unstable manifold of E2 is first attracted to a vicinity of C0r in a very fast speed, and follow C0r slowly until it arrives nearby the fold point (xM , F (xM )), then it passes the fold point. After that, it moves fast to a neighborhood of C0l , and is attracted to the stable node E1 , and forms a heteroclinic orbit. By the dynamics of the fast and slow subsystems, see Fig. 5.8 (a), it follows easily that the unstable branch Ll is also heteroclinic to E1 . Using similar arguments, one can prove that all orbits emanating from the center manifolds W u , except the unique stable one, of system (2.2) at E2 positively flight to E1 , and form infinitely many heteroclinic orbits. This proves the theorem. 2

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3419

Fig. 5.9. (a) The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for x1 = xm < x2 < xM . (b) A schematic diagram of the proof of Theorem 5.6: Cεl , Cεm and Cεr are slow manifolds. l , m and r are three horizontal sections.

Next, we give a numerical example illustrating Theorem 5.5. Example 5.3. Set ε = 0.001, K = 10, b = −9/5 and h = 2.7061. Then system (2.2) has two positive equilibrium: a local stable node E1 (x1 , y1 ) = (0.7369, 0.2723) and a degenerate equilibrium E2 (3.6839, 1.3613). The two contact points are (xm , ym ) = (1.0113, 0.1799) and (xM , yM ) = (5.7487, 1.7527). Fig. 5.8 (c) shows two typical heteroclinic orbits of system (2.2). 5.2.2. Case (b): x1 = xm < x2 < xM In this case, the equilibrium E1 (x1 , h−1 x1 ) coincides with the contact point (xm , ym ) and the equilibrium E2 (x2 , h−1 x2 ) lies on C0m . The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) are shown in Fig. 5.9 (a). By [25], E1 is a canard point. Similar to the proof for (a), E2 is a saddle-node provided H1 = 0. Here we also assume that H1 = 0. Next we show the existence of homoclinic connections to the saddle-node E2 . Theorem 5.6. Assume (K, b, h) ∈ S3 , H1 = 0 and 0√< x1 = xm < x2 < xM < K. Then for 0 < ε  1, there exists a continuous function h = hc ( ε) having the expression (4.8) such that system √ (2.2) has a unique orbit homoclinic to the saddle-node E2 (x2 , y2 ) if and only if h = hc ( ε). Furthermore, if it is the case, system (2.2) has a limit cycle ε , which approaches either to the canard-type slow-fast cycle with head Fig. 1.1 (c) or to the common slow-fast cycle Fig. 1.1 (d). Proof. Similar to the proof of Theorem 5.5, one has: C0l , C0m and C0r perturb to the nearby slow submanifolds Cεl , Cεm and Cεr respectively, and Cεr can pass the jump point (xM , yM ) and follow approximately a layer of the fast subsystem (2.3) to a neighborhood of C0l . Let γ be one of the orbits emanating from the infinitely many center manifolds of the saddle-node E2 . It will finally reach the O(ε) neighborhood of C0l at an exponential rate O(e−1/ε ). By Fenichel’s theory and the structure of E2 , one gets that the unique stable branch of the center manifolds of E2 must coincide with one of the slow manifolds, saying Cεm∗ , bifurcating from C0m and with E2 located on it. Fig. 5.9 (b) shows this information. Thus, there is unique way that the homoclinic orbit (if exists) return to the saddle-node. However, it can leave the saddle-node along any one (priori) of the infinitely many orbits composing the unstable set W u of E2 .

3420

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Now turn to any chosen orbit γ emanating from E2 , since it is rapidly attracted to a O(ε) neighborhood of C0l , and then follow C0l positively, during this last period γ is a slow manifold bifurcating from the critical manifold C0l , and becomes a Cεl∗ . Next we analyze the conditions such that Cεl∗ and Cεm∗ could connect. For j = l, m, let j = {(x, ρ 2 + ym ) : x ∈ Ij } be the horizontal sections illustrated in Fig. 5.9 j (b) for sufficiently small ρ and suitably small intervals Ij for which j is centered at C0 . Let 1 be the transition map from l to m defined by following the flow of system (2.2). Set j∗ qj,ε = j ∩ Cε . Theorem 3.2 in [27] shows that a connection from Cεl∗ to Cεm∗ exists, i.e. √ 1 (ql,ε ) = qm,ε , if and only if there exists a continuous function hc ( ε) which is given by (4.8). √ This proves that under the condition h = hc ( ε) given by (4.8), the orbit γ is homoclinic to E2 . Clearly E1 is located in the region limited by γ . The uniqueness of such a homoclinic orbit follows clearly from the fact that there is a unique orbit attending to E2 positively as mentioned before. This proves the first part of the theorem. For proving the second part, we define the Poincaré map on l in Fig. 5.9 (b) for system (2.2). For any orbit starting on l and located outside the homoclinic orbit γ , it first follows C0l until arriving at a small neighborhood of the canard point (xm , ym ), and passes it in the way either flighting to the vicinity of C0r or following C0m a while and then flighting to C0r (this fact holds because of the existence of the homoclinic orbit γ at E2 ), then it follows C0r until arriving at a small vicinity of the jump point (xM , yM ), and jumps from this point to the vicinity of C0l , next it follows C0l again and arrives at l . Theorem 2.1 of [25] shows that the Poincaré map is contract, and so it has a unique fixed point on l . This fixed point provides a limit cycle of system (2.2) for ε > 0 sufficiently small, which is uniqueness for system (2.2) in the first quadrant. It completes the proof. 2 5.2.3. Case (c): xm < x1 < x2 < xM The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) are shown in Fig. 5.10 (a). For 0 < ε  1, we have the following results. Theorem 5.7. Assume that (K, b, h) ∈ S3 and x1 and x2 lie on the normally hyperbolic repelling critical submanifold C0m . Then for each fixed ε > 0 sufficiently small, system (2.2) has a unique limit cycle, say 2ε , which is a relaxation oscillation and is located in a small tubular neighborhood of the common slow-fast cycle 20 . Furthermore, 2ε converges to 20 in the Hausdorff distance as ε → 0. Proof. Now (xm , ym ) and (xM , yM ) are both jump points. We define a transversal section as shown in Fig. 5.10 (b), and a Poincaré map on it. This map is contract, and it has a unique fixed point on , which provides a limit cycle of system (2.2), a relaxation oscillation, when 0 < ε  1. Moreover, this limit cycle is a unique periodic orbit of system (2.2) in the first quadrant. 2 5.2.4. Case (d): xm < x1 < x2 < xM , case (e): xm < x1 < x2 = xM and case (f): x m < x1 < xM < x2 < K The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) are shown in Fig. 5.10 (c), (d) and (e). The same proofs as those of Theorems 5.5, 5.6 and 5.7 can yield results similar to Theorem 5.7 for xm < x1 < x2 < xM , to Theorem 5.6 for xm < x1 < x2 = xM and

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3421

Fig. 5.10. (a) The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for xm < x1 < x2 < xM . (b) The sketch of the relaxation oscillation of system (2.2) for xm < x1 < x2 < xM . (c), (d) and (e) show the dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for xm < x1 < x2 < xM , xm < x1 < x2 = xM and xm < x1 < xM < x2 , respectively. (f) Phase portrait of system (2.2) when ε = 0.001, K = 2.9153, b = −1.8135 and h = 9.2796. (g) Behavior of the prey and predator in (f) starting at initial point (1.251, 0.1352) for t = 0. (h) A heteroclinic orbit connecting equilibria E1 and E2 . (i) Magnification of the heteroclinic orbit in (h) near the equilibrium E1 .

to Theorem 5.5 for xm < x1 < xM < x2 . The details are omitted. Next, we give two numerical examples. Example 5.4. Set ε = 0.001, K = 2.9153, b = −1.8135 and h = 9.2796. System (2.2) has two positive equilibria: E1 (x1 , y1 ) = (1.2073, 0.1301) and E2 (2, 0.2155), and two contact points: (xm , ym ) = (1.0551, 0.1208) and (xM , yM ) = (2, 0.2155). Figs. 5.10 (f) and (g) show that system (2.2) exhibits Hopf bifurcation phenomenon near equilibrium E2 when h is varied. Example 5.5. Set ε = 0.001, K = 10, b = −19/10 and h = 11.4674. System (2.2) has two positive equilibria: E1 (x1 , y1 ) = (1.0595, 0.0924) and E2 (8.9090, 0.7769), and two contact points: (xm , ym ) = (1.0056, 0.0900) and (xM , yM ) = (5.8014, 1.7104). Figs. 5.10 (h) and (i) show a heteroclinic orbit of system (2.2).

3422

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Fig. 5.11. (a) The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for x1 < xm < x2 < xM < x3 . (b) A schematic diagram of heteroclinic orbits connecting equilibria E1 and E2 , and E2 and E3 of system (2.2) for x1 < xm < x2 < xM < x3 . (c) A numerical example shows the existence of the heteroclinic orbits.

5.3. The system has three positive equilibria In this case the parameters (K, b, h) belong to the set S4 := {(K, b, h)|(b, K) ∈ U, < √ 0, K/ h < K − b − 3(1 − Kb) and − 2 < b < 1/K}, and system (2.2) has the three positive equilibria E1 (x1 , h−1 x1 ), E2 (x2 , h−1 x2 ) and E3 (x3 , h−1 x3 ). According to the relative position of the prey isocline y = F (x) and the predator isocline y = h−1 x, there are nine different cases to be considered, which are: (a) x1 < xm < x2 < xM < x3 , (b) x1 < xm < x2 < x3 < xM , (c) xm < x1 < x2 < xM < x3 , (d) xm < x1 < x2 < x3 < xM , (e) x1 = xm < x2 < x3 < xM , (f) xm < x1 < x2 < x3 = xM , (g) x1 = xm < x2 < xM < x3 , (h) x1 < xm < x2 < x3 = xM and (i) x1 = xm < x2 < x3 = xM . 5.3.1. Case (a): x1 < xm < x2 < xM < x3 The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for x1 < xm < x2 < xM < x3 are shown in Fig. 5.11 (a). By Fenichel’s theory and Fig. 5.11 (a) one gets that system (2.2) has also E1 and E3 as stable nodes, and E2 as a saddle. Theorem 5.8. Assume (K, b, h) ∈ S3 , and x1 < xm < x2 < xM < x3 . For 0 < ε  1, system (2.2) has heteroclinic orbits connecting the equilibria E1 and E2 , and E2 and E3 , and has no periodic orbits in the first quadrant. Proof. Note that E1 and E3 are the only two attractors of system (2.2) in the first quadrant. Next proof follows from Fenichel’s theory and some easy calculations. 2 Next, we give a numerical example. Example 5.6. Set ε = 0.001, K = 10, b = −9/5 and h = 4. System (2.2) has three elementary positive equilibria: E1 (x1 , y1 ) = (0.8417, 0.2104), E3 (6.6797, 1.6699) and E2 (1.7786, 0.4447), and two contact points: (xm , ym ) = (1.0113, 0.1799) and (xM , yM ) = (5.7487, 1.7527). Fig. 5.11 (c) shows that heteroclinic orbits connecting the equilibria E1 and E2 , and E2 and E3 . 5.3.2. Case (b): x1 < xm < x2 < x3 < xM The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) when x1 < xm < x2 < x3 < xM are shown in Fig. 5.12 (a). By Fenichel’s theory and Fig. 5.12 (a) one gets that

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3423

Fig. 5.12. (a) The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for x1 < xm < x2 < x3 < xM . (b) A schematic diagram of heteroclinic orbits of system (2.2): two heteroclinic orbits (solid red) connecting equilibria E1 and E2 and three heteroclinic orbits (solid blue) connecting equilibria E1 and E3 . (c) A numerical example for x1 < xm < x2 < x3 < xM . (d) The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for xm < x1 < x2 < xM < x3 .

system (2.2) has a stable node E1 , a saddle E2 and a unstable node E3 . Furthermore, we have the following results. Theorem 5.9. Assume (K, b, h) ∈ S3 , and x1 < xm < x2 < x3 < xM . For 0 < ε  1, system (2.2) has two heteroclinic orbits connecting the equilibria E1 and E2 , has infinitely many heteroclinic orbits connecting E1 and E3 , and has no periodic orbits in the first quadrant. Proof. Note that the critical manifold contains neither canard points nor slow-fast cycles. Next proof is easy and is omitted. Fig. 5.12 (b) illustrates the heteroclinic orbits. 2 Next, we give a numerical example. Example 5.7. Set ε = 0.001, K = 10, b = −9/5 and h = 3. System (2.2) has the positive equilibria: E1 (x1 , y1 ) = (0.7614, 0.2538), E3 (5.1602, 1.7201) and E2 (2.545, 0.8483), and the two contact points: (xm , ym ) = (1.0113, 0.1799) and (xM , yM ) = (5.7487, 1.7527). Fig. 5.12 (c) shows that heteroclinic orbits connecting E1 and E2 , and E2 and E3 . 5.3.3. Case (c): xm < x1 < x2 < xM < x3 This case is similar to the case (b), and so the detail statements are omitted.

3424

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Fig. 5.13. (a) The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for xm < x1 < x2 < x3 < xM . (b) A schematic diagram of the relaxation oscillation for xm < x1 < x2 < x3 < xM . (c) A relaxation oscillation for xm < x1 < x2 < x3 < xM generated by a numerical example.

5.3.4. Case (d): xm < x1 < x2 < x3 < xM For xm < x1 < x2 < x3 < xM the dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) are shown in Fig. 5.13 (a). By Fenichel’s theory we can obtain that the equilibria E1 (x1 , h−1 x1 ) and E3 (x3 , h−1 x3 ) are both unstable nodes of system (2.2) and E2 (x2 , h−1 x2 ) is a saddle of system (2.2). On limit cycles we have the next conclusion. Theorem 5.10. Let U be a small tubular neighborhood of the common slow-fast cycle 0 . Then for each fixed ε > 0 sufficiently small, there exists a unique relaxation oscillation 3ε ⊂ U , which is the only limit cycle of system (2.2) and converges to 0 in the Hausdorff distance as ε → 0. Proof. The two contact points are both jump points. The geometrical description of the proof is shown in Fig. 5.13 (b), and the details are omitted. 2 Next, we give a numerical example. Example 5.8. Set ε = 0.001, K = 2.9153, b = −1.8135 and h = 9. Fig. 5.13 (c) shows the existence of a relaxation oscillation for system (2.2). 5.3.5. Case (e): x1 = xm < x2 < x3 < xM For xm = x1 < x2 < x3 < xM the dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) are shown in Fig. 5.14 (a). By Fenichel’s theory we can obtain that E2 (x2 , h−1 x2 ) is a saddle and E3 (x3 , h−1 x3 ) is an unstable node of system (2.2). Next result shows the existence of homoclinic orbits “small” or “big”, where the “small” homoclinic orbit is the one passing only one contact point; whereas the “big” homoclinic orbit is the one passing the two contact points. Theorem 5.11. Assume (K, b, h) ∈ S4 , and 0 < x1 = xm < x2 < x3 < xM < K. Then for 0 < ε  1, the following statements hold. √ (a) There exists a continuous function hc ( ε) having the expression (4.8) such that system (2.2) has an orbit √ homoclinic to the saddle E2 (x2 , y2 ), which√is either small or big, if and only if h = hc ( ε). Depending on the different choice of hc ( ε), both of the small and big ones could appear but not simultaneously.

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3425

Fig. 5.14. (a) The dynamics of fast subsystem (2.3) and of slow subsystem (2.5) for xm = x1 < x2 < x3 < xM . (b) A schematic diagram for the proof of Theorem 5.11. (c) A “small” homoclinic orbit. (d) A “big” homoclinic orbit.

√ (b) When varying h from hc ( ε) sufficiently small, the homoclinic orbit bifurcates a unique canard cycle without head when it is small, or a unique canard cycle with √ head if it is big. Furthermore, the unique canard cycle is unstable and exists for 0 < hc ( ε) − h  1. Proof. The existence of unique homoclinic orbit follows that of Theorem 5.6. Here the homoclinic orbit could be small or big, as shown in Fig. 5.14 (c) or (d). Next, we prove that a unique canard cycle without head bifurcates from the small homoclinic orbit. For j = l, m, let 2j = {(x, ρ 2 + ym ) : x ∈ Ij } be the horizontal sections for sufficiently j

small ρ > 0 and suitably small intervals Ij such that 2j centered at its intersection with C0 . Let 2 be the transition map from 2l to 2r following the flow of system (2.2), and set qj,ε = √ j

2j ∩ Cε . By Theorem 3.2 in [27] there exists a smooth function hc ( ε) of the form (4.6) such √ that 2 (ql,ε ) = qm,ε if and only if h = hc ( ε). Since the saddle quantity of E2 is   σ0 = −x2−2 K −1 (2x23 + (b − K)x22 + K) + ε(x22 + bx2 + 1) > 0, The split function can be written as d := 2 (ql,ε − qm,ε ).

0 < ε  1.

3426

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Fig. 5.15. (a) The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for xm < x1 < x2 < x3 = xM . (b) The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for x1 = xm < x2 < xM < x3 . (c) A small homoclinic orbit. (d) The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for x1 < xm < x2 < x3 = xM .

  √ According to Theorem 3.1 in [25], one has ∂d ∂h h=hc ( ε) > 0. Thus, Andronov-Leontovich theorem (see p. 200 in [29]) implies that a √ unique limit cycle bifurcates from the homoclinic orbit under small variation of h from h = hc ( ε).√Furthermore, since the saddle quantity σ0 > 0, the limit cycle is unstable and exists for 0 < hc ( ε) − h  1. Moreover, this limit cycle is a canard one without head. Similarly when the homoclinic orbit is the big one, it will bifurcate a unique canard cycle with √ head for 0 < hc ( ε) − h  1. The details are omitted. 2 Note that ε = 0.001, K = 2.9153, b = −1.8135 and h = 8.7332 realizes the condition on locations of the equilibria in Theorem 5.11. 5.3.6. Case (f ): xm < x1 < x2 < x3 = xM The results are similar to Theorem 5.11 for xm = x1 < x2 < x3 < xM , and are omitted. 5.3.7. Case (g): x1 = xm < x2 < xM < x3 The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for 0 < x1 = xm < x2 < xM < x3 are shown in Fig. 5.15 (b). Now system (2.2) for 0 < ε  1 has E2 as a saddle and E3 as a stable node respectively, whereas E1 is a canard point. We have the following results.

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3427

Fig. 5.16. (a) The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) for x1 = xm < x2 < x3 = xM . (b) A small homoclinic orbit. (c) A big homoclinic orbit.

Theorem 5.12. Assume (K, b, h) ∈ S4 , and 0 < x1 = xm < x2 < xM < x3 < K. Then for 0 < ε  1, the following statements hold. √ (a) There exists a continuous function hc ( ε) having the expression (4.8) such √ that system (2.2) has an orbit homoclinic to√the saddle E2 (x2 , y2 ) if and only if h = hc ( ε). (b) When varying h from hc ( ε) sufficiently small, the homoclinic orbit bifurcates to a unique canard √ cycle without head. Furthermore, the unique canard cycle is unstable and exists for 0 < hc ( ε) − h  1. Proof. The proofs are similar to those of Theorem 5.6, and Theorem 5.11, and are omitted.

2

5.3.8. Case (h): x1 < xm < x2 < x3 = xM The results are similar to Theorem 5.11, and omitted. 5.3.9. Case (i): x1 = xm < x2 < x3 = xM The dynamics of the fast subsystem (2.3) and of the slow subsystem (2.5) are shown in Fig. 5.16 (a). Now E2 is a saddle for system (2.2), and E1 and E3 are both located at the contact points of the critical manifold. Theorem 5.13. Assume (K, b, h) ∈ S4 and x1 = xm < x2 < xM = x3 . Then for 0 < ε  1, the following statements hold. √ (a) There exists a continuous function hc ( √ ε) of form (4.8) such that system (2.2) has a small homoclinic orbit if and only if h = hc ( ε). Furthermore, a unique canard cycle√without head bifurcates from the homoclinic orbit when h varies slightly from h = hc ( ε). The √ unique canard cycle is unstable and exists for 0 < hc ( ε) − h  1. (b) Assume that E yM ) when h varies slightly √3 always coincides with the contact point (xM ,√ from h = hc ( ε). Then there exists a continuous function hc (√ ε) of form (4.8) such that system (2.2) has a big homoclinic orbit if and only if h = hc ( ε). A unique canard cycle √ with head bifurcates from the big homoclinic orbit when λ√ varies slightly from h = hc ( ε). The unique canard cycle is unstable and exists for 0 < hc ( ε) − h  1. Proof. (a) Similar to that of Theorem 5.6. A small homoclinic orbit is shown√in Fig. 5.16 (b). (b) Note that E3 is always a canard point when h varies slightly from h = hc ( ε). Next analysis is similar to that of Theorem 5.6, we skip the details. See Fig. 5.16 (c) for a geometrical descrip-

3428

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Fig. 6.1. Slow-fast cycles when system (6.1) (or system (2.2)) has one positive equilibrium (black dot).

Fig. 6.2. Slow-fast cycles when system (6.1) (or system (2.2)) has two positive equilibria (black dots) where one equilibrium is a saddle-node.

tion of the big homoclinic orbit. Here the big homoclinic orbit with (xM , yM ) a canard point is different from that in case x1 = xm < x2 < x3 < xM with (xM , yM ) a jump point. 2 6. Cyclicity of the slow-fast cycles In this section we turn to study the maximal number of canard cycles bifurcated from slow-fast cycles. The slow divergence integral (SDI) is the main tool, see, e.g., [11,12,5,6,13,7,9]. 6.1. The existence and types of slow-fast cycles for ε = 0 Recall that in the interior of the first quadrant system (2.2) is smoothly equivalent to dx 1 x 2 = 1− (x + bx + 1) − y ≡ f (x, y, η), dt x K dy 1  y = εy 2 1 − h (x 2 + bx + 1) ≡ εg(x, y, η), dt x x

(6.1)

with η = (b, K, h) and (b, K) ∈ U , defined in Lemma 2.2. As before, let F (x) be the function defining the critical manifold of the system, see (2.6), and let (xm , F (xm )) and (xM , F (xM )) with xm < xM be the two contact points of the critical manifold. We have clearly the next result. Lemma 6.1. Assume that (xj , F (xj )) is a canard point of system (6.1), j ∈ {m, M}. Then for any (b, K) ∈ U system (6.1) for ε = 0 could have a slow-fast cycle passing (xj , F (xj )). Next, according to the analysis on the number, types and relative locations of positive equilibria of the slow subsystem of system (6.1) presented in Section 5, we can provide a complete list of slow-fast cycles, which are illustrated in Figs. 6.1, 6.2 and 6.3 as a schematic view. f f Let s be the set of all slow-fast cycles shown in Figs. 6.1, 6.2 and 6.3. We can express s as the disjoint union of the next subsets: (a) Small cycles: • := {(I-1), (I-6), (II-1), (II-7), (III-1), (III-7), (III-12), (III-15), (III-18)}. (b) Slow-fast cycles without head: m U := {(I-2), (II-2), (II-3), (III-2), (III-3), (III-13), (III-14), (III-19), (III-20)}, M := {(I-7), (II-8), (II-9), (III-8), (III-9), (III-16), (III-17), (III-24), U (III-25)}.

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3429

Fig. 6.3. Slow-fast cycles when system (6.1) (or system (2.2)) has three positive equilibria (black dots) where equilibrium E2 is a saddle.

M (c) Slow-fast cycles with head: m S := {(I-4), (II-4), (II-5), (III-4), (III-5)}, S := {(I-9), (II-10), (II-11), (III-10), (III-11)}, S :={(III-21), (III-22), (III-23)}. (d) Transitory slow-fast cycles: TU :={(I-3), (I-8)}. (e) Common slow-fast cycles: C :={(I-5), (II-6), (III-6)}.

Let s0 = yM − ym . For s ∈ (0, s0 ), let xl (s) < xm (s) < xr (s) be the three distinct roots of F (x) = ym + s as shown in Fig. 6.4 (b). Set xl (s0 ) = xl , xl (0) = xm (0) = xm , xm (s0 ) = xr (s0 ) = xM and xr (0) = xr . Then F (xl ) = F (xM ) = yM and F (xr ) = F (xm ) = ym . Next we give explicit expressions of the different types of slow-fast cycles γ (s) for s ∈ [0, s0 ] mentioned above: (a) If γ (s) ∈ m U , s ∈ (0, s0 ), then γ (s) = {(x, F (x)) : x ∈ [xl (s), xm (s)]} ∪ {(x, ym + s) : x ∈ [xl (s), xm (s)]}.

(6.2)

(b) If γ (s) ∈ M U , s ∈ (0, s0 ), then γ (s) = {(x, F (x)) : x ∈ [xm (s), xr (s)]} ∪ {(x, ym + s) : x ∈ [xm (s), xr (s)]}.

(6.3)

(c) If γ (s) ∈ m S , s ∈ (0, s0 ), then γ (s) ={(x, F (x)) : x ∈ [xl , xm (s)]} ∪ {(x, ym + s) : x ∈ [xm (s), xr (s)]} ∪ {(x, F (x)) : x ∈ [xM , xr (s)]} ∪ {(x, yM ) : x ∈ [xl , xM ]}.

(6.4)

3430

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

(d) If γ (s) ∈ M S , s ∈ (0, s0 ), then γ (s) ={(x, F (x)) : x ∈ [xl (s), xm ]} ∪ {(x, ym ) : x ∈ [xm , xr ]} ∪ {(x, F (x)) : x ∈ [xm (s), xr ]} ∪ {(x, ym + s) : x ∈ [xl (s), xm (s)]}, s ∈ (0, s0 ). (6.5) (e) If γ (s1 , s2 ) ∈ S , s1 , s2 ∈ (0, s0 ) and s1 < s2 , then γ (s1 , s2 ) ={(x, F (x)) : x ∈ [xl (s2 ), xm (s1 )]} ∪ {(x, ym + s1 ) : x ∈ [xm (s1 ), xr (s1 )]} ∪ {(x, F (x)) : x ∈ [xm (s2 ), xr (s1 )]} ∪ {(x, ym + s2 ) : x ∈ [xl (s2 ), xm (s2 )]}, s1 , s2 ∈ (0, s0 ) and s1 < s2 . (f ) If γ (s) ∈ TU , then γ (s) = {(x, F (x)) : x ∈ [xl , xM ]} ∪ {(x, yM ) : x ∈ [xl , xM ]} for (I-3), and γ (s) = {(x, F (x)) : x ∈ [xm , xr ]} ∪ {(x, ym ) : x ∈ [xm , xr ]} for (I-8). (g) If γ (s) ∈ C , then γ (s) = {(x, F (x)) : x ∈ [xl , xm ]} ∪ {(x, ym ) : x ∈ [xm , xr ]} ∪ {(x, F (x)) : x ∈ [xM , xr ]} ∪ {(x, yM ) : x ∈ [xl , xM ]}. 6.2. The slow divergence integral m M M T By the symmetry in m U and in U , in S and in S , and (I-3) and (I-8) in U we only study the limit cycles bifurcated from the formers in each of pairs, except the ones which are different. Denote by Xε,η the vector field associated to system (6.1), with η = (K, b, h). Set η0 := (K, b, h0 ) with h0 = xymm , where (xm , ym ) is the local minimum of the critical manim fold. Then for η = η0 , and in case that there exists a slow-fast cycle in m U or in S , the contact point (xm , ym ) becomes a non-degenerate canard point and the other contact point (xM , yM ) is a non-degenerate jump point. Hence there exists a solution x0(t) of the slow subsystem (2.5) connecting C0l and C0m , which is defined on (−∞, T ), where T > 0 is such that x0 (0) = xm , and either T is finite and limt→T − x0 (t) = xM if the reduced system has no an equilibrium on C0m , or T = ∞ and limt→∞ x0 (t) is an equilibrium most close to (xm , ym ). Similarly, there is a solution xˆ0 (t) of the slow subsystem (2.5) defined on (−∞, T ) which corresponds to the slow flow on C0r and satisfies limt→T − xˆ0 (t) = xM . Let tl (s) ≤ 0 ≤ tm (s) be such that

F (x0 (tl (s))) = F (x0 (tm (s))) = ym + s, and let tr (s) be defined by F (xˆ0 (tr (s))) = ym + s. Next, we defined the so-called slow divergence integral, see for instance [6–9,12]. For the slow-fast cycle in m U ∪ {(I-3)}, define the function I (s, η0 ) as t m (s)

I (s, η0 ) = tl (s)

∂f (x0 (t), F (x0 (t)), η0 )dt, s ∈ [0, s0 ], ∂x

(6.6)

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3431

Fig. 6.4. Illustration for the slow divergence integral when system (2.2) has a unique positive equilibrium (xm , ym ). (a) A slow-fast cycle γ (s, η0 ) without head (blue curve). (b) A slow-fast cycle γ (s, η0 ) with head (blue curve).

where s0 = yM − ym is defined above. For the slow-fast cycle in m S , define the function I (s, η0 ) as t m (s)

I (s, η0 ) =

∂f (x0 (t), F (x0 (t)), η0 )dt ∂x

tl (s0 )

T +

(6.7) ∂f (xˆ0 (t), F (xˆ0 (t)), η0 )dt, s ∈ [0, s0 ]. ∂x

tr (s)

The above two integrals I (s, η0 ) are called slow divergence integrals (SDIs). Making the change of integration variable from t to x in (6.6) and (6.7) yields xm (s)

I (s, η0 ) =

∂f F  (x) (x, F (x), η0 ) dx, ∂x g(x, F (x), η0 )

s ∈ [0, s0 ],

(6.8)

xl (s)

and xm (s)

I (s, η0 ) = xl

∂f F  (x) (x, F (x), η0 ) dx ∂x g(x, F (x), η0 )

xM +

(6.9) F  (x)

∂f (x, F (x), η0 ) dx, s ∈ [0, s0 ]. ∂x g(x, F (x), η0 )

xr (s)

If the canard-type slow-fast cycle has no a head (see Fig. 6.4 (a)), for x ∈ [xl , xm ], as in [31] we define σ (x) by F (x) = F (σ (x)). Hence for x ∈ [xl , xm ), we have

3432

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

σ  (x) =

F  (x) < 0. F  (σ (x))

If the canard-type slow-fast cycle has a head (see Fig. 6.4 (b)), for each x ∈ [xl , xm ], we define σ1 (x) ∈ [xm , xM ] and σ2 (x) ∈ [xM , xr ] by F (x) = F (σj (x)), j = 1, 2. For x = xm , x = xM one has σ1 (x) =

F  (x) F  (x) < 0, σ2 (x) =  > 0. F (σ2 (x)) 1 (x))

F  (σ

Let h(x) =

∂f ∂x (x, F (x), η0 )

(6.10)

g(x, F (x), η0 )

and let x = F −1 (y) be the single-valued inverse function of y = F (x) for x ∈ [xl , xm ]. The following results hold. Lemma 6.2. For s ∈ [0, s0 ], the slow divergence integrals (6.8) and (6.9) can be written respectively as    (h(σ (x)) − h(x)) 

ym +s

I (s, η0 ) = ym

dy

(6.11)

x=F −1 (y)

and    (h(σ1 (x)) − h(x)) 

ym +s

I (s, η0 ) = ym

yM dy + x=F −1 (y)

ym +s

   (h(σ2 (x)) − h(x)) 

Proof. It follows from direct computations, see also [31].

dy. (6.12) x=F −1 (y)

2

Now we recall some results from [12,6,7,31,8,9], which will be used later on. Lemma 6.3 ([9, Theorem 1]). When I (s, η0 ) = 0 for the transitory slow-fast cycle γ (s, η0 ) as shown in Fig. 1.1 (b), system (2.2) has at most one periodic orbit Hausdorff close to γ (s, η0 ) for 0 < ε  1. When I (s, η0 ) = 0, there are at most two periodic orbits Hausdorff close to γ (s, η0 ). Definition 6.1. A slow-fast cycle γ (s, η0 ) is called non-degenerate, if it is not transitory and all of its contact points (fold points) are non-degenerate. Here a transitory slow-fast cycle is the one given in Fig. 1.1 (b).

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3433

Lemma 6.4 ([31, Lemma 2.2]). If a canard-type slow-fast cycle γ (s, η0 ) is non-degenerate as shown in Fig. 6.4 (a) or (b), then the following statements hold. (a) If I (s, η0 ) = 0, then Cycl(Xε,η , γ (s), (0, η0 )) ≤ 1. Moreover, if I (s, η0 ) < 0 (or > 0) then the perturbed canard cycle from γ (s, η0 ) is stable (or unstable). ∂I (b) If I (s, η0 ) = 0 and (s, η0 ) = 0, then Cycl(Xε,η , γ (s), (0, η0 )) ≤ 2. ∂s ∂I (c) If I (s, η0 ) = 0 and (s, η0 ) is a zero point of with multiplicity m, then Cycl(Xε,η , γ (s), ∂s (0, η0 )) ≤ 2 + m. We remark that Lemma 2.2 of [31] was obtained by combining a series of results of De Maesschalck, Dumortier and Roussarie [11,12,6]. Lemma 6.5 ([7, Theorem 2]). When γ (s, η0 ) is a common slow-fast cycle as shown in Fig. 1.1 (d), there is a unique hyperbolic periodic orbit Hausdorff close to γ (s, η0 ) for sufficiently small ε > 0. 6.3. Cyclicity f

m M Section 6.1 provides a classification of the slow-fast cycles, i.e., s = • ∪ m U ∪ U ∪ S ∪ m m T M M T M S ∪ S ∪ U ∪ C . Here we confine ourselves to the cases U , U , S , S , U and C for parameters in U and b ≥ −1. The cases −2 < b < −1 and S cannot be solved here. In fact, the slow-fast cycle in S contains two canard points, at the moment there are no results in treating them.

6.3.1. Cyclicity of the small slow-fast cycle Note that each small slow-fast cycle in • is a canard point. A small limit cycles born from it belongs Hopf bifurcation. The number of limit cycles created via Hopf bifurcation is determined by the multiplicity of a weak focus [34] or correspondingly the codimension of the Hopf bifurcation. It has been shown in [21] that system (1.1) can exhibit two limit cycles due to degenerate Hopf bifurcation. We remark that for ε > 0 sufficiently small, by the canard explosion described in Section 5.1.2 it is relatively easy to prove the coexistence of two limit cycles when system (2.2) has a unique positive equilibrium, where one is from Hopf bifurcation and another from relaxation oscillation. Fig. 6.5 illustrates the coexistence of the small and big cycles for b < −1. Next analysis shows that this coexistence could happen for b ≥ −1. Recall from Theorem 5.2 that system (2.2) exhibits subcritical singular Hopf bifurcation when Aˆ 1 > 0. Note that Krupa and Szmolyan had obtained the bifurcation diagram for the standard slow-fast normal form (4.4), which is shown in Figure 7 (b) in [27]. Then, for system (2.2) we can similarly draw the bifurcation diagram for fixed ε corresponding to the subcritical singular Hopf bifurcation and the canard explosion, as shown in Fig. 6.5 (a). In the next results, S1 and S2 are defined at the beginning of subsection 5.1, and hH (ε) and hc (ε) are defined in (4.7) and (4.8), respectively.

3434

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Fig. 6.5. (a) Sketch of the bifurcation diagram corresponding to the subcritical singular Hopf bifurcation for system (2.2). Amp denotes the amplitude of the limit cycles. hs and hr marking the beginning and the ending of canard explosion. (b) The coexistence of a canard cycle without head and a canard cycle with head when K = 10, b = −11/10 in system (2.2).

Proposition 6.6. Assume that (K, b, h) ∈ S1 ∪ S2 , and that the unique positive equilibrium of the slow subsystem of system (2.2) coincides with the contact point (xm , ym ). Let Aˆ 1 defined in (4.5) be positive. Then for 0 < ε  1 the following statements hold. (a) If h > hH (ε), system (2.2) has a unique stable relaxation oscillation surrounding an unstable equilibrium. (b) If h = hH (ε), system (2.2) exhibits a subcritical singular Hopf bifurcation. (c) If hc (ε) < h < hH (ε), system (2.2) has two limit cycles surrounding a stable equilibrium, where the inner cycle is unstable and the outer cycle is stable. More precisely, the two limit cycles either are the small Hopf cycle and the big relaxation oscillation or the canard cycle without head and the canard cycle with head. Here we do not consider the cases that there are two or three equilibria on the critical manifold. In these cases we need further study the codimension of the Hopf bifurcation. In what follows we denote σ (x), σ1 (x) and σ2 (x) in Section 6.2 by σ , σ1 and σ2 , respectively. 6.3.2. Cyclicity of the slow-fast cycle without head Here and next subsections we assume b ≥ −1. First we discuss the cyclicity of the slow-fast cycle γum (s) ∈ m U , which is defined in (6.2). For system (6.1), the function h(x) defined in (6.10) along the slow-fast cycle γum (s) has the expression

h(x) = −

Kx 3 (2x 3 + (b − K)x 2 + K) , x ∈ (xl (s), xm ). (K − x)(x 2 + bx + 1)2 (Kx 2 − h0 (K − x)(x 2 + bx + 1))

According to Lemma 6.2, we need to consider the following SDI    (h(σ (x)) − h(x)) 

ym +s

I (s, η0 ) = ym

dy, x=F −1 (y)

(b, K) ∈ U,

(6.13)

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3435

where σ (x) ∈ (xm , xm (s)) is defined by F (x) = F (σ (x)) for x ∈ (xl (s), xm ) as in Section 6.2. x = F −1 (y) is the single-valued inverse function of y = F (x) for x ∈ (xl (s), xm ). Some calculations show that h(σ (x)) − h(x) =

K N m (x, σ (x)), x ∈ (xl (s), xm ), Xum (σ (x))Xum (x) u

(6.14)

where   Xum (x) := (K − x)(x 2 + bx + 1)2 Kx 2 − h0 (K − x)(x 2 + bx + 1) ,   Xum (σ (x)) := (K − σ )(σ 2 + bσ + 1)2 Kσ 2 − h0 (K − σ )(σ 2 + bσ + 1) , Num (x, σ (x)) := x 3 (2x 3 + (b − K)x 2 + K)Xum (σ (x)) − σ 3 (2σ 3 + (b − K)σ 2 + K)Xum (x). For getting the properties of I (s, η0 ), we need to discuss the sign of formula (6.14). Considering the slow-fast cycle γum (s) defined in (6.2) (note that h0 = xymm at the moment), we have

< x1 (1 − Kx )(x 2 + bx + 1) for x ∈ (xl (s), xm ) and h10 x > x1 (1 − Kx )(x 2 + bx + 1) for x ∈ (xm , xm (s)). Therefore, Xum (σ (x)) > 0 and Xum (x) < 0 for x ∈ (xl (s), xm ) and (b, K) ∈ U . Now Num (x, σ (x)) can be written in 1 h0 x

Num (x, σ (x)) =x 3 (2x 3 + (b − K)x 2 + K)(K − σ )Yum (x) − σ 3 (2σ 3 + (b − K)σ 2 + K)(K − x)Yum (σ (x)),

(6.15)

where Yum (x) = ((b + 1)x(Kx 2 − h0 (K − x)(x 2 + bx + 1))2 + 1)2 and Yum (σ (x)) = ((b + 1)σ (Kσ 2 − h0 (K − σ )(σ 2 + bσ + 1))2 + 1)2 . Obviously, Yum (x) ≥ 1 and Yum (σ (x)) ≥ 1 for (b, K) ∈ U and b ≥ −1. Note from Fig. 2.1 (b) that 2x 3 + (b − K)x 2 + K > 0 for x ∈ [xl (s), xm ) ⊂ (0, xm ) and 2σ 3 + (b − K)σ 2 + K < 0 for σ ∈ (xm , xm (s)] ⊂ (xm , xM ). Therefore, Num (x, σ (x)) > 0 for (b, K) ∈ U and b ≥ −1. Summarizing the above analysis gives h(σ (x)) − h(x) < 0 for (b, K) ∈ U and b ≥ −1. It verifies that the slow divergence integral (6.13) is negative for (b, K) ∈ U and b ≥ −1. For the slow-fast cycle γuM (s) ∈ M U defined in (6.3), its cyclicity can be similarly considered. In this case, the SDI can be written in    (h(σ1 (x)) − h(σ2 (x))) 

ym +s

I (s, η0 ) = yM

dy,

(b, K) ∈ U,

(6.16)

x=F −1 (y)

where σ1 (x) ∈ (xm (s), xM ) ⊂ (xm , xM ) and σ2 (x) ∈ (xM , xr (s)) ⊂ (xM , xr ) are defined by F (x) = F (σ1 (x)) = F (σ2 (x)) for x ∈ (xl , xl (s)). x = F −1 (y) is the single-valued inverse function of y = F (x) for x ∈ (xl , xl (s)). The similar arguments as did for (6.13) show that the slow divergence integral defined in (6.16) is positive for (b, K) ∈ U and b ≥ −1. Applying Lemma 6.4 to the two slow divergence integrals defined in (6.13) and (6.16) together with the facts that the limit cycle bifurcated from the γum (s) (resp. γuM (s)) is always stable (resp. unstable) and that the two contact points are respectively canard and jump points, we get the next results.

3436

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Proposition 6.7. For system (2.2) with 0 < ε  1, if (b, K) ∈ U and b ≥ −1, then the cyclicity M of the slow-fast cycles without head contained in m U ∪ U is at most one. Furthermore, the M canard cycles bifurcated from the slow-fast cycles belonging to m U (resp. U ) are stable (resp. unstable). 6.3.3. Cyclicity of the slow-fast cycle with head M In this subsection, we consider the cyclicity of the slow-fast cycles in m S ∪ S with (b, K) ∈ U and b ≥ −1, which contain the two contact points (xm , ym ) and (xM , yM ). We first study γsm (s) ∈ m S defined in (6.4). According to Lemma 6.2, we need to consider the following SDI    (h(σ1 (x)) − h(x)) 

ym +s

I (s, η0 ) = ym

yM dy + ym +s

x=F −1 (y)

   (h(σ2 (x)) − h(x)) 

dy, x=F −1 (y)

where σ1 (x) ∈ (xm , xm (s)) ⊂ (xm , xM ) and σ2 (x) ∈ (xM , xr (s)) ⊂ (xM , xr ) are defined by F (x) = F (σ1 (x)) for x ∈ (xl (s), xm ) and F (x) = F (σ2 (x)) for x ∈ (xl , xl (s)). x = F −1 (y) is the single-valued inverse function of y = F (x) for x ∈ (xl , xm ). We now consider the sign of h(σ1 (x)) − h(x) for (b, K) ∈ U . From (6.14) and (6.15) one has h(σ1 (x)) − h(x) =

K Xsm (σ1 (x))Xsm (x)

Nsm (x, σ1 (x)),

where Xsm (x) :=(K − x)(x 2 + bx + 1)2 (Kx 2 − h0 (K − x)(x 2 + bx + 1)) < 0, Xsm (σ1 (x)) :=(K − σ1 )(σ12 + bσ1 + 1)2 (Kσ12 − h0 (K − σ1 )(σ12 + bσ1 + 1)) > 0, Nsm (x, σ1 (x)) :=x 3 (2x 3 + (b − K)x 2 + K)(K − σ1 )Ysm (x) − σ13 (2σ13 + (b − K)σ12 + K)(K − x)Ysm (σ1 (x)), with Ysm (x) = ((b + 1)x(Kx 2 − h0 (K − x)(x 2 + bx + 1))2 + 1)2 and Ysm (σ1 (x)) = ((b + 1)σ1 (Kσ12 − h0 (K − σ1 )(σ12 + bσ1 + 1))2 + 1)2 . Obviously, Ysm (x) ≥ 1 and Ysm (σ1 (x)) ≥ 1 for (b, K) ∈ U and b ≥ −1. We can see from Fig. 2.1 (b) that 2x 3 + (b − K)x 2 + K > 0 and 2σ13 + (b − K)σ12 + K < 0. Hence we have Nsm (x, σ1 (x)) > 0 for (b, K) ∈ U and b ≥ −1. Consequently h(σ1 (x)) − h(x) < 0 for (b, K) ∈ U and b ≥ −1. Similar to the analysis for (6.14), we have h(σ2 (x)) − h(x) =

K N m (x, σ2 (x)), Xsm (σ2 (x))Xsm (x) s

where   Xsm (x) :=(K − x)(x 2 + bx + 1)2 Kx 2 − h0 (K − x)(x 2 + bx + 1) < 0,   Xsm (σ2 (x)) :=(K − σ2 )(σ22 + bσ2 + 1)2 Kσ22 − h0 (K − σ2 )(σ22 + bσ2 + 1) > 0, Nsm (x, σ2 (x)) :=x 3 (2x 3 + (b − K)x 2 + K)Xsm (σ2 (x)) − σ23 (2σ23 + (b − K)σ22 + K)Xsm (x).

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3437

By Fig. 2.1 (b), one has 2x 3 + (b − K)x 2 + K > 0 and 2σ23 + (b − K)σ22 + K > 0 for (b, K) ∈ U . It follows that Nsm (x, σ2 (x)) > 0 for (b, K) ∈ U . Hence, h(σ2 (x)) − h(x) < 0 for (b, K) ∈ U . Putting the above analysis together, we have the conclusion:   ym +s   yM    I (s, η0 ) = (h(σ1 (x)) − h(x)) dy + (h(σ2 (x)) − h(x)) dy   ym

ym +s

x=F −1 (y)

x=F −1 (y)

is negative for (b, K) ∈ U and b ≥ −1. Next we consider γsM (s) ⊂ M S defined by (6.5), and its associated SDI    (h(σ2 (x)) − h(x)) 

ym +s

I (s, η0 ) = ym

yM dy + x=F −1 (y)

ym +s

   (h(σ2 (x)) − h(σ1 (x))) 

dy, x=F −1 (y)

(6.17) where σ1 (x) ∈ (xm (s), xM ) ⊂ (xm , xM ) and σ2 (x) ∈ (xM , xr ) are defined by F (x) = F (σ1 (x)) = F (σ2 (x)) for x ∈ (xl , xl (s)) and F (x) = F (σ2 (x)) for x ∈ (xl (s), xm ), respectively. x = F −1 (y) is the single-valued inverse function of y = F (x) for x ∈ (xl , xm ). Note that the integration limits for h(σ2 (x)) − h(x) should be taken from ym to ym + s and for h(σ2 (x)) − h(σ1 (x)) should be from ym + s to yM . Similar to the analysis for γsm (s), one gets that h(σ2 (x)) − h(x) < 0 for (b, K) ∈ U and h(σ2 (x)) − h(σ1 (x)) > 0 for (b, K) ∈ U and b ≥ −1. It follows that    (h(σ2 (x)) − h(x)) 

ym +s

ym

dy < 0, for (b, K) ∈ U

(6.18)

dy > 0, for (b, K) ∈ U and b ≥ −1.

(6.19)

x=F −1 (y)

and yM ym +s

   (h(σ2 (x)) − h(σ1 (x))) 

x=F −1 (y)

To discuss the sign of (6.17) for (b, K) ∈ U and b ≥ 1. By (6.18) and (6.19) we have    (h(σ2 (x)) − h(x)) 

ym +s

lim I (s, η0 ) = lim

s→0

s→0 ym

yM + lim

s→0

ym +s

   (h(σ2 (x)) − h(σ1 (x))) 

  yM  = (h(σ2 (x)) − h(σ1 (x)))  ym

dy x=F −1 (y)

dy, x=F −1 (y)

dy > 0, x=F −1 (y)

(6.20)

3438

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

and    (h(σ2 (x)) − h(x)) 

ym +s

lim

s→(yM −ym )

I (s, η0 ) =

lim

s→(yM −ym )

ym

yM +

lim

s→(yM −ym )

ym +s

dy x=F −1 (y)

   (h(σ2 (x)) − h(σ1 (x))) 

   (h(σ2 (x)) − h(x)) 

dy,

(6.21)

x=F −1 (y)

yM = ym

dy < 0. x=F −1 (y)

Furthermore, for (b, K) ∈ U and b ≥ 1 we have dI (s, η0 ) = (h(σ2 (x)) − h(x))|x=F −1 (y) − (h(σ2 (x)) − h(σ1 (x)))|x=F −1 (y) < 0, y = ym + s. ds (6.22) From (6.20), (6.21) and (6.22), by continuity of I (s, η0 ) with respect to its variables, we get from the intermediate value theorem that for (b, K) ∈ U and b ≥ −1 there exists a unique s ∗ ∈ (0, yM − ym ) such that the function given in (6.17) satisfies I (s ∗ , η0 ) = 0, and I (s, η0 ) > 0 for s ∈ (0, s ∗ ), and I (s, η0 ) < 0 for s ∈ (s ∗ , yM − ym ). Based on the above discussion and Lemma 6.4, we obtain the following results. Proposition 6.8. For system (2.2) with 0 < ε  1, and (b, K) ∈ U and b ≥ −1, the following statements hold. (a) The cyclicity of the slow-fast cycles with head in m S is at most one. Furthermore, any canard cycle bifurcated from a slow-fast cycle belonging to m S is stable. , there exists a unique slow-fast cycle, denoted by γsM (s ∗ ), whose cyclicity is at (b) Inside M S M ∗ most two. All the other slow-fast cycles belonging to M S \ {γs (s )} has cyclicity at most M ∗ one. Furthermore, any canard cycle bifurcated from γs (s) ∈ M S with s ∈ (0, s )(resp. s ∈ ∗ (s , yM − ym )) is unstable (resp. stable). 6.3.4. Cyclicity of the transitory and common slow-fast cycles The case for cyclicity of the transitory slow-fast cycle is similar to that in subsection 6.3.2 via careful computations, we skip them here, because of the length of this paper. For cyclicity of the common slow-fast cycle in TC , as did in the proof of Theorem 5.6 we get exactly one limit cycle, which is a relaxation oscillation. Summarizing all the results obtained above we achieve cyclicity of the slow-fast cycles.

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3439

Theorem 6.9. For system (2.2) with 0 < ε  1, and (b, K) ∈ U and b ≥ −1, the following statements hold. M (a) The cyclicity of the slow-fast cycles without head in m U ∪ U and of the transitory slow-fast T cycles in U is at most one. Furthermore, the canard cycles bifurcated from the slow-fast M cycles belonging to m U ∪ {(I-3)} (resp. U ∪ {(I-8)}) are stable (resp. unstable). (b) The cyclicity of the slow-fast cycles with head in m S is at most one, and the canard cycles bifurcated from the slow-fast cycles belonging to m S are stable. (c) For the slow-fast cycles with head belonging to M S , there exists a unique slow-fast cycle, denoted by γsM (s ∗ ), whose cyclicity is at most two, and the cyclicity of slow-fast cycles M ∗ belonging to M S \ {γs (s )} is at most one. Furthermore, the canard cycles bifurcated from M M γs (s) ∈ S with s ∈ (0, s ∗ )(resp. s ∈ (s ∗ , yM − ym )) is unstable (resp. stable). (d) The cyclicity of the common slow-fast cycle in C is one. Furthermore, the unique periodic orbit bifurcated from the common slow-fast cycle is hyperbolic and stable. (e) If one of the two contact points is a canard point and the other is a jump point, system (2.2) has at most two limit cycles which can be bifurcated from slow-fast cycles.

Proof. Statements (a)–(d) were proved previously. Statement (e) follows from stability of the canard cycles bifurcated from the slow-fast cycles via the statements (a)–(d). This last fact forces that the maximum number, 2, of limit cycles (if exist) can only be nested, with the outside one being stable and the inner one unstable. Moreover, the two limit cycles could either both have M head and are bifurcated from M S , or have one with head from S and the other without head M from U , or have one from Hopf bifurcation at the maximum and the other being relaxation oscillation or canard cycle with head at the minimum. The reasons are the following: • By statement (a) and the assumption on contact points, system (2.2) does not have two canard cycles without head. • In case there is a canard cycle with the canard point at the minimum (with head or not), it must be stable. Then system (2.2) cannot have other canard cycles and relaxation oscillation, because of the stability of the canard cycle. • If there exists a canard cycle without head and with canard point at the maximum, it must be unstable. Then there has at most one another limit cycle, which is a stable canard cycle with head at the minimum and height s > s ∗ or a stable relaxation oscillation. • If there exists a canard cycle with head at the minimum, which is stable (resp. unstable), the head of the cycle must have height s > s ∗ (resp. s < s ∗ ). If system (2.2) has another limit cycle, by stability of all the cycles obtained before, this second limit cycle must be unstable (resp. stable) and have height s < s ∗ (resp. s > s ∗ ). We remark that if the slow-fast cycle γsm (s ∗ ) can bifurcate the maximum two limit cycles, then by (6.20), (6.21) and (6.22) the inner one must be unstable and the outer one stable. 2 7. Biological interpretation of the relaxation oscillation for a numerical example There were a huge number of studies devoted to the biological interpretation of prey-predator models, whereas there were a few on slow-fast motions. Here by a concrete numerical example we interpret biological phenomenon of the relaxation oscillation. Set ε = 0.00001, K = 225/32,

3440

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

Fig. 7.1. Relaxation oscillation for system (2.2) with given values of the parameters: (a) phase diagram; (b) the behavior of the prey and predator depending on time.

b = 29/32 and h = 14400/17549, system (2.2) has the phase diagram and the behavior of prey and predator plotted in Fig. 7.1, which show the existence of a relaxation oscillation. As shown in Fig. 7.1, the relaxation oscillation is a global stable limit cycle for system (2.2) consisting of alternating fast and slow segments. The biological interpretation of this interesting feature is as follows: the existence of the stable limit cycle indicates that the prey and predator can coexist in this concrete system. When the predator density is lower than the minimum value of the critical curve, we have a prey outbreak in a very short time. After the prey density arrives at some level which is enough to support the reproduction of the predators, the predator density begins to grow slowly for a long period until they reach the density which is higher than the maximum value of the critical curve. As a consequence, the prey density decreases dramatically. After a short time, the predator density declines slowly due to the less food. After some time, the prey density eventually recovers, starting a new cycle. We can conclude that when preys are numerous, their predators increase in numbers, and that when reducing the prey population, it will in turn causes predator number to decline. Acknowledgments The authors sincerely appreciate the anonymous referee for his/her nice comments and suggestions, which greatly improve this paper both in mathematics and presentations. The both authors are partially supported by NNSF of China grant number 11671254, and the second author is also partially supported by NNSF of China grant number 11871334. References [1] V.I. Arnol’d, Dynamical Systems V: Bifurcation Theory and Catastrophe Theory, Encyclopaedia Math. Sci., vol. 5, Springer, Berlin, 1994. [2] E. Benoit, équations différentielles: relation entrée-sortie, C. R. Acad. Sci., Sér. 1 Math. 293 (5) (1981) 293–296. [3] B. Braaksma, Singular Hopf bifurcation in systems with fast and slow variables, J. Nonlinear Sci. 8 (5) (1998) 457–490. [4] J.B. Collings, The effects of the functional response on the bifurcation behavior of a mite predator-prey interaction model, J. Math. Biol. 36 (2) (1997) 149–168. [5] P. De Maesschalck, F. Dumortier, Canard solutions at non-generic turning points, Trans. Am. Math. Soc. 358 (5) (2006) 2291–2334. [6] P. De Maesschalck, F. Dumortier, Canard cycles in the presence of slow dynamics with singularities, Proc. R. Soc. Edinb., Sect. A 138 (2) (2008) 265–299.

C. Wang, X. Zhang / J. Differential Equations 267 (2019) 3397–3441

3441

[7] P. De Maesschalck, F. Dumortier, R. Roussarie, Cyclicity of common slow-fast cycles, Indag. Math. (N.S.) 22 (3–4) (2011) 165–206. [8] P. De Maesschalck, F. Dumortier, R. Roussarie, Canard-cycle transition at a fast-fast passage through a jump point, C. R. Math. Acad. Sci. Paris 352 (1) (2014) 27–30. [9] P. De Maesschalck, F. Dumortier, R. Roussarie, Canard cycle transition at a slow-fast passage through a jump point, C. R. Math. Acad. Sci. Paris 352 (4) (2014) 317–320. [10] F. Dumortier, J. Llibre, J.C. Artés, Qualitative Theory of Planar Differential Systems, Universitext, Springer, Berlin, 2006. [11] F. Dumortier, R. Roussarie, Canard cycles and center manifolds, Mem. Am. Math. Soc. 121 (577) (1996), x+100. [12] F. Dumortier, R. Roussarie, Multiple canard cycles in generalized Liénard equations, J. Differ. Equ. 174 (1) (2001) 1–29. [13] F. Dumortier, R. Roussarie, Birth of canard cycles, Discrete Contin. Dyn. Syst., Ser. S 2 (4) (2009) 723–781. [14] N. Fenichel, Persistence and smoothness of invariant manifolds for flows, Indiana Univ. Math. J. 21 (1971/1972) 193–226. [15] N. Fenichel, Asymptotic stability with rate conditions, Indiana Univ. Math. J. 23 (1973/1974) 1109–1137. [16] N. Fenichel, Asymptotic stability with rate conditions. II, Indiana Univ. Math. J. 26 (1) (1977) 81–93. [17] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations, J. Differ. Equ. 31 (1) (1979) 53–98. [18] H.I. Freedman, R.M. Mathsen, Persistence in predator-prey systems with ratio-dependent predator influence, Bull. Math. Biol. 55 (4) (1993) 817–827. [19] G. Hek, Geometric singular perturbation theory in biological practice, J. Math. Biol. 60 (3) (2010) 347–386. [20] S.B. Hsu, T.W. Huang, Global stability for a class of predator-prey systems, SIAM J. Appl. Math. 55 (3) (1995) 763–783. [21] J. Huang, S. Ruan, J. Song, Bifurcations in a predator-prey system of Leslie type with generalized Holling type III functional response, J. Differ. Equ. 257 (6) (2014) 1721–1752. [22] C.K.R.T. Jones, Geometric singular perturbation theory, in: Dynamical Systems, in: Lecture Notes in Math., vol. 1609, Springer, Berlin, 1995, pp. 44–118. [23] C.K.R.T. Jones, S.-K. Tin, Generalized exchange lemmas and orbits heteroclinic to invariant manifolds, Discrete Contin. Dyn. Syst., Ser. S 2 (4) (2009) 967–1023. [24] T.J. Kaper, C.K.R.T. Jones, A primer on the exchange lemma for fast-slow systems, in: Multiple-Time-Scale Dynamical Systems, in: IMA Vol. Math. Appl., vol. 122, Springer, New York, 2001, pp. 65–87. [25] M. Krupa, P. Szmolyan, Extending geometric singular perturbation theory to nonhyperbolic points-fold and canard points in two dimensions, SIAM J. Math. Anal. 33 (2) (2001) 286–314. [26] M. Krupa, P. Szmolyan, Extending slow manifolds near transcritical and pitchfork singularities, Nonlinearity 14 (6) (2001) 1473–1491. [27] M. Krupa, P. Szmolyan, Relaxation oscillation and canard explosion, J. Differ. Equ. 174 (2) (2001) 312–368. [28] M. Krupa, M. Wechselberger, Local analysis near a folded saddle-node singularity, J. Differ. Equ. 248 (12) (2010) 2841–2888. [29] Y.A. Kuznetsov, Elements of Applied Bifurcation Theory, 3rd edition, Appl. Math. Sci., vol. 112, Springer, New York, 2004. [30] P.H. Leslie, Some further notes on the use of matrices in population mathematics, Biometrika 35 (1948) 213–245. [31] C. Li, H. Zhu, Canard cycles for predator-prey systems with Holling types of functional response, J. Differ. Equ. 254 (2) (2013) 879–910. [32] W. Liu, Geometric singular perturbations for multiple turning points: invariant manifolds and exchange lemmas, J. Dyn. Differ. Equ. 18 (3) (2006) 667–691. [33] D. Ludwig, D.D. Jones, C.S. Holling, Qualitative analysis of insect outbreak systems: the spruce budworm and forest, J. Anim. Ecol. 47 (1) (1978) 315–332. [34] L. Perko, Differential Equations and Dynamical Systems, 3rd edition, Texts in Appl. Math., vol. 7, Springer, New York, 2001. [35] S. Schecter, Exchange lemmas. I. Deng’s lemma, J. Differ. Equ. 245 (2) (2008) 392–410. [36] S. Schecter, Exchange lemmas. II. General exchange lemma, J. Differ. Equ. 245 (2) (2008) 411–441.