SS symmetry Article Integrability and Limit Cycles via First Integrals Jaume Llibre Departament de Matemàtiques, Universitat Autònoma de Barcelona, Bellaterra, 08193 Barcelona, Spain; jllibre@mat.uab.cat Abstract: In many problems appearing in applied mathematics in the nonlinear ordinary differential systems, as in physics, chemist, economics, etc., if we have a differential system on a manifold of dimension, two of them having a first integral, then its phase portrait is completely determined. While the existence of first integrals for differential systems on manifolds of a dimension higher than two allows to reduce the dimension of the space in as many dimensions as independent first integrals we have. Hence, to know first integrals is important, but the following question appears: Given a differential system, how to know if it has a first integral? The symmetries of many differential systems force the existence of first integrals. This paper has two main objectives. First, we study how to compute first integrals for polynomial differential systems using the so-called Darboux theory of integrability. Furthermore, second, we show how to use the existence of first integrals for finding limit cycles in piecewise differential systems. Keywords: limit cycles; Darboux theory of integrability; first integrals  1. Introduction to the Darboux Theory of Integrability  For a differential system on a two dimensional manifold, its phase portrait is deter- Citation: Llibre, J. Integrability and mined by the existence of a first integral. The Hamiltonian differential systems are the Limit Cycles via First Integrals. easiest differential systems having a first integral. Symmetry 2021, 13, 1736. https:// A differential system of the form: doi.org/10.3390/sym13091736 ∂H ∂H Academic Editors: Danny Arrigo and ẋ = , ẏ = − , ∂y ∂x Dumitru Baleanu where H : R2 → R is a C2 function, is a Hamiltonian differential system or a simple Hamilto- Received: 10 August 2021 nian system in R2 . Accepted: 16 September 2021 The integrable planar differential systems different from the Hamiltonian ones, in Published: 18 September 2021 general, are not easy to find. First, we stated the basic results of the Darbouxian theory of integrability for finding first integrals for planar polynomial differential systems. The Publisher’s Note: MDPI stays neutral Darbouxian theory of integrability connects the integrability of polynomial differential with regard to jurisdictional claims in systems with the invariant algebraic curves that those systems have. published maps and institutional affil- iations. 1.1. Polynomial Differential Systems Let P and Q be real polynomials in the real variables x and y. Then, a differential system: Copyright: © 2021 by the author. dx dy = ẋ = P( x, y), = ẏ = Q( x, y), (1) Licensee MDPI, Basel, Switzerland. dt dt This article is an open access article distributed under the terms and is a two dimensional planar polynomial differential system or simply a polynomial system. As conditions of the Creative Commons usual, m = max{deg P, deg Q} is the degree of the polynomial system. In this paper, we Attribution (CC BY) license (https:// supposed that the polynomials P and Q were coprime in the ring of real polynomials in creativecommons.org/licenses/by/ the variables x and y. 4.0/). Symmetry 2021, 13, 1736. https://doi.org/10.3390/sym13091736 https://www.mdpi.com/journal/symmetry Symmetry 2021, 13, 1736 2 of 21 The Darboux theory of integrability illustrates for a polynomial differential system the relationships between the existence of exact algebraic solutions (an algebraic phenomenon) and the integrability (a topological phenomenon). 1.2. First Integrals As usual, the vector field: ∂ ∂ X=P +Q ∂x ∂y is associated to the differential system (1). Let U be an open subset of R2 . We stated that the polynomial differential system (1) was integrable in U if there was a non-constant analytic function H : U → R that was constant on all orbits ( x (t), y(t)) of the system (1) contained in U. Such a function is called a first integral. It was clear that H is a first integral of system (1) in U if, and only if, the following equality holds: dH = Hx ẋ + Hy ẏ = Hx P + Hy Q = XH = 0 in U. dt Of course, the curves H ( x, y) = constant in U were formed by the orbits of the differential system (1). Example 1. Consider the Hamiltonian system: ∂H ∂H ẋ = = P, ẏ = − = Q. ∂y ∂x Then, the function H is called the Hamiltonian of this system, and H is a first integral of it, because: ∂H ∂H XH = P +Q = Hy Hx − Hx Hy = 0. ∂x ∂y 1.3. Integrating Factors Again, U denotes an open subset of R2 . An analytic function R : U → R non- identically zero is an integrating factor of the polynomial differential system (1) in U if one of the next three equivalent conditions is satisfied in U: ∂( RP) ∂( RQ) XR = − R div( P, Q), div( RP, RQ) = 0, + = 0. ∂x ∂y As usual, we defined the divergence of the vector field X = ( P, Q) by: ∂P ∂Q div( X ) = + = div( P, Q). ∂x ∂y Doing the change of time dt = rds, the differential system (1) becomes: ∂H ∂H ẋ = RP = , ẏ = RQ = − . (2) ∂y ∂x Note that, since system (2) has a divergence of zero, it is Hamiltonian. To each integrating factor R, there is a first integral H associated, given by: Z H ( x, y) = R( x, y) P( x, y) dy + h( x ), (3) and the function h must satisfy ∂H = − RQ. ∂x Symmetry 2021, 13, 1736 3 of 21 Example 2. The quadratic polynomial differential system: ẋ = −y − b( x2 + y2 ) = P, ẏ = x = Q, (4) has the integrating factor R = 1/( x2 + y2 ). Indeed, ∂( RP) ∂( RQ) ∂(−y − b( x2 + y2 ))/( x2 + y2 ) ∂x/( x2 + y2 ) + = + = 0. ∂x ∂y ∂x ∂y In order to find a first integral of system (4), we follow (3), i.e.,: − y − b ( x 2 + y2 ) 1 Z Z H= RPdy + h( x ) = dy + h( x ) = −by − log( x2 + y2 ) + h( x ). x 2 + y2 2 Then, ∂H + RQ = h0 ( x ) = 0. ∂x Therefore, h( x ) is a constant and we can omit it from the first integral H; hence: 1 H = −by − log( x2 + y2 ), 2 or, as the first integral, we can take: e−by F = eH = p . x 2 + y2 The result of the next proposition is well known; for instance, see [1]. Since its proof is very short, we provided it. Proposition 1. Let R1 and R2 be two integrating factors of the polynomial differential system (1) in the open subset U of R2 ; then, the function R1 /R2 is a first integral in the open set U \ { R2 = 0}, if R1 /R2 is non-constant. Proof. We have XRi = − Ri div( P, Q) because Ri is an integrating factor for i = 1, 2. Therefore, from: ( XR1 ) R2 − R1 ( XR2 )   R1 X = = 0, R2 R22 the proposition follows. 1.4. Invariant Algebraic Curves A non-constant polynomial f ∈ R[ x, y] defines the invariant algebraic curve of the polynomial differential system (1) if there is a polynomial K ∈ R[ x, y], such that: ∂f ∂f Xf = P +Q = Kf. (5) ∂x ∂y We stated that the invariant algebraic curve f = 0 has, as a cofactor, the polynomial K. Of course, the cofactor has a degree of at most m − 1, because m is the degree of the polynomial differential system. From (5), it follows that the gradient (∂ f /∂x, ∂ f /∂y) of the invariant algebraic curve f = 0 is orthogonal to the vector field X = ( P, Q) at the points of the curve f = 0. Therefore, at each point of the curve f = 0, the vector X is tangent to the curve f = 0; therefore, the curve f = 0 is formed by orbits of the vector field X. Therefore, since the curve f = 0 is invariant under the flow defined by X, this curve is called an “invariant algebraic curve”. Symmetry 2021, 13, 1736 4 of 21 If the polynomial f was irreducible in the ring R[ x, y], then we stated that the invariant algebraic curve f = 0 was irreducible. Example 3. We claim that f 1 = ay + b = 0 with cofactor K1 = ax and f 2 = x2 + y2 − 1 = 0 with cofactor K2 = −2x are two invariant algebraic curves of the quadratic polynomial differential system: ẋ = −y( ay + b) − ( x2 + y2 − 1) = P, ẏ = x ( ay + b) = Q, (6) a 6= 0. Indeed, ∂ay + b ∂ay + b X f1 = P +Q = x ( ay + b) a = K1 f 1 , ∂x ∂y and ∂ f2 ∂f X f2 = P + Q 2 = (−y( ay + b) − ( x2 + y2 − 1))2x + x ( ay + b)2y ∂x ∂y = −( x2 + y2 − 1)2x = K2 f 2 . 1.5. Exponential Factors Other objects playing a similar role as the invariant algebraic curves in order to obtain first integrals of a polynomial differential system (1) are the exponential factors. Let h, g ∈ R[ x, y] either be coprime polynomials in the ring R[ x, y], or h ≡ 1. Then, for the polynomial differential system (1), the function exp( g/h) is an exponential factor if there is a polynomial K ∈ R[ x, y] of a degree of at most m − 1, such that:   g  g X exp = K exp . (7) h h Then, for the exponential factor exp( g/h) we stated that K was its cofactor. Since the exponential factor cannot vanish, it does not define invariant curves of the polynomial system (1). Example 4. Consider the polynomial differential quadratic system: ẋ = x (y + a) = P, ẏ = y = Q. Such a system has the exponential factor ey . Indeed: ∂ey ∂ey Xey = P +Q = yey = Key , ∂x ∂y with cofactor K = y. 1.6. The Method of Darboux For polynomial differential systems (1), we summarized the Darboux theory of inte- grability in the next theorem. Theorem 1. Let f i = 0 be irreducible invariant algebraic curves with cofactors Ki for i = 1, . . . , p, and let exp( g j /h j ) be exponential factors with cofactors L j for j = 1, . . . , q for a polynomial differential system (1) of degree m. (i) The function: µ1  µq gq    λ λp g1 f1 1 . . . f p exp . . . exp (8) h1 hq is a first integral of system (1) if, and only if, there are λi , µ j ∈ R, not all zero, such that: p q ∑ λi Ki + ∑ µ j L j = 0. i =1 j =1 Symmetry 2021, 13, 1736 5 of 21 p (ii) When p + q ≥ m(m + 1)/2 + 1, there are λi , µ j ∈ R, not all zero, satisfying ∑ λi Ki + i =1 q ∑ µ j L j = 0. j =1 (iii) If p + q ≥ m(m + 1)/2 + 2, all orbits of the differential system (1) are contained in invariant algebraic curves, because the system has a rational first integral. (iv) Function (8) is an integrating factor of system (1) if, and only if, there are λi , µ j ∈ R, not all p q zero, such that ∑ λi Ki + ∑ µ j L j = −div( P, Q). i =1 j =1 (v) We defined Fj = exp( g j /h j ). Since f i = 0 are invariant algebraic curves with cofactors, and Ki and Fj are exponential factors with cofactors L j , we have X f i = Ki f i , and XFj = L j Fj . Therefore, the statement (i) follows from the equality:  λp µ µq  λ X f 1 1 . . . f p F1 1 . . . Fq =  p q !  X fi XFj ∑ i fi ∑ µ j Fj = λ λp µ µq f 1 1 . . . f p F1 1 . . . Fq λ + i =1 j =1 p q !   f 1 . . . f p F1 . . . Fq , ∑ λi Ki + ∑ µ j L j = 0. λ1 λ p µ1 µq i =1 j =1 Example 5. Assume that a 6= 0 in the quadratic system: ẋ = −y( ay + b) − ( x2 + y2 − 1), ẏ = x ( ay + b). (9) Then, this system has f 2 = x2 + y2 − 1 = 0 with cofactor K2 = −2x, and f 1 = ay + b = 0 with cofactor K1 = ax, as invariant algebraic curves. Due to the fact 2K1 + aK2 = 0, by Theorem 1 (i) We obtained that H = ( ay + b)2 ( x2 + y2 − 1) a is a first integral of system (9). (ii) Since the degree of the polynomial cofactors Ki and L j is at most m − 1, we obtained that Ki , L j ∈ Rm−1 [ x, y], the space of all polynomials of R[ x, y] of a degree of at most m − 1. We observed that the vector space Cm−1 [ x, y] over C has a dimension m(m + 1)/2. Since all the polynomials Ki and L j belong to the vector space Cm−1 [ x, y] of the dimension (m + 1)/2], and we have p + q polynomials Ki and L j with p + q > m(m + 1)/2, and we obtained that the p + q polynomials must be linearly dependent in Cm−1 [ x, y]. Therefore, there are p q λi , µ j ∈ C, not all zero, such that ∑ λi Ki + ∑ µ j L j = 0. Hence, statement (ii) was proved. i =1 j =1 Example 6. Consider the real quadratic system: ẋ = x ( ax + c), ẏ = y(2ax + by + c), (10) with abc 6= 0. Then, this system has the following five invariant straight lines: f 1 = x = 0, f 2 = ax + c = 0, f 3 = y = 0, f 4 = ax + by = 0, f 5 = ax + by + c = 0. Therefore, from λ λ λ λ Theorem 1(ii) we obtained that system (10) has the first integral H = f 1 1 f 2λ2 f 3 3 f 4 4 f 5 5 with 5 λi ∈ R, satisfying ∑ λi Ki = 0, where Ki is the cofactor of f i . It is easy to check that K1 = ax + c, i =1 5 K2 = ax, K3 = 2ax + by + c, K4 = ax + by + c, K5 = ax + by. Hence, a solution of ∑ λi Ki = 0 i =1 is λ1 = λ5 = −1, λ2 = λ4 = 1 and λ3 = 0. Therefore a first integral of system (10) is: ( ax + c)( ax + by) H= . x ( ax + by + c) Symmetry 2021, 13, 1736 6 of 21 (iii) Under the hypotheses of this statement, we applied statement (ii) to the two following subsets of p + q − 1 > 0 functions formed by the invariant algebraic curves and the exponential factors. In this way, we obtained two linear dependencies between their corresponding cofactors, with which, after some relabeling and linear algebra, we could obtain: M1 + α3 M3 + . . . + α p+q−1 M p+q = 0, M2 + β 3 M3 + . . . + β p+q−1 M p+q = 0, where Ml are the cofactors L j and Ki , and αl , β l ∈ R. Therefore, by statement (i), we obtained that the two first integrals: α α p+q β β p+q G1 G3 3 . . . G p+q−1 , G2 G3 3 . . . G p+q−1 , of the differential system (1), Gl is either an invariant algebraic curve or an exponential factor with cofactor Ml for l = 1, . . . , p + q. Therefore, taking logarithms of these two, these previous first integrals, we obtained that: H1 = log( G1 ) + α3 log( G3 ) + . . . + α p+q log( G p+q ), H2 = log( G2 ) + β 3 log( G3 ) + . . . + β p+q log( G p+q ), are also first integrals of the differential system (1) where they are defined. Each one of these first integrals has associated an integrating factor Ri satisfying: ∂Hi ∂Hi P = Ri , Q = − Ri . ∂y ∂x Therefore, we have that: R1 ∂H2 ∂H1 = / . R2 ∂x ∂x Due to the fact that the functions Gl are exponentials of a quotient of polynomials, or poly- nomials, we obtained that ∂Hi /∂x are rational functions for i = 1, 2. Therefore, we obtained that R1 /R2 is a rational function, and by Proposition 1 we know that R1 /R2 is a rational first integral. Hence, statement (iii) was proved. In Example 5, we observed that, since this system has five invariant algebraic curves, by Theorem 1(iii) it has a rational first integral. The one that we obtained. p q (iv) Since equality ∑ λi Ki + ∑ µ j L j = −div( P, Q) is equivalent to the equality: i =1 j =1  λp µ µq  λ X f 1 1 . . . f p F1 1 . . . Fq =  p q !  ∑ λ i Ki + ∑ µ j L j = λ λp µ µq f 1 1 . . . f p F1 1 . . . Fq i =1 j =1  λp µ µq  λ − f1 1 ... f p F1 1 . . . Fq div( P, Q). Therefore, statement (iv) was proved. Example 7. The algebraic curve f 1 = x2 + y2 is invariant for the differential system: ẋ = −y − b( x2 + y2 ) = P, ẏ = x = Q, (11) with cofactor K1 = −2bx. Since K1 = div( P, Q), we obtained from Theorem 1(iv) that f 1−1 is an integrating factor of this differential system. Using this integrating factor, we could compute the first integral H = exp(−2by)( x2 + y2 ) of system (11). Symmetry 2021, 13, 1736 7 of 21 The Darboux theory of integrability, here presented for polynomial differential systems in R2 , was extended to polynomial differential differential systems in Rn and Cn . For more information on the Darboux theory of integrability, see Chapter 8 of [2]. We must mention that there is another nice theory for finding first integrals of the ordinary differential equations using the Lie symmetries; for instance, see book [3]. 2. Limit Cycles in Piecewise Differential Systems via First Integrals The study of limit cycles is one of the most important objectives in the qualitative theory of the planar ordinary differential equations. We remark that, to obtain an upper bound for the maximum number of limit cycles for a given differential system in the plane R2 , in general, is a very difficult problem. The study of the discontinuous piecewise differential systems, more recently also called Filippov systems, has attracted the attention of mathematicians during these past decades due to their applications. These piecewise differential systems in the plane are formed by different differential systems defined in distinct regions separated by a curve. A pioneering work on this subject was due to Andronov, Vitt and Khaikin in the 1920s, and later on, Filippov, in 1988, provided the theoretical bases for these kinds of differential systems. Nowadays, a vast literature on these differential systems is available; for instance, see the books of [4–7] and the survey by [8]. As for the smooth differential systems, the study of the existence and location of limit cycles in the piecewise differential systems is also of great importance. The main tools for computing analytical limit cycles of differential systems are based on the averaging theory, the Melnikov integral, the Poincaré map, and the Poincaré map together with the Newton–Kantorovich Theorem or the Poincaré–Miranda theorem. To these tools, in the particular case of the piecewise differential systems, we had to add the use of the first integrals of the differential systems forming the piecewise differential systems for computing their limit cycles. To show how to use the first integrals for computing the limit cycles and the periodic orbits of the piecewise differential systems is the objective of this second part of this paper. Of course, this tool is restricted to the piecewise differential systems such that all their differential systems be integrable, in the sense that we know for each of them a further invariant st integral. 2.1. Discontinuous Piecewise Differential Systems A discontinuous piecewise differential system on R2 is a pair of Cr (with r ≥ 1) differential systems in R2 separated by a smooth curve Σ. The line of discontinuity Σ of the discontinuous piecewise differential system is given by Σ = h−1 (0), where h : R2 −→ R is a C1 function having 0 as a regular value. Observe that Σ is the boundary between the regions Σ+ = {( x, y) ∈ R2 | h( x, y) > 0} and Σ− = {( x, y) ∈ R2 | h( x, y) < 0}. Hence: ( X ( x, y), if h( x, y) ≥ 0, Z ( x, y) = (12) Y ( x, y), if h( x, y) ≤ 0, is the vector field corresponding to a piecewise differential system with a line of discontinuity Σ. When the vector fields X and Y coincide on the line Σ, we obtain a continuous piecewise differential system on R2 , that, in general, will not be smooth on Σ. The vector field (12) is usually denoted by Z = ( X, Y, Σ) or simply by Z = ( X, Y ), if the separation line Σ is known. In order to establish a definition for the trajectories of Z, we had to have a criterion for the transition of the trajectories between Σ+ and Σ− across the curve of discontinuity Σ. The contact between the curve of discontinuity Σ and the vector field X (or Y) is described by the directional derivative of h with respect to the vector field X, i.e.,: Xh( p) = h∇h( p), X ( p)i. Symmetry 2021, 13, 1736 8 of 21 Here, h., .i denotes the usual inner product of the plane R2 . Filippov in [9] stated the main results of the discontinuous piecewise differential systems. The curve of discontinuity Σ is divided into the three following sets: (a) Σc : { p ∈ Σ : Xh(x) · Yh(x) > 0}, the Crossing set. (b) Σe : { p ∈ Σ : Xh(x) > 0 and Yh(x) < 0}, the Escaping set. (c) Σs : { p ∈ Σ : Xh(x) < 0 and Yh(x) > 0}, the Sliding set. The points of Σ, where both vector fields X and Y simultaneously point outwards or inwards, define the escaping Σe or sliding Σs regions, while the interior in Σ of their complement defines the crossing region Σc (see Figure 1). The points of Σ, which are not in Σc ∪ Σe ∪ Σs , are the tangency points between X or Y and Σ. Figure 1. Crossing, sliding, and escaping regions, respectively. There are many papers studying the limit cycles of continuous and discontinuous piecewise differential systems in R2 ; for instance, see [8,10–23]. 2.2. Limit Cycles of a Piecewise Differential System Formed by a Linear Differential System and a Quadratic Polynomial Differential System Separated by the Straight Line x = 0 In what follows, we wanted to study the limit cycles of the following discontinuous piecewise differential system separated by the straight line x = 0. In the half-plane x ≥ 0, there is the linear differential system: ẋ = 2 + 2x − 2y, ẏ = 6 − 2y, (13) and in the half-plane x ≤ 0, there is the quadratic polynomial differential system: ẋ = −9x + 15y + 4x2 + 8xy − 28y2 , (14) ẏ = −6x + 9y − 4x2 + 32xy − 44y2 . Theorem 2. The discontinuous piecewise differential systems (13) and (14) have a unique limit cycle shown in Figure 2. Proof. Since we wanted to compute the limit cycles of this piecewise differential system using their first integrals, we had to find such first integrals. It is known that all the linear differential systems in Rn for n ≥ 2 are Darboux integrable (see [24]), so, in particular, the differential system (13) must have a first integral. Clearly, the differential system (13) is Hamiltonian with Hamiltonian: H1 ( x, y) = 6x − 2y − 2xy + y2 . Therefore, we obtained a first integral of system (13). System (14) has a center at the origin of coordinates, because the eigenvalues of the linear part of the system at the origin are ±3i, so the origin is either a weak focus or a center, but when computing their Lyapunov constants, we saw that all of them are zero, so it is a center; for more details, see Chapter 5 of [2]. Moreover, it is well known that all quadratic polynomial differential systems having a center are Darboux integrable; for more details see the proof of Theorem 8.15 of [2] or paper [25]. Now, in order to find a first integral for the differential system (14) we applied the Darboux theory of integrability. We started looking for their invariant algebraic curves of degrees one and two. Symmetry 2021, 13, 1736 9 of 21 Figure 2. The unique limit cycle that exists for the piecewise differential systems (13) and (14). The limit cycle is travels in counterclockwise sense. First, we looked for invariant straight lines f = f ( x, y) = ax + by + c = 0, and since the polynomial differential system has a degree of two, their cofactors must be polynomials of a degree of at most one, so they must be of the form K = k0 + k1 x + k2 y, and f and K must satisfy Equation (5). Passing the right hand side of this equation to the left, we obtained the polynomial: −ck0 + (−9a − 6b − ak0 − ck1 ) x + (15a + 9b − bk0 − ck2 )y + (4a − 4b − ak1 ) x2 + (8a + 32b − bk1 − ak2 ) xy + (−28a − 44b − bk2 )y2 = 0. Therefore, we had to solve the system: ck0 = 0, 9a + 6b + ak0 + ck1 = 0, 15a + 9b − bk0 − ck2 = 0, 4a − 4b − ak1 = 0, 8a + 32b − bk1 − ak2 = 0, 28a + 44b + bk2 = 0, in the unknowns a, b, c, k0 , k1 , and k2 for obtaining the possible invariant straight lines of the differential system (14). This system has a unique solution with ( a, b, c) = (0, 0, 0). Namely, b = − a, c = −3a/8, k0 = 0, k1 = 8, and k2 = −16. Therefore, we obtained: 1 f = a(8x − 8y − 3) and K = 8( x − 2y). 8 Therefore, without a loss of generality, we could assume that the invariant straight line is f 1 = 8x − 8y − 3 = 0 with the cofactor K1 = 8( x − 2y). Now, we looked for possible invariant algebraic curves of degree two. Therefore, we had to solve Equation (5) with f = a0 + a1 x + a2 y + a3 x2 + a4 xy + a5 y2 and K = k0 + k1 x + k2 y. Solving it, in a similar way for obtaining the invariant algebraic curves of degree one, we obtained only two invariant algebraic curves of degrees two for one, namely: 1 1 − a (8x − 8y − 3)2 = 0, and − a (−9 + 96x − 96y + 128x2 − 768xy + 896y2 ). 128 4 768 4 Symmetry 2021, 13, 1736 10 of 21 The first solution really is an invariant straight line. Therefore, essentially, there is a unique invariant algebraic curve of degree two, which we could take: f 2 = −9 + 96x − 96y + 128x2 − 768xy + 896y2 = 0 with cofactor K2 = 32( x − 2y). We knew that the equation λ1 K1 + λ2 K2 = 0 is satisfied with λ1 = −4 and λ2 = 1; therefore, by Theorem 2(ii) we obtained the first integral: −9 + 96x − 96y + 128x2 − 768xy + 896y2 H2 ( x, y) = . (8x − 8y − 3)4 Clearly, in the half-planes x ≥ 0 and x ≤ 0, the differential systems (13) and (14) have no limit cycles, because the polynomial H1 and rational first integral H2 prevent their existence, respectively. Therefore, the piecewise differential system formed with the systems (13) and (14) separated by the straight line x = 0 has limit cycles, which must cross line x = 0 in exactly two points, denoted by (0, y) and (0, Y ) with y < Y. These two points must be crossing points and satisfy the system: e1 = H1 (0, y) − H1 (0, Y ) = 0, e2 = H2 (0, y) − H2 (0, Y ) = 0. The unique solution of this system satisfying y < Y is: 1 q √ q √  (y, Y ) = 1616 − 2222(1013 − 9 1257), 1616 + 2222(1013 − 9 1257) . (15) 1616 This solution provides the limit cycle of Figure 2. 2.3. Limit Cycles of Piecewise Differential Systems Formed by Three Linear Centers These last years, since the piecewise linear differential systems, had many relevant applications to physical phenomena; the interest for studying them has increased strongly. As in the smooth differential systems also in the piecewise linear differential systems, the study of their limit cycles plays a main role. Almost all papers studying the limit cycles of the piecewise linear differential systems consider piecewise linear differential systems formed only by two pieces. In this subsection, we studied piecewise linear differential systems formed with three pieces. In [26], we studied the limit cycles in R2 of the discontinuous piecewise linear differ- ential systems separated by the line of discontinuity Σ = {( x, y) ∈ R2 : y = 0 or x = 0 and y ≥ 0}, and formed by three arbitrary linear centers. Such discontinuous piecewise linear differen- tial systems can exhibit, at most, three limit cycles, three being the maximum number of limit cycles that they can exhibit. In particular, it was proved that there are such piecewise linear differential systems with three limit cycles, each limit cycle having a unique point in each branch of the three branches of Σ \ {(0, 0)}. The three components of R2 \ Σ are Q1 = {( x, y) ∈ R2 : x > 0 and y > 0}, Q2 = {( x, y) ∈ R2 : x < 0 and y > 0}, and H = {( x, y) ∈ R2 : y < 0}. The objective of this section was to study the limit cycles of the discontinuous piece- wise linear differential system defined by: 2 379 237 ẋ = y+ , ẏ = −2x + , in Q1 , 1565 1565 313 √ 4 11566 4 4430533 ẋ = y+ , ẏ = −8x − , in Q2 , (16) 1565 10955 2191 √  2 4430533 − 1299 ẋ = 2y, ẏ = −8x − , in H. 2191 Symmetry 2021, 13, 1736 11 of 21 Theorem 3. The discontinuous piecewise differential system (16) has three limit cycles intersecting each branch of the three branches of Σ in one point. These limit cycles travel in clockwise sense, see Figure 3. Figure 3. The three limit cycles of the discontinuous piecewise linear differential system (16). These limit cycles travel in counterclockwise sense. Proof. Of course, we wanted to prove this objective using the first integrals of the three linear differential centers. It is well known that all the linear differential centers can be obtained upon conducting an affine transformation of the linear differential center ẋ = −y, ẏ = x, which has the first integral H = x2 + y2 . Therefore, all the linear differential centers are Hamiltonian systems, and, then, by computing the Hamiltonians of the linear centers in Q1 , Q2 and H, we obtained the first integrals: H1 ( x, y) = 1565x2 + y2 − 1185x + 379y, √ H2 ( x, y) = 21910x2 + 7y2 + 10 4430533x + 5783y, √ H3 ( x, y) = 8764x2 + 2191y2 + (2 4430533 − 2598) x, for each one of these three systems, respectively. These limit cycles must intersect each branch of Σ \ {(0, 0)} in one point. These points are ( x+ , 0) with x+ > 0, (0, y+ ) with y+ > 0, and ( x− , 0) with x− < 0. These three points are crossing points. Then, the first integrals H1 , H2 , and H3 satisfy the following three equations: H1 ( x+ , 0) − H1 (0, y+ ) = 0, H2 (0, y+ ) − H2 ( x− , 0) = 0, (17) H3 ( x− , 0) − H3 ( x+ , 0) = 0, or, equivalently: 1565x+2 − y2 − 1185x − 379y = 0, + + + 2 2 √ 21910x− − 7y+ + 10 4430533x− − 5783y+ = 0, (18)  √  ( x+ − x− ) 4382x+ + 4382x− + 4430533 − 1299 = 0. Symmetry 2021, 13, 1736 12 of 21 We looked for the solutions ( x+ , y+ , x− ) of these three equations satisfying x+ > 0, x− < 0 and y+ > 0, and these solutions are: √ ! 1 1 1 −3083 − 4430533 ( x+ , y+ , x− ) = 1, ,1 , 4382 √ ! 2 2 2 −7465 − 4430533 ( x+ , y+ , x− ) = 2, , 10 , 4382 √ ! 3 , y3 , x3 ) = 3, − 11847 − 4430533 ( x+ + − , 26 . 4382 The solution of system (16): ! ! √ r r 1 2 2 x1 ( t ) = sin 2 t 3130(14u + 5783) cos 2 t 21910 1565 1565 !! √ r 2 −10 4430533 sin 2 t , 1565   r ! r r ! 5783 2 1 22152665 2 y1 ( t ) = u+ cos 4 t − sin 4 t 14 1565 7 626 1565 5783 − , 14 satisfies the initial conditions x1 (0) = u and y1 (0) = 0. The solution of system (16):       237 2t 379 2t 237 x2 ( t ) = v − cos √ + √ sin √ + , 626 1565 2 1565 1565 626 r     1 5 2t 379 2t 379 y2 ( t ) = (237 − 626v) sin √ + cos √ − , 2 313 1565 2 1565 2 satisfies the initial conditions x3 (0) = v and y3 (0) = 0. The solution of system (16): 1  √   x3 ( t ) = 8764w cos(4t) + 4430533 − 1299 (cos(4t) − 1) , 8764  √  1299 − 8764w − 4430533 sin(4t) y3 ( t ) = , 4382 satisfies the initial conditions x3 (0) = w and y3 (0) = 0. Now, we considered the solution ( xk1 (t), y1k (t)) for k = 1, 2, 3 of the discontinuous piecewise linear differential system (16) given by the solution ( x+ 1 , y1 , x1 ) of system (18). + − 1 1 Then, the time that the solution ( x1 (t), y1 (t)) in Q1 needs to reach point (0, v) is t1 = 0.785398163397448. The time that the solution ( x21 (t), y12 (t)) in Q2 needs to reach point (w, 0) is t2 = 4.10363864680248. Finally, t3 = 1.11762450719575. is the time that the solution ( x31 (t), y13 (t)) in H needs to reach point (u, 0). Let ( xk2 (t), y2k (t)) for k = 1, 2, 3 be the solution of the discontinuous piecewise linear differential system (16) given by the solution ( x+ 2 , y2 , x2 ) of system (18). Then the time that + − 2 2 the solution ( x1 (t), y1 (t)) in Q1 needs to reach the point (0, v) is r1 = 0.785398163397448. The time that the solution ( x22 (t), y22 (t)) in Q2 needs to reach the point (w, 0) is r2 = 7.93799227264621. The time that solution ( x32 (t), y23 (t)) in H needs to reach point (u, 0) is r3 = 2.02943545009903. Symmetry 2021, 13, 1736 13 of 21 Let ( xk3 (t), y3k (t)) for k = 1, 2, 3 be the solution of the discontinuous piecewise linear differential system (16), given by solution ( x+ 3 , y3 , x3 ) = of system (18). Then, the time that + − 3 3 solution ( x1 (t), y1 (t)) in Q1 needs to reach point (0, v) is s1 = 0.785398163397448. The time that solution ( x23 (t), y32 (t)) in Q2 needs to reach point (w, 0) is s2 = 11.27688306691738. The time that solution ( x33 (t), y33 (t)) in H needs to reach point (u, 0) is s3 = 2.88219547492608. j j Drawing the three orbits ( xk (t), yk (t)) for j = 1, 2, 3 and for the times t ∈ [0, tk ], t ∈ [0, rk ], and t ∈ [0, sk ] for k = 1, 2, 3, respectively, we obtained the three limit cycles of Figure 3, which are travel in a clockwise sense. 2.4. Periodic Orbits of a Relay System in R3 Consider the discontinuous piecewise linear differential system: ẋ = y, ẏ = z, ż = −sign( x )y. (19) Then, the goal of this subsection was to analytically study the periodic orbits of this differential system using their first integrals. The sign function is given by:   −1 if x < 0, sign( x ) = 0 if x = 0, 1 if x > 0.  Therefore, x = 0 is the plane of discontinuity of system (19). System (19) was studied in [27]. Many discontinuous piecewise differential systems (19) appear in a natural way from the control theory. In fact, system (19) is a particular relay system of the ones studied in [28]. In general, a difficult problem is to analytically find the periodic solutions of a differ- ential system; usually, this problem cannot be solved. We wanted to compute the periodic orbits of continuous or discontinuous piecewise differential systems such that they are completely integrable in each piece. The objective of this subsection was to prove the next result. Theorem 4. The following statements hold: (a) Assume that the discontinuous piecewise linear differential system (19) has a periodic orbit γ such that intersects x = 0 in two points. Then, these two points are (0, y, z) and (0, −y, z) with z > 0 and y2 − z2 < 0. See Figure 4. (b) For every pair of points (0, y, z) and (0, −y, z) with z > 0 and y2 − z2 < 0, the discontinuous piecewise linear differential system (19) has a periodic orbit γ intersecting x = 0 in these two points. Proof. The discontinuous piecewise linear differential system (19) in x > 0 is given by: ẋ = y, ẏ = z, ż = −y, (20) and in x < 0 by: ẋ = y, ẏ = z, ż = y. (21) Observe that the points of the x–axis are equilibrium points. We could obtain two independent first integrals for the systems (20) and (21) because these differential systems are linear. Linear differential systems are completely Darboux integrable; for more details, see [24]. Therefore, for system’s (20) two independent first integrals are: H1 = x + z, and H2 = y2 + z2 . Symmetry 2021, 13, 1736 14 of 21 Hence, the trajectories of system (20) are contained in: γh1 h2 = { H1 = h1 } ∩ { H2 = h2 } ∩ { x > 0}, for all (h1 , h2 ) ∈ R2 . For system (21), the two independent first integrals are F1 = x − z, and F2 = y2 − z2 . Hence, the trajectories of system (21) are contained in: γ f1 f2 = { F1 = f 1 } ∩ { F2 = f 2 } ∩ { x < 0}, for all ( f 1 , f 2 ) ∈ R2 . The set γh1 h2 with h2 > 0 is composed in x > 0 by the intersection of the plane H1 = h1 with the cylinder H2 = h2 . Therefore, the set γh1 h2 is an arc without equilibria, and it is a trajectory of system (20). If h2 = 0, then γh1 h2 is an equilibrium point. Hence, γh1 h2 is a unique orbit of system (20). Figure 4. The periodic orbit of Theorem 4 intersecting x = 0 in the points (0, −1, 2) and (0, 1, 2). If f 2 6= 0, the set γ f1 f2 in x < 0 is one or two arcs of the intersection of the plane F1 = f 1 with the hyperboloid cylinder F2 = f 2 , these arcs have no equilibria. Therefore, γ f1 f2 has one or two trayectories of system (21). If f 2 = 0, then γ f1 f2 in x < 0 is the intersection of the plane F1 = f 1 with the two planes F2 = 0. The intersection { F1 = f 1 } ∩ { F2 = 0} in x < 0 is two straight lines which intersect at the equilibrium point ( f 1 , 0, 0). Such intersections are formed by 5 or 2 trajectories, respectively. We analyzed when a trajectory of γh1 h2 and a trajectory of γ f1 f2 gave place to a periodic orbit of the system (19). From above, if h2 > 0 and f 2 6= 0, then the trajectories of systems (20) and (21) could connect, forming a periodic orbit. In the plane of discontinuity x = 0, we considered the point (0, y0 , z0 ). Let h1 = z0 , h2 = y20 + z20 , f 1 = −z0 and f 2 = y20 − z20 be the four values of the first integrals H1 , H2 , F1 , and F2 at this point, respectively. Now, we studied the points of the trajectory Symmetry 2021, 13, 1736 15 of 21 { H1 = h1 } ∩ { H2 = h2 } ∩ { x ≥ 0} contained in x = 0. This was conducted by solving the system: H1 = h1 , H2 = h2 , x = 0. This system provides the two points (0, ±y0 , z0 ). The points of the trajectory { F1 = f 1 } ∩ { F2 = f 2 } ∩ { x ≤ 0} in x = 0, are studied solving the system F1 = f 1 , F2 = f 2 , x = 0. This system provides the two points (0, ±y0 , z0 ). When these pair of points are in the same trajectory { F1 = f 1 } ∩ { F2 = f 2 } ∩ { x ≤ 0}, we obtained a periodic orbit of the system (19). In x ≥ 0, using the variable x, we parameterized the trajectory { H1 = h1 } ∩ { H2 = h2 }, obtaining the arc: q q {( x, ± y20 + z20 − (z0 − x )2 , z0 − x ) : 0 ≤ x ≤ z0 + y20 + z20 }. (22) This trajectory in x ≥ 0 is symmetric with respect to the y–axis with endpoints (0, ±y0 , z0 ) in x = 0. In x ≤ 0, using the variable x, we parameterized the curve { F1 = f 1 } ∩ { F2 = f 2 }, which is formed by the two trajectories: q • {( x, ± x2 + y20 , x ) : x ≤ 0} if z0 = 0, each trajectory has one endpoint in x = 0; q • { x, ± ( x + z0 )2 + y20 − z20 , z0 + x )} : x ≤ 0} if either z0 < 0 or z0 > 0 and y20 − z20 > 0, each trajectory has one endpoint in x = 0; • q q { x, ± ( x + z0 )2 + y20 − z20 , z0 + x )} : −z0 + z20 − y20 ≤ x ≤ 0} (23) and q q { x, ± ( x + z0 )2 + y20 − z20 , z0 + x )} : x ≤ −z0 − z20 − y20 }, if z0 > 0 and y20 − z20 < 0, the first trajectory has its two endpoints at the points (0, ±y0 , z0 ) in x = 0, and the second trajectory has its endpoints at infinity. We note that y20 − z20 6= 0; otherwise, f 2 = 0. Summarizing, if z0 > 0 and y20 − z20 < 0, the trajectory (22) of system (20) with the trajectory (23) of system (21) provides a periodic orbit of system (19), and this periodic orbit intersects x = 0 in the points (0, ±y0 , z0 ). Hence, Theorem 4 was proved. 2.5. Limit Cycles of a Class of Piecewise Differential Systems Separated by a Parabola The goal of this subsection was to analyze the limit cycles of discontinuous piece- wise differential systems separated by the parabola y = x2 and formed by two linear Hamiltonian systems without equilibrium points. Easy computations show that a linear Hamiltonian system without equilibrium points must be of the form: Xi ( x, y) = (−λi bi x + bi y + γi , −λ2i bi x + λi bi y + δi ), δi 6= λi γi and bi 6= 0, with i = 1 . . . 4, and its corresponding Hamiltonian function is: Hi ( x, y) = (−λ2i bi /2) x2 + λi bi xy − (bi /2)y2 + δi x − γi y. We wanted to prove the following result, which came from [29]. Theorem 5. Generically, the maximum number of limit cycles of the piecewise differential systems separated by the parabola y = x2 and formed by two linear Hamiltonian systems without equilibrium points is two, and this maximum is reached, see Figure 5. Symmetry 2021, 13, 1736 16 of 21 1.0 0.5 0.0 -0.5 -1.0 -1.5 -2.0 -1.0 -0.5 0.0 0.5 1.0 Figure 5. Two limit cycles of a piecewise differential system separated by the parabola y = x2 and formed by two linear Hamiltonian systems without equilibrium points. Both limit cycles travel in counterclockwise sense. Proof. In region R1 = {( x, y) : y − x2 ≥ 0}, we considered the linear Hamiltonian system without equilibrium points: ẋ = −λ1 b1 x + b1 y + γ1 , ẏ = −λ21 b1 x + λ1 b1 y + δ1 , (24) with b1 6= 0 and δ1 6= λ1 γ1 . Its corresponding Hamiltonian function is: H1 ( x, y) = −(λ21 b1 /2) x2 + λ1 b1 xy − (b1 /2)y2 + δ1 x − γ1 y. (25) In region R2 = {( x, y) : y − x2 ≤ 0}, we considered another linear Hamiltonian system without equilibrium points: ẋ = −λ2 b2 x + b2 y + γ2 , ẏ = −λ22 b2 x + λ2 b2 y + δ2 , (26) with b2 6= 0 and δ2 6= λ2 γ2 . Its corresponding Hamiltonian function is: H2 ( x, y) = −(λ22 b2 /2) x2 + λ2 b2 xy − (b2 /2)y2 + δ2 x − γ2 y. (27) In order to have a crossing limit cycle which intersects the parabola y − x2 = 0 in the points ( xi , xi2 ) and ( xk , xk2 ), these points must satisfy the following system: H1 ( xi , xi2 ) − H1 ( xk , xk2 ) = 0, (28) H2 ( xi , xi2 ) − H2 ( xk , xk2 ) = 0, We supposed that the two systems (24) and (26) have three crossing limit cycles, and we arrived to a contradiction. Then, system (28) must have three pairs of points as solutions, namely, pi = (ri , ri2 ) and qi = (si , s2i ), with i = 1, 2, 3. Since the points p1 = (r1 , r12 ) and q1 = (s1 , s21 ) satisfy system (28), we obtained that the parameters γ1 and γ2 must be: 1 γ1 = (−r1 r13 − b1 r12 s1 − b1 r1 s21 − b1 s31 + 2δ1 + 2b1 r12 λ1 + 2b1 r1 s1 λ1 2(r1 + s1 ) +2b1 s21 λ1 − b1 r1 λ21 − b1 s1 λ21 ), and γ2 has the same expression that γ1 changes (b1 , λ1 , δ1 ) by (b2 , λ2 , δ2 ). Symmetry 2021, 13, 1736 17 of 21 If the second points p2 = (r2 , r22 ) and q2 = (s2 , s22 ) satisfy system (28), then the parameters δ1 and δ2 must be: b1 δ1 = (−r13 r2 − r1 r23 + r12 r2 s1 − r23 s1 + r1 r2 s21 + r2 s31 + r13 s2 2(r1 − r2 + s1 − s2 ) −r1 r22 s2 + r12 s1 s2 − r22 s1 s2 + r1 s21 s2 + s31 s2 − r1 r2 s22 − r2 s1 s22 − r1 s32 − s1 s32 −2r12 r2 λ1 + 2r1 r22 λ1 − 2r1 r2 s1 λ1 + 2r22 s1 λ1 − 2r2 s21 λ1 − 2r12 s2 λ1 + 2r1 r2 s2 λ1 −2r1 s1 s2 λ1 + 2r2 s1 s2 λ1 − 2s21 s2 λ1 + 2r1 s22 λ1 + 2s1 s22 λ1 ), and δ2 has the same expression that δ1 changes (b1 , λ1 ) by (b2 , λ2 ). Finally, we supposed that the points p3 = (r3 , r32 ) and q3 = (s3 , s23 ) satisfy system (28); then, the parameters λ1 and λ2 must be λ1 = A/B, where: A= r13 (r2 − r3 + s2 − s3 ) + r12 s1 (r2 − r3 + s2 − s3 ) + r23 (r3 − s1 + s3 ) + r22 s2 (r3 − s1 +s3 ) + r1 (−r23 + r33 − r3 s21 − r22 s2 + s21 s2 − s32 + r2 (s21 − s22 ) + r32 s3 − s21 s3 + r3 s23 +s33 ) + (s1 − s2 )(r33 + r32 s3 + (s1 − s3 )(s2 − s3 )(s1 + s2 + s3 ) − r3 (s21 + s1 s2 + s22 −s23 )) − r2 (r33 − s31 + s1 s22 + r32 s3 − s22 s3 + s33 + r3 (−s22 + s23 )), B= 2((s1 − s2 )(r32 + (s1 − s3 )(s2 − s3 ) − r3 (s1 + s2 − s3 )) + r12 (r2 − r3 + s2 − s3 ) +r22 (r3 − s1 + s3 ) + r1 (−r22 + r32 − r3 s1 + r2 (s1 − s2 ) + s1 s2 − s22 + r3 s3 − s1 s3 +s23 ) − r2 (r32 + r3 (−s2 + s3 ) − (s1 − s3 )(s1 − s2 + s3 )). Furthermore, λ2 has the same expression that λ1 changes b1 by b2 . We replaced γ1 , λ1 and δ1 in the expression of H1 ( x, y), and γ2 , λ2 and δ2 in the expression of H2 ( x, y), and we obtained H1 ( x, y) = H2 ( x, y). Therefore, the two linear differential systems forming the piecewise system coincide. Therefore the piecewise system has no limit cycles. Consequently, two is the maximum number of limit cycles. To complete the proof of the theorem, we presented a discontinuous piecewise differ- ential system satisfying the assumptions of the theorem, having two limit cycles. Let the parabola y = x2 be the discontinuity line of the piecewise differential system formed by the following two linear Hamiltonian systems without equilibria: ẋ = 5.5x − 0.5y + 3, ẏ = 60.5x − 5.5y + 0.2, (29) in the region R1 its Hamiltonian is: H1 ( x, y) = 30.25x2 − 5.5xy + 0.2x + 0.25y2 − 3y. The second system is: ẋ = 0.2x − 0.1y − 0.778814, ẏ = 0.4x − 0.2y + 0.00727332, (30) in the region R2 , its Hamiltonian is: H2 ( x, y) = 0.2x2 − 0.2xy + 0.00727332x + 0.05y2 + 0.778814y. This piecewise differential system has the limit cycles shown in Figure 5. Therefore, Theorem 5 was proved. 2.6. Piecewise Differential System with a Non-Regular Discontinuity Line The extended 16th Hilbert problem consists of finding, for a given class of differential systems, an upper bound for the maximum number of their limit cycles. This is, in general, a very difficult problem, in general unsolved. Only for very few classes of differential systems this problem has been solved. Now, we studied the extended 16th Hilbert problem for the piecewise differential systems separated by the non-regular line R formed by the two positive half-axes x and y, and formed by two linear centers. Symmetry 2021, 13, 1736 18 of 21 We denoted by R1 the open positive quadrant of R2 , and by R2 the interior of R2 \ R1 . It is known that an arbitrary linear center can be written as: ẋ = − Ax − ( A2 + Ω2 )y + B, for ( x, y) ∈ R1 , (31) ẏ = x + Ay + C, and ẋ = − ax − ( a2 + ω 2 )y + b, for ( x, y) ∈ R2 , (32) ẏ = x + ay + c, with Ω, ω > 0, A, B, C, a, b, c ∈ R, and A, a 6= 0. Each system (31) and (32) have, respectively, the first integrals: H1 ( x, y) = ( x + Ay)2 + 2(Cx − By) + y2 Ω2 , (33) H2 ( x, y) = ( x + ay)2 + 2(cx − by) + y2 ω 2 . The next result appears in [30]. Theorem 6. Consider the two arbitrary linear differential centers (31) and (32) forming the discontinuous piecewise differential systems separated by the non-regular line R. Then, the maximum number of limit cycles of these piecewise systems intersecting R in two points is two. Moreover there exist systems with exactly two limit cycles of this type, see Figure 6. 5 4 3 2 1 0 -1 -8 -6 -4 -2 0 Figure 6. Both limit cycles are travelled in counter-clockwise sense. Proof. We considered the discontinuous planar linear differential systems (31) and (32). If there exists a crossing limit cycle intersecting the non-regular separation curve R, in two points of the forms ( x, 0) and (0, y) both are different from the origin. Since the functions H1 and H2 defined in (33) are first integrals of the systems (31) and (32), respectively, these points must satisfy the equations: e1 := H1 ( x, 0) − H1 (0, y) = 2Cx + x2 + 2By − A2 y2 − y2 Ω2 = 0, (34) e2 := H2 ( x, 0) − H2 (0, y) = 2cx + x2 + 2by − a2 y2 − y2 ω 2 = 0. Symmetry 2021, 13, 1736 19 of 21 By using the Bézout theorem, this system can have, at most, four isolated solutions, one of them being the origin. In order to obtain crossing limit cycles, equations e1 = 0 and e2 = 0 must have isolated solutions ( xi , yi ) with xi , yi > 0. Therefore, there are at most three crossing limit cycles of system e1 = e2 = 0. In order for these three solutions ( xi , yi ) to produce limit cycles, it is necessary that: 0 < x1 < x2 < x3 and 0 < y1 < y2 < y3 . We claimed that there are at most two solutions ( x1 , y1 ) and ( x2 , y2 ) providing limit cycles, that is, satisfying: 0 < x1 < x2 and 0 < y1 < y2 . (35) Now, we proved the claim. If a2 + ω 2 − A2 − Ω2 = 0, then the piecewise system has at most one limit cycle. Indeed, the resultant of the polynomial e1 and e2 with respect to the variable y is: 4x a2 + ω 2 ( x a2 c2 − 2a2 cC + a2 C2 − b2 + 2bB − B2 + c2 ω 2 − 2cCω 2 + C2 ω 2   −2(b − B)(bC − Bc)). Therefore, at most, one positive solution of x, consequently, is at most one limit cycle. Assume now that a2 + ω 2 − A2 − Ω2 6= 0. Equations e1 = e2 = 0 are equivalent to equations E1 = e1 − e2 = 0 and E2 = e1 ( a2 + ω 2 ) − e2 ( A2 + Ω2 ) = 0, i.e.,: E1 =2(C − c) x + 2( B − b)y − ( A2 + Ω2 − a2 − ω 2 ))y2 = 0, E2 =2 ( a2 + ω 2 )C − ( A2 + Ω2 )c x + 2 ( a2 + ω 2 ) B − ( A2 + Ω2 )b y   (36) + ( a2 + ω 2 − A2 − Ω2 ) x2 = 0. If C = c, then E1 = 0 reduces to either one horizontal straight line, or two horizontal parallel straight lines passing one of these two straight lines through the origin. The equation E2 = 0 is either a parabola symmetric with respect to some vertical straight line, one vertical straight line, or two vertical parallel straight lines passing one of these two straight lines through the origin. Since E1 = E2 = 0 pass through the origin, there are at most two intersection points satisfying (35) and so, at most, two limit cycles. Assume now that C 6= c. In this case, E1 = 0 is a parabola symmetric with respect to some horizontal straight line and E2 = 0 is a parabola symmetric with respect to some vertical straight line. Since both parabolas intersect at the origin, there are at most two intersection points satisfying (35), and so, at most, two limit cycles. This proves the claim and, consequently, the theorem once we provided an example with two limit cycles. Now, we gave a discontinuous piecewise linear differential system (31)–(32) having exactly two limit cycles intersecting in two points the discontinuity line R. In region R1 , we considered the linear differential center: 3 43 ẋ = −2x − 8y − , ẏ = x + 2y + , (37) 2 4 with the first integral:   2 43 3 H1 ( x, y) = 4y + 2 x + y + ( x + 2y)2 ; 4 2 and in region R2 , we considered the linear differential center: 7 ẋ = − x − 2y, ẏ = x + y + , (38) 4 Symmetry 2021, 13, 1736 20 of 21 with the first integral: 7 H2 ( x, y) = y2 + x + ( x + y)2 . 2 In this case, the two solutions of Equation (36) are:     1 3 ( x1 , y1 ) = ,1 , ( x2 , y2 ) = 1, , 2 2 and the corresponding limit cycles are shown in Figure 2. 3. Discussion In this paper, we summarized the main results on the Darboux theory of integrability for finding first integrals. We illustrated, with some relevant examples, the different main ingredients of this theory, as the invariant algebraic curves, the exponential factors, the integrating factors, and the first integrals. After, we used the first integrals of distinct classes of piecewise differential systems for studying the limit cycles of these differential systems in R2 and R3 . Thus, first we studied the limit cycles of a discontinuous piecewise differential system in R2 with two zones separated by a straight line and formed by a linear differential center and a quadratic polynomial differential system. After, we computed the three limit cycles of a piecewise differential system in R2 with three zones separated by the non-regular line : {( x, y) ∈ R2 : y = 0 or x = 0 and y ≥ 0}, and each zone having an arbitrary linear differential system. We also studied the periodic orbits of a relay system in R3 . We analyzed the limit cycles of discontinuous piecewise differential systems in R2 with two zones separated by a parabola and each zone having a Hamiltonian system without equilibrium points. Finally, we proved that the maximum number of limit cycles of a piecewise differential system in R2 with two zones separated by the non-regular line formed by the positive x and y half-axes and having, in each zone, an arbitrary linear differential system was two. We also provided an example of these differential systems having exactly two limit cycles. 4. Conclusions We illustrated how to compute first integrals of the polynomial differential systems via the Darboux theory of integrability, and we also illustrated how to compute periodic orbits and limit cycles of different classes of piecewise differential systems using their first integrals. Funding: The author was supported by the Agencia Estatal de Investigación grant PID2019-104658GB- I00, the Agència de Gestió d’Ajuts Universitaris i de Recerca grant 2017SGR1617, and the H2020 European Research Council grant MSCA-RISE-2017-777911. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Conflicts of Interest: The author declares no conflict of interest. References 1. Whittaker, E.T. A Treatise on the Analytical Dynamics of Particles and Rigid Bodies; Cambridge Mathematical Library, Cambridge University Press: Cambridge, UK, 1988. 2. Dumortier, F.; Llibre, J.; Artés, J.C. Qualitative Theory of Planar Differential Systems; UniversiText; Springer: New York, NY, USA, 2006. 3. Olver, P.J. Applications of Lie Groups to Differential Equations, 2nd ed.; Graduate Texts in Mathematics; Springer: New York, NY, USA, 1993; Volume 107. Symmetry 2021, 13, 1736 21 of 21 4. di Bernardo, M.; Budd, C.J.; Champneys, A.R.; Kowalczyk, P. Piecewise-Smooth Dynamical Systems: Theory and Applications; Applied Mathematical Sciences; Springer: London, UK, 2008; Volume 163. 5. Leine, R.I.; Nijmeijer, H. Dynamics and Bifurcations of Non-Smooth Mechanical Systems; Lecture Notes in Applied and Computational Mechanics; Springer: Berlin/Heidelberg, Germany, 2004; Volume 18. 6. Liberzon, D. Switching in Systems and Control: Foundations and Applications; Birkhäuse: Boston, MA, USA, 2003. 7. Simpson, D.J.W. Bifurcations in Piecewise-Smooth Continuous Systems; World Scientific Series on Nonlinear Science Series A; World Scientific: Singapore, 2010; Volume 69. 8. Makarenkov, O.; Lamb, J.S.W. Dynamics and bifurcations of nonsmooth systems: A survey. Physica D 2012, 241, 1826–1844. [CrossRef] 9. Filippov, A.F. Differential Equations with Discontinuous Right-Hand Sides; Kluwer Academic: Dordrecht, The Netherlands, 1988. 10. Lum, R.; Chua, L.O. Global propierties of continuous piecewise-linear vector fields. Part I: Simplest case in R2 . Int. J. Circuit Theory Appl. 1991, 19, 251–307. [CrossRef] 11. Huan, S.M.; Yang, X.S. Existence of limit cycles in general planar piecewise linear systems of saddle–saddle dynamics. Nonlinear Anal. 2013, 92, 82–95. [CrossRef] 12. Braga, D.C.; Mello, L.F. Limit cycles in a family of discontinuous piecewise linear differential systems with two zones in the plane. Nonlinear Dynam. 2013, 73, 1283–1288. [CrossRef] 13. Buzzi, C.; Gasull, A.; Torregrosa, J. Algebraic limit cycles in piecewise linear differential systems. Int. J. Bifurc. Chaos 2018, 28, 1850039. [CrossRef] 14. Buzzi, C.; Pessoa, C.; Torregrosa, J. Piecewise linear perturbations of a linear center. Discret. Contin. Dyn. Syst. 2013, 33, 3915–3936. [CrossRef] 15. Freire, E.; Ponce, E.; Rodrigo, F.; Torres, F. Bifurcation sets of continuous piecewise linear systems with two zones. Int. J. Bifurc. Chaos 1998, 8, 2073–2097. [CrossRef] 16. Freire, E.; Ponce, E.; Torres, F. Canonical Discontinuous Planar Piecewise Linear Systems. SIAM J. Appl. Dyn. Syst. 2012, 11, 181–211. [CrossRef] 17. Freire, E.; Ponce, E.; Torres, F. A general mechanism to generate three limit cycles in planar Filippov systems with two zones. Nonlinear Dyn. 2014, 78, 251–263. [CrossRef] 18. Giannakopoulos, F.; Pliete, K. Planar systems of piecewise linear differential equations with a line of discontinuity. Nonlinearity 2001, 14, 1611–1632. [CrossRef] 19. Han, M.; Zhang, W. On Hopf bifurcation in non—Smooth planar systems. J. Differ. Equ. 2010, 248, 2399–2416. [CrossRef] 20. Huan, S.M.; Yang, X.S. On the number of limit cycles in general planar piecewise systems. Discret. Cont. Dyn. Syst. Ser. A 2012, 32, 2147–2164. [CrossRef] 21. Li, L. Three crossing limit cycles in planar piecewise linear systems with saddle-focus type. Electron. J. Qual. Theory Differ. Equ. 2014, 70, 1–14. [CrossRef] 22. Wang, J.; Huang, C.; Huang, L. Discontinuity-induced limit cycles in a general planar piecewise linear system of saddle-focus type. Nonlinear Anal. Hybrid Syst. 2019, 33, 162–178. [CrossRef] 23. Chen, H.; Li, D.; Xie, J.; Yue, Y. Limit cycles in planar continuous piecewise linear systems. Commun. Nonlinear Sci. Numer. Simul. 2017, 47, 438–454. [CrossRef] 24. Falconi, M.; Llibre, J. n − 1 independent first integrals for Linear differential systems in Rn and Cn . Qual. Theory Dyn. Syst. 2004, 4, 233–254. [CrossRef] 25. Schlomiuk, D. Algebraic particular integrals, integrability and the problem of the center. Trans. Am. Math. Soc. 1993, 338, 799–841. [CrossRef] 26. Llibre, J.; Zhang, X. Limit cycles created by piecewise linear centers. Chaos 2019, 29, 053116. [CrossRef] [PubMed] 27. Llibre, J.; Teixeira, M.A. Periodic orbits of continuous and discontinuous piecewise linear differential systems via first integrals. Sao Paulo J. Math. Sci. 2018, 12, 121–135. [CrossRef] 28. Anosov, D.V. Stability of the equilibrium positions in relay systems. Avtomatika i Telemehanika 1959, 20, 135–149. 29. Benterki, R.; Llibre, J. Crossing Limit cycles of planar piecewise linear Hamiltonian systems without equilibrium points. Mathematics 2020, 8, 755. [CrossRef] 30. Esteban, M.; Llibre, J.; Valls, C. The extended 16-th Hilbert problem for discontinuous piecewise isochronous centers of degree one or two separated by a straight line. Int. J. Bifurc. Chaos 2021, 31, 043112. [CrossRef]