On Polynomial Solutions of Linear Operator Equations

Sergei A. Abramov^{∗}

Computer Center of the Russian Academy of Science, Vavilova 40, Moscow 117967, Russia.

abramov@sms.ccas.msk.su

Manuel Bronstein^{†}

Institute for Scientific Computation, ETH Zentrum IFW, CH-8092 Zürich, Switzerland.

            bronstein@inf.ethz.ch

Marko Petkovšek^{‡}

Department of Mathematics, University of Ljubljana, Jadranska 19, 61111 Ljubljana, Slovenia.

marko.petkovsek@uni-lj.si

1 Introduction

Let K be a field of characteristic 0 and L : K[x] → K[x] an endomorphism of the K-linear space of univariate polynomials over K. We consider the following computational tasks concerning L:

T1. Homogeneous equation Ly = 0: Compute a basis of Ker L in K[x].

T2. Inhomogeneous equation Ly = f: Given f ∈ K[x], compute a basis of the affine space L^{−1}(f) in K[x].

T3. Parametric inhomogeneous equation Ly = ∑_{i=1}^{m} λ_{i}f_{i}: Given f_{1}, f_{2}, …, f_{m} ∈ K[x], compute a basis of Ker L′ where L′ : (K[x] ⊕ K^{m}) → K[x] and L′ : (y, λ) ↦ Ly − ∑_{i=1}^{m} λ_{i}f_{i}, for y ∈ K[x], λ ∈ K^{m}.

Many problems and algorithms in differential and difference algebra contain these tasks as subproblems which, however conceptually simple, often account for a fair share of the overall computing time. For instance, the algorithms for

• finding all the rational solutions of differential and (q-)difference equations [Abr89b, Abr95],

• finding Liouvillian solutions of differential equations [Sin91],

• finding (q-)hypergeometric solutions of (q-)difference equations [Pet92, Abr&Pet95],

• factoring linear differential and difference operators with rational coefficients [Bro&Pet94],

• indefinite hypergeometric summation (Gosper’s algorithm) [Gos78],

• definite hypergeometric summation (Zeilberger’s “creative telescoping” algorithm) [Zei91], all need to find polynomial solutions of operator equations at some stage. For instance, task T3 is crucial for the algorithms presented in [Sin91] and [Zei91]. Clearly, it is important to use as efficient algorithms as possible for such time-critical subproblems.

We consider three types of linear operators: differential, difference, and q-difference operators, all with polynomial coefficients. Throughout the paper, L is the operator, r its order and d the maximum degree of its coefficients. To find polynomial solutions of an equation of the form Ly = f, it is natural to proceed in two steps:

1. Compute a degree bound N for polynomial solutions of the equation.

2. Given N, compute polynomial solution(s) of degree(s) ≤ N.

Degree bounds required in step 1 are given in [Abr89b] for the differential and difference cases. Here we improve slightly upon these bounds and also give them for the qdifference case. Our main contribution, however, is to step 2 of the above procedure. Usually, it is carried out by means of undetermined coefficients: expand y(x) with respect to the usual power basis (x^{n})_{n=0}^{∞} and substitute this expression into the equation, obtaining a triangular system of linear algebraic equations for c_{i}’s. Then there is a one-to-one correspondence between solutions of this system and those of the original equation.

The problem with the method of undetermined coefficients is that the number of unknown coefficients, N + 1 (respectively N + m + 1 in the parametric case), is often much larger than the order of the equation, r. We show that by using appropriate polynomial bases, the resulting linear system has a band-diagonal matrix which represents a recurrence for the unknown coefficients. After taking full advantage of this recurrence, the remaining linear system that needs to be solved is of size corresponding to the bandwidth of that matrix.

In particular, for differential, difference and q-difference equations, we reduce step 2 of the above procedure to a system of r+d linear algebraic equations with only r unknowns (respectively r + m in the parametric case). Thus, in order to minimize the number of unknowns in the final linear system, one should use the method proposed here when N ≥ r, and the method of undetermined coefficients when N < r.

∗Supported in part by grant 95-01-01138a of the Russian Fund for Fundamental Research.

†Research supported in part by CATHODE (ESPRIT WG 7213).

‡Supported in part by grant J2-6193-0101-94 of the Ministry of Science and Technology of Slovenia.

In Section 2 we define formal series expansions with respect to distinct-degree bases of K[x], and use them to find polynomial solutions of operator equations. Section 3 gives algorithms to solve problems T1, T2 and T3. Given a recurrence we show how to find all sequences satisfying it on a specified interval, on which its leading coefficient may vanish at some points. Section 4 applies these results to the cases when L is a linear differential, difference, or q-difference operator. In Section 5 we present some empirical evidence suggesting that our method can be much faster than the commonly used method of undetermined coefficients. Finally we show in Section 6 how our method can also be used to compute formal series solutions of the equations we are considering.

Throughout the paper, the set of nonnegative integers is denoted by IN. We take the degree of the zero polynomial to be −∞.

2 Formal series expansions in K[x]

Let (P_{n}(x))_{n=0}^{∞} be a sequence of polynomials from K[x] such that

P1. deg P_{n} = n for n ≥ 0,

P2. P_{n} | P_{m} for 0 ≤ n < m.

Then {P_{0}, P_{1}, …} is a basis of the K-linear space K[x] and every polynomial from K[x] has a unique expansion in terms of polynomials P_{n}(x). Let l_{n} : K[x] → K be linear functionals such that l_{n}(P_{m}) = δ_{mn}. Then

p(x) = ∑_{n=0}^{deg p} l_{n}(p) P_{n}(x)  (1)

for all p ∈ K[x].

Example 1 Take P_{n}(x) = (x − a)^{n}/n! where a ∈ K is arbitrary. Then l_{n}(p) = p^{(n)}(a), and (1) is Taylor’s formula for p(x).

Example 2 Take P_{n}(x) = \binom{x−a}{n} where a ∈ K is arbitrary. Then l_{n}(p) = Δ^{n}p(a), and (1) is Newton’s interpolation formula for p(x).

Example 3 Take P_{n}(x) = x^{n}/(n_{q})! where (n_{q})! = 1_{q}2_{q} … n_{q} and m_{q} = 1+q+…+q^{m−1}. Then l_{n}(p) = D_{q}^{n}p(0), where D_{q}p(x) = (p(qx) − p(x))/(qx − x), and (1) is a qanalogue of Taylor’s formula for p(x).

Let K[[(P_{n})_{n=0}^{∞}]] denote the algebra of formal series of the form

s(x) = ∑_{n=0}^{∞} c_{n}P_{n}(x)

where c_{n} ∈ K, and let l̃_{n} : K[[(P_{n})_{n=0}^{∞}]] → K be linear functionals such that l̃_{n}(s(x)) = c_{n}. Operations in K[[(P_{n})_{n=0}^{∞}]] are defined by setting

l̃_{n}(s + t) = l̃_{n}(s) + l̃_{n}(t),  (2)
l̃_{n}(s t) = ∑_{k,m: k+m≥n} l̃_{k}(s)l̃_{m}(t) l_{n}(P_{k}P_{m}), (3)

for all s, t ∈ K[[(P_{n})_{n=0}^{∞}]]. From property P2 it follows that l_{n}(P_{k}P_{m}) = 0 when n < max{k, m}, hence the sum in (3) is finite. By virtue of (1), (2) and (3), K[x] is a subalgebra of K[[(P_{n})_{n=0}^{∞}]], and l̃_{n} agrees with l_{n} on K[x].

Example 4 If P_{n}(x) = x^{n} then K[[(P_{n})_{n=0}^{∞}]] = K[[x]], the usual algebra of formal power series over K.

In order to extend L to an endomorphism L̃ of the K-linear space K[[(P_{n})_{n=0}^{∞}]], we assume additionally that

P3. there are A, B ∈ IN, and elements α_{i}(n) ∈ K for n = 0, 1, 2, … and i = −A, 1 − A, …, B, such that

LP_{n} = ∑_{i=−A}^{B} α_{i}(n) P_{n+i}  (4)

where α_{−A}(n) is not identically zero, and P_{k} is taken to be 0 when k < 0.

Define L̃ on s ∈ K[[(P_{n})_{n=0}^{∞}]] by

L̃s = ∑_{n=0}^{∞} ( ∑_{i=−A}^{B} α_{i}(n − i) l̃_{n−i}(s)) P_{n}  (5)

where l̃_{k} is taken to be 0 when k < 0.

The operator L̃ has been defined in (5) so as to agree with L on K[x]. This enables us to replace the equation Ly = f in K[x] with L̃y = f in K[[(P_{n})_{n=0}^{∞}]]. By (5), the latter equation is equivalent to

∑_{i=−B}^{A} α_{−i}(n + i) l̃_{n+i}(y) = l_{n}(f),  (6)

for all n ≥ 0. This is a recurrence for the unknown sequence ⟨l̃_{n}(y)⟩_{n=0}^{∞} ∈ K, of effective order at most A + B (in fact A + b where b is given below in (8)), with leading coefficient α_{−A}(n + A), and homogeneous for n > deg f.

Theorem 1 Let Ly = f where y is a polynomial. Then deg y ≤ N where

N = max{−b − 1, deg f − b} ∪ {n ∈ IN; α_{b}(n) = 0)}, (7)
b = max{k ∈ ZZ; α_{k}(n) ≢ 0}  (8)

and α_{k} is as in (4).

Proof: Let y be a polynomial of degree N satisfying Ly = f. Clearly b ≤ B; it can even happen that b < 0. We distinguish two cases: either N + b < 0 and hence N < −b, or N + b ≥ 0. In the latter case, equation (6) holds for n ≥ N + b. When n > N + b, the left-hand side of (6) is zero. For n = N + b it equals α_{b}(N)l_{N}(y), as l_{N+b−i}(y) = 0 when i < b. Now either α_{b}(N) ≠ 0 and then N + b = deg f, or else N is a root of α_{b}(x). In summary: either N < −b, or N = deg f − b, or α_{b}(N) = 0. □ Our algorithm is based on the following simple observa-

tion.

Theorem 2 Let N and b be as in (7) and (8), respectively, with N ≥ 0. Then for every y ∈ K[[(P_{n})_{n=0}^{∞}]] the following are equivalent:

(i) y is a polynomial which satisfies Ly = f,

(ii) the sequence ⟨l̃_{n}(y)⟩_{n=0}^{∞} satisfies (6) for n ≤ N +b, and l̃_{n}(y) = 0 for n > N.

Proof: By (5) and since N is an upper bound for degree of polynomial solutions of Ly = f, it is clear that (i) implies (ii).

Conversely, to see that (ii) implies (i), let n > N + b. Then n − b > N, hence all terms appearing on the left of (6) are zero. Also, n > N + b ≥ deg f − b + b = deg f, thus the right side of (6) is zero as well. It follows that (6) is satisfied for all n ≥ 0. Then by (5), y is a polynomial which satisfies Ly = f. □

Thus to find all polynomial solutions y ∈ K[[(P_{n})_{n=0}^{∞}]] of Ly = f, it suffices to compute all vectors (l̃_{−B}(y), l̃_{1−B}(y), …, l̃_{N+A+B}(y)) which satisfy (6) for 0 ≤ n ≤ N +b, then select those with l̃_{n}(y) = 0 for n < 0 and for n > N and set y = ∑_{n=0}^{N} l̃_{n}(y)P_{n}(x). This can be done by using recurrence (6) in either the forward or backward direction, taking conditions at one end as initial conditions and those at the other end as constraints on the free parameters. The backward direction approach in the case of differential and difference equations and the power basis P_{n}(x) = x^{n} has been essentially used already in [Abr89a, Abr&Kva91]. Here we describe the forward direction approach, which can be implemented more efficiently, specially in the nonsingular case where computing with indeterminates can be completely avoided. In the differential and difference case, it is always possible to choose a basis for which the forward direction will be nonsingular (see section 4).

3 The algorithm

3.1 Forward solutions of recurrences

Here we show how to find a generating set for the affine space of vectors v = (v_{0}, v_{1}, …, v_{N+A+b}) ∈ K^{N+A+b+1} which satisfy

∑_{i=−b}^{A} α_{−i}(n + i) v_{n+i} = l_{n}(f),  (9)

for 0 ≤ n ≤ N + b (take v_{n} = 0 for n < 0). Let

S = {n ∈ IN; α_{−A}(n) = 0}

be the set of singularities of (9), and N = {0, 1, …, A−1}∪ S. Denote by σ the cardinality of the set {n ∈ S; A ≤ n ≤ N + A + b}.

As we compute from (9), we maintain a list of vectors V, a list of indeterminates I, a list of equations E, and an additional vector g. Throughout the course of the algorithm, V and I will have equal length whose current value will be denoted by t. The only operations performed on these lists are appending a new element, and extending all vectors in V by one component. We denote the current elements of V and I by V = (v^{(1)}, v^{(2)}, …, v^{(t)}) and I = (c_{i_{1}}, c_{i_{2}}, …, c_{i_{t}}), respectively. Note that we start indexing elements of lists with 1 and components of vectors with 0.

Initially, the lists V, I, E are empty, and g is the empty vector. Next, for n = 0, 1, …, N+A+b, repeat the following: If n ∉ N then set

v_{n} = − (∑_{i=1}^{A+b} α_{i−A}(n − i) v_{n−i}) / α_{−A}(n)  (10)

for all v ∈ V, and

g_{n} = (l_{n−A}(f) − ∑_{i=1}^{A+b} α_{i−A}(n − i) g_{n−i}) / α_{−A}(n) (11)

while leaving E and I unchanged.

Otherwise (n ∈ N) set v_{n} = 0 for all v ∈ V, g_{n} = 0, and append a new element to each of the lists V, I, E:

V := V + (e^{(n+1)}),
I := I + (c_{n}),
E := E + (E_{n−A}), provided that n ≥ A.

Here e^{(k)} denotes the vector of length k with a 1 in the last position and zeros everywhere else, + denotes concatenation of lists, and E_{n} is Equation (9). Before appending E_{n−A} to E, replace each indeterminate c_{k} with ∑_{j=1}^{t−1} c_{i_{j}} v_{k}^{(j)}. This ensures that only indeterminates from I appear in E. Should the equation be 0 = 0 after this replacement, we do not need to add it to E.

After termination of this extension loop, t = A + σ, E contains σ equations, and the vectors v^{(j)} and g have length N +A+b+1. Denote by s_{j}(L̃) and h(L̃, f) the corresponding polynomials

s_{j}(L̃) = ∑_{n=0}^{N+A+b} v_{n}^{(j)}P_{n} for j = 1, 2, …, t ,

and

h(L̃, f) = ∑_{n=0}^{N+A+b} g_{n}P_{n} .

Denote by E(L̃, f) the final list of equations. Our algorithm proceeds from here on in slightly different ways for tasks T1, T2, and T3. Those are covered in the next 3 subsections. In the nonsingular case (σ = 0), this algorithm reduces

to starting with

V := (e^{(1)}, e^{(2)}, …, e^{(A)}),
I := (c_{0}, c_{1}, …, c_{A−1}),
E := (),
g := (0, 0, …, 0) ∈ K^{A}

and simply using (10) and (11) to extend g and the vectors in V for n = A, A + 1, …, N + A + b. This makes I and E unnecessary in that case.

Note that in the course of repeatedly applying formulas (10) and (11), we have to evaluate the polynomials α_{i} at successive integer values of n. Through either the method of finite differences or chains of recurrences [Zima84, Bach&al94], this can be done efficiently, using only deg(α_{i}) additions in K for each evaluation.

A detailed example of the singular case of this algorithm is given in Section 4.3.

3.2 Homogeneous equation

Compute s_{j} = s_{j}(L̃) for j = 1, 2, …, t and find the general solution (c_{i_{1}}, c_{i_{2}}, …, c_{i_{t}}) of the homogeneous linear system composed of E(L̃, 0) and of

∑_{j=1}^{t} c_{i_{j}} l_{n}(s_{j}) = 0, for N < n ≤ N + A + b.

Then y = ∑_{j=1}^{t} c_{i_{j}} s_{j} is the general polynomial solution of Ly = 0 over K.

3.3 Inhomogeneous equation

Compute s_{j} = s_{j}(L̃) for j = 1, 2, …, t and h = h(L̃, f). Then find the general solution (c_{i_{1}}, c_{i_{2}}, …, c_{i_{t}}) of the inhomogeneous linear system composed of E(L̃, f) and of

∑_{j=1}^{t} c_{i_{j}} l_{n}(s_{j}) = −l_{n}(h), for N < n ≤ N + A + b.

Then y = ∑_{j=1}^{t} c_{i_{j}} s_{j} + h is the general polynomial solution of Ly = f over K. Clearly, any particular solution of the above linear system yields a particular polynomial solution of Ly = f.

3.4 Parametric inhomogeneous equation

Let N_{i} denote the degree bound corresponding to f_{i}, and N = max_{1≤i≤m} N_{i}. Compute s_{j} = s_{j}(L̃) for j = 1, 2, …, t and h_{i} = h(L̃, f_{i}), for i = 1, 2, …, m. Let E(L̃, f_{i}) = (u^{(1)} · c = r_{i1}, u^{(2)} · c = r_{i2}, …, u^{(σ)} · c = r_{iσ}), where c = (c_{i_{1}}, c_{i_{2}}, …, c_{i_{t}}) is the vector of indeterminates, u^{(k)} are constant vectors of length t, and · denotes inner product. Now find the general solution (c_{i_{1}}, c_{i_{2}}, …, c_{i_{t}}, λ_{1}, λ_{2}, …, λ_{m}) of the homogeneous linear system composed of

u^{(k)} · c = ∑_{i=1}^{m} λ_{i}r_{ik}, for 1 ≤ k ≤ σ,

and

∑_{j=1}^{t} c_{i_{j}} l_{n}(s_{j}) + ∑_{i=1}^{m} λ_{i}l_{n}(h_{i}) = 0, for N < n ≤ N + A + b.

Then y = ∑_{j=1}^{t} c_{i_{j}} s_{j} +∑_{i=1}^{m} λ_{i}h_{i} is the general polynomial solution of Ly = ∑_{i=1}^{m} λ_{i}f_{i} over K.

3.5 The size of the linear system

The final linear system has at most σ + A + b equations and σ + A unknowns (resp. σ + A + m unknowns, in the parametric case).

4 Choosing the polynomial basis

4.1 Differential equations

Let L = ∑_{k=0}^{r} p_{k}(x)D^{k} where p_{k}(x) ∈ K[x], p_{r} ≠ 0, and Dp(x) = dp(x)/dx. Choose a ∈ K such that p_{r}(a) ≠ 0. The polynomials P_{n}(x) = (x − a)^{n}/n! clearly satisfy properties P1 and P2 from Section 2. Let d = max ∪_{k=0}^{r}{j; l_{j}(p_{k}) ≠ 0}. Since P_{m}P_{n} = \binom{m+n}{m}P_{m+n} and D^{k}P_{n} = P_{n−k}, using (1) on p_{k} we obtain

LP_{n} = ∑_{k=0}^{r} p_{k}P_{n−k} = ∑_{i=−r}^{d} P_{n+i} ∑_{j=0}^{d} \binom{n + i}{j}l_{j}(p_{j−i})

where p_{k} = 0 when k < 0 or k > r. Comparing this with (4) we read off

A = r, B = d,

α_{i}(n) = ∑_{j=0}^{d} \binom{n + i}{j}l_{j}(p_{j−i}),

thus P3 is satisfied. Obviously, α_{i}(n) is a polynomial in n. Since α_{−A}(n + A) = ∑_{j=0}^{d} \binom{n}{j}l_{j}(p_{j+r}) = l_{0}(p_{r}) = p_{r}(a) ≠ 0 by the choice of a, recurrence (9) has no singularities and σ = 0. Hence the final linear system has at most r + d equations and r unknowns.

The value of b required for the degree bound is the maximum i such that l_{j}(p_{j−i}) ≠ 0, for some j. This is the same as the maximum i such that deg p_{j−i} = j. Writing k = j −i, we obtain b = max_{0≤k≤r}(deg p_{k} − k).

Alternatively, we can work with the power basis P_{n}(x) = x^{n} which also satisfies P1, P2 and P3, but here recurrence (9) can have singularities.

4.2 Difference equations

Let L = ∑_{k=0}^{r} p_{k}(x)Δ^{k} where p_{k}(x) ∈ K[x], p_{r} ≠ 0, and Δp(x) = p(x + 1) − p(x). Choose a ∈ K such that p_{r}(n + a) ≠ 0 for all nonnegative integer n. For example, take a = 0 if p_{r}(x) has no nonnegative integer zero, else take a = max{x ∈ IN; p_{r}(x) = 0} + 1. The polynomials P_{n}(x) = \binom{x−a}{n} clearly satisfy properties P1 and P2 from Section 2. Let d = max ∪_{k=0}^{r}{i; l_{i}(p_{k}) ≠ 0} as before. Multiplying a variant of the well-known Chu-Vandermonde identity

\binom{x}{n} = ∑_{i} \binom{x − m}{i − m}\binom{m}{m + n − i}

by \binom{x}{m}, revising binomials and replacing x by x − a, we get

P_{m}P_{n} = ∑_{i} \binom{i}{m}\binom{m}{i − n} P_{i}.  (12)

Since Δ^{k}P_{n} = P_{n−k}, using (1) on p_{k} and then (12) we find

LP_{n} = ∑_{k=0}^{r} p_{k}P_{n−k}

= ∑_{i=−r}^{d} P_{n+i} ∑_{k=0}^{r} ∑_{j=0}^{d} \binom{n + i}{j}\binom{j}{i + k}l_{j}(p_{k})

where p_{k} = 0 when k < 0 or k > r. Comparing this with (4) we read off

A = r, B = d,

α_{i}(n) = ∑_{k=0}^{r} ∑_{j=0}^{d} \binom{n + i}{j}\binom{j}{i + k}l_{j}(p_{k}),

thus P3 is satisfied. Obviously, α_{i}(n) is a polynomial in n. Since α_{−A}(n + A) = ∑_{k=0}^{r} ∑_{j=0}^{d} \binom{n}{j}\binom{j}{k−r}l_{j}(p_{k}) =

∑_{j=0}^{d} \binom{n}{j}l_{j}(p_{r}) = p_{r}(n + a) ≠ 0 for all nonnegative integer n, by the choice of a, recurrence (9) has no singularities and σ = 0. Hence the final linear system has at most r + d equations and r unknowns.

The value of b required for the degree bound can be shown to be b = max_{0≤k≤r}(deg p_{k} − k) as in the differential case.

If we work with the basis P_{n}(x) = \binom{x}{n} conditions P1, P2, P3 are satisfied but recurrence (9) may have singularities. The power basis P_{n}(x) = x^{n} does not satisfy P3.

4.3 q-Difference equations

Here we assume that q ∈ K is not zero and not a root of unity. Let L = ∑_{k=0}^{r} p_{k}(x)Q^{k} where Qp(x) = p(qx), p_{k}(x) ∈ K[x], p_{r} ≠ 0, and not all l_{0}(p_{k}) are zero. The polynomials P_{n}(x) = x^{n} clearly satisfy properties P1 and P2 from Section 2. Let d = max ∪_{k=0}^{r}{j; l_{j}(p_{k}) ≠ 0}. Since P_{m}P_{n} = P_{m+n} and Q^{k}P_{n} = q^{nk}P_{n}, using (1) on p_{k} we obtain

LP_{n} = ∑_{k=0}^{r} p_{k}q^{nk}P_{n} = ∑_{i=0}^{d} P_{n+i} ∑_{k=0}^{r} q^{nk}l_{i}(p_{k})

where p_{k} = 0 when k < 0 or k > r. Comparing this with (4) we read off

A = 0, B = d,

α_{i}(n) = ∑_{k=0}^{r} q^{nk}l_{i}(p_{k}),  (13)

thus P3 is satisfied. Obviously, α_{i}(n) is a polynomial in q^{n}. The value of b required for the degree bound equals d. However, as α_{−A}(n + A) = α_{0}(n) = ∑_{k=0}^{r} q^{nk}l_{0}(p_{k}) can

have nonnegative integer roots, recurrence (9) may be singular with σ ≤ r. Hence the final linear system has at most r + d equations and r unknowns.

Unlike the cases of difference and differential equations, there seems to be no obvious choice of a polynomial basis P_{n} which would guarantee nonsingularity of (9) in the qdifference case.

Example 5 Consider the equation

(1 − q^{10} − (q − q^{10}) x) y(q^{2}x) −

(1 − q^{20} − (q^{2} − q^{20}) x) y(qx) +

q^{10}(1 − q^{10} − (q^{2} − q^{11}) x) y(x) =

(q^{21} − q^{20} − q^{12} + q^{10} + q^{2} − q) x (14)

Here r = 2, d = 1 and deg f = 1. From (13) we find

α_{0}(n) = (1 − q^{10})(q^{n} − 1)(q^{n} − q^{10}),  (15)
α_{1}(n) = (q^{10} − q)(q^{n} − q)(q^{n} − q^{10}).

Since b = d = 1, degrees of polynomial solutions either equal deg f − b = 0 or are roots of α_{d}(n) = α_{1}(n), therefore they belong to the set {0, 1, 10} and we take N = 10. From (15) we find N = S = {n ∈ IN; α_{0}(n) = 0} = {0, 10}. Recurrence (9) in this case is

α_{0}(n)v_{n} + α_{1}(n − 1)v_{n−1} = l_{n}(f), for 0 ≤ n ≤ 11.

We start with V = I = E = () and g = (). At n = 0 we set g = (0) and put the vector (1) into V, c_{0} into I and equation α_{0}(0)c_{0} = l_{0}(f) into E. This equation turns out to be the identity 0 = 0, so E remains empty. We perform next the extension step for n = 1, 2, …, 11. Except at n = 10 we use the formulas v_{n} = −α_{1}(n − 1)v_{n−1}/α_{0}(n) and g_{n} = (l_{n}(f)−α_{1}(n−1)g_{n−1})/α_{0}(n), obtained from (10) and (11). At n = 9 we have V = ((1, −1, 0, 0, 0, 0, 0, 0, 0, 0)) and g = (0, 1, 0, 0, 0, 0, 0, 0, 0, 0). At n = 10 we extend both vectors with 0 and append (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1) to V, c_{10} to I and equation α_{1}(9)c_{9} = 0 to E. However, since c_{9} is not in I we replace it by v_{9}^{(1)}c_{0}. But v_{9}^{(1)} = 0, hence the new equation is 0 = 0 again, and E remains empty. After one more step (which extends all three vectors by 0 thanks to α_{1}(10) = 0) the extension loop terminates with V = ((1, −1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0)), g = (0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), I = (c_{0}, c_{10}) and E = (). The corresponding polynomials are s_{1} = 1 − x, s_{2} = x^{10} and h = x. The final linear system consists of a single equation c_{0}l_{11}(s_{1}) + c_{10}l_{11}(s_{2}) = −l_{11}(h) which holds identically. Hence the general polynomial solution of (14) is y(x) = x + C_{1}(1 − x) + C_{2}x^{10}. As it happens, there exist solutions for all three possible degrees (0, 1, and 10).

5 Some time comparisons

Table 1 shows the times that were needed to find the rational solutions of the Legendre equation

L_{m}(y) = (1 − x^{2}) y′′ − 2x y′ + m(m + 1) y = 0

where n is a fixed positive integer. This equation has a onedimensional space of rational solutions generated by the m^{th} Legendre polynomial which has degree m. For that equation and the basis P_{n}(x) = x^{n}/n! we get r = 2, d = 2, A = 2, B = 2, p_{0} = m(m + 1), p_{1} = −2x, p_{2} = 1 − x^{2}, α_{−2} = 1, α_{−1} = α_{1} = α_{2} = 0, α_{0}(n) = m(m + 1) − n(n + 1), b = 0 and N = max{0, m} = m, so recurrence (9) is:

v_{n+2} + (m(m + 1) − n(n + 1)) v_{n} = 0

for n ≥ 0.

m  Ratlode Series Orthopoly
20  0.25  0.18  0.06
30  0.35  0.20  0.13
40  0.60  0.23  0.23
50  0.83  0.28  0.43
100  1.98  0.35  3.62
200  6.35  0.67  56.17
500  40.05  1.52  2736.72
1000 161.67  5.70  •
10000  •  398.92  •
15000  •  935.70  •
19000  •  1612.50  •

Table 1: Maple CPU seconds to solve L_{m}(y) = 0 (SPARC 10/41).

Ratlode is the ratlode package of the Maple share library, which uses a straightforward implementation of the method of undetermined coefficients; Series is a straightforward implementation of the method presented here, and Orthopoly is the built-in orthopoly[P] function of Maple which computes Legendre polynomials. A dot indicates that the computation did not terminate after 24 hours.

Table 2 shows the times that were needed to find polynomial solutions of the equation

B̃_{m}(y) = y(x + 1) − y(x) = mx^{m−1}

where m is a fixed positive integer. This equation has general solution of the form y(x) = B_{m}(x) + C where B_{m}(x) is the Bernoulli polynomial of degree m, and C is an arbitrary constant. For that equation and the basis P_{n}(x) = \binom{x}{n} we get r = 1, d = 0, A = 1, B = 0, p_{0} = 0, p_{1} = 1, α_{−1} = 1, α_{0} = 0, b = −1 and N = max{0, m} = m, so recurrence (9) is simply the recurrence of order 0

v_{n+1} = Δ^{n}f(0)

for n ≥ 0 where f(x) = mx^{m−1}.

n  Poly  Series BernoulliB
20  2.48  0.43  0.25
30  7.90  0.72  1.82
40  20.40  1.12  7.30
50  45.12  1.48  20.90
100 635.83  6.07  554.43
200 9471.98 36.70  14768.7
300 49663.60 114.57  •
400  •  282.83  •
700  •  2084.08  •

Table 2: Mathematica CPU seconds to solve B̃_{n}(y) = nx^{n−1} (SPARC 10/41).

Poly is a straightforward implementation of the method of undetermined coefficients; Series is a straightforward implementation of the method presented here, and BernoulliB is the built-in BernoulliB[n, x] function of Mathematica which computes Bernoulli polynomials. As the first two implementations return solutions satisfying y(0) = 0, the nth Bernoulli number (BernoulliB[n] of Mathematica) has been added to their respective outputs in order to obtain identical results in all three cases. A dot indicates that the computation did not terminate after 24 hours.

6 Formal series solutions

In conclusion, we remark that our algorithm can also be used to compute formal series solutions in K[[(P_{n})_{n=0}^{∞}]] of equations of the form Ly = f: let M = max(S). Using our extension loop for n = 0, 1, …, M, we can construct the set P of all the polynomials p ∈ K[x] satisfying deg(p) ≤ M and

Lp(x) ≡ f (mod P_{M+1}(x)) .  (16)

The set P together with recurrence (9) describes all the solutions in K[[(P_{n})_{n=0}^{∞}]] of Ly = f, in the sense that any such solution must be an element of P prolongated by (9).

Example 6 Consider the equation

(q^{6}x^{3} + 1) y(q^{2}x) − q^{14}y(x) = 0

with the usual power basis P_{n}(x) = x^{n}. The recurrence (9) in this case is

((q^{n})^{2} − q^{14}) v_{n} + (q^{n})^{2}v_{n−3} = 0 for n ≥ 0  (17)

so N = S = {7} and M = 7. After performing the extension step for n = 0, 1, …, 7, solving (16) yields P = {Cx^{7}} where C is an arbitrary constant. P together with (17) describes all the power series solutions of the original equation, namely

y(x) = C (x^{7} + \frac{q^{6}}{1 − q^{6}} x^{10} + \frac{q^{6}}{1 − q^{6}} \frac{q^{12}}{1 − q^{12}} x^{13} + ⋯)

= C ∑_{k=0}^{∞} \frac{q^{3k(k+1)}}{(q^{6}; q^{6})_{k}} x^{3k+7}

where (a; q)_{k} = (1 − a)(1 − aq) ⋯ (1 − aq^{k−1}).

In the differential and (q-)difference cases, we can compute [Abr89b, Abr95] a universal denominator d(x) ∈ K[x] such that any solution of Ly = f of the form

\frac{z(x)}{p(x)}  (18)

where z(x) ∈ K[[(P_{n})_{n=0}^{∞}]] and p ∈ K[x] can be written as z̃(x)/d(x) where z̃(x) ∈ K[[(P_{n})_{n=0}^{∞}]]. So we can apply our algorithm after doing the change of variable z(x) = d(x)y(x), obtaining all solutions of the form (18). This leads to formal Laurent series in the differential case.

Acknowledgements

The third author wishes to express his thanks to Peter Paule who pointed out that (12) is just the Chu-Vandermonde identity in disguise, and also provided its q-analogue.

References

[Abr89a] S. A. Abramov (1989): Problems of computer algebra involved in the search for polynomial solutions of linear differential and difference equations, Moscow Univ. Comput. Math. and Cybernet., no. 3, 63–68.

[Abr89b] S. A. Abramov (1989): Rational solutions of linear difference and differential equations with polynomial coefficients, USSR Comput. Maths. Math. Phys., vol. 29, 7–12.

[Abr95] S.A. Abramov (1995): Rational solutions of linear difference and q-difference equations with polynomial coefficients, these proceedings.

[Abr&Kva91] S. A. Abramov, K. Yu. Kvashenko (1991): Fast algorithm to search for the rational solutions of linear differential equations, Proc. ISSAC ’91, 267–270.

[Abr&Pet95] S. A. Abramov, M. Petkovšek (1995): Finding all q-hypergeometric solutions of q-difference equations, Proc. FPSAC ’95, in press.

[Bach&al94] O. Bachmann, P. Wang, E. V. Zima (1994): Chains of Recurrences – a method to expedite the evaluation of closed-form functions, Proc. ISSAC ’94, 242– 249.

[Bro&Pet94] M. Bronstein, M. Petkovšek (1994): On Ore rings, linear operators and factorisation, Programmirovanie 20, N.1, 27–45 (also Research Report 200, Informatik, ETH Zürich).

[Gos78] R. W. Gosper, Jr. (1978): Decision procedure for indefinite hypergeometric summation, Proc. Natl. Acad. Sci. USA 75, 40–42.

[Pet92] M. Petkovšek (1992): Hypergeometric solutions of linear recurrences with polynomial coefficients, J. Symb. Comp. 14, 243–264.

[Sin91] M. F. Singer (1991): Liouvillian solutions of linear differential equations with Liouvillian coefficients, J. Symb. Comp., vol. 11, 251–273.

[Zei91] D. Zeilberger (1991): The method of creative telescoping, J. Symb. Comput. 11, 195–204.

[Zima84] E. V. Zima (1984): Automatic Construction of Systems of Recurrence Relations, USSR Comput. Maths. Math. Phys., vol. 24, 193–197.