Solving linear ordinary differential equations over C(x, e^{∫ f(x)dx})

Manuel Bronstein INRIA – Projet Café 2004 Route des Lucioles, BP 93 F-06902 Sophia Antipolis Cedex

                        Manuel.Bronstein@sophia.inria.fr

Anne Fredet Laboratoire GAGE Ecole Polytechnique F-91128 Palaiseau Cedex

                                                                    fredet@ariana.polytechnique.fr

Abstract

We describe a new algorithm for computing the solutions in

F = C(x, e^{∫ f(x)dx}) of linear ordinary differential equations with coefficients in F. Compared to the general algorithm of [9], our algorithm avoids the computation of exponential solutions of equations with coefficients in C(x), as well as the solving of linear differential systems over C(x). Our method is effective and has been implemented.

1 Introduction

We recall that a differential field is a field K together with a map ′ : K → K satisfying (a + b)′ = a′ + b′ and (ab)′ = a′b + ab′ for any a, b ∈ K. In that case the constant field of K is defined as

Const(K) = {a ∈ K s.t. a′ = 0} .

All fields in this paper are of characteristic 0.

Definition 1 Let K be a differential field and L = ∑_{i=0}^{n} a_{i}∂^{i} be a linear ordinary differential operator with coefficients in K. A rational solution of L is any element y ∈ K such that Ly = 0, and an exponential solution of L

is any element u ∈ K such that Le^{∫ u} = 0.

Although our main goal is to compute the rational solutions of such equations, it turns out to be more useful to solve the following, slightly more general problem, instances of which are generated by our algorithm even when starting with a homogeneous equation.

Definition 2 We say that we can effectively solve parameterized linear differential equations over K if given any linear ordinary differential operator L = ∑_{i=0}^{n} a_{i}∂^{i} with coefficients in K and g_{1}, …, g_{m} in K, we can effectively find h_{1}, …, h_{r} in K and a matrix M with r + m columns and coefficients in Const(K) such that Ly = ∑_{j=1}^{m} c_{j}g_{j} for y ∈ K and c_{1}, …, c_{m} ∈ Const(K) if and only if y = ∑_{k=1}^{r} y_{k}h_{k} where y_{1}, …, y_{r} ∈ Const(K) and M(y_{1}, …, y_{r}, c_{1}, …, c_{m})^{T} = 0.

In [9], Singer gave an algorithm that reduced the problem of solving parameterized linear differential equations over a Liouvillian extension E = K(t) of K to solving such equations over K and, in the case where t is exponential over K, to compute all the exponential solutions of linear ordinary differential equations over K. His method, although algorithmic, suffers from several difficulties, mainly the requirement to compute exponential solutions in Liouvillian towers, and the reduction to linear differential systems rather than scalar equations once the order and degree of (Laurent)– polynomial solutions are bounded. To our knowledge, that method has not been implemented, even for the case where K = ℚ(x). In this paper, we describe an alternative algorithm that does not require computing exponential solutions

over C(x) when E is of the form C(x, e^{∫ f(x)dx}), ′ = d/dx, and E is not an algebraic function field. Furthermore, for transcendental exponential extensions of an arbitrary differential field K, our algorithm directly produces scalar equations over K (without increasing the order) rather than differential systems.

Let now C be an arbitrary field^{1}, x be transcendental over C and E = C(x, t) be a differential extension of (C(x), d/dx) where dt/dx = ft for some f ∈ C(x),

i.e. t = e^{∫ f(x)dx}. We write′ for d/dx throughout this paper. There are three separate steps in solving the parameterized equation Ly = ∑_{j=1}^{m} c_{j}g_{j} for y ∈ E and c_{1}, …, c_{m} ∈ C:

1. compute the “denominator” of y, i.e. compute q ∈ C(x)[t] such that qy ∈ C(x)[t, t^{−1}] for any solution y. Replacing y by Y/q in the parameterized equation, this reduces the problem to solving a similar equation for Y ∈ C(x)[t, t^{−1}];

2. bound the degree and order of Y, i.e. compute integers δ ≤ γ such that any solution Y ∈ C(x)[t, t^{−1}] must be of the form

Y = y_{δ}t^{δ} + … + y_{γ}t^{γ} where y_{δ}, …, y_{γ} ∈ C(x) ,

3. compute the solution Y knowing δ and γ.

In the case when t is transcendental over C(x) and the constant subfield of E is C (we will eventually ensure those conditions), then the first step is done using either Lemma 3.2 of [9] or Theorem 5.5 of [3], so we address only the last two steps in this paper. The third step is done in [9] by replacing Y by ∑_{i=δ}^{γ} y_{i}t^{i} in the equation and computing the rational solutions of the resulting linear differential system for the y_{i}’s.

1For the algorithm to be complete, we eventually require that a vector space basis for C over ℚ be known.

The rest of this paper is organized as follows: some preliminary results about exponential extensions are given in section 2, then a specialized bounding algorithm for exponential extensions of C(x) is described in section 3, and a new algorithm for the third step in arbitrary exponential extensions is described in section 4. Notation: for any field k, Z transcendental over k, and any f ∈ k(Z)^{∗}, let ∑_{i≥0} a_{i}Z^{d+i} and ∑_{j≥0} b_{j}Z^{e−j} be the Laurent expansions of f at Z = 0 and Z = ∞ respectively where a_{0} ≠ 0 ≠ b_{0}. We write d = ν_{Z}(f) (the order of f), e = deg_{Z}(f) (the degree of f) and b_{0} = lc_{Z}(f) (the leading coefficient of f). For any nonzero linear ordinary differential operator with coefficients in k(Z), L = f_{n}∂^{n}+f_{n−1}∂^{n−1}+…+f_{0}, we write deg_{Z}(L) = max_{0≤i≤n}(deg_{Z}(f_{i})) and ν_{Z}(L) = min_{0≤i≤n}(ν_{Z}(f_{i})). We omit the subscripts when the context is clear (to prevent confusion, we never use deg(L) for the order of L).

2 Preliminaries

We group in this section some known results about exponential extensions that will be used by our algorithm. Let (K,′) be an arbitrary differential field.

Lemma 1 For any a_{0}, …, a_{n}, a, b ∈ K with b ≠ 0,

(∑_{i=0}^{n} a_{i}∂^{i})(ab) = b(∑_{i=0}^{n} a_{i}(∂ + \frac{b′}{b} )^{i})a .

Proof. The formula is trivial for n = 0, so assume it holds for operators of order up to n − 1. We have

b(∑_{i=0}^{n} a_{i}(∂ + \frac{b′}{b} )^{i})a

= ba_{0}a + b(∑_{i=0}^{n−1} a_{i+1}(∂ + \frac{b′}{b} )^{i})(a′ + \frac{b′}{b} a)

= a_{0}ab + (∑_{i=0}^{n−1} a_{i+1}∂^{i})(a′b) + (∑_{i=0}^{n−1} a_{i+1}∂^{i})(ab′)

= (∑_{i=0}^{n} a_{i}∂^{i})(ab) . □

Lemma 2 Let E be a differential extension of K and t ∈ E^{∗} be such that t′/t ∈ K. Then, for any linear ordinary differential operator L with coefficients in K,

\frac{L(at^{m})}{t^{m}} ∈ K for any integer m and any a ∈ K .

Proof. Write L = ∑_{i=0}^{n} a_{i}∂^{i} where the a_{i}’s are in K. By Lemma 1,

L(at^{m}) = (∑_{i=0}^{n} a_{i}∂^{i})(at^{m})

= t^{m} [(∑_{i=0}^{n} a_{i}(∂ + m \frac{t′}{t} )^{i})a] = t^{m}b

where b ∈ K since t′/t ∈ K. □

We now introduce the construction of “projected” operators: let E be a differential extension of K and t ∈ E be transcendental over K. We can uniquely decompose any linear ordinary differential operator L with coefficients in K[t] as

L = ∑_{k=ν(L)}^{deg(L)} t^{k}L_{k}  (1)

where the L_{k}’s are linear ordinary differential operators with coefficients in K (the “projections” of L with respect to t). Note that L_{ν(L)} ≠ 0 ≠ L_{deg(L)}.

The following, whose proof is present in [9], is the basis of all the algorithms for bounding the order and degree of (Laurent)–polynomial solutions.

Proposition 1 Let E be a differential extension of K, t ∈ E be transcendental and exponential over K, L be a linear ordinary differential operator with coefficients in K[t] and b_{1}, …, b_{m} ∈ K(t). If there are integers δ ≤ γ, y_{δ}, …, y_{γ} ∈ K and c_{1}, …, c_{m} ∈ Const(K) such that y_{δ} ≠ 0 ≠ y_{γ} and

L(y_{δ}t^{δ} + … + y_{γ}t^{γ}) = c_{1}b_{1} + … + c_{m}b_{m} ,  (2)

then,

(i) either δ ≥ min_{1≤j≤m}(ν(b_{j})) − ν(L) or L_{ν(L)}(y_{δ}t^{δ}) = 0;

(ii) either γ ≤ max_{1≤j≤m}(deg(b_{j})) − deg(L) or L_{deg(L)}(y_{γ}t^{γ}) = 0;

where the minimum and maximum are taken over the nonzero b_{j}’s.

Proof. Using the decomposition (1) we have

L(y_{δ}t^{δ} + … + y_{γ}t^{γ}) = ∑_{k=ν(L)}^{deg(L)} t^{k}L_{k}(∑_{i=δ}^{γ} y_{i}t^{i})

= ∑_{k=ν(L)}^{deg(L)} ∑_{i=δ}^{γ} z_{ki}t^{k+i}  (3)

where z_{ki} = L_{k}(y_{i}t^{i})/t^{i} is in K by Lemma 2.

(i) Suppose that δ < min_{1≤j≤m}(ν(b_{j}))−ν(L). It then follows from (2) and (3) that z_{ν(L)δ} = 0, hence that L_{ν(L)}(y_{δ}t^{δ}) = 0. (ii) Suppose that γ > max_{1≤j≤m}(deg(b_{j})) − deg(L). It then follows from (2) and (3) that z_{deg(L)δ} = 0, hence that L_{deg(L)}(y_{γ}t^{γ}) = 0. □

This reduces the problem to bounding the integers m such that L(at^{m}) = 0 for a ∈ K^{∗} and L an operator with coefficients in K. This problem is solved for the case K = C(x) and ′ = d/dx in the next section.

3 Bounding the order and degree

We now consider the case E = C(x, t) where ′ = d/dx and

t′/t ∈ C(x), i.e. t = e^{∫ f(x)dx}. Since we do not consider algebraic extensions in this paper, we make the additional assumption that t′/t is not the logarithmic derivative of a C(x)-radical, i.e. t′/t cannot be written as n^{−1}f′/f for f ∈ C(x) and n a nonzero integer. By [4, Theorem 5.1.2], this implies that t is transcendental over C(x) and that the

constant field of E is C. Let u ∈ C(x) and y = e^{∫ udx}. Then, y^{(i)} = P_{i}(u, u′, …, u^{(i−1)})y where the P_{i}’s are polynomials with integer coefficients given by

P_{0} = 1 and P_{i} = P′_{i−1} + uP_{i−1} for i > 0 . (4)

Let L = ∑_{i=0}^{n} a_{i}∂^{i} be a linear ordinary differential operator

with coefficients in C(x). If Le^{∫ udx} = 0, then

a_{n}P_{n}(u, …, u^{(n−1)})+a_{n−1}P_{n−1}(u, …, u^{(n−2)})+…+a_{0} = 0 ,

and the above equation is called the Riccati equation associated with L. The following two Lemmas are special cases of Lemma 2.2 of [9] and their proofs can be found there. For the second one, we use that P_{i}(c, …, c^{(i−1)}) = c^{i} for c ∈ C.

Lemma 3 ([9]) Let p ∈ C[x] be irreducible, and u ∈ C(x) have p-adic expansion of the form u = u_{γ}/p^{γ} + … where γ > 0, u_{γ} ≠ 0 and deg u_{γ} < deg p. Then, for i ≥ 1 the padic expansion of P_{i}(u, …, u^{(i−1)}) has the form v_{iγ}/p^{iγ} +… where

\{ \begin{matrix}v_{iγ} ≡ u_{γ}^{i} & (mod p) & if γ > 1 , \\ v_{iγ} ≡ ∏_{j=0}^{i−1}(u_{1} − jp′) & (mod p) & if γ = 1 . \end{matrix}

Lemma 4 ([9]) Let u ∈ C(x) have a Laurent expansion at x = ∞ of the form u = ∑_{i≥0} u_{i}x^{γ−i} + … where γ ≥ 0 and u_{0} ≠ 0. Then, for i ≥ 1 the Laurent expansion at x = ∞ of P_{i}(u, …, u^{(i−1)}) has the form u_{0}^{i}x^{iγ} + . . ..

During the bounding process, we may have to compute the integer solutions of some algebraic equations. The following Lemma shows that the systems that arise have finitely many such solutions.

Lemma 5 Let K ⊂ F be fields and p ∈ F [Z] be irreducible with deg(p) = m. Let a_{0}, …, a_{n} ∈ F [Z] and q, u ∈ F [Z] be linearly independent over K and such that deg(q) < m and deg(u) < m. If p/| gcd(a_{0}, …, a_{n}), then the set

S = {(γ, α) ∈ K^{2} such that ∑_{i=0}^{n} a_{i}(γu+αq)^{i} ≡ 0 (mod p)}

is finite with at most n elements.

Proof. Let π : F [Z] → F [Z]/(p) be the canonical projection and Q = ∑_{i=0}^{n} π(a_{i})T^{i} ∈ (F [Z]/(p))[T]. Since p/| gcd(a_{0}, …, a_{n}), Q is not identically 0 and has degree at most n, so let β_{1}, …, β_{k} with k ≤ n be its distinct roots in F [Z]/(p), and let p_{1}, …, p_{k} ∈ F [Z] be such that deg(p_{j}) < m and π(p_{j}) = β_{j} for each j. For any (γ, α) ∈ S, we have Q(γπ(u)+απ(q)) = 0, which implies that π(γu+αq) = π(p_{j}) for some j, hence that γu + αq = p_{j} since the degrees of u, q and p_{j} are strictly less than m. But for each j, there is at most one pair (a, b) ∈ K^{2} such that au + bq = p_{j}, since the difference of two such relations with distinct pairs would yield a linear dependence over K for u and q. Therefore, S has at most n elements. □

We now return to bounding the integers h such that L(at^{h}) = 0 for a ∈ C(x)^{∗} and any operator L with coefficients in C(x). For exponential extensions of C(x), we replace the search for exponential solutions by the following adaptation of the algorithm of Proposition 2.3 of [9]. Our modification computes only the potential degrees of “monomial” solutions rather than all the exponential solutions.

Proposition 2 Suppose that E = C(x, t) where ′ = d/dx, t′/t ∈ C(x) is not the logarithmic derivative of a C(x)- radical, that we know a vector space basis for C over ℚ, that we can factor polynomials in C[Z] and compute their integer roots. Then, for any nonzero linear ordinary differential operator L with coefficients in C(x), we can compute a finite set β(L) ⊂ ℤ such that if L(at^{γ}) = 0 for an integer γ and a ∈ C(x)^{∗}, then γ ∈ β(L).

Proof. Write

L = A_{n}∂^{n} + … + A_{0} where A_{i} ∈ C[x], A_{n} ≠ 0 ,

and

\frac{t′}{t} = G + \frac{A}{D}

where G, A, D ∈ C[x], D ≠ 0 and deg(A) < deg(D). Suppose that L(at^{γ}) = 0 for some a ∈ C(x)^{∗} and an integer γ. By adding 0 to β(L), we can assume without loss of generality that γ ≠ 0. We have

∑_{i=0}^{n} A_{i}P_{i}(u, u′, …, u^{(i−1)}) = 0  (5)

where u = a′/a + γt′/t and the P_{i}’s are given by (4). Note that for any irreducible p ∈ C[x], the p-adic expansion of a′/a is a′/a = νp′p^{−1} + … where ν ∈ ℤ is the order of a at p. We distinguish several cases.

1. If G ≠ 0, then the Laurent expansion of t′/t at x = ∞ has the form lc(G)x^{M} + … where M = deg(G), so the expansion of u is γlc(G)x^{M} +…, and by Lemma 4, the expansion of P_{i}(u, u′, …, u^{(i−1)}) is γ^{i}lc(G)^{i}x^{iM} + . . .. Therefore, the expansion of the left-hand side of (5) is

( ∑_{i s.t. deg(A_{i})+iM=μ} lc(A_{i})γ^{i}lc(G)^{i})x^{μ} + …

where μ = max_{i}(deg(A_{i})+iM). It follows that Q(γ) = 0 where

Q = ∑_{i s.t. deg(A_{i})+iM=μ} lc(A_{i})lc(G)^{i}Z^{i} ∈ C[Z]

and β(L) is the set of all the integer roots of Q.

2. If there is an irreducible p ∈ C[x] such that p^{2} | D, then let t′/t = u_{β}/p^{β} + … be the p-adic expansion of t′/t where β ≥ 2 and u_{β} ∈ C[x] is nonzero and has degree strictly less than deg(p). The p-adic expansion of u is then γu_{β}/p^{β} + …, so by Lemma 3, the p-adic expansion of P_{i}(u, u′, …, u^{(i−1)}) is v_{iβ}/p^{iβ} + … where v_{iβ} ≡ γ^{i}u_{β}^{i} (mod p). Let

A_{i} = a_{α_{i}}p^{α_{i}} + …  (6)

with α_{i} ≥ 0 be the p-adic expansion of each nonzero A_{i} and μ = min_{i}(α_{i} − iβ). Then, the p-adic expansion of the left-hand side of (5) is

( ∑_{i s.t. α_{i}−iβ=μ} w_{i})p^{μ} + …

where w_{i} ≡ a_{α_{i}}γ^{i}u_{β}^{i} (mod p). It follows that Q(γ) ≡ 0 (mod p) where

Q = ∑_{i s.t. α_{i}−iβ=μ} a_{α_{i}}u_{β}^{i}Z^{i} ∈ C[x][Z] .

Computing the remainder of Q by p in C[Z][x] (the leading coefficient of p being a unit in C[Z]) and equating each coefficient of x^{j} to 0, we get a system of polynomial equations Q_{0}(γ) = … = Q_{deg(p)−1}(γ) = 0, hence, H(γ) = 0 where H = gcd(Q_{0}, …, Q_{deg(p)−1}) ∈ C[Z]. Note that H is not identically 0, since the leading coefficient of Q is a_{α_{I}} u_{β}^{I} ∈ (C[x]/(p))^{∗} where I = max{i s.t. α_{i} − iβ = μ}. As earlier, β(L) is the set of all the integer roots of Q.

3. If there is an irreducible p ∈ C[x] such that p | D, p^{2} /|D and u_{1} and p′ are linearly independent over ℚ, where u_{1} ∈ C[x] \ {0} is given by deg(u_{1}) < deg(p) and t′/t = u_{1}p^{−1} + … is the p-adic expansion of t′/t. Then, the p-adic expansion of u is (γu_{1} + νp′)p^{−1} + …, where ν ∈ ℤ (note that ν can be 0 if p does not divide the numerator or denominator of a). By Lemma 3, the p-adic expansion of P_{i}(u, u′, …, u^{(i−1)}) is v_{i}/p^{i} + … where v_{i} ≡ ∏_{j=0}^{i−1}(γu_{1} +νp′ −jp′) (mod p). Using (6), the p-adic expansion of the left-hand side of (5) is

( ∑_{i s.t. α_{i}−i=μ} w_{i})p^{μ} + …

where μ = min_{i}(α_{i} − i) and w_{i} ≡ a_{α_{i}} ∏_{j=0}^{i−1}(γu_{1} + νp′ − jp′) (mod p). It follows that Q(γu_{1} + νp′) ≡ 0 (mod p) where

Q = ∑_{i s.t. α_{i}−i=μ} a_{α_{i}} ∏_{j=0}^{i−1}(Z − jp′) ∈ C[x][Z] .

As in the previous case, computing the remainder of Q by p in C[Z][x] and equating it to 0, we get a system of polynomial equations Q_{0}(γ, ν) = … = Q_{deg(p)−1}(γ, ν) = 0 with coefficients in C. If that system has dimension 0, then we can compute (via resultants or Gröbner bases) a nontrivial univariate equation H(γ) = 0 where H ∈ C[Z] is in the ideal generated by Q_{0}, …, Q_{deg(p)−1}. In that case, β(L) is as earlier the set of all the integer roots of H. If the system has a strictly positive dimension^{2}, we then write Q = ∑_{i=0}^{n} Q_{i}Z^{i} where Q_{0}, …, Q_{n} ∈ C[x]. Since the leading coefficient of Q is a_{α_{I}} where I = max{i s.t. α_{i}−i = μ} and deg(a_{α_{I}}) < deg(p), it follows that p/| gcd(Q_{0}, …, Q_{n}), so by Lemma 5, the equation Q(γu_{1} + νp′) ≡ 0 (mod p) has a finite set of rational solutions, which we can compute in the following way: since we can factor elements of C[Z], we can also factor elements of (C[x]/(p))[Z] [7, §8.8], so the roots β_{1}, …, β_{k} of Q in C[x]/(p) can be computed from its linear factors over C[x]/(p). For each β_{j}, equating the powers of x on both sides of the equation γu_{1}+αp′ = β_{j} yields a linear system A_{j}(γ, α)^{T} = v_{j} with coefficients in C for γ and α. Projecting that system on our basis for C over ℚ, we get a matrix B_{j} and a vector w_{j}, both with coefficients in ℚ such that the rational solutions of A_{j}(γ, α)^{T} = v_{j} and B_{j}(γ, α)^{T} = w_{j} are the same. If B_{j} is nonsingular, then that system has a unique solution, which we compute. Otherwise, that system has no solution since its set of rational solutions is finite. Selecting the integer solutions from the rational ones yields the set β(L).

Suppose that none of the above cases is applicable, then t′/t = A/D where deg(A) < deg(D) and every irreducible factor p of D must satisfy the following conditions: p^{2} /|D, and {u_{1}, p′} are linearly dependent over ℚ, where u_{1} is defined as in the third case. This linear dependence implies that for each such p, there is a rational number r_{p} ∈ ℚ such that p does not divide the denominator of t′/t − r_{p}p′/p. Therefore, η = t′/t − ∑_{p|D} r_{p}p′/p ∈ C[x], where the sum is taken over all the irreducible factors of D. But the polynomial part of η is 0, so η = 0. Writing r_{p} = n_{p}/m for a common denominator m, we get t′/t = m^{−1}f′/f where f = ∏_{p|D} p^{n_{p}} ∈ C(x)^{∗}, in contradiction with t′/t not being the logarithmic derivative of a C(x)-radical. □

Theorem 1 Suppose that E = C(x, t) where ′ = d/dx, t′/t ∈ C(x) is not the logarithmic derivative of a C(x)- radical, that we know a vector space basis for C over ℚ, that we can factor polynomials in C[Z] and compute their integer roots. Then, for any nonzero linear ordinary differential operator L with coefficients in E and b_{1}, …, b_{m} ∈ E, we can compute finite sets Δ ⊂ ℤ and Γ ⊂ ℤ such that if there are integers δ ≤ γ, y_{δ}, …, y_{γ} ∈ C(x) and c_{1}, …, c_{m} ∈ C such that y_{δ} ≠ 0 ≠ y_{γ} and

L(y_{δ}t^{δ} + … + y_{γ}t^{γ}) = c_{1}b_{1} + … + c_{m}b_{m} ,

then δ ∈ Δ and γ ∈ Γ.

Proof. This follows immediately from Propositions 1 and 2. □

4 Computing the solutions

We now turn back to the third step, namely computing the solutions once bounds on their degree and order are known. Throughout this section, let K be an arbitrary differential field and t be transcendental and exponential over K. Suppose that Const(K(t)) = Const(K) and that we can effectively solve parameterized linear differential equations over K. Let L be a linear ordinary differential operator with coefficients in K[t], b_{1}, …, b_{m} ∈ K(t) and integers δ ≤ γ be given. We look for y_{δ}, …, y_{γ} ∈ K and c_{1}, …, c_{m} ∈ Const(K) such that L(∑_{i=δ}^{γ} y_{i}t^{i}) = ∑_{j=1}^{m} c_{j}b_{j}. Let z = t^{−δ}y. We have z ∈ K[t], deg(z) ≤ N = γ − δ and, by Lemma 1, Ly = L(t^{δ}z) = t^{δ}L̃z where L̃ is also an operator with coefficients in K[t]. It is thus sufficient to look for solutions p ∈ K[t] of degree at most N ≥ 0 of Lp = ∑_{j=1}^{m} c_{j}b_{j}. Dividing both sides by a power of t if necessary, we can assume that ν(L) = 0. Let d ∈ K[t] be a common denominator for the b_{i}’s and write b_{i} = f_{i} +g_{i}/d where f_{i}, g_{i} ∈ K[t] and deg(g_{i}) < deg(d). If Lp = ∑_{j=1}^{m} c_{j}b_{j} for p ∈ K[t], then

\frac{∑_{j=1}^{m} c_{j}g_{j}}{d} = Lp − ∑_{j=1}^{m} c_{j}f_{j} ∈ K[t]

2According to Lemma 5 this can happen only when u_{1} and p′ are linearly dependent over C.

which implies that Lp = ∑_{j=1}^{m} c_{j}f_{j} and that ∑_{j=1}^{m} c_{j}g_{j} = 0. Equating the coefficients of each t^{i} to 0 in that last condition yields a matrix A with coefficients in K such that A(c_{1}, …, c_{m})^{T} = 0. Using Lemma 7.1.2 of [4] and the associated algorithm, we construct a matrix B with coefficients in Const(K) such that the constant solutions of Ax = 0 and Bx = 0 are the same.

We are reduced to solving Lp = ∑_{j=1}^{m} c_{j}f_{j} where p, f_{1}, …, f_{m} ∈ K[t] and deg(p) ≤ N for a given N.

If N = 0, then p ∈ K, so using the decomposition (1) our equation becomes

∑_{k=0}^{deg(L)} t^{k}L_{k}p = ∑_{j=1}^{m} c_{j}f_{j}

and equating the coefficients of each t^{k} on both sides we get either new linear constraints for the c_{j}’s (if k > deg(L) or L_{k} = 0), or parameterized linear differential equations over K for p (if L_{k} ≠ 0), which we can solve using the induction hypothesis.

If N > 0, write p = p_{0} + tq where p_{0} = p(0) ∈ K and q ∈ K[t] with deg(q) < N. Since t′/t ∈ K, the ideal (t) of K[t] is closed under ′, so the canonical projection π : K[t] → K = K[t]/(t), which sends h ∈ K[t] to h(0) is a differential homomorphism [4, Theorem 4.2.1]. Applying it to our equation we get L_{0}p_{0} = ∑_{j=1}^{m} c_{j}f_{j}(0) where L_{0} = L_{ν(L)} is given by the decomposition (1). It is a nonzero linear ordinary differential operator with coefficients in K, so using the induction hypothesis we compute h_{1}, …, h_{r} ∈ K and a matrix M with coefficients in Const(K) such that for any solution p_{0} ∈ K we must have p_{0} = ∑_{i=1}^{r} e_{i}h_{i} and M(e_{1}, …, e_{r}, c_{1}, …, c_{m})^{T} = 0. Replacing p by ∑_{i=1}^{r} e_{i}h_{i} + tq in the equation, we get

∑_{j=1}^{m} c_{j}f_{j} = Lp = Lp_{0} + L(tq)

= ∑_{j=1}^{m} c_{j}f_{j}(0) + (L − L_{0})(∑_{i=1}^{r} e_{i}h_{i}) + tL̃q

where L̃ is given by Lemma 1. Therefore,

L̃q = ∑_{j=1}^{m} c_{j} \frac{f_{j} − f_{j}(0)}{t} − ∑_{i=1}^{r} e_{i} \frac{1}{t} (L − L_{0})h_{i} .

The above is a parameterized equation with right–hand side in K[t] and we solve it with the bound N −1 for deg(q) using the same method. This process eventually terminates since the degree bound is reduced at each step.

We remark that our algorithm becomes simpler when L = ∑_{k} a_{k}∂^{k} has its coefficients a_{k} ∈ K. In that case any solution y = ∑_{i=δ}^{γ} y_{i}t^{i} of Ly = ∑_{j=1}^{m} c_{j}f_{j} where f_{j} ∈ K[t, t^{−1}] satisfies

L(∑_{i=δ}^{γ} y_{i}t^{i}) = ∑_{i=δ}^{γ} z_{i}t^{i} = ∑_{j=1}^{m} c_{j}f_{j}

where z_{i} = L(y_{i}t^{i})/t^{i} ∈ K by Lemma 2. Writing f_{j} = ∑_{i} f_{ji}t^{i} for f_{ji} ∈ K, we get the equations

(∑_{k} a_{k}(∂ + i \frac{t′}{t} )^{k})y_{i} = ∑_{j=1}^{m} c_{j}f_{ji} for δ ≤ i ≤ γ (7)

as well as the linear constraints

∑_{j=1}^{m} c_{j}f_{ji} = 0 for i < δ or i > γ .

Those constraints can be reduced to linear constraints with constant coefficients as earlier, and the equations (7) can be solved independently of each other. In the homogeneous case and when the coefficients of L are in K, then any solution y ∈ K(t) of Ly = 0 must be in K[t, t^{−1}] [3, Theorem 5.5], and the equations (7) become

(∑_{k} a_{k}(∂ + i \frac{t′}{t} )^{k})y_{i} = 0

so parametric equations are not introduced, unlike when t appears in the coefficients of L. Furthermore, if K = C(x) and ′ = d/dx, then the set β(L) of Proposition 2 contains all the integers for which the above equations can have solutions in K^{∗}.

Together with the algorithms of [3, 9] for the denominator up to a power of t, the bounding methods of the previous sections, and the known algorithms for solving parameterized linear differential equations over C(x) [1, 2, 3, 9], the above algorithm proves the following result:

Theorem 2 Let C be any field, and E = C(x, t) where t is transcendental and exponential over (C(x), d/dx). If we know a vector space basis for C over ℚ and if we can factor polynomials in C[Z] and compute their integer roots, then we can effectively solve parameterized linear differential equations over E.

Note that the requirement of factoring over an algebraic extension of C appears only in one branch of the bounding method, and only when the nonlinear system considered for γ and ν has strictly positive dimension. In particular, this cannot happen when C = ℚ.

5 Examples

Consider first ℚ(x, e^{∫ x+1/x}) = ℚ(x, t) where t′/t = x+1/x, which is not the logarithmic derivative of a ℚ(x)-radical, and look for solutions y ∈ ℚ(x)[t, t^{−1}] of Ly = 0 where

L = ((x^{4} + 2x^{2})t^{2} − x^{3} − x)∂^{2}
+(x^{4} + 3x^{2} − (3x^{5} + 11x^{3} + 6x)t^{2})∂
+(2x^{6} + 9x^{4} + 11x^{2} + 6)t^{2} .

• Bounding the degree and order of y:

– Bounding the degree: deg_{t}(L) = 2 so L_{2}(at^{deg_{t}(y)}) = 0 for some a ∈ ℚ(x) where

L_{2} = (x^{4} + 2x^{2})∂^{2} − (3x^{5} + 11x^{3} + 6x)∂
+2x^{6} + 9x^{4} + 11x^{2} + 6 .

We have t′/t = x + 1/x whose polynomial part has degree 1, so Q(deg_{t}(y)) = 0 where

Q = Z^{2} − 3Z + 2 = (Z − 1)(Z − 2)

so deg_{t}(y) ≤ 2.

– Bounding the order: ν_{t}(L) = 0 so L_{0}(at^{ν_{t}(y)}) = 0 for some a ∈ ℚ(x) where

L_{0} = −(x + x^{3})∂^{2} + (x^{4} + 3x^{2})∂ .  (8)

As before, the polynomial part of t′/t is x, so

Q = −Z^{2} + Z = Z(1 − Z)

so ν_{t}(y) ≥ 0.

• Computing the solution: from the bounding step, we conclude that any solution y ∈ ℚ(x)[t, t^{−1}] of Ly = 0 must be of the form y = y_{0} + y_{1}t + y_{2}t^{2}.

– Reducing modulo t, we have L_{0}y_{0} = 0 where L_{0} is given by (8). Recursively solving that equation over ℚ(x) we find that y_{0} = c_{1} · 1 where c_{1} is an arbitrary constant. We have

\frac{1}{t} (L − L_{0})1 = (2x^{6} + 9x^{4} + 11x^{2} + 6)t ,

so, replacing y by c_{1} + tz where z = y_{1} + y_{2}t we get the new equation

L̃z = −c_{1}(2x^{6} + 9x^{4} + 11x^{2} + 6)t

where

L̃ = ((x^{4} + 2x^{2})t^{2} − x^{3} − x)(∂ + x + \frac{1}{x})^{2}

+(x^{4} + 3x^{2} − (3x^{5} + 11x^{3} + 6x)t^{2})(∂ + x + \frac{1}{x})

+(2x^{6} + 9x^{4} + 11x^{2} + 6)t^{2}

= ((x^{4} + 2x^{2})t^{2} − x^{3} − x)∂^{2}

−((x^{5} + 5x^{3} + 2x)t^{2} + x^{4} + x^{2} + 2)∂

– Reducing modulo t, we have L̃_{0}y_{1} = 0, where

L̃_{0} = −(x^{3} + x)∂^{2} − (x^{4} + x^{2} + 2)∂ .

Recursively solving that equation over ℚ(x) we find that y_{1} = c_{2} · 1 where c_{2} is an arbitrary constant. We have 1/t(L̃ − L̃_{0})1 = 0, so, replacing z by c_{2} + ty_{2} we get the new equation

(t^{2}\overline{L}_{2} + \overline{L}_{0})y_{2} = −c_{1}(2x^{6} + 9x^{4} + 11x^{2} + 6) (9)

where

\overline{L}_{2} = (2x^{2} +x^{4})∂^{2} +(x^{5} +x^{3} +2x)∂ −x^{4} −x^{2} −2 ,

and

\overline{L}_{0} = −(x^{3} + x)∂^{2} − (3x^{4} + 5x^{2} + 4)∂

− \frac{2x^{6} + 6x^{4} + 6x^{2} + 2}{x} .

– Equating the coefficients of t^{2} and t^{0} on both sides of (9) we get the equations \overline{L}_{2}y_{2} = 0 and

\overline{L}_{0}y_{2} = −c_{1}(2x^{6} + 9x^{4} + 11x^{2} + 6) .

Recursively solving those equations over ℚ(x) we find that y_{2} = c_{3}x with the linear constraint c_{3} − c_{1} = 0. Therefore, the general solution of Ly = 0 is

y = c_{1}(1 + xe^{2∫ x+1/x}) + c_{2}e^{∫ x+1/x}

where c_{1} and c_{2} are arbitrary constants.

We illustrate in the second example the simpler version of the algorithm mentioned at the end of section 4: consider the differential equation

\frac{d^{2}y}{dx^{2}} − xy = (4x^{3} − x^{2} + 6x)e^{x^{2}+c}  (10)

where c is an arbitrary parameter. Since the homonogeneous equation has no Liouvillian solutions [8], the inhomogeneous one has Liouvillian solutions if and only it has a solution y ∈ ℚ(c)(x, e^{x^{2}+c}) = C(x, t) where C = ℚ(c) and t′/t = 2x [6]. Since the coefficients of the operator and the denominator of the right–hand side do not involve t, we only need to look for solutions y ∈ C(x)[t, t^{−1}].

• Bounding the degree and order of y: we get Q(γ) = 0 where Q = 4Z^{2}, which implies that β(L) = {0}. Proposition 1 then implies that ν_{t}(y) ≥ 0 and deg_{t}(y) ≤ 1.

• Computing the solution: from the bounding step, we conclude that any solution y ∈ C(x)[t, t^{−1}] must be of the form y = y_{0} + y_{1}t. The equations (7) become d^{2}y_{0}/dx^{2} −xy_{0} = 0, which has no solution in C(x), and

((\frac{d}{dx} + 2x)^{2} − x)y_{1} = 4x^{3} − x^{2} + 6x .

Recursively solving that equation over C(x), we find the unique solution y_{1} = x, so the only Liouvillian particular solution of (10) is y = xe^{x^{2}+c}.

6 Conclusion

We have presented a variant of Singer’s algorithm for the “base case”: simple exponential extensions of (C(x), d/dx). Our method is effective and has been implemented^{3}. Since the steps for finding the denominator and finding the solutions with known bounds work over arbitrary exponential extensions, the only step specific to the base case in our algorithm is the bounding method of Proposition 2. It turns out that if we replace C(x) by a Liouvillian (or monomial) transcendental extension K(θ) of a differential field K (i.e. E = K(θ, t)), then the algorithm of Proposition 2 can still be applied if we replace the hypothesis t′/t − rf′/f ≠ 0 for r ∈ ℚ and f ∈ C(x)^{∗} by the stronger one t′/t−rf′/f ∉ K for r ∈ ℚ and f ∈ K(θ)^{∗}. Work on the details is under way and will be presented in a subsequent paper. We remark that together with the results of [5, 10], our algorithm yields a Liouvillian solver for completely reducible second order linear ordinary differential equations

with coefficients in C(x, e^{∫ f(x)dx}).

3http://www.inria.fr/cafe/Manuel.Bronstein/bernina.html

Ackowledgements

We are thankful to the UMS Medicis for the use of their computing facilities during the programming and testing of this algorithm. In addition, we would like to thank Jacques– Arthur Weil and Jean Moulin–Ollagnier for many useful discussions during the preparation of this article.

References

[1] Abramov, S., Bronstein, M., and Petkovšek, M. On polynomial solutions of linear operator equations. In Proceedings of ISSAC’95 (1995), ACM Press, pp. 290–296.

[2] Abramov, S. A., and Kvashenko, K. Y. Fast algorithms to search for the rational solutions of linear differential equations with polynomial coefficients. In Proceedings of ISSAC’91 (1991), S. Watt, Ed., ACM Press, pp. 267–270.

[3] Bronstein, M. On solutions of linear ordinary differential equations in their coefficient field. Journal of Symbolic Computation 13, 4 (April 1992), 413–440.

[4] Bronstein, M. Symbolic Integration I – Transcendental Functions. Springer, Heidelberg, 1997.

[5] Bronstein, M., Mulders, T., and Weil, J.-A. On symmetric powers of differential operators. In Proceedings of ISSAC’97 (1997), ACM Press, pp. 156–163.

[6] Davenport, J., and Singer, M. Elementary and Liouvillian Solutions of Linear Differential Equations. Journal of Symbolic Computation 2, 3 (September 1986), 237–260.

[7] Geddes, K., Czapor, S., and Labahn, G. Algorithms for Computer Algebra. Kluwer Academic Publishers, Boston, 1992.

[8] Kovacic, J. An Algorithm for Solving Second Order Linear Homogeneous Differential Equations. Journal of Symbolic Computation 2, 1 (March 1986), 3–43.

[9] Singer, M. Liouvillian solutions of linear differential equations with liouvillian coefficients. Journal of Symbolic Computation 11 (1991), 251–273.

[10] Ulmer, F., and Weil, J.-A. Note on Kovacic’s algorithm. Journal of Symbolic Computation 22, 2 (August 1996), 179–200.