Integration of Solutions of Linear Functional Equations

Sergei A. Abramov^{∗}

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

                    abramov@ccas.ru
                   sabramov@cs.msu.su

Mark van Hoeij

Department of mathematics Florida State University, Tallahassee, FL 32306-3027, USA

                                                      hoeij@math.fsu.edu

Abstract We introduce the notion of the adjoint Ore ring and give a definition of adjoint polynomial, operator

and equation. We apply this for integrating solutions of Ore equations.

1 Introduction

The goal of this paper is integration (in the difference case: summation) of solutions of Ore equations. For this purpose we first define an adjoint for an Ore ring, similar to the well-known adjoint for differential operators, and also similar to ideas in [10]. The use of Ore rings allows to handle the case of differential, difference and q-difference equations simultaneously.

An application is integration of special functions, like Bessel functions or hypergeometric functions. If a function satisfies a differential operator L for which the operator L̃ that we will compute has the same order then we have an easy way to integrate. This will be illustrated with an example. The integrals of special functions that we obtain this way are often much less complicated than the integrals given by computer algebra systems.

Another situation where solutions of linear differential equations need to be integrated is the following. For solving linear differential equations one often applies “reduction of order” in case one of the solutions was found. Reduction of order leads to the problem of integrating solutions of a differential equation. In this paper we give a simple and easy to implement method for this problem. Given an operator L, our algorithm computes an operator L̃ of minimal order such that the derivatives of the solutions of L̃ are the solutions of L. In the case that the order of L̃ equals the order of L, this effectively removes, at low computational cost, one integration symbol from the symbolic solutions of the original differential equation. The use of Ore rings makes our algorithm more general, so that it can be applied to the cases of difference and q-difference equations as well.

A preliminary version of this paper appeared as [4].

2 Integrating factors and adjoints

Let k be a field and let K be a ring that contains k. We consider Ore rings k[Δ] and K[Δ] for two different types of Δ:

• Case 1: Δ is a derivation on K, i.e. Δ(ab) = bΔ(a) + aΔ(b) for all a, b ∈ K, and it is also a derivation on k (so Δ(a) ∈ k for a ∈ k).

• Case 2: Δ = σ − 1 where σ is an automorphism of k and of K.

In both cases we will assume that the set of constants Const = {a ∈ K|Δ(a) = 0} is a subfield of k. Let k[Δ] be the ring of all operators ∑_{i=0}^{n} a_{i}Δ^{i}. Similarly define K[Δ]. An element operator L ∈ K[Δ] can be viewed as a Const-linear map from K to K, L(y) = ∑_{i=0}^{n} a_{i}Δ^{i}(y) ∈ K. We will assume that ∑ a_{i}Δ^{i} ∈ K[Δ] acts as the zero map on K if and only if all a_{i} are zero. A common situation for such Ore rings is that one is given an equation L(y) = 0 for some L ∈ k[Δ] and one is interested in finding solutions y of this equation in some (field or ring) extension K of k. By “rational solutions” of L we mean solutions y ∈ k.

∗Supported in part by the RFBR and INTAS under Grant 95-IN-RU-412 and by RFBR under Grant 98-01-00860.

Now k[Δ] and K[Δ] are rings. The multiplication in these rings corresponds to composition of operators. Using the relation

• Case 1: Δ ∘ a = aΔ + Δ(a)

• Case 2: Δ ∘ a = σ(a)Δ + Δ(a)

any product of elements in K[Δ] can be written in a standard form ∑_{i=0}^{n} a_{i}Δ^{i}. We define Δ^{∗} as follows:

• Case 1: Δ^{∗} = −Δ

• Case 2: Δ^{∗} = σ^{−1} − 1

and we can define the adjoint ring of k[Δ] as k[Δ^{∗}]. Note that in case 1 we have k[Δ] = k[Δ^{∗}] but in case 2 k[Δ] needs not be equal to its adjoint ring. Now we can define the adjoint map from

ad : k[Δ] → k[Δ^{∗}]

for an operator L = ∑_{i} a_{i}Δ^{i} as follows:

ad L = ∑_{i}(Δ^{∗})^{i} ∘ a_{i}.

This can be rewritten to the standard form ad L = ∑_{i} b_{i}(Δ^{∗})^{i} for some b_{i} ∈ k. For brevity we will often write L^{∗} instead of ad L.

Now one can verify that the adjoint is a Const-linear bijective map and that

(L ∘ M)^{∗} = M^{∗} ∘ L^{∗}

for all L, M ∈ k[Δ].

Proposition 1

Δ(f) = 0 ⟺ f ∈ Const ⟺ Δ^{∗}(f) = 0.

Additionally for any L ∈ K[Δ]

L(1) = 0 ⟺ ∃_{M} L = M ∘ Δ

and

L^{∗}(1) = 0 ⟺ ∃_{M} L = Δ ∘ M.

Proof: The first statement follows from the definitions. Write L = ∑_{i} a_{i}Δ^{i} for some a_{i} ∈ K. Now L = M ∘Δ for some M if and only if a_{0} = 0. Now the second statement follows because a_{0} = L(1). For the third statement, write L^{∗} ∈ K[Δ^{∗}] (note that in general L^{∗} needs not be an element of K[Δ]) as L^{∗} = ∑_{i} a_{i}(Δ^{∗})^{i} for some a_{i} ∈ K. Now L = ∑_{i} Δ^{i}a_{i} and L = Δ ∘ M for some M iff a_{0} = L^{∗}(1) = 0. □

An element l ∈ K is an integrating factor for L ∈ K[Δ] if lL = Δ∘M, for some M ∈ K[Δ]. The following proposition shows that in the general case the adjoint equation has an important feature which is well-known in the differential case.

Proposition 2 l ∈ K is an integrating factor for L iff L^{∗}(l) = 0.

Proof: By proposition 1 we have lL = Δ∘M for some M iff (lL)^{∗}(1) = 0. Since l ∈ K we have l^{∗}(1) = l(1) = l and so (lL)^{∗}(1) = L^{∗}(l). □

3 Accurate integration

An element g ∈ K is a primitive of f ∈ K if Δ(g) = f. Consider the following problem:

Let f ∈ K and the minimal annihilating operator L ∈ k[Δ] for f be given. So n = ord L is minimal with the property that L ∈ k[Δ] and L(f) = 0. Decide whether there exists a primitive g of f such that the minimal annihilating operator L̃ for g has order n. If so, then construct all such g together with their minimal annihilating operators. We show that this problem (the problem of the accurate integration) can be solved with the help of finding integrating factors.

Let g be any primitive of f and L̃ ∈ k[Δ] be the minimal annihilating operator for g. Now L ∘ Δ(g) = L(f) = 0 hence by the minimality of L̃ (and by the fact that k[Δ] is a Euclidean ring, c.f. [13]) it follows that L̃ is a right-hand factor of L ∘ Δ. Hence

ord L̃ ≤ ord L ∘ Δ = n + 1 and if ord L̃ = n + 1 then L̃ = L ∘ Δ.

Consider the least common left multiple (LCLM) of L̃ and Δ presented in the form

LCLM(L̃, Δ) = L_{1} ∘ Δ,  (1)

L_{1} ∈ k[Δ], ord L_{1} ≤ ord L̃. We have L̃(g) = 0, so L_{1} ∘ Δ(g) = 0, hence L_{1}(f) = 0 and so ord L_{1} ≥ n by the minimality of L. So ord LCLM(L̃, Δ) ≥ n + 1 and hence

ord L̃ ≥ n and if ord L̃ = n then GCRD(L̃, Δ) = 1  (2)

where GCRD stands for greatest common right divisor.

Thus there are two alternatives for ord L̃: n or n + 1. The questions are: when is ord L̃ = n and what is L̃ in this case?

If ord L̃ = n then from equation (2) and the extended Euclidean algorithm it follows that

r ∘ Δ + l̃∘ L̃ = 1  (3)

for some l̃, r ∈ k[Δ] with ord r < ord L̃ = n and ord l̃< ord Δ = 1. Applying equation (3) on g results in

r(f) = g.

Applying Δ on this equation yields Δ ∘ r(f) = f so (1 − Δ ∘ r)(f) = 0. By the minimality of L it follows that 1 − Δ ∘ r = l ∘ L for some operator l; hence

Δ ∘ r + l ∘ L = 1.  (4)

Conversely, if equations (3),(4), ord r < n and ord l̃< 1 hold then one can easily verify that ord L̃ = n, that L̃ is the minimal annihilating operator for r(f) and that Δ(r(f)) = f. Hence equations (3),(4) with the conditions on ord r and ord l̃ are equivalent to the problem of accurate integration.

The inequality ord l̃< 1 implies ord l = ord l̃= 0, i.e. l, l̃∈ k. Both sides of (4) are operators and if we take the adjoints we get

r^{∗} ∘ Δ^{∗} + L^{∗} ∘ l^{∗} = 1.  (5)

Applying the left- and the right-hand sides of (5) to the constant function 1 we obtain

L^{∗}(l) = 1.  (6)

For each solution l ∈ k of (6) we have (1 − lL)^{∗}(1) = 1 − L^{∗}(l) = 0 and so by proposition 1 it follows that equation (4) allows a unique solution r. The minimal annihilating operator of g is defined up to a left-hand factor in k. Therefore we can take l̃= 1 and

L̃ = 1 − r ∘ Δ.  (7)

This operator annihilates one-unique primitive r(f) of f. If operators r_{0} and r_{1} correspond to different solutions l_{0} and l_{1} of (6) then the primitives r_{0}(f) and r_{1}(f) of f are also different (otherwise the operator r_{0}−r_{1} of order < n annihilates f). Since primitives are determined up to constants it follows that (r_{0}−r_{1})(f) must be a constant.

L̃ maps the primitive r(f) of f to 0. Furthermore it maps any constant to itself. Hence it maps any primitive of f to a constant.

The preceding can be formulated as the following

Proposition 3 Let L ∈ k[Δ] be the minimal annihilating operator for f ∈ K and L^{∗}(l) = 1, l ∈ k. Then the equality Δ ∘ r + l ∘ L = 1 uniquely determines r. In turns r lets find the operator L̃ ∈ k[Δ] (up to a factor in k) annihilating the primitive

g = r(f)  (8)

of f. If formula (7) is used to construct L̃ then L̃(g_{1}) ∈ Const for any primitive g_{1} of f. □

If (6) has no solution in k then no primitive of f has a minimal annihilating operator over k of order n. If (6) has a unique solution in k then a primitive and its minimal annihilating operator can be defined uniquely by (8),(7).

Proposition 4 Let M be the set of all solutions of (6) in k. Then M is empty, or M has only one element, or M has the form

M = {l_{0} + Ch|C ∈ Const}  (9)

where l_{0}, h ∈ k, h ≠ 0. In the last case any primitive of f has a minimal annihilating operator of order n.

Proof: Suppose there exists a solution l_{0} of (6). Then the solution space of (6) is of the form l_{0} + V where V is the solution space of L^{∗}(l) = 0. The map l ↦ r(f) (r depends on l by (4)) is an injective (here we use that L is minimal) linear map from l_{0} + V to the set of primitives of f. Since the set of primitives is an affine space of dimension 1, V must have dimension ≤ 1. □

Note that if the solution space of L^{∗}(l) = 0 has dimension > 1 then the map l ↦ r(f) can not be injective because the image of this map has dimension ≤ 1. The fact that the map is not injective means that there exists an r, ord r < ord L, with r(f) = 0 which contradicts our assumption that L is minimal.

Let now M have the form (9). Denote by l_{C} the solution l_{0} + Ch of (6) and by r_{C} and L̃_{C} the operators which are found starting with l_{C}. Since h is an integrating factor for L we have hL = Δ ∘ M, ord M = n − 1. Now from (4) and (7) we obtain

r_{C} = r_{0} − CM  (10)
L̃_{C} = L̃_{0} + CM ∘ Δ  (11)

where r_{0} and L̃_{0} correspond to the solution l_{0} of (6). The operator L̃_{C} is the minimal annihilating operator for the primitive

g_{C} = r_{C}(f)  (12)

of f.

Let g be a primitive of f and C ∈ Const. Then L̃_{C}(g) = L̃_{0}(g) + CM(Δg) = L̃_{0}(g) + CM(f) and M(Δg) ∈ Const because L̃_{C}(g), L̃_{0}(g) ∈ Const. Additionally M(f) ≠ 0 because ord M < ord L. Taking

C = −\frac{L̃_{0}(g)}{M(f)}  (13)

we obtain the value of C such that g = r_{C}(f).

The price which we pay for solving the problem of the accurate integration is finding solutions in k of the equation L^{∗}(y) = 1. If k is the rational function field, then the last problem can be solved effectively in all cases mentioned in the examples below (c.f. [1, 2, 3]).

An implementation called integrate sols is available in the DEtools package in Maple V release 5. As one can see below the algorithm is very short and easy to implement.

Procedure IntegrateSolutions Input: L ∈ k[Δ]

L^{∗} := adjoint(L)

Compute the rational solutions of L^{∗}(y) = 1 if there exists a rational solution then Let l be a rational solution. r := LeftQuotient(1 − lL, Δ) L̃ := 1 − r ∘ Δ

else L̃ := L ∘ Δ The integration operator r does not exist.

end if Output: L̃ and, if it exists, r as well.

4 Examples

Example 1. Let k = C(x), Δ = D = \frac{d}{dx}. Applying the described approach to f = ln x,

L = xD^{2} + D  (14)

gives L^{∗} = xD^{2} + D, and the general rational solution of the equation L^{∗}(y) = 1 is l_{C} = x + C. Therefore l_{0} = x, h = 1. Any primitive of ln x is annihilated by a second order operator. We obtain

L̃_{C} = (x^{2} + Cx)D^{2} − xD + 1, r_{C} = (−x^{2} − Cx)D + x.

It obviously holds for any function f(x) whose minimal annihilating operator has the form (14). For f(x) = ln x we have r_{C}(f(x)) = x ln x − x − C.

Example 2. The algorithm proposed above lets in some cases integrate special functions.

a) The minimal annihilating operator for Bessel function J_{1} is x^{2}D^{2} + xD + (x^{2} − 1). Now L^{∗} = x^{2}D^{2} + 3xD + x^{2}, and L^{∗}(y) = 1 has a unique rational solution \frac{1}{x^{2}} . We obtain

L̃ = D^{2} + \frac{1}{x}D + 1, r = −D − \frac{1}{x}.

Thus we get a primitive of J_{1} in the form

r(J_{1}) = (−D − \frac{1}{x})(J_{1}),

with minimal annihilating operator L̃. Other primitives (r(J_{1}) plus a constant) are annihilated by the operator L ∘ D of order 3.

b) Another example of integration of special functions is the following: L = xD^{2} + (C_{1} + C_{2}x)D + C_{3}. The solutions of this operator L can be expressed in terms of Whittaker functions. Our algorithm produces the operator

r = \frac{x}{C_{2} − C_{3}} D + \frac{C_{1} + C_{2}x − 1}{C_{2} − C_{3}},

so the solutions y of L can be integrated by our method ∫ y dx = r(y).

Example 3. Let k = Q(n), Δ = E −1 where E(n) = n+1, E^{∗} = E^{−1}. Let u_{0}, u_{1}, ... be Fibonacci numbers. Apply the described approach to u_{n}^{2}, L = E^{3} − 2E^{2} − 2E + 1. We obtain L^{∗} = E^{−3} − 2E^{−2} − 2E^{−1} + 1 (note that E = Δ + 1 and therefore E^{∗} = Δ^{∗} + 1 = E^{−1} − 1 + 1 = E^{−1}, it lets one work with linear operators described in terms of E. In general in case 2 one can consider operators from k[σ], setting σ^{∗} = σ^{−1}). The equation L^{∗}(y) = 1 has the unique rational solution −\frac{1}{2}. It shows that one unique primitive of u_{n}^{2} can be annihilated by an operator of order 3: L̃ = −\frac{1}{2}L, while r = \frac{1}{2}(E^{2} − E − 3). This primitive is

\frac{1}{2}(u_{n+2}^{2} − u_{n+1}^{2} − 3u_{n}^{2}).

The inverse of Δ is the summation operator ∑_{i=0}^{n−1}, up to a constant which is \frac{1}{2}(u_{2}^{2} − u_{1}^{2} − 3u_{0}^{2}) = 0 in this example. So

∑_{i=0}^{n−1} u_{i}^{2} = \frac{1}{2}(u_{n+2}^{2} − u_{n+1}^{2} − 3u_{n}^{2}).

There are several methods for proving such formulas. Our algorithm does more in this example, it also finds this formula. Certainly, we need the minimal annihilator for u_{n}^{2}. Since u_{n}^{2} is a d’Alembertian sequence this annihilator can easily be constructed, for example, by algorithm [6]. Remark that algorithm [6] itself uses the accurate integrating algorithm.

Example 4. This is an example of integration of algebraic functions. If the minimal polynomial f(x, y) ∈ C(x)[y] for an algebraic function α(x) is given then the minimal annihilating differential operator L ∈ C(x)[D] for α(x) can be constructed as follows. The algebraic function α(x) and its derivatives are elements of C(x, α), which is a C(x)-vector space of dimension deg_{y} f(x, y). One can compute α(x), α′(x), α′′(x), … in this vector space. Take ρ the minimal integer such that α(x), α′(x), α′′(x), …, α^{(ρ)}(x) are C(x)-linearly dependent. This linear relation gives L.

Now again k = C(x), Δ = D = \frac{d}{dx}. Let f(x, y) = y^{3} + xy + x^{2} and let α ∈ \overline{k} be a root of f(x, y) as a polynomial in y. One obtains the following annihilating operator for α

L = D^{2} − 2\frac{1}{x(4 + 27x)}D + 2\frac{3x + 1}{x^{2}(4 + 27x)} ∈ k[D]

The equation L^{∗}(l) = 1 has a unique rational solution −1/45 − x/12 + (9x^{2})/20. This yields

r = (\frac{1}{45} + \frac{x}{12} − \frac{9x^{2}}{20} )D + \frac{162x^{2} − 9x − 2}{180x}

and so

∫ α dx = r(α) = (\frac{1}{45} + \frac{x}{12} − \frac{9x^{2}}{20} )α′ + \frac{162x^{2} − 9x − 2}{180x}α.

We obtain a unique primitive, despite the fact that primitives are not unique, but only unique up to a constant. Note that our method produces a primitive of α if and only if there is a primitive which is a k-linear combination of α and its derivatives. So in the hardest case in elementary integration (the case when logarithmic extensions are needed) our algorithm will not produce a primitive (so then L̃ must be L ∘ D).

Example 5. Let k = C(x) and v(x) = 1/x. The minimal annihilating differential operator L over k for v(x) is L = xD + 1. The equation L^{∗}(l) = 1 has no solutions in k. So every primitive of 1/x is only annihilated by operators of order ≥ 2. The primitives can not be obtained by applying linear differential operators over k to v(x).

Example 6. Given a first order (i.e. hypergeometric) sequence over C(n), the well-known Gosper’s algorithm ([11]) decides whether there exists another sequence of such a kind that is a primitive for the given sequence. The algorithm in this paper generalizes Gosper’s algorithm in two ways: it solves the analogous problem for a wider class of equations, and for any order n instead of only n = 1. Using (8) we can express the mentioned primitive explicitly in terms of the given sequence (function). We will illustrate this in the difference and q-difference cases.

a) Hypergeometric case. Let k = Q(n), Δ = E − 1. Let s_{n} = \binom{2n}{n}/4^{n}, n = 0, 1, … This sequence is hypergeometric; s_{n+1}/s_{n} = (2n + 1)/(2n + 2). We have L = 2(n + 1)E − (2n + 1). We obtain L^{∗} = 2nE^{−1} − (2n + 1), and the equation L^{∗}(y) = 1 has the unique rational solution −1. It shows that one-unique primitive of s_{n} can be annihilated by an operator of order 1, and r = 2n. This primitive is

2ns_{n} = \frac{2n\binom{2n}{n}}{4^{n}}

and in particular

∑_{n=0}^{N−1} \frac{\binom{2n}{n}}{4^{n}} = \frac{2N\binom{2N}{N}}{4^{N}}.

b) q-Hypergeometric case. Let q be a new variable. A sequence {h_{n}} is q-hypergeometric if h_{n+1}/h_{n} is a rational function of q, q^{n}. The sequence

(a; q)_{n} = \{ \begin{matrix}1, & if n = 0, \\ (1 − a)(1 − aq) ⋯ (1 − aq^{n−1}), & if n > 0, \end{matrix}

where a is a parameter is a q-hypergeometric for any fixed value of a due to (a; q)_{n+1}/(a; q)_{n} = 1 − aq^{n−1} = 1 − a\frac{q^{n}}{q}. We will consider the sequence h_{n} = q^{n}(q; q)_{n}. Here h_{n+1}/h_{n} = q(1 − q · q^{n}). Write x for q^{n}. Then h_{n+1}/h_{n} = q(1 − qx).

Let k = Q(q)(x), Δ = Q − 1 where Q(x) = qx, Q^{∗} = Q^{−1}. We have L = Q − q(1 − qx). We obtain L^{∗} = Q^{−1} − q(1 − qx), and the equation L^{∗}(y) = 1 has the unique rational solution 1/(q^{2}x). It shows that one-unique primitive of h_{n} can be annihilated by an operator of order 1, and r = −1/(qx). This primitive is

−\frac{1}{q · q^{n}} h_{n} = −\frac{(q; q)_{n}}{q}

and in particular

∑_{n=0}^{N−1} q^{n}(q; q)_{n} = \frac{1 − (q; q)_{N}}{q} .

Another way to obtain the last formula was demonstrated in [5].

5 Conclusion

It follows from the results in this paper that the following 3 problems are equivalent.

• Find the solutions l ∈ k of L^{∗}(l) = 1.

• Let f ∈ K. Let the minimal annihilating operator L ∈ k[Δ] for f be given. Decide whether there exists r ∈ k[Δ] such that r(f) is a primitive of f. If so, then construct such r.

• The problem of accurate integration.

• Computing solutions (r, l) of equation (4). Note that according to section 3.1 in [12] this is equivalent to computing a complement of Const in the solution space of L ∘ Δ.

Similar problems, but in a more general situation, are studied in [8, 9]. Our approach is less general but it has the advantage of simplicity, it only uses an adjoint and rational solutions, which are quite efficient.

Acknowledgement

We would like to thank Marko Petkovšek who provided us with useful comments on an earlier draft.

References

[1] S. A. Abramov (1989): Rational solutions of linear difference and differential equations with polynomial coefficients, USSR Comput. Maths. Math. Phys. 29, 7 – 12. Transl. from Zl. Vychislit. matem. mat. fiz. 29, 1611 – 1620.

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

[3] S. A. Abramov (1995): Rational solutions of linear difference and q-difference equations with polynomial coefficients, Programming and Comput. Software 21, No 6, 273 – 278. Transl. from Programmirovanie, No 6, 3 – 11.

[4] S. A. Abramov, M. van Hoeij (1997): A method for the Integration of Solutions of Ore Equations, Proc. ISSAC’97, 172 – 175.

[5] S. A. Abramov, P. Paule, M. Petkovšek (1998): q-Hypergeometric solutions of q-difference equations, Discrete Math. 180, 3–22.

[6] S. A. Abramov, E. V. Zima (1997): Minimal completely factorable annihilators, Proc. ISSAC’97, 290 – 297.

[7] M. Bronstein, M. Petkovšek (1995): An introduction to pseudo-linear algebra, Theoretical Computer Science 157, 3 – 33.

[8] F. Chyzak (1997): An extension of Zeilberger’s fast algorithm to general holonomic functions, Proc. FPSAC’97 , 172–183.

[9] F. Chyzak, B. Salvy (1996): Non-commutative elimination in Ore algebra proves multivariate holonomic identities, INRIA Research Report, No 2799.

[10] P. M. Cohn (1995): Skew Fields, Theory of General Division Rings, Encyclopedia of Mathematics and its Applications, Cambridge University Press, 57.

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

[12] M. van Hoeij (1996): Rational Solutions of the Mixed Differential Equation and its Application to Factorization of Differential Operators, Proc. ISSAC’96, 219–225.

[13] O. Ore (1933): The theory of non-commutative polynomials, Ann. Maths. 34, 480 – 508.