A solution algorithm for the Abel Inverse Riccati ODE

Austin D. Roche The Centre for Experimental and Constructive Mathematics Simon Fraser University, Burnaby, BC, Canada

Abstract

A solution is presented for the three-parameter class of ordinary differential equation called Abel Inverse Riccati. To match this general family, a new theory for equivalence of Abel equations is developed, based on the idea of finding rational equivalence transformations. The concepts of minimal representative equations and invariant component polynomials for a class are introduced. A set of algorithms is then described which can solve the general class of Abel Inverse Riccati equations.

1 Introduction

1.1 The equivalence problem for Abel equations

The ordinary differential equation attributed to Abel following his pioneering work in [1] takes many forms, the most general being known as Abel’s equation of second kind:

(g_{1}y + g_{0})y′ = f_{3}y^{3} + f_{2}y^{2} + f_{1}y + f_{0}, (1)

where y ≡ y(x), and g_{0}, g_{1}, f_{0}, f_{1}, f_{2}, and f_{3} are analytic functions of x. Abel himself considered the special case where f_{3} = 0; the further specialization g_{1} = 0 yields a Riccati equation, which itself is equivalent to a linear equation of second order. However, it appears that there is no way to linearize the general Abel equation.

As with most ordinary differential equations, one of the best solution techniques for Abel equations is to change variables, looking for an equation of simpler type. Any solution to the latter will generate a solution to the former by the reverse transformation. Accordingly, we start by considering those transformations which do not make the equation’s type more complicated. The most general point transformation which preserves the form of Eq.(1) is

(x = F (x̃), y = \frac{P (x̃)ỹ + Q(x̃)}{R(x̃)ỹ + S(x̃)}) , (2)

where x̃ and ỹ are the new independent and dependent variables, and F, P, Q, R, S are arbitrary functions of x̃ satisfying F′ ≠ 0 and P S ≠ QR. Two Abel equations are said to be equivalent if one can be obtained from the other by means of such a transformation, which are consequently known as equivalence transformations. The set of all equations equivalent to one given is its equivalence class.

A comprehensive long term solving strategy involves two phases:

1. Determination of solvable classes. Generating of a database S of Abel classes, each represented by an equation with known solution;

The first step means not only generating solvable equations [1–5,9,11], but also collecting and classifying them and removing redundancies [4–6,8,11–14].

2. The equivalence problem. Developing an algorithm to match a given input equation with its equivalent, if it exists, among the solvable set S, and to calculate the associated equivalence transformation.

The main algorithm will involve many sub-algorithms, each specialized to the nature of the class it is designed to match. Because there are infinitely many known solvable classes, it is necessary to collect related ones into parameterized super-classes, which themselves may compose S. By removing parameters through transformation, a representative for each super-class should be found which carries the minimal number of parameters, essentially the dimension of the super-class. The matching algorithms for parameterized classes must find not only the equivalence transformations but also the values of the parameters identifying the given class. Concerning the invariantbased algorithm for matching zero-parameter classes presented by Liouville in [11], a computer implementation was developed and presented, along with an extension to one-parameter classes, in [5]. A different algorithm, based on finding symmetries, was developed for the two-parameter Abel Inverse Linear (AIL) class described in [6]. However, these algorithms were do not generalize to provide a practical method for the case of a multi-parameter class, and there has been no attempt to provide a general matching strategy for classes of Abel equations which depend on multiple parameters.

A more general three parameter super-class known as Abel Inverse Riccati (AIR) was also discussed in [6]. Together with AIL these classes were found to contain as special cases almost all the solvable Abel equations found in the literature. Unfortunately, there was no solution known for the full class.

In this work I will present both a new solution for the complete three-parameter AIR super-class, as well as a set of practical algorithms for solving all member equations of the class. New routines implementing these algorithms have been implemented in Maple. The implementation is not yet complete, and will not be described in this work, except to say that it represents a significant improvement in speed over the aforementioned routines described in [5].

In the following section we will review and summarize the existing theory on classifying and solving Abel equations, including canonical forms, invariants, and the standard equivalence method due to Liouville. A new and more useful set of differential invariants will be described. We will then in Section 3 describe a new theory for solving the equivalence problem, based on the idea of rational function decomposition. The concept of a minimal representative for an Abel class will be seen as the key, allowing us to restrict our search to equivalence functions rational in the independent variable. In the fourth section we will introduce and show how to compute invariant polynomials, the polynomial components of the rational invariant functions which characterize any class. This leads us to the description of a new general method for Abel equations. Finally in the last section we apply these ideas to the problem of the Abel Inverse Riccati class of equations, developing a set of minimal equations and corresponding solution algorithms that together represent a new and unique solving strategy for the complete class.

2 Canonical forms, invariants, and Liouville’s equivalence method

Classifying an Abel equation involves generating its canonical form, calculating a set of invariants defined for that form, and then comparing the relations among those invariants with those of known classes. A preliminary form for Eq.(1) is the Abel equation of first kind,

y′ = f_{3}y^{3} + f_{2}y^{2} + f_{1}y + f_{0}, f_{3} ≠ 0, (3)

obtained using the transformation^{1} y = 1/ỹ−g_{0}/g_{1} when g_{1} ≠ 0. Starting from Eq.(3), a further simplification is obtained through y = ỹ− \frac{1}{3}f_{2}/f_{3}, which makes the coefficient of ỹ^{2} vanish. Finally, assuming f_{2} = 0 (and f_{0} ≠ 0), the transformation y = ỹf_{0} yields rational normal form:^{2}

Ω_{RNF} : y′ = c_{3}(x)y^{3} + c_{1}(x)y + 1, c_{3} ≠ 0. (4)

Other canonical forms, notably those suggested respectively by Appell [4] and Liouville,^{3}

Ω_{A} : Y′(x) = (Y (x)^{3} + J(x))X′(x), X′ ≠ 0, (5)

1 As a convention, we shall use y and x for the dependent and independent variables of an equation respectively, ỹ and x̃ for the new variables defined in a transformation. After the transformation, we will drop the tildes in order to keep the notation simple. 2 This form was suggested in [15]. If f_{2} = f_{0} = 0, Eq.(3) is of Bernoulli type and readily solvable. 3 The equation which appears in [11], \frac{dz_{1}}{dt} + \frac{1}{T} (z_{1}^{3} + \frac{1−T}{3t} z_{1} + \frac{1}{t}) = 0, is in error.

Ω_{L} : z′_{1}(x) + \frac{t′(x)}{T (x)} (−\frac{z_{1}(x)^{3}}{t(x)^{2}} + \frac{1 − T (x)}{3t(x)}z_{1}(x) + \frac{1}{t(x)}) = 0, t′, T ≠ 0. (6)

are also obtained from Eq.(3) by changing the dependent variable: y = P (x)ỹ + Q(x).

The group of transformations which preserves each form is known as its structure invariance group. For the most general Abel equation Eq.(1) this group consists of all the equivalence transformations of the form Eq.(2); those for the canonical forms are more restricted. Structure invariance transformations for Abel equations of first kind may depend only on F, P and Q, with R = 0 and S = 1; those for the rational normal form, (x = F (x̃), y(x) = F′(x̃)ỹ), and those for the canonical forms of Appell and Liouville (x = F (x̃), y(x) = ỹ), depend only on F. Assuming that both input and representative equations take one of these canonical forms, the equivalence problem for Abel equations reduces to determining this function F. We shall call F the equivalence function for the transformation Eq.(2).

Generally, a function J taking ODEs of a given form F = {Ω_{i}(x, y, y′, ...) = 0} to functions of the variables x, y, y′, ... is an invariant if

J (g(Ω)) = φ_{J} (g) g(J (Ω)) (7)

for some multiplier φ_{J}, all equations Ω ∈ F, and all transformations g in the structure invariance group of F. The invariant J is absolute when φ_{J} is identically 1, relative otherwise. Related invariants J_{i} have weights i respectively if there is a function φ̃ satisfying φ_{J_{i}} ≡ φ̃^{i} for all i. See for example the development by Halphen [7].

The functions J, X′, t, T of Eqs. (5,6) are examples of absolute invariants. Whereas J and X′ can be computed from Eq.(1) only by taking roots of integrals, t and T can be computed with just derivatives. We will review their derivation and show that they they form part of a sequence.

With repect to Abel equations of first kind Eq.(3), Liouville defined the relative invariants s_{2n+3} of weight 2n + 3 and multiplier φ̃ = F′P:

s_{3} ≡ f_{0}f_{3}^{2} + \frac{2f_{2}^{3}}{27} + \frac{f_{3}f_{2}′ − f_{2}f_{3}′ − f_{1}f_{2}f_{3}}{3} , s_{n+2} ≡ f_{3}s_{n}′ − ns_{n} (f_{3}′ + f_{1}f_{3} − \frac{f_{2}^{2}}{3} ) .

For equations in rational normal form, the corresponding invariants have multiplier φ̃ = F′^{2} and can be written simply^{4}

s_{3} = c_{3}^{2}, s_{2n+3} = D̃^{n}s_{3}, (8)

where the differential operator D̃ defined by

D̃(z_{n}) ≡ c_{3}z_{n}′ − nz_{n} (c_{3}′ + c_{1}c_{3})

takes invariants z_{n} of weight n to those of weight n + 2. An especially convenient sequence of absolute invariants is:

L_{1} ≡ \begin{matrix}s_{5} & 3 \\ s_{3} & 5 \end{matrix} , L_{2} ≡ \frac{DL_{1}}{L_{1}} = \frac{t_{10}}{s_{5}^{2}} , L_{3} ≡ \frac{Dt_{10}}{t_{10}} = \frac{s_{3}t_{12}}{s_{5}t_{10}} , ... L_{n} ≡ \frac{Dt_{2n+4}}{t_{2n+4}} . (9)

where t_{10} ≡ 3s_{3}s_{7} − 5s_{5}^{2}, t_{2n+10} ≡ D̃^{n}t_{10}, and^{5}

D ≡ \frac{s_{3}}{s_{5}} D̃. (10)

L_{1}, L_{2} are Liouville’s invariants t, T respectively of Eq.(6). Note that the operator D preserves the weight of an invariant. We will find this operator very useful below.

Integration strategies for the case where L_{1} is constant were first provided by Abel [1] and Liouville [9]. Appell showed in [3] a simple equivalence transformation bringing the equation into a form with constant coefficients.

When L_{1} is not constant the equivalence problem is more complicated. The equivalence method suggested by Liouville [10] for matching a given Abel equation Ω̃(x̃, ỹ, ỹ′) = 0 with the representative Ω(x, y, y′) = 0 of a known class requires equating their respective calculated invariants:

L̃_{1}(x̃) = L_{1}(x), L̃_{2}(x̃) = L_{2}(x), ... (11)

4 With respect to rational normal form equations, a sequence of relative invariants r_{n} with lesser weights n and a simpler multiplier φ̃ = F′ can be defined via

r_{3} ≡ c_{3}, r_{n+1} ≡ c_{3}^{n−1}D̃(r_{n}c_{3}^{−n}) = n c_{1}r_{n} + \frac{d}{dx}r_{n}, n ≥ 3.

However, these do not seem to lead to simpler absolute invariants in practice. This is probably because the solvable equations under consideration are of the AIA class and so the coefficients f_{i} of Eq.(3) are simpler than the c_{i} of Eq.(4). In the AIR case for example the f_{i} and hence the s_{i} will be polynomial. See [6]. 5 Using the notation in [4], we have D = \frac{J}{J′(X)} \frac{d}{dX} = \frac{d}{d ln(J)}.

and searching for a common solution x = F (x̃). Because the rational normal form of an equation can be expressed in terms of the first two Liouville invariants (cf. Eq.(6)), a sufficient condition for equivalence is that there be a solution to the first two equations above. When there are no parameters, two such equations will also suffice to determine a set of candidates for F; another equation must be used to eliminate each extra parameter characterizing a super-class.

As we are about to see, the typical and relevant case is when both the representative and input equations have rational invariants. If so, a GCD calculation will determine a set of candidates for F; resultant calculations are required in advance to remove any extra parameters. A computer implementation of this method was developed as part of the Maple system and presented in [5].

The complexity of these operations depends on the degrees of both sets of invariants.^{6} Consequently it is important to determine a set of invariants which minimize the degrees of the invariants for typical equations; the invariants defined in Eq.(9) seem to be optimal in this regard. It is also important to determine representative equations for each class for which the associated invariants are minimal in degree.

Finally, the complexity of this approach is also roughly exponential in the number of parameters. In practice, the algorithm present in Maple is capable of matching one-parameter classes; analogous algorithms for multi-parameter classes are still impractical. We will develop an alternative approach to overcome this problem below. In preparation for this new algorithm, we must show that there is no loss in generality in assuming that the equivalence function F is rational.

3 Class representatives and rational equivalence transformations

It happens that there are only a few classes among those with non-constant invariant L_{1} for which solutions are known. As shown in [6], the integrable Abel classes found in the literature all have representatives from a very specific family: those obtained from linear, Riccati, or other Abel equations under the inverse transformation (x = ỹ, y = x̃), which switches the roles of the dependent and independent variables. In particular, each of these classes has a representative equation with rational coefficients.

Let us call an Abel equation rational when it has coefficients rational in the independent variable, and an Abel class rational when it contains a rational representative. Each rational class admits a certain hierarchical structure related to the complexity of the invariants of its member equations. In particular, there are some members within each class which can be considered distinguished.

6 By the degree (or rational degree) of a rational function we mean the maximum of the degrees of its numerator and denominator.

Definition 3.1 A rational Abel equation with non-constant invariant L_{1} is a minimal representative for its class if all equivalence transformations Eq.(2) which bring it to another rational equation involve a rational equivalence function.

Theorem 3.1 Every rational Abel class C with non-constant invariant L_{1} has a minimal representative. Suppose Ω(x, y) ∈ C is rational with Liouville invariants L_{1}(x), L_{2}(x). The following conditions are equivalent:

(i) Ω is a minimal representative for C; (ii) among rational equations in C, the degrees of L_{1} and L_{2} are minimal for Ω; (iii) the system

L_{1}(x) = L_{1}(t),  L_{2}(x) = L_{2}(t)  (12)

admits the unique solution x = t; (iv) the only equivalence transformation which leaves Ω unchanged is the identity.

We will first prove the implications (iii) ⇒ (i), (ii), then use this result to construct a minimal representative Ω̃ for C, and finish with the reverse implications and the equivalence (iii) ⇔ (iv). Let us begin with a lemma.

Lemma 3.2. If the quotient of two polynomials A and B is not constant, then x − t divides A(x)B(t) − A(t)B(x) once only.

Proof: Let A(x) = ∑_{j=0}^{d} a_{j}x^{j}, B(x) = ∑_{j=0}^{d} b_{j}x^{j} ≠ 0, and consider the polynomial

R(x, t) ≡ \frac{A(x)B(t) − A(t)B(x)}{x − t}.

By solving sequentially the coefficients of R(x, x), an inductive argument shows that if R(x, x) = 0 then a_{j}b_{k} = b_{j}a_{k} for all j, k; hence A/B is constant.

Proof of Theorem 3.1: Let L_{i} = N_{i}/M_{i} for relatively prime polynomials M_{i}, N_{i} ∈ K[x] for some field K and d_{i} be the rational degree of L_{i}, i = 1, 2.

(iii) ⇒ (i), (ii): Suppose Ω̃ ∈ C has rational coefficients and is obtained from Ω under the transformation Eq.(2) with equivalence function x = t(x̃). By construction the invariants L̃_{i}(x̃) ≡ L_{i}(t(x̃)) are rational in x̃; thus we can write L_{i}(t(x̃)) = Ñ_{i}(x̃)/M̃_{i}(x̃) for relatively prime polynomials Ñ_{i}, M̃_{i}, i = 1, 2.

The solution to the system (12) may be obtained by taking the greatest common divisor of the two polynomials N_{i}(x)M̃_{i}(x̃) − M_{i}(x)Ñ_{i}(x̃), i = 1, 2. Since this system has a unique solution for x, the GCD must be linear in x; by construction it will be polynomial in x̃. Solving, we find that x = t(x̃) is the quotient of two polynomials in x̃, of rational degree d_{t}. Hence Ω is minimal, and the degree of L̃_{i} is d_{i}d_{t} ≥ d_{i}, i = 1, 2.

Existence of a minimal equation Ω̃: Suppose Eq.(12) has m distinct solutions x = X_{j}(t), j = 0 … m − 1. Without loss of generality we may assume that X_{0}(t) = t. Let Q(x, t) be the monic GCD in K(t)[x] of the two polynomials N_{i}(x) − M_{i}(x)L_{i}(t), i = 1, 2. Any solution x = X(t) of Q(x, t) = 0 must satisfy L_{i}(t) = L_{i}(X(t)) for each i = 1, 2. Thus X = X_{j} for some j. Hence Q(x, t) is a product of factors of the form x − X_{j}(t). Each such factor must have the same multiplicity, since Q(x, t) by construction is unchanged under each substitution t → X_{j}(t). Note however, that Q(x, t) is also the monic GCD in K(t)[x] of the polynomials N_{i}(x)M_{i}(t) − N_{i}(t)M_{i}(x), i = 1, 2, at least the first of which is divisible by x − t exactly once, by the lemma. Therefore Q(x, t) is squarefree, and more precisely,

Q(x, t) = ∏_{j=1}^{m} x − X_{j}(t).

Finally, let P (x, t) be the GCD in K[x, t] of the polynomials N_{i}(x)M_{i}(t)−N_{i}(t)M_{i}(x),i = 1, 2, and V (t) be the leading coefficient of P (x, t) with respect to x. P (x, t) must be Q(x, t)V (t), so that P has degree m in x, and by symmetry, in t as well. We may write P (x, t) = ∑_{j=0}^{m} V_{j}(t)x^{j}, where V_{m} = V, and each V_{j} has degree at most m. The coefficients V_{j}(t)/V (t) of Q have rational degree at most m and are unchanged under each of the m substitutions t → X_{j}(t). Hence they must be either constant or have degree exactly m. Since P (x, x) = 0, P (x, t) cannot be separable in x and t. Thus at least one of these coefficients V_{k}(t)/V (t) must depend on t; as a result, V_{k} and V must be relatively prime. For each j we may choose A_{j} so that V_{j}(t) + A_{j}V_{k}(t) has a factor in common with V (t). For example, if V (r) = 0, A_{j} ≡ −V_{j}(r)/V_{k}(r) suffices. The function C_{j} ≡ (V_{j}(t) + A_{j}V_{k}(t))/V (t) must have rational degree less than m, and as it is invariant under each of the m substitutions t → X_{j}(t), must be constant. Let A(x) = ∑_{j=0}^{m} A_{j}x^{j}, C(x) = ∑_{j=0}^{m} C_{j}x^{j}. By eliminating V_{j} we obtain

\frac{P (x, t)}{V (x)V (t)} = \frac{C(x)}{V (x)} − \frac{V_{k}(t)}{V (t)} \frac{A(x)}{V (x)}.

Antisymmetry of the mixed partial derivative \frac{∂^{2}}{∂x∂t} \frac{P (x,t)}{V (x)V (t)} implies that

\frac{d}{dx} (\frac{A(x)}{V (x)}) /\frac{d}{dx} (\frac{V_{k}(x)}{V (x)}) = −\frac{d}{dt} (\frac{A(t)}{V (t)}) /\frac{d}{dt} (\frac{V_{k}(t)}{V (t)}) .

As the two sides are equal they depend neither on x nor t; as they are antisymmetric in x and t, they must be 0. Hence α ≡ A/V must be constant, nonzero because A_{k} is -1. Using P (x, x) = 0 we determine that C = αV_{k}:

\frac{P (x, t)}{V (x)V (t)} = α (\frac{V_{k}(x)}{V (x)} − \frac{V_{k}(t)}{V (t)}) .

The system Eq.(12) has a solution precisely when V_{k}(x)/V (x) = V_{k}(t)/V (t). Indeed, we will show that L_{i} = L̃_{i} ∘ (V_{k}/V), i = 1, 2 for some rational functions L̃_{i}.

Let x = F (x̃) be a root of the equation V_{k}(x)/V (x) = x̃. Note that V_{k}/V = C/A and because V_{k} and V are relatively prime, it must be that A_{m} = 0, C_{m} = 1. Thus we may apply the substitution x^{m} = ∑_{j=0}^{m−1}(A_{j}x̃ − C_{j})x^{j} repeatedly to simplify the expression L_{i}(x), yielding an expression rational in x and x̃ and of rational degree at most m − 1 in x. As this expression is unchanged under each of the m substitutions x → X_{j}(x), it cannot depend on x. It is therefore L_{i}(F (x̃)).

Consider the equation Ω̃(x̃, ỹ, ỹ′) obtained from the Liouville canonical form of Ω under the transformation x = F (x̃). Its Liouville invariants are L̃_{i}(x̃) ≡ L_{i}(F (x̃)), which we have just shown to be rational in x̃. The coefficients of Ω̃ must therefore also be rational by Eq.(6). Now let t̃ = V_{k}(t)/V (t). The system Eq.(12) can therefore also be written L̃_{i}(x̃) = L̃_{i}(t̃), i = 1, 2, and we have seen it to be solved only when x̃ = t̃. Hence Ω̃ is minimal.

(i) ⇒ (iii): If Ω is minimal then since Ω̃ is rational, F (x̃) must be rational, of degree n say. Since the composition V_{k}(F (x̃))/V (F (x̃)) = x̃ has rational degree mn = 1, we must have n = m = 1; in other words, the system Eq.(12) has the unique solution x = t.

(ii) ⇒ (iii): Considering the string of equalities,

L_{i}(x) = L_{i}(F (x̃)) = L̃_{i}(x̃) = L̃_{i}(V_{k}(x)/V (x)).

we see that L_{i} = L̃_{i} ∘ (V_{k}/V). By hypothesis the degree of L_{1} is at most that of L̃_{1}. Hence the degree of V_{k}/V, namely m, must be 1, and the solution x = t is unique.

(iii) ⇔ (iv): Let g be the unique equivalence transformation y = (P (x)ỹ+Q(x))/(R(x)ỹ+ S(x)) sending Ω to the canonical form of Liouvillle, Eq.(6). Then x = F (t) is a nontrivial solution of Eq.(12) if and only if g^{−1} ∘ (x → F (x̃)) ∘ g is a nontrivial equivalence transformation preserving Ω.

Notwithstanding the clear advantage that their invariants are of low degree, using minimal representatives for a given class will allow us to simplify the equivalence problem considerably. In order to match all equations having rational coefficients in a given class, we need only consider transformations involving rational equivalence functions. A method devoted to finding these types of transformations will be presented in the following section.

We can show whether or not an Abel equation Ω(x, y, y′) is a minimal representative for its class by determining the number of solutions to the system Eq.(12). The obvious method is to compute Q(x, t) as in the proof and check whether it is equal to x − t.

We may also consider a representative of a super-class which may be minimal when its parameters are general but not when they take specific values. In this case determining minimal representatives for all the subclasses is more involved. We leave the details for another work.

If an equation Ω is not minimal, we can use it to find a minimal equivalent Ω̃ by the construction outlined in the proof. By this process we can ensure that each solvable class in S is represented by a minimal equation.

Calculating minimal equations could be used as part of a new method for solving the equivalence problem. Since the theorem provides a constructive method of computing them, minimal representatives can be considered as just another canonical form for Abel equations. Furthermore, it is clear from the definition that any two minimal representatives for a given Abel class are related through a transformation wherein the equivalence function F (x̃) is rational of first degree, that is, a fractional linear or Möbius transformation of the independent variable. Starting with equations in these minimal canonical forms, the equivalence problem would be reduced to finding the three constant parameters of this Möbius transformation, a comparably trivial operation.

The main calculation required to reduce to a minimal equation involves taking the GCD of the numerators of L_{1}(x) − L_{1}(t) and L_{2}(x) − L_{2}(t). A secondary calculation is simplifying the invariants L_{i}(x) to obtain L̃_{i}(x̃). This should be the less significant. For non-parameterized classes, these are precisely the operations used in the equivalence method of Liouville outlined in the previous section. On the other hand, for multiparameter classes, there is a significant cost improvement, since only one new variable t need be introduced - the parameters need not be solved for at this stage. One could use a similar trick as that described in [5] - giving values to various symbolic constants in the input equations in order to reduce the complexity of the calculation, and then interpolating.

This method suffers however in one respect. All four functions L_{1}(x), L_{2}(x), L_{1}(t), L_{2}(t) may have high degree; whereas with Liouville’s method we can ensure that L_{1}(x) and L_{2}(x) are the invariants of minimal representatives, and therefore of relatively low degree. In any case, this method still involves generating huge expressions by equating two invariants in terms of different variables and solving. We will keep it as our standard for generating minimal equations for the database S, but for matching there is a more efficient method yet.

4 Invariant component polynomials

Apart from being of low degree in practice, the Liouville invariants of Eq.(9) have a further special property: they can be factored into products of various polynomials, many of which are common to more than one invariant. As such, these invariants can be analyzed in terms of their constituent polynomials. We will develop the technique in this section and employ it in the next to solve the equivalence problem for Abel Inverse Riccati equations.

4.1 Definitions and properties

This proposition was given by Zippel in [16], but we supply an alternate proof:

Proposition 4.1. Let K be a field. Suppose P, Q ∈ K[x, y] are relatively prime and both homogeneous, and X, Y ∈ K[t] are relatively prime. Then P (X(t), Y (t)) and Q(X(t), Y (t))) are relatively prime in K[t].

Proof: It is clear from the definition that resultant of P and Q with respect to x will be homogeneous in y, and vice versa:

R(y) ≡ res_{x}(P, Q) = ry^{n}, S(x) ≡ res_{y}(P, Q) = sx^{m}

for some r, s ∈ K and n, m ∈ N. Since

gcd(X(t), Y (t)) = gcd(rY (t)^{n}, sX(t)^{m}) = 1,

it follows that we can write 1 as a homogeneous linear combination of R(Y (t)) and S(Y (t)), or, consequently, of P (X(t), Y (t)) and Q(X(t), Y (t)).

Let Ω(x, y, y′) be an rational Abel equation with coefficients in K(x) for some field K and A = {I_{i}(x) ∈ K(x) : i = 1, 2, ...} be an absolute invariant sequence for Ω. Consider the homogenized invariants I_{i}(x_{1}/x_{2}) ∈ K(x_{1}, x_{2}).

Definition 4.1. With respect to (Ω, A), the invariant multiplicity M of a nonconstant polynomial P_{M}(x_{1}, x_{2}) is the sequence M = (m_{i}, i = 1, 2, ...), where m_{i} is the number of times P_{M} divides the numerator of the homogenized invariant I_{i}(x_{1}/x_{2}) minus the number of times it divides the denominator. A multiplicity is prime if no prime number divides all its entries.

Definition 4.2. The set of invariant component polynomials of order r for (Ω, A) is the unique collection of monic, relatively prime, homogeneous polynomials in K[x_{1}, x_{2}] with distinct, prime invariant multiplicities which together constitute all the non-constant factors of the first r homogenized invariants. The invariant signature of the equation Ω is the corresponding set of multiplicities, usually denoted s.

We can therefore write the invariants as a product of their component polynomials, up to a constant:

I_{i} (\frac{x_{1}}{x_{2}} ) = A_{i} ∏_{M∈s} P_{M}(x_{1}, x_{2})^{m_{i}}, i = 1..r; (P_{M}, P_{N}) = 1, M, N ∈ s.

Calculating the invariant polynomials involves two steps: factoring the invariants over their common field of coefficients and then multiplying together all factors of the invariants which have proportional multiplicities.

These polynomials are not only components of the invariants, but are also in a sense invariant in their own right.

Proposition 4.2. Let P_{M}(x_{1}, x_{2}), M ∈ s be the invariant component polynomials of order r for an Abel equation Ω with signature s; let G and H be relatively prime in K[x̃]. The invariant component polynomials P̃_{M} of order r for the equation Ω̃ obtained from Ω by the transformation x = G(x̃)/H(x̃) are given by:

P̃_{M}(x̃_{1}, x̃_{2}) = α_{M}x̃_{2}^{d_{M}} P_{M} (G (\frac{x̃_{1}}{x̃_{2}} ) , H (\frac{x̃_{1}}{x̃_{2}} ))

for some nonzero constants α_{M}; d_{M} is the degree with respect to x̃ of P_{M}(G(x̃), H(x̃)).

Proof: A common factor x̃_{1} − kx̃_{2} of any two of these polynomials would imply the common factor x̃ − k among {P_{M}(G(x̃), H(x̃))}; these, however, are relatively prime by Propostion 4.1. The invariant multiplicity of P̃_{M}(x̃_{1}, x̃_{2}) is the same as that of P_{M}(x_{1}, x_{2}), namely M.

This proposition is not enough to make the invariant component polynomials invariants in the sense of Eq.(7). We can, however, use them to generate invariants. For example, if p_{M} and p_{N} have the same degree, then

\frac{p̃_{M}(x̃)}{p̃_{N}(x̃)} = \frac{α_{M}}{α_{N}} \frac{p_{M}(F (x̃))}{p_{N}(F (x̃))},

where p_{M}(x) = P_{M}(x, 1), etc. and x = F (x̃) is any rational equivalence function. The logarithmic derivative of p_{M}/p_{N} is therefore a relative invariant.

For future reference, we will give a special name to the invariant polynomials in one variable. Notationally, we will use lower-case.

Definition 4.3. The inhomogeneous invariant polynomial corresponding to an invariant polynomial P_{M}(x_{1}, x_{2}) is p_{M}(x) ≡ P_{M}(x, 1).

Corollary 4.3. The invariant signature of order r for an Abel equation is preserved under rational equivalence transformations, and is therefore a unique and easily computable characteristic of any given class.

This is an extremely important corollary. If we want to determine which of a set of known solvable equations {Ω_{i}} may be equivalent to an input equation Ω, it is enough to consider the subset of equations which each share the same signature as Ω. This can reduce the time of computation dramatically. The converse unfortunately does not hold. Many different classes, particularly parameterized classes, may share the same signature.

4.2 A general method for solving Abel equations based on invariant polynomials

We can make use of invariant component polynomials and other invariants generated from them to enhance Liouville’s equivalence method. The resulting algorithm can be outlined as follows.

Algorithm 1. Solve-Abel.

Input: An Abel equation Ω̃ and a set of S of super-classes with known solutions, each represented by an equation in rational normal form, a set of invariants and invariant polynomials, and a corresponding signature for that equation.

Output: The solution of Ω̃, if Ω̃ is a member of a class in S.

1.1. Put the input equation Ω̃ in rational normal form (cf. Eq.(4)).

1.2. Calculate its differential invariants using Eq.(4) and the preceding.

1.3. Calculate the signature s and corresponding invariant component polynomials P̃_{M}, M ∈ s as explained above.

1.4. If the signature s matches that of a super-class in S, let each invariant polynomial of the representative Ω of that class be proportional, with ratios a_{M}, to the inhomogeneous polynomial of the input equation having the same multiplicity M:

P_{M}(x_{1}, x_{2}) = a_{M} p̃_{M}(x̃), M ∈ s. (13)

Note that each equation in Eq.(11) yields an additional relation among the unknowns α_{M}. We can use the corresponding relations among the a_{M} to simplify the system (13).

1.5. Solve the system for the two polynomials x_{1} = G(x̃), x_{2} = H(x̃), the constants {a_{M} = α_{M}, M ∈ s}, and any parameters defining the class of Ω̃ within the superclass found in S. Any method may be used, such as by generating new invariants as mentioned above.

1.6. If a match is found, substitute x = G(x̃)/H(x̃) and the appropriate parameter values into the known solution for the representative equation of the class which was matched.

1.7. Reverse the dependent variable transformation which brought the equation to rational normal form, obtaining the solution of the input equation Ω̃.

In the following section we will apply this new method to matching the AIR class. In particular, we develop an appropriate method for step 5. above - solving the particular system (13) in the case of AIR equations.

5 The Abel Inverse Riccati equation

5.1 The intrinsic parameter dependence of the Abel Inverse Riccati equation

In the paper [6] the authors defined the Abel Inverse Linear (AIL) and Abel Inverse Riccati (AIR) classes of ODEs. These classes are characterized by having representatives of the form

y′ = \frac{a_{3}y^{3} + a_{2}y^{2} + a_{1}y + a_{0}}{(s_{2}x^{2} + s_{1}x + s_{0})y + r_{2}x^{2} + r_{1}x + r_{0}} , (14)

with AIL being the special case s_{2} = r_{2} = 0. Furthermore, it was shown that in general the AIL class depends on 2 parameters intrinsically, and it is straightforward to show that a two parameter representative for AIL is

y′ = −\frac{dy^{3} + ay + b}{y − x}.

We will show that the AIR class has, in general, a similar three parameter representative.

Proposition 5.1. Every equation of the form Eq.(14) is equivalent to either a separable or inverse-linear equation, or an equation of the form:

AIR_{d,a,b,c} : y′ = −\frac{dy^{3} + ay^{2} + by + c}{y − x^{2}}. (15)

Proof: The first step is to arrive at a representative with r_{2} = r_{1} = 0:

y′ = \frac{a_{3}y^{3} + a_{2}y^{2} + a_{1}y + a_{0}}{(s_{2}x^{2} + s_{1}x + s_{0})y + r_{0}} , (16)

We start by applying to Eq.(14) a linear transformation in y, y = ỹ + A:

y′ = \frac{a_{3}y^{3} + ã_{2}y^{2} + ã_{1}y + ã_{0}}{(s_{2}x^{2} + s_{1}x + s_{0})y + (r_{2} + As_{2})x^{2} + (r_{1} + As_{1})x + (r_{0} + As_{0})}, (17)

The discriminant of (r_{2} + As_{2})x^{2} + (r_{1} + As_{1})x + (r_{0} + As_{0}) with respect to x can be solved for A unless

s_{1}^{2} = 4s_{2}s_{0}, r_{1}s_{1} = 2(s_{2}r_{0} + r_{2}s_{0}).

This condition implies that Eq.(14) is either separable (s_{2} = s_{1} = s_{0} = 0) or a member of AIL: if s_{2} = s_{1} = r_{2} = 0 this is apparent; if s_{1} = s_{0} = r_{0} = 0 we need to apply x = 1/x̃; otherwise s_{1} ≠ 0 and we use x = 1/x̃ − 2s_{0}/s_{1}. When A does solve the discriminant, Eq.(17) takes on either the form Eq.(16), or else

y′ = \frac{a_{3}y^{3} + a_{2}y^{2} + a_{1}y + a_{0}}{(s_{2}x^{2} + s_{1}x + s_{0})y + r_{2}(x − r_{0})^{2}} . (18)

This latter takes the form Eq.(16) under x = 1/x̃ + r_{0}. If s_{2} = 0 Eq.(16) represents an inverse-linear equation; otherwise, the transformation y = 1/ỹ followed by appropriate linear transformations in x and y yields Eq.(15).

5.2 Solution of the Abel Inverse Riccati equation

The first step in obtaining a solution for the general AIR equation is to ensure that the corresponding Riccati equation takes the simplest possible form. We present only the most general case, where the cubic

q_{3}(y) ≡ dy^{3} + ay^{2} + by + c

of Eq.(15) has three distinct roots, α, β, γ. The equation to be solved is:

y′ = −\frac{d(y − β)(y − γ)(y − δ)}{y − x^{2}}.

If we then apply the transformation x = x̃−\sqrt{γ}, y = δ+(γ−δ)(1+ỹ(β−γ)(δ−β)^{−1})^{−1}, we arrive at an equation of the form:

y′ = \frac{(1 − y)y}{ϵ(x − ξ)(x − η)y + x(ϕx + κ)(1 − y)}

for some constants ϵ, ξ, η, ϕ, κ with ϵ, ξ ≠ 0; κ ≠ −ξϕ. The further transformation x = ξ(x̃ − κ)/(x̃ + ξϕ) yields an equation of the form:

y′ = \frac{y(y − 1)}{(x − σ)(x − λ)y − x(x − κ)}.

This Abel equation admits a solution in terms of hypergeometric functions

0 = (κ − 1)(κ − x) _{2}F_{1}([λ − κ, σ − κ], [1 − κ], y)
+(κ − σ)(κ − λ)y _{2}F_{1}([λ − κ + 1, σ − κ + 1], [2 − κ], y)
+Cy^{κ} (λσy _{2}F_{1}([λ + 1, σ + 1], [κ + 2], y) + x(κ + 1) _{2}F_{1}([λ, σ], [κ + 1], y)) .

for general values of σ, κ, λ. The other cases are similar.

5.3 Minimal representatives for the class AIR

It is easy to check that by a scaling of the variables x → ξx̃, y → ηỹ, and if necessary with an extra transformation of the form x → −1/x̃, y → 1/ỹ, Eq.(15) can be reduced to one of two canonical forms:

AIR[d, b, c] : y′ = −\frac{dy^{3} + y^{2} + by + c}{y − x^{2}}; (19)

AIR[c] : y′ = −\frac{y^{3} + c}{y − x^{2}} . (20)

Note that when b ≠ 0, the transformation

x → −b/x̃, y → b^{2}/ỹ (21)

manifests an equivalence between AIR[d, b, c] and AIR[cb^{−3}, b, db^{3}]. In particular, this transformation leaves AIR[d, b, db^{3}] unchanged. By Theorem 3.1, AIR[d, b, db^{3}] is not minimal; if this is the only symmetry of AIR[d, b, db^{3}] we can obtain a minimal equation by taking x̃ = x − b/x.

Similarly AIR[c] is unchanged by the transformation x → ζx̃, y → ζ^{2}ỹ, where ζ^{3} = 1. When c ≠ 0, AIR[c] is also unchanged by x → −z/x̃, y → z^{2}/ỹ, where z^{3} = c. If these are the only symmetries of AIR[c], we obtain a minimal equation with x̃ = x^{3} − c/x^{3}.

Let us now explicitly develop minimal equations for these two singular classes, each in a form which will facilitate their matching. First, starting from AIR[d, b, db^{3}], with b, d ≠ 0, we introduce new parameters A = 1/bd, C = b^{3}d^{2}, apply the transformation

y = \frac{A^{3}C(x^{3}ỹ − A^{3}C^{2})}{x(A^{3}Cỹ + x^{3})},

and then take x̃ = x/A − AC/x to be the new independent variable. We obtain the equation

(x^{2} + 4C)(x^{2} + C)^{2}y′ = (Cx^{2} + C^{2}(A + 1)) y^{3}
+ (CAx^{3} − 2Cx^{2} + 4C^{2}Ax − 2C^{2}) y^{2}
+ (x^{5} + CAx^{4} + 4Cx^{3} + 3C^{2}(A + 1)x^{2} + 3C^{2}x − C^{3}(A − 3)) y
+ Cx^{5} + 4C^{2}x^{3} − 4C^{2}x^{2} − C^{3}(A − 3)x − 4C^{3}, A, C ≠ 0,  (22)

which we will show to be minimal except when (A, C) = (−5, −1/64) or (−5, 1/16). When (A,C) = (-5,-1/64), we start from Eq.(22), take x̃ = 4x + 1/4x to be the new independent variable, and put the resulting equation in rational normal form, obtaining the minimal representative

y′ = −\frac{25(5x − 8)^{2}}{729x^{3}(x + 2)^{2}(x − 2)^{3}} y^{3} + \frac{60x^{3} − 89x^{2} − 128x + 256}{6x(x − 2)(5x − 8)(x + 2)}y + 1. (23)

Similarly, when (A,C) = (-5,1/16), we take x̃ = 2x + 1/2x to be the new independent variable, and put the resulting equation in rational normal form, obtaining the minimal equation

y′ = \frac{50}{729x^{3}(x + 2)^{3}} y^{3} + \frac{2(3x^{2} − 7x − 7)}{3x(x − 2)(x + 2)}y + 1. (24)

For AIR[0, b, 0] with b ≠ 0 we instead take x̃ = x − b/x as the new independent variable, put the resulting equation in rational normal form, and then apply y → 9ỹ(x̃^{2} + 4b)/(2x̃^{3} + 6x̃^{2} + (9b + 6)x̃ + 18b − 16), arriving at:

y′ = \frac{y^{3} + (−3x^{2} + 3x − 9b − 12)y + 2x^{3} + 6x^{2} + (9b + 6)x + 18b − 16}{9x^{2} + 36b}, b ≠ 0. (25)

Finally, for AIR[c], the transformation y → (2x^{3}ỹ − c)/(2ỹ + x^{3})/x yields a minimal equation for the class, in terms of the new variable x̃ = x^{3} − c/x^{3}:

ỹ′ = \frac{8ỹ^{3} − 8ỹ^{2} + (2x̃ + 6c)ỹ + cx̃ − 4c}{6x̃^{2} + 24c}. (26)

We have proved:

Theorem 5.2. Every Abel Inverse Riccati equation is equivalent to either a separable or inverse-linear equation, Eq.(19) with either c ≠ db^{3} or b = 0, (22) with A, C ≠ 0 and (A, C) ∉ \{(−5, \frac{1}{16}), (−5, −\frac{1}{64})\}, or one of Eqs.(23-26).

5.4 Absolute invariants and invariant polynomials for class AIR

We have established that any given AIR equation is equivalent to a representative equation of the form Eq.(15). For each choice of parameters d, a, b, c we obtain the representative equation of a specific Abel class. As we have seen, not all of these choices yield minimal representatives. However, the invariant polynomials for any of these equations, or for their minimal equivalents, can be obtained from those of the general equation (15).

The next step is therefore to determine the rational normal form, invariants, and invariant polynomials for this general equation. We follow the prescription of the previous sections. The expressions involved are complicated, so for convenience we will abbreviate many of the polynomials we shall encounter, beginning with q_{3}(x) ≡ dx^{3} +ax^{2} +bx+c. As a convention, the subscript of a polynomial will denote its degree (for general values of the parameters).

We start by applying the transformation y = 1/ỹ + x^{2} to Eq.(15), yielding

y′ = q_{3} (\frac{1}{y} + x^{2}) y^{3} + 2xy^{2},

or, upon expanding q_{3} in series about 1/y,

y′ = q_{6}y^{3} + 3j_{4}y^{2} + j_{2}y + d,

where q_{6}(x) ≡ q_{3}(x^{2}), j_{4}(x) ≡ \frac{1}{3}(2x + q_{3}′(x^{2})), j_{2}(x) ≡ \frac{1}{2}q_{3}′′(x^{2}). Applying the further transformation y = ỹ − j_{4}/q_{6}, we obtain

y′ = q_{6}y^{3} + \frac{j_{5}}{q_{6}} y + \frac{p_{6}}{q_{6}^{2}} ,

where

j_{5} ≡ j_{2}q_{6} − 3j_{4}^{2}, p_{6} ≡ j_{4}′q_{6} − j_{4}q′_{6} − j_{2}q_{6}j_{4} + q_{6}^{2} + 2j_{4}^{3}. (27)

A final transformation, y = p_{6}/q_{6}^{2} ỹ, brings the equation into rational normal form:

y′ = \begin{matrix}p_{6} & 2 \\ q_{6} & 3 \end{matrix} y^{3} − \frac{p_{6}′q_{6} − j_{5}p_{6} − 2q_{6}′p_{6}}{p_{6}q_{6}} y + 1.

We use Eqs. (8, 9) to calculate the relative and absolute invariants for Eq.(15):

s_{3} = \begin{matrix}p_{6} & 4 \\ q_{6} & 6 \end{matrix} , s_{5} = −\frac{p_{10}p_{6}^{5}}{q_{6}^{10}} , t_{10} = −\frac{p_{8}p_{6}^{10}}{q_{6}^{18}} , t_{12} = −\frac{p_{12}p_{6}^{12}}{q_{6}^{22}} ;

L_{1} = −\begin{matrix}p_{10} & 3 \\ p_{6} & 5 \end{matrix} , L_{2} = \frac{p_{8}q_{6}^{2}}{p_{10}^{2}} , L_{3} = \frac{p_{12}p_{6}}{p_{8}p_{10}} ; (28)

where we have defined ^{7}

p_{10} ≡ 3j_{5}p_{6} + 3p_{6}q_{6}′ − p_{6}′q_{6}, (29)

p_{8} ≡ (5p_{6}′p_{10} − 3p_{6}p_{10}′)/q_{6}, p_{12} ≡ 10p_{8}j_{5} + 8p_{8}q_{6}′ − p′_{8}q_{6}. (30)

7 Note that p_{8} is indeed a polynomial - one can check that q_{6} divides 5p_{6}′p_{10} − 3p_{6}p_{10}′ by substituting the expressions for p_{10}, p_{6}, and j_{5}.

If the parameters are general, the inhomogeneous invariant component polynomials of order 3 for Eq.(15) and invariants L_{i} are, up to constant factors, given by

p_{−5,0,1} = p_{6}, p_{0,1,−1} = p_{8}, p_{3,−2,−1} = p_{10}, p_{0,0,1} = p_{12}, p_{0,1,0} = q_{6}^{2}. (31)

The corresponding homogeneous polynomials are pairwise relatively prime and are the only nontrivial factors of the homogenized invariants.

5.5 The matching algorithms for class AIR

Given an input Abel equation Ω̃ of the form Eq.(1) with coefficients g_{0}, g_{1}, f_{0}, f_{1}, f_{2}, f_{3} polynomial in x, we would like to determine whether it is a member of the class AIR, and if so, determine both the equivalence transformation Eq.(2) and the values of the parameters which brings it to one of the forms (19) with c = db^{3} or b = 0, (22) with A, C ≠ 0 and (A, C) ∉ \{(−5, \frac{1}{16}), (−5, −\frac{1}{64})\}, or (23-26). We will assume for now that the parameters are such that each of these equations is minimal. By Definition 3.1 then, the equivalence function F must be rational.

We start by describing the most complicated case, and follow with the others, which are similar.

Algorithm 2. Match-AIR[d, b, c], c ≠ db^{3}, or b = 0.

Input: An Abel equation Ω̃.

Output: [F, d, b, c], or FALSE.

At each step there are various conditions to check which are necessary for an equivalence. If any of these conditions fail, we know to return FALSE.

Step 2.1: Find q_{6}(F)^{2}H^{12}, p_{6}(F)H^{6}, p_{10}(F)H^{10} (up to multiplication by a constant). We follow steps 1.1-1.3 of Algorithm 1. Using Propostion 2 and Eq.(31), we note that in the general case, the polynomials sought are, respectively, the inhomogeneous invariant polynomials p̃_{0,1,0}, p̃_{−5,0,1}, p̃_{3,−2,−1} for the input equation which were calculated in Step 1.3. For certain exceptional values of the parameters d, b, c, this is not the case and we must compute these expressions somehow from the invariant polynomials. We will consider these special cases in the following subsections.

Returning to the general case, there is one more simplification and a convenient abuse of notation that can be achieved. Inspired by Eq.(31), we will define

p̃_{6} ≡ p̃_{−5,0,1}, p̃_{10} ≡ p̃_{3,−2,−1}, q̃_{6} ≡ (p̃_{0,1,0})^{\frac{1}{2}} .

Also, the equation L_{1}(x) = L̃_{1}(x̃) tells us that α_{−5,0,1}^{5} = α_{3,−2,−1}^{3}; we let α = α_{−5,0,1}^{1/3} = α_{3,−2,−1}^{1/5}. By redefining appropriately, we can therefore assume that we have computed ^{8}

q̃_{6} := γq_{6}(F)H^{6}, p̃_{6} := α^{3}p_{6}(F)H^{6}, p̃_{10} := α^{5}p_{10}(F)H^{10},

where γ, α are constants to be determined. These are essentially the first three equations of system (13), and will be all that are required to determine G, H, and the parameters d, b, c. The remaining steps in Algorithm 2 are a detailed explanation of Step 1.5 of Algorithm 1, explaining how to solve this system.

Step 2.2: Find GH (up to a constant). As suggested in the discusion following Proposition 4.2, the invariant component polynomials can lead to new invariants. Since q_{6} and p_{6} are both of degree 6, D ln(q̃_{6}/p̃_{6}) is an absolute invariant (cf. (37)). After adding 1, the result happens to have a nontrivial factor depending on F^{2}:

w_{2}(F^{2}) ≡ (3db − 1)F^{4} + (9dc − b − 4)F^{2} + 3c − b^{2}.

The numerator of w_{2}(F^{2}) can be computed from the information we already have:

e_{1} := \frac{p̃_{10}}{p̃_{6}} (1 + D ln (\frac{q̃_{6}}{p̃_{6}} )) = α^{2}(3j_{5}(F) + 2q_{6}′(F))H^{4} = α^{2}w_{2}(F^{2})H^{4}. (32)

Necessary conditions are that e_{1} be polynomial in x̃ and not zero. Including q_{6}(F) = q_{3}(F^{2}), we now have two distinct polynomials in F^{2} (up to a factor of a power of H). We use this information to obtain an expression rational in F^{2}, and a differential operator on C(F):

e_{3} := \frac{e_{1}^{3}}{q̃_{6}^{2}} = \frac{α^{6}}{γ^{2}} \frac{w_{2}(F^{2})^{3}}{q_{3}(F^{2})^{2}} , D_{F} := −\frac{p̃_{10}}{e_{1}p̃_{6}} D = \frac{q_{3}(F^{2})}{w_{2}(F^{2})} \frac{d}{dF} . (33)

The chain rule conveniently gives us two expressions having F as a factor:

e_{4} := \frac{D_{F} e_{3}}{2e_{3}} = F \frac{w_{3}(F^{2})}{w_{2}(F^{2})^{2}} ,

e_{5} := e_{3}D_{F} e_{4} = \frac{α^{6}}{γ^{2}} \frac{w_{5}(F^{2})}{q_{3}(F^{2})w_{2}(F^{2})}, e_{6} := \frac{e_{3}D_{F} e_{5}}{2} = \frac{α^{12}}{γ^{4}} F \frac{w_{8}(F^{2})}{q_{3}(F^{2})^{3}} ,

where

where there is no ambiguity.

8 We will employ three notations for equality: ‘:=’ denotes the assignment of a variable as the result of a computation. ‘≡’ denotes the mathematical definition of an unknown quantity in terms of others. ‘=’ denotes an assertion of mathematical equality between quantities which may or may not have been calculated. We will also drop the functional dependence on x̃

w_{3}(x) ≡ 3w_{2}′(x)q_{3}(x) − 2w_{2}(x)q_{3}′(x),  (34)
w_{5}(x) ≡ w_{3}(x)w_{2}(x) + 2x(w′_{3}(x)w_{2}(x) − 2w_{3}(x)w′_{2}(x)),  (35)
w_{8}(x) ≡ w′_{5}(x)q_{3}(x)w_{2}(x) − w_{5}(x)q_{3}′(x)w_{2}(x) − w_{5}(x)q_{3}(x)w′_{2}(x).  (36)

We would like to obtain F from e_{4} and e_{6} as some kind of greatest common divisor. However, F itself is not polynomial. Instead we determine two expressions, necessarily polynomial in x̃ and having G(x̃)H(x̃) as their greatest common divisor:

e_{4}e_{1}^{2} = α^{4}GHw_{3}(F^{2})H^{6}, e_{6}q̃_{6}^{3} = \frac{α^{12}}{γ} GHw_{8}(F^{2})H^{16}.

e_{7} := gcd(e_{4}e_{1}^{2}, e_{6}q̃_{6}^{3}) = \frac{βγ}{α^{2}} GH.

This computation is valid as long as the polynomials w_{3} and w_{8} defined in (34,36) are relatively prime, as they are in general. We will consider the exceptional case where they have a common factor in the next subsection. β is a non-zero constant to be determined.

Step 2.3: Find β^{2}. The expression e_{7} leads us to a simpler differential operator on C(F):

D := \frac{e_{7}e_{1}}{q̃_{6}} D_{F} = βF \frac{d}{dF} . (37)

Proceeding, we compute

e_{8} := \frac{e_{1}}{e_{7}^{2}} = \frac{α^{6}}{γ^{2}β^{2}} ((3db − 1)F^{2} + \frac{3c − b^{2}}{F^{2}} + (9dc − b − 4)) ,

D^{n}e_{8} = \frac{β^{n−2}2^{n}α^{6}}{γ^{2}} ((3db − 1)F^{2} + (−1)^{n} \frac{3c − b^{2}}{F^{2}} ) , n > 0.

The condition c ≠ db^{3} or b = 0 implies that D^{n}e_{8} never vanishes. We thus obtain more necessary conditions: for all n, D^{n}e_{8} is not constant and the expression

β^{2} := \frac{D^{n+2}e_{8}}{4D^{n}e_{8}} , n = 1, 2

is constant and does not depend on n.

Step 2.4: Find F, d, b, c. Applying the operator D to e_{9} ≡ q̃_{6}/e_{7}^{3}, we obtain:

D^{n}e_{9} = β^{n−3} \frac{α^{6}}{γ^{2}} (3^{n}dF^{3} + F + (−1)^{n} \frac{b}{F} + (−3)^{n} \frac{c}{F^{3}} ) , n = 0 .. 3.

We can solve this system for the powers of F if we know β. So we choose a square root of β^{2}, and call it β̃. If β̃ = β, we obtain:

v := \frac{1}{48} ( \begin{matrix}−3 & − 1 & 3 & 1 \\ 27 & 27 & −3 & −3 \\ 27 & −27 & −3 & 3 \\ −3 & 1 & 3 & −1 \end{matrix} ) ( \begin{matrix}β̃^{3}e_{9} \\ β̃^{2}De_{9} \\ β̃D^{2}e_{9} \\ D^{3}e_{9} \end{matrix} ) = \frac{α^{6}}{γ^{2}} ( \begin{matrix}dF^{3} \\ F \\ bF^{−1} \\ cF^{−3} \end{matrix} ) . (38)

If β̃ = −β, the calculated v takes the value −α^{6}γ^{−2}(cF^{−3}, bF^{−1}, aF, dF^{3}). The additional condition b = 0 would imply v_{2} = 0; if this is found to be the case we should choose the other root β̃ := −β̃ and recalculate v as above. Otherwise b ≠ 0 and so on account of the transformation Eq.(21), the function F̃ ≡ −b/F is an equally valid equivalence function. If we redefine F := F̃ as the new equivalence function, the definition of β also must change so that Eq.(37) still holds. Since F̃ \frac{d}{dF̃} = −F \frac{d}{dF}, the required redefinition of β is to multiply it by -1. We now have β̃ = β after all.

We finally compute

\frac{α^{6}}{γ^{2}} := \frac{2}{β̃De_{8}} ((3v_{1}v_{3} − v_{2}^{2}) − (3v_{2}v_{4} − v_{3}^{2})) .

Using Eq.(38), the equivalence function F and the parameters d, b, c follow directly. Necessary conditions are that d, b, c and \frac{α^{6}}{γ^{2}} are constant and that F depends on x.

Step 2.5: Verification. These computed values for F and the parameters are as yet only candidates, in part because they have been calculated using only three of the five invariant polynomials. We can compare the invariants of the input equation with those of the candidate, using the values of F and d, b, c computed. A 3sufficient condition for a match is that the first two invariants agree, since the rational normal form can be expressed in terms of these first two invariants, as in Eq.(6).

Algorithm 3. Match-AIR[d, b, db^{3}], d, b ≠ 0.

Input: An Abel equation Ω̃.

Output: [f, d, b], or FALSE.

As we have seen, Eq.(19) fails to be minimal when c = db^{3}; a minimal representative for this class is given as Eq.(22). The corresponding invariants L_{i}[19] and L_{i}[22] satisfy

L_{i}[22](g/h) = L_{i}[19](G/H)|_{c=db^{3},b=CA^{2},d=1/CA^{3}},

where g = G^{2} −A^{2}CH^{2}, h = AGH are relatively prime polynomials and f ≡ g/h is the equivalence function to be determined. In general, the signature of Eq.(22) is the same as for Eq.(19). By Proposition 4.2 the corresponding invariant polynomials satisfy

α_{M}P_{M}[22](g, h) = P_{M}[19](G, H)|_{c=db^{3},b=CA^{2},d=1/CA^{3}}, M ∈ s.

Exceptions arise when any pair of the invariant polynomials have common factors, that is, when the resultant between such a pair vanishes. These cases will be discussed in the following sections.

We now follow the computations of Step 2.2 above, noticing that the polynomials w_{3} and w_{8} have a second common factor G^{2} − A^{2}CH^{2}, so that for some constant δ, we have, in general,

e_{7} = δγα^{−1}A^{−2}gh. (39)

It can be checked that in this situation (e_{4}e_{1}^{2}/gh)^{2} and (e_{6}q̃_{6}^{3}/gh)^{2} are homogeneous polynomials in g^{2} and h^{2}. The four exceptional cases where they have a common factor are determined by setting to 0 their resultant with respect to g^{2}:

A = 3; C = \frac{4}{(A − 3)^{2}} ; C = \frac{4}{(A + 1)(A − 3)^{2}} ;

C = \frac{k}{(A − 3)^{2}} , A = −\frac{3(9k + 16)(k − 4)}{k(27k + 20)}.

In each case there is a simple way to determine the offending factor and thereby determine the correct value for e_{7} given above. We omit the details.

We now use equations (32,33,39) to compute the expressions Ψ, Θ, Φ and operators D_{f}, D_{f}, introducing another unknown constant β:

Ψ ≡ \frac{e_{1}}{e_{7}} = −\frac{α^{3}}{δγ} (A(A − 3)f + \frac{3(A + 1)(A − 3)C + 4}{f}) ,

D_{f} ≡ \frac{Ψe_{7}^{3}}{q̃_{6}^{2}} D_{F} = \frac{α^{2}e_{7}^{2}}{γq̃_{6}H^{2}} \frac{d}{dF} = \frac{δ^{2}Cf^{2}}{f^{2} + (A + 1)C} \frac{d}{df} ,

Θ ≡ \frac{denom(D_{f}Ψ)}{e_{7}} = \frac{β}{δ^{2}C} (f + \frac{C(1 + A)}{f}) ,

Φ ≡ \frac{q̃_{6}^{2}}{Θ^{2}e_{7}^{3}} = \frac{α^{3}δ}{β^{2}γ} (f + \frac{4C}{f} ) ,

D_{f} ≡ ΘD_{f} = βf \frac{d}{df} .

We are finally able to determine a series of constants:

β^{2} := Φ^{−1}D_{f}^{2}Φ,

λ ≡ ΦD_{f}Ψ − ΨD_{f}Φ = \frac{2α^{6}}{βγ^{2}} (4 − (A − 3)^{2}C),
μ ≡ ΦD_{f}Θ − ΘD_{f}Φ = \frac{−2α^{3}}{δγ} (A − 3),
ν ≡ ΨD_{f}Θ − ΘD_{f}Ψ = \frac{2β^{2}α^{3}}{δ^{3}γ} \frac{(A + 1)(A − 3)^{2}C − 4}{C},
ψ ≡ β^{2}Ψ^{2} − (D_{f}Ψ)^{2} = \frac{4α^{6}β^{2}}{δ^{2}γ^{2}} A(A − 3)(3(A + 1)(A − 3)C + 4),
θ ≡ β^{2}Θ^{2} − (D_{f}Θ)^{2} = \frac{4β^{4}}{δ^{4}} \frac{A + 1}{C},
ϕ ≡ β^{2}Φ^{2} − (D_{f}Φ)^{2} = \frac{16δ^{2}α^{6}}{β^{2}γ^{2}} C.

The case A = −1 is identified by the condition θ = 0; we obtain C = β^{2}μ^{4}/4ϕν^{2}. Otherwise, we form an overdetermined system for A and C:

\frac{16β^{2}μ^{2}}{ϕθ} = \frac{(A − 3)^{2}}{A + 1} , \frac{16ψ}{ϕθ} = \frac{A(A − 3)(3(A + 1)(A − 3)C + 4)}{A + 1},

\frac{λ^{2}β^{2}}{ϕ^{2}θ} = \frac{((A − 3)^{2}C − 4)^{2}}{256C(A + 1)} , \frac{64β^{2}ν^{2}}{ϕθ^{2}} = \frac{((A + 1)(A − 3)^{2}C − 4)^{2}}{(A + 1)^{2}C} .

This system admits a unique solution unless C = −1/(A − 3)^{2} or C = −16/((A − 3)^{2}(A + 1)), in which case it is invariant under the substitution A → (15 − A)/(A + 1). The duplicity of solutions occurs because the classes of the corresponding equations Eq.(22) with C = −1/(A − 3)^{2} or with C = −16/(A + 1)/(A − 3)^{2} happen to be equivalent to each other. Consequently, either solution for (A, C) is valid.

Given A and C, we obatain f in three different manners, depending on the case:

A ∉ \{3, −5\} : β := \frac{μ(5 + A)}{(β^{−2}D_{f}ΦD_{f}Θ − ΦΘ)(A − 3)}, f := \frac{Ψ + β^{−1}D_{f}Ψ}{μA};

A = −5, C ∉ \{\frac{1}{16}, −\frac{1}{64}\} : f := \frac{16C − 1}{4Ψμ^{−1} + 5Φ(64C + 1)θν^{−1}} ;

A = 3 : f := \frac{θD_{f}Φ}{4D_{f}ΘD_{f}Ψ}.

In the remaining cases A = −5, C ∈ \{\frac{1}{16}, −\frac{1}{64}\} it happens that the equation Eq.(22) is not minimal; minimal equivalents are provided in (23,24) and can be matched by the standard procedure for non-parameterized classes.

Algorithm 4. Match-AIR[0, b, 0], b ≠ 0.

Input: An Abel equation Ω̃.

Output: [f, b], or FALSE.

Algorithm 3 showed how to match the class AIR[d, b, db^{3}] when d ≠ 0. The situation is similar when d = 0, even though the minimal equation for the class, namely (25), is not obtained from Eq.(22) by setting d = 0.

The equivalence function f = g/h to be determined is related to that of AIR[d, b, c], namely F = G/H, by means of g = G^{2} − bH^{2}, h = GH. We first compute the invariant polynomials, noting that in general they are as expected: q̃_{6}, p̃_{6}, p̃_{10}, ..., and exceptions arise only when the various resultants between these invariant polynomials vanish. The corresponding values of b are:

b = −1, −4, −36, −16;

the roots of the following irreducible polynomials,

b^{2} + 3b + 1, b^{2} − 4b + 16, 25b^{2} − 184b + 400, 2b^{3} + 3b^{2} − 12b + 12;

and the roots of four other irreducible polynomials of degree at least six. We find it unlikely that an input equation will have as a b-value one of the roots of these four large-degree polynomials, and thus will ignore these cases as unimportant. For the other exceptional values of b, it is straightforward to determine a priori what are the cancellations, and match the class by recognizing these cases, putting back the cancelled factors and then continuing as explained below.

Having computed the correct polynomials q̃_{6}, p̃_{6}, p̃_{10}, we follow step 2.2 to compute e_{7}, and proceed by finding Ψ, Θ, Φ, λ, μ, ν, ψ, θ, ϕ, A, C, and f as above. The only exceptional cases where e_{7} has an extra factor occur with b either 4 or -20/27. As in the previous section, there are simple methods to determine the offending factor in these cases, and thereby determine the correct value of e_{7} and so continue as normal.

We should find θ = 0, A = −1, and C = 1/16b; we will have calculated a value for f but it will not be the correct equivalence function. If d were not zero, the condition A = −1 would signify the class given by d = −1/b. In order to distinguish between the two classes AIR[−1/b, b, −b^{2}] and AIR[0, b, 0], we perform a test, calculating the variable

Υ ≡ \frac{16Cf + (3 + 64C)f^{−1}}{Ψ} .

Υ will be independent of x for the class AIR[0, b, 0], and will depend on x for the class AIR[−1/b, b, −b^{2}] unless b = −1/16. However, it is easy to check that this last case AIR[16, −1/16, −1/256] is equivalent to the class AIR[0, −1, 0], and so there is no harm in assuming that we are matching the case AIR[0, b, 0] when Υ does not depend on x.

To finish, we compute b := 1/16C and redefine f := −(3b + 4)/f unless b = −4/3, in which case we obtain f := 2Ψ/μ.

Algorithm 5. Match-AIR[c].

Input: An Abel equation Ω̃.

Output: [f, c], or FALSE.

A minimal representative for the class AIR[c] is given as Eq.(26). The corresponding invariants L_{i}[c] and L_{i}[26] satisfy

L_{i}[26](g/h) = L_{i}[c](G/H),

where g = G^{6} − cH^{6}, h = G^{3}H^{3} are relatively prime polynomials and f ≡ g/h is the equivalence function to be determined. For general values of c, the signature for Eq.(26),

{((−5, 0, 1), (0, 0, 1), (0, 1, −1), (0, 1, 0), (2, −1, −1), (3, −2, −1)}

differs from that for AIR[d, b, c]. We must treat separately the exceptional cases which arise for certain values of c because some factors have canceled among the invariant polynomials. These valaues can be determined by finding the roots of the resultant between any pair of the corresponding homogenized invariant polynomials; they are given by

c = 0, −2/3, −1/9, −4/9, −16/9, −4,

and the roots of the following polynomials:

(81c^{2} − 36c − 40), (81c^{2} − 36c + 16), (81c^{2} + 99c + 25),
(2025c^{2} − 1656c + 400), (162c^{2} + 63c − 40), (729c^{3} − 243c^{2} − 144c + 128).

In each case, it is straightforward to recognize the situation and match the class by first determining the factors which have canceled and then employing the algorithm below.

In the general case, four of the invariant polynomials are given by

P_{0,1,0} = q̃_{6}^{2} = γ^{2}(g^{2} + 4ch^{2}),

P_{−5,0,1} = p̃_{6} = α^{3} ((2c + \frac{16}{27}) h − (c + \frac{2}{3}) g) ,

P_{2,−1,−1} = −\frac{δ^{3}}{α^{15}} h,

P_{3,−2,−1} = −\frac{α^{15}p̃_{10}GH}{δ^{2}h} = \frac{α^{15}}{27δ^{2}} ((243c^{2} − 108c − 120)g − (810c^{2} + 144c − 64)h)

for some nonzero constant δ. We then define D_{f} and a sequence of expressions:

D_{f} ≡ −\frac{P_{3,−2,−1}P_{2,−1,−1}^{2}}{P_{−5,0,1}P_{0,1,0}} D = −\frac{3δ^{4}}{γ^{2}α^{18}} \frac{d}{df} ,

e_{1} ≡ \frac{P_{0,1,0}}{P_{2,−1,−1}^{2}} , e_{2} ≡ \frac{P_{3,−2,−1}}{P_{2,−1,−1}} , e_{3} ≡ \frac{P_{−5,0,1}}{P_{2,−1,−1}} , e_{4} ≡ \frac{D_{f}e_{3}}{D_{f}^{2}e_{1}} ,

e_{5} ≡ (D_{f}e_{1})^{2} − 2e_{1}D_{f}^{2}e_{1}, e_{6} ≡ e_{3} − e_{4}D_{f}e_{1}, e_{7} ≡ \frac{D_{f}e_{2}}{D_{f}^{2}e_{1}} .

Next, c can be found as the unique solution of the system

(3c + 2)(27c + 8)e_{5} + 34992ce_{4}e_{6} = 0, (81c^{2} − 36c − 40)e_{6} + (108c + 32)e_{7} = 0,

unless c = −1/9 or c = −16/9, in which case either solution is valid, as the two corresponding classes are equivalent. Finally f is given by

c ≠ −\frac{2}{3} : f = −\frac{(3c + 2)^{2}D_{f}e_{1}}{1944e_{4}^{2}} , c = −\frac{2}{3} : f = −\frac{25}{26244} \frac{e_{5}D_{f}e_{1}}{e_{3}^{2}} .

5.6 Identification of escaping classes

We have presented a sequence of four algorithms for matching various subclasses of the Abel Inverse Riccati class. For Algorithms 2 and 3, however, the computation of the invariant polynomials required a fundamental assumption:

Assumption 1. The parameters defining the subclass are general enough that the invariant polynomials Eq.(31) remain pairwise relatively prime upon substitution of those parameters.

Algorithm 2 required a second assumption as well:

Assumption 2. The parameters defining the subclass are general enough that the polynomials w_{3} and w_{8} of (34,36) remain relatively prime upon substitution of those parameters.

A choice of parameters for which either of these assumptions fail will correspond to an escaping class, a class of AIR equations not solved by the algorithms as explained in the previous section. It is important to note that for the majority of classes in AIR, the assumptions do not fail and the algorithms are complete. Indeed, the failure of either assumption implies a particular relation among the parameters, so that any escaping class can depend on at most two parameters, and thus represents at most an infinitesimal portion of the Abel Inverse Riccati super-class. On the other hand, many of the simpler choices for the parameters which cause the assumptions to fail at the same time yield Abel classes which have been described in the literature. It is therefore important to identify and develop alternative matching schemes for as many of the escaping classes as possible. These efforts are the subject of ongoing research and will be described in a future work.

6 Conclusions and further work:

The main results of this work include a new solution for the complete three-parameter AIR super-class, as well as a set of practical algorithms for solving all but a few member equations of the class. In the process of developing these ideas, we have had provided theoretical advances in the area of matching Abel’s differential equations. In particular, in the second section we provided an important new sequence of invariants. We then developed the idea of minimal representatives for a class, discussed their importance for the equivalence problem and showed how to test an equation for minimality and compute minimal equations by means of rational function decomposition. In the fourth section we showed how to decompose a sequence of invariants into its set of invariant component polynomials, and described how this idea leads directly to a new and more practical solving strategy for Abel equations. Finally, we presented a new set of routines which together constitute the first effective solver for the most general class of solvable Abel equations known, the Abel Inverse Riccati.

There are many prospects for future work which can be suggested. First, as many as possible of the the escaping classes discussed in the previous section should be identified and new routines should be developed for them. Secondly we expect that it would not be difficult to generalize these solving routines to include the case where the equivalence function is not rational. This will provide a small but important increase in the applicability of the solver. Another generalization would be to include in the solver integrable classes other than Abel Inverse Riccati, such as those described as Abel Inverse Linear and Abel Inverse Abel in [6]. Furthermore, we hope that the idea of decomposing a set of invariants into their polynomial components would apply equally well to matching other types of differential equations, including those of Riccati and second order linear type. Finally, the ideas presented in Sections 3 and 4 for the problem of overdetermined, parameterized, rational functional decomposition are relevant independent of their applicability to differential equations. Dedicated routines can be developed for this problem.

The initial motivation for developing the algorithms presented here was to facilitate classification. In fact, even after the first prototypes were developed, we were still unaware of a solution for the general Abel Iinverse Riccati equation. Nevertheless, it seemed important to have a tool to quickly determine whether a new equation or family of eqauations under consideration was a member of a previously well-studied class. These algorithms will go a long way towards that goal.

7 Acknowledgements

This research was supported by the Centre for Experimental and Constructive Mathematics, funded by MITACS (Mathematics of Information Technology and Complex Systems); and a research assistantship provided by Peter Borwein, funded by the Natural Sciences and Engineering Research Council of Canada. I would like to especially thank Edgardo Cheb-Terrab, who introduced the subject and has been the guide for much of the work. I have also had useful discussions with Peter Borwein, Michael Monagan and Theodore Kolokolnikov and I thank them for their advice.

References

[1] Abel, N. H. 1881 Oeuvres Complétes II. S.Lie and L.Sylow, Eds., Christiana.

[2] Appell, P. 1887 Sur les équations différentielles algébriques et homogènes par rapport à la fonction inconnue at à ses dérivées. Comptes Rendus 104, 1776-1779.

[3] Appell, P. 1887 Sur les invariants des équations différentielles. Comptes Rendus 105, 55-58.

[4] Appell, P. 1889 Sur les invariants de quelques équations différentielles. Journal de Mathématique 5, 361-423.

[Comptes Rendus, July 4 1887]

[5] Cheb-Terrab, E. S., Roche, A. D. 2000 Abel Equations: Equivalence and Integrable Classes. Computer Physics Communications 130, 204-231.

[6] Cheb-Terrab, E. S., Roche, A. D. 2003 An Abel ODE class generalizing known integrable classes. European Journal of Applied Mathematics 14 217-229.

[7] Halphen, G.H. 1887 Sur la réduction des équations différentielles linéaires aux formes intégrables. Mémoires présentés par divers savants à l’Académie des Sciences, [Savants étrangers] 28, no. 1.

[8] Kamke E. 1959 Differentialgleichungen: Lösungsmethoden und Lösungen. Chelsea Publishing Co, New York.

[9] Liouville, R. 1886 Sur certaines équations différentielles du premier ordre. Comptes Rendus 103, 476-479.

[10] Liouville, R. 1887 Sur une classe d’équations différentielles du premier ordre et sur les formations invariantes qui s’y rapportent. Comptes Rendus 105, 460-463.

[11] Liouville, R. 1902 Sur une équation différentielle du premier ordre. Acta Mathematica 27, 55-78.

[12] Murphy, G. M. 1960 Ordinary Differential Equations and their solutions. Van Nostrand, Princeton.

[13] Polyanin, A. D. & Zaitsev, V. F. 1995 Handbook of Exact Solutions for Ordinary Differential Equations. CRC Press, Boca Raton.

[14] Sachdev, P. L. 1997 A Compendium of Nonlinear Ordinary Differential Equations. John Wiley & Sons.

[15] Schwarz, F. 1998 Algorithmic Solution of Abel’s Equation. Computing 61, 39-46.

[16] Zippel, R. 1991 Rational function decomposition. In Proc. ISSAC’91, 1-6.