New Advances on the Lyapunov Constants of Some Families of Planar Differential Systems

. This note presents some advances regarding the Lyapunov constants of some families of planar polynomial diﬀerential systems, as a ﬁrst step towards the resolution of the center and cyclicity problems. Firstly, a parallelization approach is computationally implemented to achieve the 14th Lyapunov constant of the complete cubic family. Sec-ondly, a technique based on interpolating some speciﬁc quantities so as to reconstruct the structure of the Lyapunov constants is used to study a Kukles system, some ﬁfth-degree homogeneous systems and a quartic system with two invariant lines.


Introduction
Let us consider a real polynomial differential system in the plane with some parameters, λ ∈ R d , written in complex coordinates as ż = iz + Z(z, w, λ), where w =z and Z(z, w, λ), W (z, w, λ) =Z(z, w, λ) are polynomial perturbations having neither linear nor constant terms in z, w. The center problem consists in identifying whether the origin of (1) is a center or a focus, when the origin is a monodromic nondegenerate equilibrium point. This problem is related with the local cyclicity problem, which aims to determine the maximum number of small limit cycles which bifurcate from the origin when perturbing the system in a polynomial class of fixed degree. All of them are relevant studies in the 16th Hilbert problem.
To deal with this problem let us consider the Poincaré map Π, which maps a given point ρ in a section Σ transversal to the orbit γ to the first intersection Π(ρ) of Σ and γ in positive time. The Poincaré map can be analytically extended to ρ = 0, so we can consider its Taylor expansion and define the displacement map for certain values V n . Observe that the center problem is equivalent to determine whether all V n are zero or not, since periodic orbits are fixed points of the Poincaré map. If not all V n vanish, the first non-zero V n must have odd subindex. These V n with odd n ≥ 3 are known as Lyapunov constants, and they will be denoted by L (n−1)/2 := V n . According to [6], the Lyapunov constants are polynomials, whose variables are the parameters of system (1). In this case for which not all V n vanish, the origin of the system is a focus and its stability is determined by the first non-zero Lyapunov constant. As a consequence, the center problem reduces to the problem of finding and vanishing all the Lyapunov constants L k , that is solving the nonlinear system {L 1 = L 2 = · · · = 0}. They also are an essential tool in the study of the cyclicity problem, since the limit cycles correspond to the isolated zeroes of (2). In this short paper we will briefly present two methods to compute these Lyapunov constants or study some of their useful properties to deal with the center and cyclicity problems. This is a preprint of: "New advances on the Lyapunov constants of some families of planar differential systems", Iván Sánchez Sánchez, Joan

A parallelization of Lyapunov method
An algorithm to find Lyapunov constants is the so called Lyapunov method, which is based on the utilization of a first integral of system (1). The computations could be made using real values (see [8]), but we consider complex coordinates because the obtained expressions are shorter. The objective is then to find a formal first integral F = F 2 + F 3 + F 4 + · · · of system (1), with F k (z, w) := k j=0 h k−j,j z k−j w j homogeneous degree k polynomials. We aim to study whetherḞ vanishes or not. We computė The last equality is a consequence of a result which states that if there exists such F, then in suitable coordinates it is analytic on zw (see [3]). ThereforeḞ is also analytic on zw. Observe that if all L k (λ) vanish thenḞ = 0, and therefore F is a first integral, so the origin is a center. These L k (λ) are actually the Lyapunov constants, which are polynomials in the parameters λ (see [4]). For the sake of simplicity, we will denote them simply as L k . F can be found recursively, starting by imposing equation (3) and performing formal operations as follows: Here we solve the linear system obtained by equating the same degree coefficients in z, w. This way one can find the coefficients h k−j,j of each term F k of the first integral and use it to find the corresponding Lyapunov constants.
The above technique has been implemented in PARI/GP (or simply PARI) programming language, see [9]. As the computation of Lyapunov constants is a highly computationally expensive procedure, this algorithm has been optimized and improved by means of parallelization, which allows to significantly increase computation velocity. The idea is to find each of the Lyapunov constants and the coefficients h k−j,j of F k of degree k in terms of the coefficients of lower degree. This part is relatively fast computationally speaking since the manipulated expressions are not too large. Then we parallelize the substitution of those coefficients with their actual value, and here parallelization is essential because this process deals with very large expressions.
The results of this parallelization technique are amazing, and its efficiency has allowed our method to find Lyapunov constants in a relatively short time for cases which had not been solved before due to the huge amount of time and computational complexity required. In particular, we have applied this method to the complete cubic system ż = iz +r 20 z 2 +r 11 zw +r 02 w 2 +r 30 z 3 +r 21 z 2 w +r 12 zw 2 +r 03 w 3 , w = −iw +ŝ 20 w 2 +ŝ 11 wz +ŝ 02 z 2 +ŝ 30 w 3 +ŝ 21 w 2 z +ŝ 12 wz 2 +ŝ 03 z 3 .
We have observed that if time is rescaled by dividing by the imaginary unit i, computations are much more efficient and the calculation time decreases. Actually, the computations we describe here cannot be performed if this time rescaling is not done. If we denote r jk =r jk /i and s jk =ŝ jk /i, the system in the new time variable can be written as z = z + r 20 z 2 + r 11 zw + r 02 w 2 + r 30 z 3 + r 21 z 2 w + r 12 zw 2 + r 03 w 3 , w = −w + s 20 w 2 + s 11 wz + s 02 z 2 + s 30 w 3 + s 21 w 2 z + s 12 wz 2 + s 03 z 3 .
Up to our knowledge, the highest known Lyapunov constant for the above system is the 10th (see [7]). But with our parallelization technique we have been able to reach the 14th. To perform this computation we have used the computer network of our department. The parallelization has been done with the software PBala, see [5]. This server has eight nodes with Intel Xeon 2.60GHz processors. The total memory is 640GB and 96 threads can be run at the same time. The size of the new Lyapunov constants are shown in Table 1 It is well-known, see for example [2], that the solution of the center problem for the general cubic differential system needs at least 11 Lyapunov constants. Now, the obstacle to obtain a complete characterization of the cubic centers is how can we solve the nonlinear system {L 1 = L 2 = · · · = 0}, and not how to construct it because we think that we have computed enough Lyapunov constants.

Interpolation and reconstruction technique
Let us consider the ideal generated by all Lyapunov constants L 1 , L 2 , L 3 , . . . associated to the differential system (1) with some fixed degree. This is an ideal of the polynomial ring with coefficients in C. Due to the Hilbert Basis Theorem, this ideal is finitely generated, so there must exist m ∈ N such that L 1 , L 2 , L 3 , . . . = L 1 , L 2 , L 3 , . . . , L m .
To know this m would significantly simplify the problem, because by computing the first m Lyapunov constants we would obtain the center conditions. Nevertheless, as [1] states, there are no general methods to find this m and this is the reason why the center problem has been solved only for certain polynomial families. Let B k := L 1 , . . . , L k be the Bautin ideal generated by the first k Lyapunov constants. The method we suggest aims to check whether a certain Lyapunov constant L n belongs to B n−1 , and therefore it vanishes when the previous are equal to zero. It is important to remark that to apply this techique at this step we assume that we have been able to compute the first n Lyapunov constants.
Let us start by writing where A j are polynomials whose variables are the parameters of (1). Our method consists in trying to see whether we can determine these polynomials A j , since this will tell if expression (7) is possible or not. Using the notation of (5) for the parameters, let us consider a monomial M = k, r p k k s q k k , where r k and s k correspond to the coefficients of z k w of Z(z, w, λ) and W (z, w, λ), respectively. We define the quasi degree of M as k, (k + − 1)(p k + q k ) and its weight as k, (1 − k − )(p k − q k ). Then, by [4], the monomials of a Lyapunov constant L j satisfy that they have quasi-degree 2j and weight 0. Now using these properties together with the degree of L j , we can select which monomials are candidates to be part of each A j , but with undetermined coefficients. Thus, we have that A j are polynomials whose monomials have been selected and have undetermined coefficients, and these coefficients of A j are what we try to compute. Now knowing the structure of A j , we would substitute it in (7), expand the products and the sum and finally equate the coefficients of monomials with the same literal part. This gives a set of linear equations consisting of the coefficients of equality (7). If this system of linear equations is compatible, then the polynomials A j do exist and L n vanishes when L 1 , . . . , L n−1 are zero; otherwise, if the system is incompatible then the polynomials A j do not exist and L n does not belong to B n−1 . Instead of explicitly solving the system of equations, we have compared ranks of the system matrices to see if they are equal or not.
With this method we have studied three different families and we have obtained the following results. Proposition 1. Consider the Kukles differential system which in complex coordinates is written as ż = iz + r 20 z 2 + r 11 zw + r 02 w 2 + r 30 z 3 + r 21 z 2 w + r 12 zw 2 + r 03 w 3 , w = −iw − r 20 z 2 − r 11 zw − r 02 w 2 − r 30 z 3 − r 21 z 2 w − r 12 zw 2 − r 03 w 3 .
Then L 9 does not belong to B 8 . In the case r 12 = 0, L 9 again does not belong to B 8 , but L 10 does belong to B 9 since there exist A j such that L 10 = 9 j=1 A j L j .
From the above result, we can guess that only the first 9 Lyapunov constants are enough to solve the center problem for this family when r 12 = 0. This interpolation method works better than the standard approach using Groebner basis for the simplifications.
The center problem for degree 5 homogeneous perturbations of the linear oscillator is an open problem, and even the value m in (6) is unknown. The mechanism proposed here fails for the general family due to the size of the computations. The next result presents some particular cases.
Then the next properties hold: • If r 41 = s 41 = 0, then L 10 does not belong to B 9 but L 11 does belong to B 10 .
• If r 32 = s 32 = 0, then L 11 and L 12 do not belong to B 10 and B 11 , respectively.
The last considered family is a special quartic differential system with four invariant straight lines.
For this system both L 7 , L 8 , and L 9 do not belong to B 6 , B 7 , and B 8 , respectively. However, when r 11 = s 11 = 0, L 7 does belong to B 6 .