Efficient Algorithms for Computing Rational First Integrals and Darboux Polynomials of Planar Polynomial Vector Fields

Let ${\mathbb{K}}$ denote an effective field of characteristic zero, i.e, one can perform arithmetic operations and test equality of two elements (typically, ${\mathbb{K}}={\mathbb{Q}}$ or ${\mathbb{Q}}(\alpha)$, where $\alpha$ is an algebraic number). Given two polynomials $A,\,B$ in ${\mathbb{K}}[x,y]$, we consider the planar polynomial vector field

and discuss the problem of computing rational first integrals of ($\mathsf{S}$), i.e., rational functions $F\in{\mathbb{K}}(x,y)\setminus{\mathbb{K}}$ that are constant along the solutions $(x(t),y(t))$ of ($\mathsf{S}$). More precisely, the present article is concerned with the following algorithmic problem:

$({\mathcal{P}}_{N})$: given a degree bound $N\in{\mathbb{N}}$, either compute a rational first integral $F\in{\mathbb{K}}(x,y)\setminus{\mathbb{K}}$ of ($\mathsf{S}$) of total degree at most $N$, or prove that no such $F$ exists.

This old problem was already studied by Darboux in 1878 ([Dar78]) and has been the subject of numerous works ever since. The naive approach (by indeterminate coefficients) leads to a polynomial system of quadratic equations in the coefficients of $F$. Other methods use what is called nowadays Darboux polynomials, in the spirit of the celebrated Prelle-Singer’s method [PS83]; see Subsection 2.3 for a review. These methods also require solving a polynomial system of quadratic equations. Recently, Chèze [Chè11] has shown that problem $(\mathcal{P}_{N})$ can be solved in polynomial time in $N$. The importance of this result is mainly theoretical since the exponent in the polynomial complexity estimate is bigger than 10.

To improve upon this current state of affairs, our starting point is the article [FG10] of Ferragut and Giacomini. The key observation is that ($\mathsf{S}$) has a rational first integral if and only if all power series solutions in ${\mathbb{K}}[[x]]$ of the first order non-linear differential equation

are algebraic over ${\mathbb{K}}(x)$. Furthermore, minimal polynomials of these algebraic power series lead to rational first integrals.

The algorithm in [FG10] still involves solving a polynomial system of quadratic equations. Indeed, the key observation above is merely used to reduce the number of equations in the quadratic system provided by the naive approach.

In the present article, we push further the observation of Ferragut and Giacomini, so as to give fast algorithms solving Problem $({\mathcal{P}}_{N})$. In particular, we prove that this can be done by considering only systems of linear equations instead of systems of quadratic equations.

We design a probabilistic algorithm that uses $\tilde{\mathcal{O}}(N^{2\omega})$ arithmetic operations in ${\mathbb{K}}$, where $\omega\in[2,3]$ is the exponent of linear algebra over ${\mathbb{K}}$, and the soft-O notation $\tilde{\mathcal{O}}(\,)$ indicates that polylogarithmic factors are neglected. The probabilistic algorithm is then turned into a deterministic one, that solves Problem $({\mathcal{P}}_{N})$ in arithmetic complexity $\tilde{\mathcal{O}}(d^{2}\,N^{2\omega+1})$, where $d=\max(\deg(A),\deg(B))$ denotes the degree of the polynomial vector field ($\mathsf{S}$). This compares well to the previous polynomial time algorithm given in [Chè11], which uses at least $d^{\omega+1}N^{4\omega+4}$ arithmetic operations. Note that if we take $\omega=3$ (i.e., the cost of naive linear algebra), then the above means that the best previously known complexity would be in $\tilde{\mathcal{O}}(d^{4}N^{16})$ whereas our deterministic algorithm would use at most $\tilde{\mathcal{O}}(d^{2}N^{7})$ arithmetic operations, and our probabilistic one would use $\tilde{\mathcal{O}}(N^{6})$. Lastly, we sketch a heuristic method that uses $\tilde{\mathcal{O}}(N^{\omega+2})$ arithmetic operations (i.e., $\tilde{\mathcal{O}}(N^{5})$ using classical linear algebra) which is sub-cubic, given that the output has size $\mathcal{O}(N^{2})$.

We provide algorithmic details, notably precise degree bounds and complexity estimates. The algorithms developed in the article are implemented in a Maple package called RationalFirstIntegrals which is available with various examples at http://www.ensil.unilim.fr/~cluzeau/RationalFirstIntegrals.html. Using this implementation, we demonstrate the efficiency of our algorithms on some examples. Finally, we show how to apply a similar method to the computation of Darboux polynomials.

In Section 2, we recall Darboux’s approach to the integrability of polynomial vector fields, related works, and existing results about the problem $({\mathcal{P}}_{N})$. We also give useful facts on the so-called spectrum problem. We recall in Section 3 the connection between rational first integrals of the polynomial vector field ($\mathsf{S}$) and algebraic power series solutions of $\eqref{diff-eq-intro}$. We then propose a first algorithm, based on linear algebra, that solves Problem $({\mathcal{P}}_{N})$. Building on this, we develop in Section 4 an efficient probabilistic algorithm, and then turn it into an efficient deterministic algorithm. In Section 5, we study the arithmetic complexity of the algorithms developed in Section 4, and discuss several algorithmic issues. Then, in Section 6 we present our implementation and display its behavior on various examples. Finally, Section 7 shows how similar ideas can be used for computing the set of all irreducible Darboux polynomials (of a given degree) of planar polynomial vector fields.

The degree $\deg(P)$ of a bivariate polynomial $P\in{\mathbb{K}}[x,y]$ is the total degree of $P$. A rational function $P/Q$ with $P,Q\in{\mathbb{K}}[x,y]$ is said to be reduced when $P$ and $Q$ are coprime. The degree $\deg(F)$ of a reduced rational function $F=P/Q$ is the maximum of $\deg(P)$ and $\deg(Q)$.
We denote by $\overline{{\mathbb{K}}}$ an algebraic closure of the field ${\mathbb{K}}$.
We write $\dot{f}:=\frac{\partial f}{\partial t}$ for the usual formal derivative of the “function” (polynomial, or power series) $f$ with respect to the variable $t$.
For a set $\Omega$, we denote by $|\Omega|$ its cardinality.

We consider an autonomous planar polynomial vector field

where $x$ and $y$ are unknown “functions” of the time variable $t$, $A$ and $B$ are polynomials in ${\mathbb{K}}[x,y]$, and $d:=\max(\deg(A),\deg(B))$ denotes the degree of the polynomial vector field. Without any loss of generality, $A$ and $B$ will be assumed to be coprime in the remaining of the article.

To ($\mathsf{S}$) is attached the derivation

acting on the polynomial ring ${\mathbb{K}}[x,y]$. We thus view ${\mathbb{K}}[x,y]$ as a differential ring endowed with the derivation ${\mathcal{D}}$. We denote by ${\mathbb{K}}(x,y)$ its field of fractions.

A rational first integral $F$ of ($\mathsf{S}$) is thus a non-trivial constant for the derivation ${\mathcal{D}}$. Intuitively, this means that if $(x(t),y(t))$ is a pair of “functions” satisfying ($\mathsf{S}$), then $F(x(t),y(t))$ is constant when $t$ varies. We explain in Theorem 11 below why no algebraic extension of the base field is necessary in Definition 1.

A starting observation is that the rational function $F=P/Q$ is a first integral for ($\mathsf{S}$) if and only if ${\mathcal{D}}(P)\,Q={\mathcal{D}}(Q)\,P$. Therefore, if $F$ is a reduced rational first integral for ($\mathsf{S}$), then $P$ divides ${\mathcal{D}}(P)$, and $Q$ divides ${\mathcal{D}}(Q)$ in ${\mathbb{K}}[x,y]$. This motivates the following definition.

Darboux polynomials were introduced by G. Darboux in [Dar78]. These polynomials correspond to algebraic curves invariant under the vector field. The following lemma will be used in the sequel: it means that if we have a non-singular initial condition, then there is a unique irreducible invariant algebraic curve satisfying this initial condition, see [Sin92, Lemma A.1].

Darboux polynomials are sometimes called partial first integrals in the literature. The reason is that rational first integrals and Darboux polynomials are intimately related notions: as sketched above, numerators and denominators of reduced rational first integrals are Darboux polynomials. The converse is also true, see Corollary 5 below.

A fundamental property of Darboux polynomials is given in the following lemma (see, e.g., [DLA06, Lemma 8.3, p. 216]) that can be proved by a straightforward calculation.

As a corollary we get:

The previous corollary gives a relation between Darboux polynomials and rational first integrals. The next theorem shows that if we have enough Darboux polynomials, then we have a rational first integral, see [Sin92, Appendix] for a modern proof.

A useful corollary of Theorem 6 is the following, see [Sin92]:

Given a derivation ${\mathcal{D}}$, the problem of finding a bound for the degree of irreducible Darboux polynomials is known to be difficult: this is the so-called Poincaré problem. It has been deeply studied in the literature and many partial results exist ([Poi91, CLN91, Car94, Per02, Wal00, CG06, LY05] and others) though the question is not fully solved yet. The fact that the derivation ${\mathcal{D}}=n\,x\,\frac{\partial}{\partial x}+y\,\frac{\partial}{\partial y}$ with $n\in{\mathbb{N}}^{*}$ admits $x-y^{n}$ as an irreducible Darboux polynomial shows that a bound depending only on the degrees of the entries cannot exist: arithmetic conditions on the coefficients of ${\mathcal{D}}$ have to be taken into account as well.

Consequently, given a planar polynomial vector field ($\mathsf{S}$), or equivalently a derivation ${\mathcal{D}}$, two distinct problems occur when we want to compute rational first integrals:

1. (1)

   Find a bound on the degree of the numerator and denominator of a rational first integral, that is a bound on the degree of irreducible Darboux polynomials;

2. (2)

   $({\mathcal{P}}_{N})$: given a degree bound $N\in{\mathbb{N}}$, either compute a rational first integral $F\in{\mathbb{K}}(x,y)\setminus{\mathbb{K}}$ of ($\mathsf{S}$) of total degree at most $N$, or prove that no such $F$ exists.

Our aim is to give an efficient algorithm to handle the second problem $({\mathcal{P}}_{N})$.

In this article we suppose that $d\leq N$. This hypothesis is natural because if a derivation has a polynomial first integral of degree $N$, then we can show that $d\leq N-1$, see [FL07, Theorem 6] or [Poi91].

We recall here the definition of composite rational functions and what is called the spectrum of a rational function. We then use these notions to describe the kernel of the derivation ${\mathcal{D}}$ and to give some of its properties.

In [MO04, Chè12b], the authors propose different algorithms for the decomposition of rational functions using properties of Darboux polynomials and rational first integrals of the Jacobian derivation.

A more precise version of Lemma 9 is given by the next theorem which completely describes the ${\mathbb{K}}$-algebra structure of the set of all rational first integrals of ($\mathsf{S}$). This theorem seems to be a folklore result but we have not found a suitable reference. Consequently, a complete proof is provided here.

As a consequence of Theorem 10, non-composite reduced rational first integrals coincide with rational first integrals with minimal degree; they will play a key role in the remaining of this text.

In Definition 1, we have defined rational first integrals as elements of ${\mathbb{K}}(x,y)$. The next theorem explains why it is in general not necessary to consider rational first integrals in $\overline{{\mathbb{K}}}(x,y)$. To our knowledge this result is proved here for the first time. In [MM97], the authors show that if there exists a rational first integral in $\overline{{\mathbb{K}}}(x,y)$, then there also exists a rational first integral in ${\mathbb{K}}(x,y)$. We improve this result by taking into account the degrees of these rational first integrals.

Now, we introduce the spectrum of a rational function which will play a crucial role in our algorithms.

In the context of rational first integrals of polynomial vector fields, the elements of the spectrum are sometimes called remarkable values, see, e.g., [FL07]. There exists a vast bibliography about the spectrum, see for example [Rup86, Lor93, Vis93b, Vis93a, AHS03, Bod08, BC11].

The spectrum $\sigma(P,Q)$ is finite if and only if $P/Q$ is non-composite and if and only if the pencil of algebraic curves $\lambda\,P-\mu\,Q=0$ has an irreducible general element, see for instance [Jou79, Chapitre 2, Théorème 3.4.6] or [Bod08, Theorem 2.2] for detailed proofs.

To the authors’ knowledge, the first effective result on the spectrum is due to Poincaré. In [Poi91], he establishes a relation between the number of saddle points and the cardinal of the associated spectrum, in the case where all the singular points of the polynomial vector field are distinct. In particular this yields the bound $|\sigma(P,Q)|\leq d^{2}$ on the cardinality of the spectrum. This bound was improved recently in [Chè12a]:

As a consequence of Theorem 13 and Corollary 5, if $P/Q$ is a reduced non-composite rational first integral of ($\mathsf{S}$), then for all but $\mathcal{B}(d)+1$ constants $\sigma\in\overline{{\mathbb{K}}}$, the polynomial $P-\sigma\,Q$ has degree $N$ and is an irreducible Darboux polynomial for ${\mathcal{D}}$ with the same cofactor as $P$ and $Q$. This means that if ($\mathsf{S}$) has a rational first integral, there exist an infinite number of irreducible Darboux polynomials which all have the same degree (and the same cofactor).

The following lemma will be useful in Section 4 for the study of our probabilistic and deterministic algorithms.

There is a vast literature regarding the computation of rational or elementary first integrals (see for example [FG10, MM97, Chè11, PS83, SGR90, Sin92, LZ10, DDdMS01, Poi91]) and Darboux polynomials (see for example [CMS09, CMS06, CGG05, Wei95, Dar78]). Note that, among these articles, very few restrict to the specific question of rational first integrals. Surveys on computing first integrals (not restricted to planar systems) can be found for example in [Gor01, Sch93, DLA06] and [Olv93] for symmetry methods which we do not address here.

Given a degree bound $N$, the naive approach to solve $({\mathcal{P}}_{N})$ consists in using the method of undetermined coefficients. This leads to a system of polynomial (quadratic) equations in the unknown coefficients of the rational first integral, see [Chè11] for a complexity estimate of this approach.

Interest in Darboux polynomials has been revived by the appearance of the Prelle-Singer’s method, [PS83, MM97, DDdMS01]. In [CMS06, CMS09], Coutinho and Menasché Schechter give necessary conditions for the existence of Darboux polynomials. Other necessary conditions are contained in [CGGL03, CGG05] and also in works on inverse integrating factors [CGGL03, CFL10]. The bottleneck of the Prelle-Singer’s method and all of its variants is the computation by undetermined coefficients of all irreducible Darboux polynomials of bounded degree, which leads again to solving a polynomial system. This yields an exponential complexity algorithm, see [Chè11].

In [Chè11], Chèze shows that if the derivation ${\mathcal{D}}$ admits only finitely many irreducible Darboux polynomials of degree at most $N$, then it is possible to compute all of them by using the so-called ecstatic curve introduced in [Per01] within a number of binary operations that is polynomial in the bound $N$, in the degree $d=\max(\deg(A),\deg(B))$ of ${\mathcal{D}}$ and in the logarithm of the height of $A$ and $B$. A nontrivial modification of this algorithm provides a polynomial-time method to solve $({\mathcal{P}}_{N})$, see again [Chè11]. To our knowledge, this is the first algorithm solving $({\mathcal{P}}_{N})$ in polynomial-time. Unfortunately, the exponent is quite large, making the algorithm unpractical even for moderate values of $N$. This drawback is due to the fact that algorithm Rat-First-Int in [Chè11] needs to compute the irreducible factors of a bivariate polynomial of degree $\approx dN^{4}$, and the best known algorithms for solving this subtask have arithmetic complexity, e.g., $\tilde{\mathcal{O}}(d^{\omega+1}N^{4\omega+4})$, see [BLS⁺04, Lec06].

Last, we mention the article [FG10] of Ferragut and Giacomini, where the algebraicity of a generic power series solution of the differential equation $\eqref{diff-eq-intro}$ is used to improve the efficiency of the naive algorithm. Precisely, the system of quadratic equations yielding the coefficients of a rational first integral is reduced to a simpler system of (still quadratic) equations. Although this gives a good heuristic improvement on the naive method, we show in the present article how to turn it into a fast algorithm. Starting from the link between ($\mathsf{S}$) and $\eqref{diff-eq-intro}$, we reduce $({\mathcal{P}}_{N})$ to solving a system of linear equations. Furthermore, we give tight bound on the number of terms of power series solutions of $\eqref{diff-eq-intro}$ that are sufficient to detect the existence of rational first integrals and to compute one of them when it exists. This enables us to turn the heuristic in [FG10] into an algorithm with polynomial complexity, that is more efficient both in theory and in practice than all previous algorithms, see Section 6.

With the planar polynomial vector field ($\mathsf{S}$), we associate the first order non-linear differential equation:

We may assume without any loss of generality that $x$ does not divide $A$, i.e., $A(0,y)\not\equiv 0$. We will explain how we can reduce to this situation and study the complexity of this reduction in Subsection 5.1.

Then, the formal version of the Cauchy-Lipschitz theorem for non-linear (first-order) differential equations ensures that for any $c\in{\mathbb{K}}$ such that $A(0,c)\neq 0$, the equation ($\mathsf{E}$) admits a unique power series solution $y_{c}(x)\in{\mathbb{K}}[[x]]$ satisfying $y_{c}(0)=c$. Note that high-order truncations of the power series $y_{c}(x)$ can be computed efficiently using the algorithm of Brent and Kung  [BK78].

The following standard result is fundamental to both our method and the one of Ferragut and Giacomini in [FG10].

Proposition 16 shows in particular that if ($\mathsf{S}$) has a rational first integral $P/Q$, then all power series solutions of ($\mathsf{E}$) are algebraic. The next proposition which is well known (see [FG10, Wei95]) asserts that the converse is also true.

We have seen in Proposition 16 that if $P/Q$ is a reduced non-composite rational first integral of ($\mathsf{S}$), then a minimal polynomial of a power series solution of ($\mathsf{E}$) is generically of the form $\lambda\,P-\mu\,Q$ for some constants $\lambda$ and $\mu$. Thus, if we are able to compute such a minimal polynomial, we can deduce the rational first integral $P/Q$. In practice, we do not compute a power series $y_{c}(x)\in{\mathbb{K}}[[x]]$ solution of ($\mathsf{E}$) but only a truncation of $y_{c}(x)$, i.e., a finite number of terms of its expansion on the monomial basis. Given a degree bound $N$ for the rational first integral that we are searching for, the following lemma shows that computing $y_{c}(x)\mod x^{N^{2}+1}$, i.e., the first $N^{2}+1$ terms of its expansion, is enough for our purposes.

Such an analysis of the needed precision for the power series solutions of ($\mathsf{E}$) that we compute is not included in [FG10]. Note that this kind of strategy was already used in a polynomial factorization setting (see, e.g., [Kal85]) and in a differential equations setting (see [ACFG05, BCS⁺07]). The next lemma is a small improvement of [ACFG05, Lemma 2.4].

Note that in $(\star)$, the power series $\hat{y}(x)$ can be replaced by its truncation $\hat{y}(x)\mod x^{N^{2}+1}$.

In the sequel, we say that $\tilde{M}$ is a minimal solution of $(\star)$ if it is a solution of $(\star)$ with minimal degree in $y$.

We now propose a first algorithm, based on linear algebra, for solving $({\mathcal{P}}_{N})$. More efficient algorithms, based on this one, are given in Section 4. The strategy of this algorithm is then used in Section 7 for computing Darboux polynomials.

1. (1)

   For an indeterminate $c$, compute the polynomial $y_{c}\in{\mathbb{K}}(c)[x]$ of degree at most $(N^{2}+1)$ s.t. $y_{c}(0)=c$ and $\frac{dy_{c}}{dx}\equiv\frac{B(x,y_{c})}{A(x,y_{c})}\mod x^{N^{2}+1}$.
   

2. (2)

   Compute all¹¹1i.e., a basis over ${\mathbb{K}}(c)$, see Remark 19. non-trivial polynomials $\tilde{M}\in{\mathbb{K}}(c)[x,y]$ of degree $\leq N$ s.t.

   $$(\star):\tilde{M}(c,x,y_{c}(x))\equiv 0\mod x^{N^{2}+1}.$$

   If no such $\tilde{M}$ exists, then Return “None”. Else, among the solutions of $(\star)$, pick a minimal solution $\overline{M}\in{\mathbb{K}}[c][x,y]$.
   

3. (3)

   Let $M$ denote the primitive part of $\overline{M}$ relatively to $y$.
   Set $P(x,y):=M(0,x,y)$.
   Pick any $c_{1}\in{\mathbb{K}}$ s.t. $\frac{M(c_{1},x,y)}{P(x,y)}\not\in{\mathbb{K}}$ and set $Q(x,y):=M(c_{1},x,y)$.
   

4. (4)

   If ${\mathcal{D}}(P/Q)=0$, then Return $P/Q$. Else Return “None”.