Finding Liouvillian first integrals of rational ODEs of any order in finite terms

It is known, due to Mordukhai-Boltovski, Ritt, Prelle, Singer, Christopher and others, that if a given rational ODE has a Liouvillian first integral then the corresponding integrating factor of the ODE must be of a very special form of a product of powers and exponents of irreducible polynomials. These results lead to a partial algorithm for finding Liouvillian first integrals. However, there are two main complications on the way to obtaining polynomials in the integrating factor form. First of all, one has to find an upper bound for the degrees of the polynomials in the product above, an unsolved problem, and then the set of coefficients for each of the polynomials by the computationally-intensive method of undetermined parameters. As a result, this approach was implemented in CAS only for first and relatively simple second order ODEs. We propose an algebraic method for finding polynomials of the integrating factors for rational ODEs of any order, based on examination of the resultants of the polynomials in the numerator and the denominator of the right-hand side of such equation. If both the numerator and the denominator of the right-hand side of such ODE are not constants, the method can determine in finite terms an explicit expression of an integrating factor if the ODE permits integrating factors of the above mentioned form and then the Liouvillian first integral. The tests of this procedure based on the proposed method, implemented in Maple in the case of rational integrating factors, confirm the consistence and efficiency of the method.

Then (2) becomes where We conclude that the integrating factor of ODE (1) is µ = ∂ζ ∂y n−1 , where the function ζ is a first integral, a solution of the linear first-order PDE (4).
If we eliminate the function ζ from the system (3), (4), we obtain a PDE system for the integrating factor. One of the equations of the PDE system obtained by double differentiation of (4) with respect to y n−1 is as follows (n ≥ 2): Assuming that an integrating factor is known, it is well-known that we are able to obtain corresponding first integral by quadratures. If we know n independent first integrals for ODE (1), we can then find its general solution, at least in implicit form. As it is seen from (3), the algebraic structure of the integrating factor is simpler than the structure of the first integral, so in some cases finding the integrating factors is an easier problem.
The most beautiful way of finding integrating factors is the Darboux method and its refinements. If f in (1) is a rational function with respect to all variables x, y 0 , y 1 , . . . , y n−1 , we will call such an ODE a rational one. In 1878, Darboux [2] proved that if a rational ODE has at least m(m + 1)/2, where m = max(order(numerator(f ), order(denominator(f )) invariant algebraic curves P i , now known as Darboux polynomials, then it has a first integral or an integrating factor of the form i P a i i for suitable complex constants a i . So, in principle, the problem of finding an integrating factor is reduced here to the problem of finding Darboux polynomials.
In 1983, Prelle and Singer [3] proved the essential theorem that all first order ODEs, which possess elementary first integrals, have an integrating factor of the above mentioned form i P a i i . Their result unifies and generalizes a number of results originally due to Mordukhai-Boltovski [4], Ritt [5] and others. Singer in [6] refines this result for ODEs, which possess Liouvillian first integral.
So, to find an integrating factor to ODE (6) in the conditions of Theorem 1, we first should obtain the sets of polynomials P i and Q j , and if we succeed in this task then we can obtain the constants a i 's and b j 's from the full PDE system for the integrating factor.
To consider the ways to find the sets of polynomials P i and Q j , let us substitute the integrating factor in the form (7) into (5). After some rearrangement, we arrive at where L 1 is a polynomial. As A and B are relatively prime polynomials, we conclude that B divides one of the polynomials P i or Q j . Suppose that A = A 1 (y n−1 )A 2 (x, y 0 , y 1 , . . . , y n−1 ), . . , y n−2 )B 3 (x, y 0 , y 1 , . . . , y n−1 ) and Substituting these expressions for A, B and µ into (5), we can conclude, by considering divisibility of polynomials, that it is necessary that ε 3 = 1 and then such polynomials are factors of A or B, and they can be easily selected and included into the set of candidates for the integrating factor structure.
Without loss of generality, we will consider an integrating factor of the form and will focus our attention on the cases when ∂P i ∂y n−1 = 0 and D(P i ) = 0 ( ∂Q j ∂y n−1 = 0 and D(Q j ) = 0).
For the case of a k = 1, substitution of (8) into (5) leads to and similarly for the case of b k < 0 where L 2 and L 3 are polynomials. All components involved in (9) and (10) are polynomials. We assume that P i and Q j are irreducible, relatively prime polynomials, so gcd P i , ∂P i ∂y n−1 and gcd Q i , ∂Q i ∂y n−1 are constants, where greatest common divisor is denoted by gcd. We can conclude that (f | g means that f divides g) and The basic significance of relations (11) and (12) is that we can investigate each of the polynomials P i or Q j independently. Properties (11) and (12) constitute the basis of the so-called Prelle-Singer method, which originates from Darboux [2], to obtain the integrating factors and first integrals of ODEs of type (6). Supposing that and substituting P k into (11), we can in principle obtain c kij 0 ···j n−1 's and as a result P k (or similarly Q k ) in explicit form.
To do so, we first have to establish an upper bound N for the degrees of the polynomial P i (or Q k ). It is known that N exists, but there is no constructive method to evaluate it so far. It is the principal weakness of the Prelle-Singer method. In the existing implementations of the Prelle-Singer method, an upper bound N is predefined by a user of the method.
With an increase of equation order the number of undetermined variables c kij 0 ···j n−1 and, correspondingly, the dimension of the system of algebraic equations for these constants increases. As a result, solving of this system becomes problematical. This approach was implemented in CAS only for first and relatively simple second order ODEs [8,9,10,11,12,13,14,15], and as a rule the method can not be used for N > 4. In contrast, the modified method proposed in the next section does not suffer from these restrictions.
We must also mention that a rational ODE does not always have a local Liouvillian first integral. While generalizations of the Darboux integrability theory allow in some ("integrable") cases to find non-Liouvillian first integrals (see, e.g. [16,17,18,19,20]), we consider only Liouvillian first integral theory.

The method of resultants
Let us investigate closely the relation (11) or equivalently (12). Let R z (f, g) be the resultant of polynomials f and g with respect to the indeterminate z. In the sequel, we will consider only the case when both A and B are not constants and there are such z ∈ (x, y 0 , y 1 , . . . , y n−1 ) that R z (A, B) is not a constant.
The following Theorem 2 and Corollary are our main results: Theorem 2. If a polynomial P i satisfies the condition (11) then for any indeterminate z ∈ (x, y 0 , y 1 , . . . , y n−1 ) the resultant R z (P i , B) must divide the resultant R z (A, B).
Proof . We may write P i L = BD(P i ) + A ∂P i ∂y n−1 for some polynomial L. Taking resultants of both sides of this equation with B and with respect to an indeterminate z ∈ (x, y 0 , y 1 , . . . , y n−1 ), we get the following polynomial equation As P i and ∂P i ∂y n−1 have no common roots, R z (P i , B) does not divide R z ∂P i ∂y n−1 , B . Therefore, and it is obvious that we can choose R z (P i , B) as a divisor of R z (A, B) up to a constant factor.
How can we recover P i if we know the resultant R z (P i , B)? First of all, we note that the case when A/B depends on only one indeterminate is trivial. As we have mentioned above, if P i depends only on one indeterminate, then such polynomials are factors of A or B and we have already included them into the set of candidates. So, we assume that P i depends at least on two indeterminates and there exists such z that R z (P i , B) is a polynomial. We do not insist that R z (P i , B) is irreducible.
It is known from elementary facts about resultants that there exist such polynomials α i , β i ∈ K[x, y 0 , y 1 , . . . , y n−1 ] that β i P i = α i B + R z (P i , B).
As a corollary of Theorem 2, we have the following result. Assuming that P i 's (and Q j 's) are irreducible polynomials we conclude that Corollary. If polynomials P i (and respectively Q j ) satisfy (14) then the following holds where P i hyp = β i P i , β i is some polynomial, α i ∈ (−1, 0, 1) and c i are constants.
Proof . If α i = 0 identically, then for any z ∈ (x, y 0 , y 1 , . . . , y n−1 ) from where γ ∈ (−1, 1), we conclude that R z (α i , P i ) = 1/γ. So as P i depends at least on two indeterminates then α i must be a constant and α i ∈ (−1, 0, 1). Now we are able to modify the Prelle-Singer method by replacing the hypothesis (13) with (15) and finding constants α i and c i from the following equation Then P i is a irreducible factor of P i hyp which satisfies (11). Thus, we can use algebraic, not differential, relations for obtaining polynomials as in the original Darboux (and Prelle-Singer) approach. First of all, we bypass the principal weakness of the Prelle-Singer method as we do not need to establish an upper bound N for the degrees of polynomial P i here. Second, we have the advantage that finding only two undetermined constants α i and c i in the P i hyp structure is much easier than finding the set of c kij 0 ···j n−1 's.
By substitution of the irreducible factors of P 1 hyp and P 2 hyp to (11), we select the following candidates for the integrating factor structure: P 1 = (y 2 1 − 2y 1 − y 0 ) and P 2 = (y 1 + y 0 + x). So the hypothesis of µ is After substitution of this µ hyp into the PDE system for the integrating factor, we obtain that X2 = −(X1 + 2) so an integrating factor of given ODE is where X1 is an arbitrary constant.
Examinations of the cases when a k = 1 and b k > 0, are more difficult. We are not able to give here the full analysis of the problem which necessitates consideration of many subcases. For brevity, we will only demonstrate that the method of resultants in principle enables us to obtain all the needed polynomials in a finite number of steps. The main tool here is usage of repeated resultants.
For compactness, we will use the following notation for the repeated resultants g, h)).
For the case a k = 1, we obtain the following equation similar to (9): where L 4 , M 1 , M 2 are some polynomials, andP i (i = k),Q j are only such polynomials from the integrating factor structure which do not satisfy conditions (11) and (12). So for this case we can obtain the following properties ( ∂P k ∂y n−1 = 0): and further From (16) we can obtain the hypothesis for P k and then with the help of (17) all its undetermined constants.
If there is a polynomial P s in the integrating factor structure which satisfy (11), then finding P k can be simplified. Since Finding Liouvillian First Integrals of Rational ODEs of Any Order in Finite Terms Among other things we can observe from (16) and (18) is that it is very likely that the P k 's satisfy the relation (11) even if a k = 1. Also, there are relations for Q k with b k > 0, similar to (16)- (19).
The main elements of the method described above were implemented in demonstration procedures in Maple [21,22] for the case of rational integrating factors.
We have tested routines in the following area: for a given ODE order we assigned the first integral ζ of the following type are some constants. Such algebraic first integral leads to an ODE of type (1) with f = −D(ζ)/ ∂ζ k ∂y n−1 , which is potentially tractable by the proposed procedure as this ODE has a rational integrating factor µ = ∂ζ k ∂y n−1 = i P a i i . The first stage -finding candidates and hypothesis of the structure of an integrating factor -is very fast. The second stage is relatively resource consuming -finding the unknown parameters a i of the structure of an integrating factor.
We can conclude that our prototypes are workable for rational ODEs even with symbolic constant parameters of orders from n = 1 to n = 4 − 5. Although in principle the method provides ability to find an integrating factor for almost any of the above mentioned type of ODEs, in practice the procedure heavily relies on some basic CAS functions such as factor, simplify and so on, and in some cases fails. In addition, different outputs are observed in different Maple versions for the same ODE.
As a rule, Maple cannot find any integrating factor for ODEs from the considered test area while proposed procedure is usually successful. Often our procedure produces more than one integrating factor. For example, the output may contain a set of µ's or µ with arbitrary parameters, and sometimes they lead to independent first integrals.

Conclusion
We have proposed an algebraic method for finding polynomials P i and Q j of the integrating factor form (7) of rational ODE (6) for any order using resultants of the numerator and denominator of the ODE's right-hand-side.
The main advantage of the method of resultants is ability to obtain the polynomials P i and Q j without an a priori assumption of the upper bound for their respective degrees. Our method bypasses the principal weakness of the Prelle-Singer (Darboux) method.
Moreover, the number of undetermined variables involved in calculations by our method is radically smaller than in the Prelle-Singer (Darboux) method. This property allows us to obtain polynomials P i and Q j of arbitrary degrees without any complications whereas the Prelle-Singer (Darboux) method can not be useful as a rule for N > 4.