Rational Solutions of the Mixed Differential Equation and its Application

to Factorization of Differential Operators

Mark van Hoeij Department of mathematics University of Nijmegen 6525 ED Nijmegen The Netherlands e-mail: hoeij@sci.kun.nl

March 27, 1996

Abstract

The topic of this paper is a fast method to compute the rational solutions of a certain differential equation that will be called the mixed differential equation. This can be applied to speed up the factorization of completely reducible linear differential operators with rational functions coefficients.

1 Introduction

A differential equation

y^{(n)} + a_{n−1}y^{(n−1)} + … + a_{1}y′ + a_{0}y = 0

corresponds to a differential operator of order n

f = ∂^{n} + a_{n−1}∂^{n−1} + … + a_{0}∂^{0}

acting on y. In the present paper the coefficients a_{i} are elements of the differential field k(x) and ∂ is the differentiation d/dx. The field k is the field of constants. It is assumed to have characteristic 0. \overline{k} is the algebraic closure of k. The differential operator f is an element of the non-commutative ring k(x)[∂]. This is an example of an Ore ring [12]. The equation ∂x = x∂ + 1 holds in this ring. We denote the solution space of f as V (f). It is a \overline{k} vector space of dimension order(f). One way to define V (f) without ambiguity is to define it as a subset of the universal extension, cf. [7].

A factorization f = LR where L, R ∈ k(x)[∂] is useful for computing solutions of f because the solutions V (R) of the right hand factor R are solutions of f as well. Given L and R we can ask the following question: How can one compute (if it exists) an operator R_{2} ∈ k(x)[∂] such that V (f) = V (R) ⨁ V (R_{2})? In section 3.1 it is shown how this question can be reduced to the mixed equation with c = 1. Such an operator R_{2} could also be computed by one of the factorization algorithms for differential operators, cf. [4, 5, 9, 13]. However, these factorization algorithms can be quite costly so we want to have an alternative method.

The mixed equation (called gemischte Gleichung in [11]) is the following: Let L and R in k(x)[∂]. We will assume that L and R are monic (i.e. the coefficient in k(x) of the highest power of ∂ is 1). Let c ∈ k(x)[∂] with order(c) < order(L). In the applications c is either 1 (cf. section 3.1) or 0 (cf. section 3.2). Compute the set of all r ∈ k(x)[∂] such that there exists an l ∈ k(x)[∂] with

Rr + lL = c.  (1)

We call this the mixed equation. An equivalent equation is

RRem(Rr, L) = c.  (2)

Here RRem stands for the remainder after a right hand division. Note that it is sufficient to compute only solutions r of order smaller than order(L) because the other solutions are obtained from these by adding elements of k(x)[∂]L. Therefore, we will restrict ourselves to solutions r with

order(r) < order(L).

The set of solutions of the mixed equation in the case c = 0 is called E_{D}(R, L) in [14]. Singer gives a very interesting application of computing the set E_{D}(f, f), namely the following: If E_{D}(f, f) contains a non-constant element r then he shows how one can compute a non-trivial factor in \overline{k}(x)[∂] of f in an efficient way, cf. section 3.1 in [14]. We will give an example of this in section 3.2. If f is completely reducible then either E_{D}(f, f) contains a non-constant element (and so f can be factored) or E_{D}(f, f) is the set of constants and then f is irreducible. This will be the main application of the mixed equation.

Acknowledgments: I would like to thank J.A. Weil and S.P. Tsarev for useful discussions about these topics. It should also be noted that J.A. Weil has ideas (which have not yet been written down) for similar results about factoring completely reducible operators as well. Both approaches have their own advantages. The advantage of his approach is it can be applied to other problems as well. The advantage of the approach in this paper is that the algorithm is short and, even though the bound we give in proposition 1 is very technical, easier to implement (at least for the most important case L = R and c = 0). A proper comparison between the two methods can be done when both methods are implemented. This will be done somewhere in 1996.

Outline of the paper: The part of this paper that is new is section 4. Sections 2 and 3 are not new, see also [11, 14]. These sections are intended as an introduction for section 4 and to give applications. Section 4.1 can be viewed as an introduction for section 4.2 because the result of section 4.1 is re-done more generally in section 4.2. The algorithm in section 4 for solving the mixed equation consists of two parts:

1. For each singularity p ∈ P^{1}(\overline{k}) compute a bound for the valuation at p of the coefficients in k(x) of r.

2. Solve linear equations over k (a different approach is given in section 4 as well).

Computer algebra systems already have code for solving linear equations, so the only thing that needs to be implemented for solving the mixed equation is the bound in section 4.2.

An implementation of the algorithm (currently only the computation of E_{D}(f, f) is implemented, but it is not much work to adapt this code for the case of a more general mixed equation) is available from

http://www-math.sci.kun.nl/math/compalg/diffop/

2 Preliminaries

This section lists a few facts that we use about differential operators. For a more complete introduction see [14]. In section 4.2 the reader is assumed to be familiar with section 3 in [9].

The least common left multiple f = LCLM(R_{1}, R_{2}) of two operators R_{1} and R_{2} is defined as the monic operator f with minimal order such that R_{1} and R_{2} are right hand factors of f. V (f) = V (R_{1}) + V (R_{2}).

An operator f is called completely reducible if f is the LCLM of irreducible operators R_{1}, …, R_{n}. Note that an irreducible operator is completely reducible as well.

The greatest common right divisor R = GCRD(f_{1}, f_{2}) of two operators f_{1} and f_{2} is defined as the monic operator R with maximal order such that R is a right hand factor of both f_{1} and f_{2}. V (R) = V (f_{1}) ⋂ V (f_{2}). By the Euclidean algorithm (cf. [12]) one can find two operators g_{1} and g_{2} such that

R = g_{1}f_{1} + g_{2}f_{2}

and order(g_{1}) < order(f_{2}). If R = 1 then the pair g_{1}, g_{2} is uniquely determined. This can be shown as follows: If there were two different pairs then the difference h_{1}, h_{2} of these pairs would satisfy the equation h_{1}f_{1} + h_{2}f_{2} = 0 and order(h_{1}) < order(f_{2}). Then V (f_{2}) ⊂ V (h_{1}f_{1}) so f_{1}(V (f_{2})) ⊂ V (h_{1}). However, f_{1} is injective on V (f_{2}) because V (f_{1}) ⋂ V (f_{2}) = V (R) = 0 so order(h_{1}) = dim(V (h_{1})) ≥ dim(f_{1}(V (f_{2}))) = dim(V (f_{2})) = order(f_{2}) which is a contradiction.

A valuation on a ring R is a map v : R → ZZ ⋃{∞} (other additive groups than ZZ, such as Q, are allowed as well) such that v(0) = ∞, v(fg) = v(f) + v(g) for all f, g ∈ R \ {0}. Furthermore v(f + g) ≥ min(v(f), v(g)) and v(f + g) = min(v(f), v(g)) if v(f) ≠ v(g). One can define different valuations on the ring of local (i.e. power series coefficients) differential operators, cf. section 2 in [8]. The valuation of a power series a ∈ k((x)) \ {0} is defined as the smallest n for which the coefficient of x^{n} in a is non-zero. The valuation of a non-zero rational function a ∈ k(x) \ {0} at a point p ∈ P^{1}(\overline{k}) = \overline{k} ⋃\{∞\} is defined as follows: If p = ∞ then the valuation of a at p is the degree of the denominator minus the degree of the numerator. If p ∈ \overline{k} then the valuation of a at p is defined as the integer n for which b = a/(x − p)^{n} has no pole at p and b(p) ≠ 0. Suppose that B_{p} for p ∈ P^{1}(\overline{k}) are given integers for which B_{p} ≠ 0 for only finitely many p and that N ≥ 0 where N = ∑ B_{p}. Then the set of all a in \overline{k}(x) for which the valuation of a at p is ≥ −B_{p} is a \overline{k} vector space of dimension N + 1. The elements of this vector space are of the form n/D where D ∈ \overline{k}(x) can be determined from the B_{p} and where n ∈ \overline{k}[x] is a polynomial of degree ≤ N.

In section 3.1 in [14] Singer gives three methods for computing E_{D}(f, f). We will write here the first method (this method is computationally very costly, however. That is why we give an alternative method in section 4.) because we can draw two useful conclusions from this method. Let n = order(L) and m = order(R). Write r = r_{0}∂^{0} +…+r_{n−1}∂^{n−1} where the r_{i} are indeterminates. Write c = c_{0}∂^{0} +…+c_{n−1}∂^{n−1}. Now RRem(Rr, L) is a k(x) linear combination of the r_{i} and the derivatives of the r_{i}. The highest derivative of r_{i} that appears in RRem(Rr, L) is r_{i}^{(m)}. This r_{i}^{(m)} only appears in the coefficient of ∂^{i} in RRem(Rr, L) and the coefficient of r_{i}^{(m)} in that expression is 1 (we use that R is monic). So equation 2 is equivalent with the following set of differential equations

{r_{i}^{(m)} + … = c_{i}|0 ≤ i ≤ n − 1}

where the dots stand for k(x) linear expressions in r_{j}^{(m′)}, j ∈ {0, …, n − 1}, m′ ∈ {0, …, m−1}. If c ≠ 0 then this system is not homogeneous. To make it homogeneous we add one extra variable z, one equation z′ = 0, and replace the c_{i} by zc_{i}. We can convert this system to a system of first order equations by introducing new indeterminates r_{i,j} for the j’th derivative of r_{i} and adding the equations r_{i,j+1} = r′_{i,j}. This way we obtain a system of equations of the form

Ay = y′  (3)

where A is a nm + 1 by nm + 1 matrix (if c = 0 then nm instead of nm + 1 because then z is not needed) over k(x) and y is a vector consisting of the r_{i,j} and z. Such a matrix differential equation can be reduced to a single equation of order nm+1 (or nm if c = 0) by a cyclic vector computation and can then be solved, see [1] for computing rational solutions of a differential operator. We note the following:

Remark 1. If p ∈ \overline{k} is a regular point for L, R and the c_{i} have no pole at p then the matrix A has no pole at the point p. If the vector y would have a pole of order t at the point p then the pole order of y′ would be t + 1, but the pole order of Ay is ≤ t. This contradicts Ay = y′ so y (and hence the r_{i} as well) have no pole at the point p.

Remark 2. Every basis of the solutions y of equation 3 over k(x) is a basis of solutions y over \overline{k}(x) as well. Hence, for the mixed equation, to compute a basis for the solutions r ∈ \overline{k}(x)[∂] it suffices to compute a basis for the solutions r ∈ k(x)[∂].

3 Applications of the mixed equation

3.1 Computing a complement of V (R)

Let f = LR where f, L, R ∈ k(x)[∂] are monic and suppose there exists a different factorization f = L_{2}R_{2} in \overline{k}(x)[∂] such that V (f) = V (R) ⨁ V (R_{2}). We assume that the operators f, L, R, L_{2} and R_{2} are monic. The greatest common right divisor GCRD(R, R_{2}) of R and R_{2} must be 1 because R and R_{2} have no common non-zero solutions. Then, by the Euclidean algorithm, it follows that

rR + r_{2}R_{2} = 1

for some r, r_{2} ∈ \overline{k}(x)[∂] which are uniquely determined under the condition order(r) < order(R_{2}) (note that order(R_{2}) = order(L)).

The map rR + r_{2}R_{2} = 1 is the identity and the map r_{2}R_{2} acts like the zero map on V (R_{2}) ⊂ V (f). Hence the map rR acts like the projection of V (f) to V (R_{2}). So rR acts like the identity on V (R_{2}). Now R maps V (f) onto V (L), the kernel is V (R) hence R is a bijection from V (R_{2}) to V (L). So r is the inverse bijection from V (L) to V (R_{2}). Hence Rr acts like the identity on V (L), then Rr − 1 maps V (L) to 0 so L must be a right hand factor of Rr − 1. In other words

Rr + lL = 1  (4)

for some l ∈ \overline{k}(x)[∂]. R_{2} can be constructed from r using the equation V (R_{2}) = r(V (L)) as follows: Write z = r(y). We can write the derivatives of z as vectors over \overline{k}(x) on the basis y, y′, …, y^{(n−1)} (the higher order derivatives of y can be simplified using the relation L(y) = 0). Here n = order(L). The n’th derivative of z must be \overline{k}(x) linearly dependent on the lower order derivatives z, z′, …, z^{(n−1)}. Computing this linear dependence gives R_{2}.

So: computing the set of monic operators R_{2} with the property V (f) = V (R) ⨁ V (R_{2}) is equivalent with solving the mixed equation for c = 1. These monic R_{2} are in 1-1 correspondence to the solutions r of this mixed equation. Because of remark 2 in section 2 the existence of such an R_{2} ∈ \overline{k}(x)[∂] is equivalent with the existence of such an R_{2} ∈ k(x)[∂].

3.2 Singer’s Factorization method

In this section we describe a factorization method of Singer (cf. section 3.1 in [14]) and show in an example how to combine Singer’s method with the method from section 4. Note that our f and r are called L and R in [14], we call this f and r to be consistent with the rest of our paper.

Suppose the dimension of E_{D}(f, f) is more than 1. Then we can take an element r ∈ \overline{k}(x)[∂] in E_{D}(f, f) which is not a constant (we always take r of order < order(f)).

Now r is a \overline{k} linear map from V (f) to V (f). We can compute a basis of V (f) by computing formal solutions of f at a point (this is easiest at a regular point). Then compute the matrix of the map r in this basis and take an eigenvalue a ∈ \overline{k}. Then GCRD(f, r − a) ∈ \overline{k}(x)[∂] is a non-trivial factor of f. Note that for f ∈ k(x)[∂] with dim(E_{D}(f, f)) > 1 the only algebraic extension that is used over k to compute a right hand factor is the eigenvalue a. Example: Given is the following differential operator

f = ∂^{4} + \frac{6}{x}∂^{3} + \frac{2(x^{2} − 1)}{x^{4}}∂^{2} − \frac{2(3x^{2} − 1)}{x^{5}}∂ + \frac{1}{x^{8}}

which is the LCLM of two irreducible operators in \overline{Q}(x)[∂]. The problem is to find a right hand factor of f. E_{D}(f, f) is the set of operators r of order < order(f) such that

fr + lf = 0

for some l ∈ \overline{Q}(x)[∂]. In section 4 we will compute the solutions r of order < order(f) of this mixed equation. The solution space is 2 dimensional. Hence we can choose a non-constant solution:

r = −x^{5}∂^{3} − x^{4}∂^{2} + 2x^{3}∂ + x∂.

Now we need a basis for V (f). We choose a basis b_{1}, …, b_{4} ∈ Q((x − 1)) of formal solutions at the point x = 1. These b_{i} are uniquely determined by requiring that b_{i} mod (x − 1)^{4} is (x − 1)^{i−1}. The operator r acts on this basis as

( \begin{matrix}0 & 3 & −2 & −6 \\ 1 & 3 & −2 & 0 \\ −3/2 & 1/2 & 2 & −3 \\ 13/6 & −5/3 & 1/3 & 7 \end{matrix} ) .

The eigenvalues of this matrix are 3 + \sqrt{2} and 3 − \sqrt{2}. From the eigenvalue 3 + \sqrt{2} we obtain the following right hand factor of f

GCRD(f, r − 3 − \sqrt{2}) = ∂^{2} + \frac{\sqrt{2}}{x}∂ − \frac{1}{x^{4}} .

The other eigenvalue 3 − \sqrt{2} gives the conjugate over Q of this factor.

Note that dim(E_{D}(f, f)) > 1 implies that f is reducible but not that f is completely reducible. For example if f = (∂ + 1/x)∂ then r = x∂ ∈ E_{D}(f, f) (giving the right hand factor ∂) but f is not the LCLM of irreducible operators.

4 Solving the mixed equation

Write

r = ∑_{i=0}^{order(L)−1} r_{i}∂^{i}.

In section 2 we have seen that the r_{i} can not have a pole at a point p ∈ \overline{k} if R, L and c have no pole at p. In section 4.1 (lemma 1) and section 4.2 (proposition 1 and the comments after proposition 1) we show how to compute a bound for the valuation of each r_{i} in the remaining places (i.e. the point ∞ and the points in \overline{k} where R, L or c has a pole). Then we can write

r_{i} = \frac{n_{i}}{D_{i}}

with n_{i} ∈ k[x] of degree ≤ N_{i} and D_{i} ∈ k(x) where N_{i} and D_{i} are computed from the bounds. So given the bounds, we only need to determine the polynomials n_{i}. This can be done in several different ways:

• Approach 1. Write the n_{i} as a polynomials in x of degree N_{i} with undetermined coefficients. Then substitute r in the mixed equation and find linear equations for these undetermined coefficients. Solving these equations gives the solutions of the mixed equation.

• Approach 2.

Let L_{1} be an operator such that V (L_{1}) = c(V (L)). If c = 0 then take L_{1} = 1, if c = 1 then L_{1} = L. In the remaining cases L_{1} can be obtained as follows: write z = c(y) where y is a solution of L. Denote n = order(L). Using L(y) = 0 we write the derivatives of z as k(x) linear expressions in y, y′, …, y^{(n−1)}. Computing a k(x) linear dependence between z, z′, …, z^{(n)} gives L_{1}. Compute a basis for V (L) and a basis z_{1}, …, z_{m} for V (L_{1}R) (for example a basis of formal solutions at a regular point). Let y ∈ V (L). Then Rr(y) = c(y) ∈ V (L_{1}) so r(y) ∈ V (L_{1}R). Assume that all the D_{i} are equal to one polynomial D (if they are not polynomials we can take the numerator, and if they are not equal we can replace the D_{i} by the least common multiple). Then

r(y) = \frac{1}{D} ∑ n_{i}∂^{i}(y) ∈ V (L_{1}R)

so

∑ n_{i}y^{(i)} = C_{1}Dz_{1} + … + C_{m}Dz_{m}.  (5)

Here the n_{i} and C_{i} are polynomials in x with a bounded degree (the degree of the C_{i} is 0). The y^{(i)} and Dz_{i} are power series in x − p (if we computed the formal solutions at a regular point p). For each basis element y of V (L) we obtain an equation of this form. The problem of computing all polynomials n_{i} and C_{i} satisfying the given conditions on their degrees, such that equation (5) holds up to a given accuracy a ∈ N (i.e. modulo (x − p)^{a}) is handled efficiently by the Beckermann-Labahn algorithm, cf. [6, 3]. A good guess for the accuracy a that is required to obtain solutions of the mixed equation is the sum of the degrees of the n_{i} plus the number of n_{i} and C_{i}. If we took the accuracy a too small then we find too many solutions, i.e. we find a basis b_{1}, …, b_{t} such that the solutions of the mixed equation form a subspace of the vector space spanned by b_{1}, …, b_{t}. Then we can pick out the correct solutions from this vector space as follows: Substitute r = c_{1}b_{1} + … c_{t}b_{t}, where the c_{i} are variables, in the mixed equation (2) and solve linear equations to find the right c_{1}, …, c_{t}.

• Approach 3. One can compute the n_{i} from the formal solutions r_{i} ∈ k((x − p)) of the mixed equation where p is a regular point, in a way that is similar to the way that the rational solutions of a differential operator are obtained from the formal solutions in [1]. This way we have to solve linear equations over k in order(L) · order(R) variables.

4.1 The local bound problem, regular singular case

Define δ = x∂. The ring k(x)[∂] is a subring of the ring k((x))[∂] = k((x))[δ]. We can define a valuation (see also section 2 in [8])

v_{0} : k((x))[δ] → ZZ ⋃{∞}

as follows: For non-zero f = ∑_{i,j} f_{i,j}x^{i}δ^{j} the valuation v_{0}(f) is the smallest i for which f_{i,j} ≠ 0 for some j. For this f we can define the Newton polynomial for slope 0 (see also section 3.4 in [8]) using a variable T as follows

N_{0}(f) = ∑_{j} f_{v_{0}(f),j}T^{j} ∈ k[T].

The substitution map

S_{T =T +i} : k[T] → k[T]

is a k homomorphism defined by T → T + i.

We recall a few facts about the Newton polynomial N_{0}

• degree(N_{0}(f)) = order(f) if and only if f is regular singular.

• The roots of N_{0}(f) in \overline{k} are called the exponents of f.

• For all L, R in k((x))[δ] we have

N_{0}(L · R) = S_{T =T +v_{0}(R)}(N_{0}(L)) · N_{0}(R).

Note that we assumed that L and R are monic in k(x)[∂]. Then they are not monic when considered as elements of k((x))[δ].

The local bound problem is the following: Given are R, L and c in k((x))[δ] with order(c) < order(L). Give a lower bound for v_{0}(r) (or give for each coefficient of r in k((x)) a lower bound separately) for all solutions r ∈ k((x))[δ] of the mixed equation, i.e. for all r with order(r) < order(L) for which there exists an l ∈ k((x))[δ] such that Rr + lL = c.

Lemma 1 Bound for v_{0}(r) in the regular singular case. Let R, r, l, L and c in k((x))[δ] with Rr + lL = c, order(c) < order(L), order(r) < order(L), R ≠ 0, L ≠ 0 and r ≠ 0. Assume that L is regular singular. Then

c ≠ 0 and v_{0}(r) ≥ v_{0}(c) − v_{0}(R)

or

gcd(S_{T =T +v_{0}(r)}(N_{0}(R)), N_{0}(L)) ≠ 1.  (6)

Proof: If v_{0}(r) ≥ v_{0}(c) − v_{0}(R) then v_{0}(c) must be finite and so c ≠ 0. Now assume v_{0}(r) < v_{0}(c) − v_{0}(R). Then v_{0}(Rr) < v_{0}(c) = v_{0}(Rr + lL). This is only possible if the lowest power of x in Rr and lL cancel against each other. Hence

N_{0}(Rr) + N_{0}(lL) = 0.

The assumption that L is regular singular means

degree(N_{0}(L)) = order(L).

Apply the multiplication formula for the Newton polynomials

0 = N_{0}(Rr) + N_{0}(lL)
= S_{T =T +v_{0}(r)}(N_{0}(R)) · N_{0}(r) + S_{T =T +v_{0}(L)}(N_{0}(l)) · N_{0}(L).

Because degree(N_{0}(r)) < order(L) = degree(N_{0}(L)) equation 6 follows.

□

Note that equation 6 can hold for only finitely many integers v_{0}(r). So the minimum of these integers (and the integer v_{0}(c) − v_{0}(R) if c ≠ 0) is a lower bound for v_{0}(r). The bound can be computed from v_{0}(c), v_{0}(R), N_{0}(L) and N_{0}(R).

4.2 The local bound problem, general case

We can generalize the valuation v_{0} to \overline{k((x))}[δ] as follows: If f ∈ \overline{k((x))}[δ] then f is an element of \overline{k}((x^{1/n}))[δ] for some n ∈ N. Write f = ∑_{i} x^{i}f_{i} where f_{i} ∈ \overline{k}[δ] and where the sum is taken over i ∈ \frac{1}{n}ZZ. If f ≠ 0 then v_{0}(f) is defined as the smallest i for which f_{i} ≠ 0 and the Newton polynomial N_{0}(f) is defined as this f_{i} (with δ replaced by the variable T).

Define the set E = ⋃_{n} \overline{k}[x^{−1/n}]. In section 3.3 in [9] we have defined a valuation from a subset V_{∗} of the universal extension V to the set E

v : V_{∗} → E

The map

Exp : E → V_{∗}

is defined as Exp(e) = exp(∫ \frac{e}{x}dx). We have v(Exp(e)) = e.

For a definition of the canonical list see [9]. The canonical exponential parts are the elements of this list. The multiplicity of a canonical exponential part is the number of times that it appears in the canonical list. For e ∈ \overline{k((x))} the substitution map

S_{e} : \overline{k((x))}[δ] → \overline{k((x))}[δ]

is defined as the \overline{k((x))} homomorphism given by S_{e}(δ) = δ + e.

Definition 1 Let e ∈ E and f ∈ \overline{k((x))}[δ] \ \{0\}. Then the multiplicity ν_{e}(f) of the canonical exponential part e in f is defined as the multiplicity of the root 0 in N_{0}(S_{e}(f)).

Note: ν_{e}(f) is not the same as the multiplicity of the exponential part μ_{e}(f), cf. section 3 in [9]. The exponential parts are the canonical exponential parts modulo a certain equivalence, hence their multiplicity μ_{e}(f) is ≥ the multiplicity ν_{e}(f) of the canonical exponential part. The sum of ν_{e}(f) taken over all e ∈ E is the number of elements of the canonical list which is order(f).

The canonical exponential parts are a generalization of exponents. The exponents of an operator f are those canonical exponential parts which are in \overline{k}. An operator is regular singular if and only if all canonical exponential parts are exponents, i.e. if they are elements of \overline{k}.

Our approach for the general (i.e. not necessarily regular singular) case is quite technical. To explain the idea behind this approach we will first reformulate the approach of the previous section into the terminology of exponents instead of Newton polynomials. Then we can generalize by replacing the exponents by canonical exponential parts. If

Rr + lL = c and c = 0

and e is an exponent of L then e is an exponent of Rr as well. If the multiplicity of the exponent e in r is smaller than ν_{e}(L) then (the proof follows later in the more general case) v_{0}(r) + e is an exponent of R. This must happen for at least one exponent e of L because if L is regular singular then the number of exponents e (counting with multiplicity) is more than the order of r. By comparing the exponents of R and L, and taking the smallest possible integer difference, we obtain a bound for v_{0}(r). If c ≠ 0 we have to consider the possibility v_{0}(Rr) ≥ v_{0}(c) as well.

If L is not regular singular (L is irregular singular) then the number of exponents e is not necessarily larger than order(r). And so the argument that v_{0}(r) + e must be an exponent of R for some exponent e of L does not hold anymore. This problem can be fixed by using the canonical exponential parts instead of exponents. We will first relate the multiplicity ν_{e}(f) to a property degl of formal solutions. Denote

\overline{V}_{0} = \overline{k((x))}[log(x)]

as in [8, 9].

Definition 2 For an element y ∈ \overline{V}_{0} define degl(y) as follows: y can be written as y = ∑_{i,j} a_{ij}x^{i}log(x)^{j} where the sum is taken over j ∈ N and i ∈ \frac{1}{n}ZZ for some n ∈ N. Then degl(y) is the maximal j for which a_{v(y),j} ≠ 0.

Every element y ∈ V_{∗} is of the form Exp(e)z for some e ∈ E and z ∈ \overline{V}_{0}. Then degl(y) is defined as degl(z).

Lemma 2 Let e ∈ \overline{k}, f = δ − e + s with s ∈ \overline{k((x))} and v(s) > 0 and y ∈ \overline{V}_{0} \ \{0\}. If e ≠ v(y) then

v(f(y)) = v(y) and degl(f(y)) = degl(y).

If e = v(y) and degl(y) > 0 then

v(f(y)) = v(y) and degl(f(y)) = degl(y) − 1

and if e = v(y) and degl(y) = 0 then

v(f(y)) > v(y).

The proof of the lemma is easy, we skip it. Note that f(y) = 0 is only possible in the case e = v(y) and degl(y) = 0.

Lemma 3 Let f ∈ \overline{k((x))}[δ] \ \{0\}, y ∈ \overline{V}_{0} \ \{0\} and d = degl(y). Let e = v(y) ∈ Q and d′ = ν_{e}(f). If d′ ≤ d then

v(f(y)) = v_{0}(f) + v(y) and degl(f(y)) = d − d′.

If d′ > d then

v(f(y)) > v_{0}(f) + v(y)

Note that f(y) = 0 is only possible in the case d′ > d. Proof: Factor (cf. section 5 in [8]) f as f = L · (δ − e_{1} + s_{1}) ⋯ (δ − e_{n} + s_{n}) where v(s_{i}) > 0, e_{1}, …, e_{n} ∈ \overline{k} are the exponents of f and L has no regular singular factor (i.e. L has no slope 0 in the Newton polygon). Then v_{0}(f) = v_{0}(L). Denote z = (δ − e_{1} + s_{1}) ⋯ (δ − e_{n} + s_{n})(y) ∈ \overline{V}_{0}. Now either v(z) ∈ Q or z = 0. d′ is the number of i for which e_{i} = v(y).

Assume z ≠ 0. Write L = L_{0}δ^{0} + … + L_{m}δ^{m}. Now v(L_{i}) > v(L_{0}) for i > 0 (because L has no slope 0) and v(δ^{i}(z)) ≥ v(z) so v(L_{i}δ^{i}(z)) > v(L_{0} · z). Hence v(L(z)) = v(L_{0} · z) and degl(L(z)) = degl(L_{0} · z). Now f(y) = L(z) so v(f(y)) = v(L_{0} · z) = v(L_{0}) + v(z) = v_{0}(L) + v(z) = v_{0}(f) + v((δ − e_{1} + s_{1}) ⋯ (δ − e_{n} + s_{n})(y)) and degl(f(y)) = degl(L_{0} · z) = degl(z) = degl((δ − e_{1} + s_{1}) ⋯ (δ − e_{n} + s_{n})(y)). Now the lemma follows by repeated use of the previous lemma.

□

Lemma 4 Let f ∈ k((x))[δ]. Then f has a solution y in V_{∗} with degl(y) = d and v(y) = e if and only if e, e, …, e (d + 1 times) is a sublist of the canonical list of f (in other words: ν_{e}(f) > d).

The canonical list is defined up to permutations, so we can view it as a set with multiplicities. A solution in V_{∗} with degl(y) = d means that the multiplicity of v(y) in the canonical list is at least d + 1, i.e. ν_{v(y)}(f) ≥ d + 1. Note the following consequence of the lemma: If there exists a solution y with v(y) = e and degl(y) > 0 then there exists a solution z with v(z) = e and degl(z) = degl(y)−1. This can easily be shown in a different way as well, take z = y − S_{log}(y) where S_{log} is the map that replaces log(x) by log(x) + 1, cf. section 9 in [8]. Proof: Let y ∈ V_{∗} with degl(y) = d, v(y) = e and f(y) = 0. Write z = Exp(−e)y, so v(z) = 0 and z ∈ \overline{V}_{0} \ \{0\}. We have V (f) = Exp(e) · V (S_{e}(f)). So S_{e}(f)(z) = 0 and hence by lemma 3 it follows that ν_{0}(S_{e}(f)) > degl(z) = d. Since ν_{e}(f) = ν_{0}(S_{e}(f)) one part of the lemma follows.

Now suppose ν_{e}(f) > d. We must prove that f has a solution in V_{∗} with valuation e and degl d. Let R ∈ k((x))[e, δ] be the right hand factor of S_{e}(f) of maximal order which is semi-regular (cf. section 3.2 in [9]) over k((x))[e]. Now ν_{0}(R) = ν_{0}(S_{e}(f)) = ν_{e}(f) > d.

It is sufficient to prove that R has a solution y with valuation 0 and degl d, because then Exp(e)y is a solution of f with the desired property. Section 8.1 in [8] gives a recursive algorithm for computing a basis of solutions of R. This algorithm makes repeated use of integration s_{i} = ∫ \frac{a_{i}}{x} dx. In this integration process (we take the constant term in the integral equal to 0) we have v(s_{i}) = v(a_{i}). Furthermore degl(s_{i}) = degl(a_{i}) if v(a_{i}) ≠ 0 and degl(s_{i}) = degl(a_{i}) + 1 if v(a_{i}) = 0. Using these relations and induction with respect to the order of R it follows that the algorithm in section 8.1 in [8] produces a solution y with valuation e and degl(y) = j for every exponent e of R and every integer j with 0 ≤ j < ν_{e}(R).

□

Lemma 5 Let f ∈ k((x))[δ] be of order n and e ∈ E. Let v′ be 0 if e = 0 and v(e) otherwise (so v′ ∈ Q and v′ ≤ 0). Let d = v_{0}(S_{e}(f)). Then d is an integer divided by the ramification index of e. The coefficient of δ^{i} in f has valuation ≥ d + (n − i)v′.

The proof of the lemma is easy, we skip it. The ramification index of e is defined as the smallest positive integer n such that e ∈ \overline{k}((x^{1/n})).

Proposition 1 Bound for v_{0}(S_{e}(r)): Let Rr + lL = c where R, r, l, L, c in k((x))[δ] with R ≠ 0, L ≠ 0, r ≠ 0, order(c) < order(L) and order(r) < order(L).

Take an e ∈ E with ν_{e}(L) > ν_{e}(r). Take a solution y ∈ V_{∗} of L with v(y) = e and degl(y) = ν_{e}(L) − 1. Let M = ∞ if c(y) = 0 and M = v(c(y)) − e ∈ Q if c(y) ≠ 0. Then

c(y) ≠ 0 and v_{0}(S_{e}(r)) = M − v_{0}(S_{e}(R))

or v_{0}(S_{e}(r)) + e is in the canonical list of R.

Note: It is not a priori known which e we can take, i.e. which e satisfies equation (9). Also note that ν_{e}(L) > ν_{e}(r) implies that ν_{e}(L) > 0 in other words: e is in the canonical list of L.

In the two applications in section 3 we have c = 0 or c = 1. If c = 0 then M = ∞ so then the first case in the proposition can not occur. If c = 1 then M = 0. So in both applications the proposition can be used without computing a solution y with the desired property. Proof:

Rr(y) = (Rr + lL)(y) = c(y).

Denote z = Exp(−e)y and w = S_{e}(r)(z). Then v(z) = 0 and degl(z) = degl(y) = ν_{e}(L) − 1 ≥ ν_{e}(r) = ν_{0}(S_{e}(r)). Now S_{e}(R)S_{e}(r) + S_{e}(l)S_{e}(L) = S_{e}(c) and z is a solution of S_{e}(L) so

S_{e}(R)(w) = S_{e}(R)(S_{e}(r)(z)) = S_{e}(c)(z)  (7)

Now S_{e}(c)(z) = Exp(−e)c(y) so

v(S_{e}(R)(w)) = v(c(y)) − e = M.

By lemma 3 and degl(z) ≥ ν_{0}(S_{e}(r)) it follows that

v(w) = v(S_{e}(r)(z)) = v_{0}(S_{e}(r)) + v(z) = v_{0}(S_{e}(r)).  (8)

According to lemma 3 and equation 7 there are two possibilities (if M = ∞ then the first case can not occur)

v(S_{e}(R)(w)) = v_{0}(S_{e}(R)) + v(w)

or ν_{v(w)}(S_{e}(R)) > degl(w) which implies that v(w) is an element of the canonical list of S_{e}(R). The latter case implies that v(w)+e is in the canonical list of R. So v_{0}(S_{e}(r)) = v(w) = v(S_{e}(R)(w)) − v_{0}(S_{e}(R)) = M − v_{0}(S_{e}(R)) or v_{0}(S_{e}(r)) + e = v(w) + e is in the canonical list of R.

□

Note that both two cases imply a lower bound for v_{0}(S_{e}(r)). Since we do not know which of these two cases holds (unless M = ∞ then the first case can not occur) we have to take the minimum of these two bounds to obtain a lower bound for v_{0}(S_{e}(r)). In the case where v_{0}(S_{e}(r)) + e is in the canonical list of R we obtain a lower bound for v_{0}(S_{e}(r)) by taking the minimal m ∈ \frac{1}{ram(e)}ZZ for which m + e is in the canonical list of R. Here ram(e) is the ramification index of e. Then by lemma 5 we obtain a lower bound for the valuation of the coefficients of r.

The order of an operator equals the sum of the ν_{e} taken over all e ∈ E. Since order(L) > order(r) we must have

ν_{e}(L) > ν_{e}(r)  (9)

for at least one e ∈ E. Note, however, that we do not know for which e equation 9 holds. So to obtain a lower bound for the valuations of the coefficients of r we must take the minimum of these lower bounds for all canonical exponential parts e of L.

Example, continued from section 3.2: Now we will use the bound to compute the rational solutions r of the mixed equation in the example of section 3.2. We can write

r = \frac{n_{3}}{D_{3}} ∂^{3} + … + \frac{n_{0}}{D_{0}} ∂^{0}.

The only singularities of f are x = 0 and x = ∞. The point x = 0 is an irregular singularity so we must compute the canonical list:

α − \frac{1}{x}, 2 − α − \frac{1}{x}, α + \frac{1}{x}, 2 − α + \frac{1}{x}.

Here α is a root of the polynomial 1 − 4Z + 2Z^{2} (note that it is not necessary to compute all canonical exponential parts, it suffices to compute them up to conjugation). Now the smallest possible difference between canonical exponential parts which is an integer divided by the ramification index is 0. So we have v_{0}(S_{e}(r)) ≥ 0 for some e in the canonical list. Then by lemma 5 it follows that the coefficient of δ^{i} in r (here r localized at the point x = 0, in δ notation instead of ∂ notation, cf. section 3.4 in [9]) has valuation ≥ i − 3. Now we should convert this bound for the δ notation to a bound in ∂ notation. The result is that r_{i} = n_{i}/D_{i} has valuation ≥ i − 3 + i at the point x = 0. So we can take D_{i} = x^{3−2i} (D_{i} is not a polynomial if i ≥ 2, however. In these cases the notion of the degree of D_{i} is problematic. But then we can simply interpret degree(D_{i}) as −1 times the valuation of D_{i} at the point infinity).

Now we want a lower bound for the valuation of r_{i} at infinity (i.e. an upper bound for degree(n_{i})−degree(D_{i})). The operator f is regular singular at infinity and the Newton polynomial is T^{4} − 5T^{2} + 2T = T (T − 2)(T^{2} + 2T − 1). Then by lemma 1 it follows that v_{0}(l_{∞}(r)) ≥ −2. Here l_{∞}(r) is r localized at infinity, cf. section 3.4 in [9]. We have to convert this to a bound for the valuation of r_{i} at infinity. The result is that the valuation of r_{i} at infinity is ≥ −2 − 2i. This means degree(n_{i}) − degree(D_{i}) ≤ 2 + 2i. Hence degree(n_{i}) ≤ 2 + 2i + degree(D_{i}) = 5. So we can write r with 4 · (5 + 1) undetermined coefficients. Twenty-four indeterminates is not very much so approach 1 in section 4, solving linear equations, will be efficient enough to be able to handle this example. These linear equations are obtained from RRem(fr, f) = 0.

References

[1] A. Abramov, M. Bronstein, M. Petkovšek On polynomial solutions of linear operator equations. Proceedings ISSAC 95, ACM Press, 290-296.

[2] B. Beckermann, G. Labahn, A uniform approach for Hermite Padé and simultaneous Padé Approximants and their Matrix-type generalizations, Numerical Algorithms, 3 (1992), pp. 45-54

[3] B. Beckermann, G. Labahn, A uniform approach for the fast computation of Matrix-type Padé approximants, SIAM J. Matrix Analysis and Applications (1994), 804-823.

[4] E. Beke Die Irreduzibilität der homogenen linearen Differentialgleichungen, Math. Ann. 45, (1894), 278-294.

[5] M. Bronstein Linear Ordinary Differential Equations: breaking through the order two barrier. Proceedings ISSAC 92, ACM Press, 42-48.

[6] H. Derksen An algorithm to compute generalized Padé-Hermite forms Manuscript.

Available by ftp at daisy.math.unibas.ch in /pub/hderksen/pade.dvi

[7] P. Hendriks, M. v.d. Put Galois action on solutions of a differential equation, Preprint, to appear in J. Symb. Comput.

[8] M. van Hoeij, Formal Solutions and Factorization of Differential Operators with Power Series Coefficients, University of Nijmegen Report nr. 9528, submitted to J. Symb. Comput, available at

http://www-math.sci.kun.nl/math/compalg/diffop/

[9] M. van Hoeij, Factorization of Differential Operators with Rational Functions Coefficients, University of Nijmegen Report nr. 9552, submitted to J. Symb. Comput, available at

http://www-math.sci.kun.nl/math/compalg/diffop/

[10] A. Loewy, Über vollständig reduzible lineare homogene Differentialgleichungen, Math. Ann., 62, (1906), 89-117.

[11] O. Ore, Formale Theorie der linearen Differentialgleichungen (Zweiter Teil), J. für d. Reine u. angew. Math., 168, (1932), 233-252.

[12] O. Ore, Theory of non-commutative polynomial rings, Ann. of Math. 34 pp. 480-508 (1933)

[13] F. Schwarz A factorization Algorithm for Linear Ordinary Differential Equations. Proceedings of ISSAC 89, ACM Press, 17-25.

[14] M.F. Singer, Testing Reducibility of Linear Differential Operators: A Group Theoretic Perspective, To appear in J. of Appl. Alg. in Eng. Comm. and Comp.

[15] S.P. Tsarev, On the problem of factorization of linear ordinary differential operators, Programming & computer software, 1994, v. 20, 1, p. 27-29.

[16] J.A. Weil, Constantes et polynômes de Darboux en algèbre différentielle : application aux systèmes différentiels linéaires, PhD dissertation, École Polytechnique, 1995.