Averaging theory at any order for computing limit cycles of discontinuous piecewise differential systems with many zones

This work is devoted to study the existence of periodic solutions for a family of discontinuous differential systems $Z(x,y;\epsilon)$ with many zones. We show that for $\epsilon$ sufficiently small the averaged functions at any order control the existence of crossing limit cycles for systems in this family. We also provide some examples dealing with nonlinear centers.


Introduction and statement of the main results
In the qualitative theory of real planar differential system the determination of limit cycles, defined by Poincaré [21], has become one of the main problems. The second part of the 16th Hilbert problem deals with planar polynomial vector fields and proposes to find a uniform upper bound H(n) (called Hilbert's number) for the number of limit cycles that these vector fields can have depending only on the polynomial degree n. The averaging method has been used to provide lower bounds for the Hilbert number H(n) see, for instance, [13]. The interest on this topic extends to what we call discontinuous piecewise vector fields.
The increasing interest in the theory of nonsmooth vector fields has been mainly motivated by its strong relation with Physics, Engineering, Biology, Economy, and other branches of science. In fact, their associated differential systems are very useful to model phenomena presenting abrupt switches such as electronic relays, mechanical impact, and neuronal networks, see for instance [2,7,23]. The extension of the averaging theory to discontinuous piecewise vector field has been the central subject of investigation of the following works [11,12,14,17].
A piecewise vector field defined on an open bounded set U ⊂ R n is a function F : U → R n which is continuous except on a set Σ of measure 0, called the set of discontinuity of the vector field F . It is assumed that U \Σ is a finite collection of disjoint open sets U i , i = 1, 2, . . . , m, such that the restriction F i = F U i is continuous and extendable to the compact set U i . The local trajectory of F at a point p ∈ U i is given by the usual notion. However the local trajectory of F at a point p ∈ Σ needs to be given with some care. In [8], taking advantage of the theory of differential inclusion (see [1]), Filippov established some conventions for what would be a local trajectory at points of discontinuity where the set Σ is locally a codimension one embedded submanifold of R n . For a such point p ∈ Σ, we consider a sufficiently small neighborhood U p of p such that Σ splits U p \ Σ in two disjoint open sets U + p and U − p and denote F ± (p) = F U ± p (p).
In short, if the vectors F ± (p) point at the same direction then the local trajectory of F at p is given as the concatenation of the local trajectories of F ± at p. In this case we say that the trajectory crosses the set of discontinuity and that p is a crossing point. If the vectors F ± (p) point in opposite directions then the local trajectory of F at p slides on Σ. In this case we say that p is a sliding point. For more details on the Filippov conventions see [8,10].
In this paper we are interested in establishing conditions for the existence of crossing limit cycles for a class of planar discontinuous piecewise vector fields, that is limit cycles which only crosses the set of discontinuity Σ. It is worth to say that if Σ is locally described as h −1 (0), being h : U → R a smooth function and 0 a regular value, then ∇h(p), F + (p) ∇h(p), F − (p) > 0 is the condition in order that p is a crossing point.
In the sequel we introduce a short review of the averaging theory for computing isolated periodic solutions of differential equations. Then we set the class of planar discontinuous piecewise differential equations that we are interested. After that the rest of the section is devoted to the statement of our main result.
1.1. Background on the averaging theory for smooth systems. Let D be an open bounded subset of R + and consider C k+1 functions F i : R × D → R for i = 1, 2, . . . , k, and R : R × D × (−ε 0 , ε 0 ) → R. We assume that all these functions are 2π-periodic in the first variable. Now consider the following differential equation and assume that the solution ϕ(θ, z) of the unperturbed system r (θ) = F 0 (θ, r), such that ϕ(0, ρ) = ρ, is 2π-periodic for every ρ ∈ D. Here the prime denotes the derivative in the variable θ.

AVERAGING THEORY AT ANY ORDER FOR PIECEWISE DIFFERENTIAL SYSTEMS 3
A central question in the study of system (1) is to understand which periodic orbits of the unperturbed system r (θ) = F 0 (θ, r) persists for |ε| = 0 sufficiently small. In others words to provide sufficient conditions for the persistence of isolated periodic solutions. The averaging theory is one of the best tools to track this problem. Summarizing, it consists in defining a collection of functions f i : D → R, for i = 1, 2, . . . , k, called averaged functions, such that their simple zeros provide the existence of isolated periodic solutions of the differential equation (1). In [15,16] it was proved that these averaged functions are where y i : R × D → R for i = 1, 2, . . . , k, are defined recurrently by the following integral equations (3) Here ∂ L G(φ, ρ) denotes the derivative order L of a function G with respect to the variable ρ, and S l is the set of all l-tuples of non-negative integers (b 1 , b 2 , . . . , b l ) satisfying b 1 + 2b 2 + · · · + lb l = l, and L = b 1 + b 2 + · · · + b l .

1.2.
A class of planar discontinuous piecewise smooth vector fields. When one consider the above problem in the world of discontinuous piecewise differential systems it is not always true that the higher averaged functions (2) allow to study the persistence of isolated periodic solutions. In [17,12] this problem was considered for general Filippov systems when F 0 (θ, r) ≡ 0 and it was proved that the averaged function of first order can provide information in the persistence of crossing isolated periodic solutions. Furthermore the authors have found conditions on those systems in order to assure that the averaged function of second order also provides information on the existence of crossing isolated periodic solutions. When F 0 (θ, r) ≡ 0 but satisfies the condition that the solution ϕ(θ, ρ) is 2π-periodic the authors in [14] have found conditions on those systems in order to assure that the averaged function of first order provides information on the existence of crossing isolated periodic solutions.
This work is devoted to study the existence of isolated periodic solutions for an εfamily of planar discontinuous piecewise differential system (ẋ,ẏ) T = Z(x, y; ε). Here the dot denotes derivative in the variable t. In short we shall provide sufficient conditions in order to show that for |ε| = 0 sufficiently small the averaged functions (2) at any order can be used for obtaining information on the existence of crossing limit cycles for systems of this family.
We start defining the family of smooth piecewise differential systems that we shall study. The construction that we shall perform in the sequel has been done in [12] for a particular class of systems. Let n > 1 be a positive integer, α n = 2π and α = (α 0 , α 1 , . . . , α n−1 ) ∈ T n is a n-tuple of angles such that 0 = α 0 < α 1 < α 2 < · · · < α n−1 < α n = 2π and let X (x, y; ε) = (X 1 , X 2 , . . . , X n ) be a n-tuple of smooth vector fields defined on an open bounded neighborhood U ⊂ R 2 of the origin and depending on a small parameter ε in the following way For j = 1, 2, . . . , n let L j be the intersection between the domain U with the ray starting at the origin and passing through the point (cos α j , sin α j ), and take Σ = n j=1 L j . We note that Σ splits the set U \Σ ⊂ R 2 in n disjoint open sectors. We denote the sector delimited by L j and L j+1 , in counterclockwise sense, by C j , for j = 1, 2, . . . , n. Now let Z X ,α : U → R 2 be a discontinuous piecewise vector field defined as Z X ,α (x, y; ε) = X j (x, y; ε) when (x, y) ∈ C j , and consider the following planar discontinuous piecewise differential system (5) (ẋ,ẏ) T = Z X ,α (x, y; ε).
The above notation means that at each sector C j we are considering the smooth differential system (6) (ẋ,ẏ) T = X j (x, y; ε).

Standard form and main result.
The averaging theory deals with periodic nonautonomous differential systems in the standard form (1). Therefore in order to use the averaging theory for studying system (5) it has to be written in the standard form. A possible approach for doing this is to consider the polar change of variables x = r cos θ and y = r sin θ. However the appropriate change of variables may depend on the initial system (5). In general, for each j = 1, 2, . . . , n, after a suitable change of variables system (6) reads functions depending on the vector fields X j i , and they are 2π-periodic in the first variable, being D an open bounded interval of R + and S 1 ≡ R/(2πZ). Furthermore system (5) becomes where the characteristic function χ A (θ) of an interval A is defined as X System (8) is now a nonautonomous periodic discontinuous piecewise differential system having its set of discontinuity formed by Σ Denote by ϕ(θ, ρ) the solution of the system r (θ) = F 0 (θ, r) such that ϕ(0, ρ) = ρ. From now on this last system will be called unperturbed system. We assume the following hypothesis: (H1) For each z ∈ D the solution ϕ(θ, ρ) is defined for every θ ∈ S 1 , it reaches Σ only at crossing points, and it is 2π-periodic.
In what follows we state our main result.
The assumption D ⊂ R is not restrictive. In fact, if one consider D as being an open subset of R n the conclusion of Theorem 1 still holds by assuming that the Jacobian matrix Jf l (ρ * ) is nonsingular, that is det(Jf l (ρ * )) = 0. In this case the derivative ∂ L G(φ, ρ) is a symmetric L-multilinear map which is applied to a "product" of L vectors of R n , denoted as L j=1 y j ∈ R nL (see [15]). For the particular class of systems (8) Theorem 1 generalizes the main results of [12,14,17], increasing the order of the averaging theory. It also generalizes the main results of [11,22] dealing now with nonvanishing unperturbed systems and allowing more zones of continuity. This paper is organized as follows. In section 2 we provide, explicitly, the formulae of the averaged functions (2) for nonsmooth systems in the standard form (8). In section Thus we have the next result.
Note that when F 0 = 0 the recurrence defined in (10) is actually an integral equation. Moreover in order to implement an algorithm to compute the averaged function, it may be easier to write each z j i in terms of the partial Bell polynomials, which are already implemented in algebraic manipulators as Mathematica and Maple. For each pair of nonnegative integers (p, q), the partial Bell polynomial is defined as where S p,q is the set of all (p − q + 1)-tuple of nonnegative integers (b 1 , b 2 , . . . , b p−q+1 ) satisfying b 1 + 2b 2 + · · · + (p − q + 1)b p−q+1 = p, and b 1 + b 2 + · · · + b p−q+1 = q. In the next proposition, following [20], we solve the integral equation (10) to provide the explicit recurrence formula for z j i in terms of the Bell polynomials.
Solving the above linear differential equation we get Now for i = 2, . . . , k and j = 1 the recurrence (10) can be written in terms of the partial Bell polynomials as (for more details, see [20]) We note that the function z 1 i appears in the right hand side of (15) only if l = i and m = 1. In this case B i,1 (z 1 1 , z 1 2 , . . . , z 1 i ) = z 1 i for every i ≥ 1. So we can rewriting (15) as the following integral equation which is equivalent to the following Cauchy problem: Solving the above linear differential equation we obtain the expressions of z 1 i (θ, ρ), for i = 2, . . . , k, given in the statement of the proposition.

Proof of the main result
In this section we shall present the proof of Theorem 1. This proof is based on a preliminary result (see Lemma 4) which expands the solutions of the discontinuous differential equation (8) in powers of ε.
In the sequel we shall expand the right hand side of the above equality in Taylor series in ε around ε = 0. To do that we first recall the Faá di Bruno's Formula about the l-th derivative of a composite function. Let g and h be sufficiently smooth functions then where S l is the set of all l-tuples of non-negative integers (b 1 , b 2 , · · · , b l ) satisfying b 1 + 2b 2 + · · · + lb l = l, and L = b 1 + b 2 + · · · + b l . So expanding F j i (φ, r j (φ, ρ, ε)) in Taylor series in ε around ε = 0 we obtain From the Faá di Bruno's Formula we compute Substituting (19) in (18) we have for i = 0, 1, ..., k − 1. Moreover for i = k we have that . Substituting (20) and (21) in (17) we get (22) After some transformations of the indexes i and l we obtain Therefore from (22) and (23) we have where for i = 0, . . . , k and j = 1, 2, . . . , n we are taking (25) Note that for i = 1, . . . , k and j = 2, . . . , n the following recurrence holds with the initial condition Putting (26) and (27) together we obtain for i = 1, 2, . . . , k and j = 1, 2, . . . , n.

Examples
In this section we present three applications of our main result (Theorem 1). In the first two examples (subsections 4.1 and 4.2) we use the averaged functions (11) up to order 7 to provide lower bounds for the maximum number of limit cycles admitted by some piecewise linear systems with four zones. The first system is a piecewise linear perturbation of the linear center (ẋ,ẏ) = (−y, x), and the second one is a piecewise linear perturbation of a discontinuous piecewise constant center. As usual, the expressions of the higher order averaged functions are extensive (see [11,15]), so we shall omit them here. We emphasize that our goal in these first two examples, by taking particular classes of perturbations, is to illustrate the using of the higher order averaged functions.
In the third example we study the quadratic isochronous center (ẋ,ẏ) = (−y +x 2 , x+ xy) perturbed inside a particular family of piecewise quadratic system with n zones. Using the first order averaged function (11) we provide lower bounds, depending on n, for the maximum number of limit cycles admitted by this system. We emphasize that our goal in this last example, again by taking a particular class of perturbation, is to illustrate the using of Theorem 1 to study discontinuous piecewise nonlinear system with many zones.
The next proposition, proved in [6], is needed to deal with our examples.  such that h j | I has constant sign, it is possible to get an h given by (31), such that it has at least n − 1 simple zeroes in I.

4.1.
Nonsmooth perturbation of the linear center. The bifurcation of limit cycles from smooth and nonsmooth perturbations of the linear center (ẋ,ẏ) = (−y, x) is a fairly studied problem in the literature, see for instance [3,4,9,18,19]. Here we apply our main result (Theorem 1) to study these limit cycles when the linear center is perturbed inside a particular of piecewise linear system with 4 zones. Following the notation introduced in subsection 1.2 we take (32) X j 0 (x, y) = − y, x , for j = 1, . . . , n, and X j i (x, y) = a ij x + b j , 0 , for j = 1, . . . , n, and i = 1, . . . , k.
First of all, in order to apply our main result (Theorem 1) to study the limit cycles of (ẋ,ẏ) T = Z X ,α (x, y; ε), we shall write it into the standard form (8). To do that we consider the polar coordinates x = r cos θ, y = r sin θ. So the set of discontinuity becomes Σ = {θ = 0} ∪ {θ = α 1 } ∪ {θ = α 2 } ∪ {θ = α 3 } and in each sector C j (see (6)), j = 1, 2, 3, 4, the differential system (ẋ,ẏ) T = Z X ,α (x, y; ε) readṡ Note thatθ(t) = 0 for |ε| sufficiently small, thus we can take θ as the new independent time variable by doing r (θ) =ṙ(t)/θ(t). Then where F j i is the coefficient of ε i in the Taylor series in ε ofṙ(t)/θ(t) around ε = 0. From here we shall use the averaged functions (11) up to order 7 to study the isolated periodic solutions of the piecewise differential equation defined by (33) or, equivalently, the limit cycles of the piecewise differential system (ẋ,ẏ) T = Z X ,α (x, y; ε) defined by (32). As we have said before, due to the complexity of the expressions of the higher order averaged functions we shall not provided them explicitly. So we first describe the methodology to obtain lower bounds for the number of their zeros, and consequently for the number of limit cycles of (32).
Assume that one have computed the list of averaged functions f i , i = 1, . . . , k, and that they are polynomials. The first step is to established a lower bound for the number of zeros that f 1 can have. To do that, one can build a vector M 1 where each entry s of M 1 is given by the coefficient of r s of the function f 1 . Clearly M 1 is a function on the parameter variable v 1 = {a 1j : j = 1, . . . , 4} ∪ {b 1j : j = 1, . . . , 4}. So taking the derivative D v 1 M 1 , a lower bound for the number of zeros of f i will be given by the rank of the matrix D v 1 M 1 decreased by 1. For instance, in our first example system (33), the averaged function f 1 reads r(a 11 + a 12 + a 13 + a 14 ) From here we shall use the averaged functions (11) up to order 7 to study the isolated periodic solutions of the piecewise differential equation defined by (35) or, equivalently, the limit cycles of the piecewise differential system (34). Following the same methodology described in subsection 4.1, we provide a table showing the lower bound N (l), l = 1, . . . , 7, for the maximum number of limit cycles of (34) obtained by studying the averaged function of order l. In this section we consider the quadratic isochronous center (ẋ,ẏ) = (−y + x 2 , x + xy) perturbed inside a class of piecewise quadratic system with n zones. Following the notation introduced in subsection 1.2 we take X 1 0 (x, y) = − y + x 2 , x + xy , for j = 1, . . . , n, and where a j , b j and c j are real numbers for all j ∈ {1, 2, . . . , n}. We consider the discontinuous piecewise differential system (ẋ,ẏ) T = Z X ,α (x, y; ε) (see (5)) where X = X 1 , . . . , X n ) (see (4)) and α = (α j ) n−1 j=0 = (2jπ/n) n−1 j=0 . As before, in order to apply our main result (Theorem 1) to study the limit cycles of (ẋ,ẏ) T = Z X ,α (x, y; ε), we shall write it into the standard form (8). To do that we consider a first change of coordinates x = −u/(v − 1), y = −v/(v − 1) (see [5]). Note that this change keeps fixed all straight lines passing through the origin and therefore does not change the set of discontinuity. In each sector C j (see (6)), j = 1, 2, 3, 4, the differential system (ẋ,ẏ) T = Z X ,α (x, y; ε) reads Now, as a second change of variables, we consider the polar coordinates u = r cos θ and v = r sin θ. Taking θ as the new independent time variable by doing r (θ) =ṙ(t)/θ(t), system (36) becomes where F j (θ, r) = cos θ c j + r − c j sin θ + cos θ b j + a j r cos θ 1 − r sin θ .
It is easy to see that this combination is linearly independent. Regarding the functions h j 's we have the following properties (1) Let j ∈ {2, 3, . . . , n}. Then h j (r) ≡ 0 if and only if n is even and j = 1 + n/2.
(2) Let j 1 , j 2 ∈ {2, 3, . . . , n}. Then h j 1 (r) ≡ h j 2 (r) if and only if n is even and (j 1 + j 2 − 2) ∈ {n/2, 3n/2}. From the above properties we first conclude that if n is odd then the function f 1 is a linearly independent combination of n + 1 linearly independent functions. From Proposition 5 we can find parameters such that f 1 has n simple zeros.