THE METHOD OF DIFFERENTIATING UNDER THE INTEGRAL SIGN

Gert Almkvist * and Doron Zeilberger**

Appeared in J. Symbolic Computation 10 (1990), 571-591

”I could never resist a definite integral” (G. H. Hardy [Co])

”One thing I never did learn was contour integration. I had learned to do integrals by various methods shown in a book that my high-school physics teacher Mr. Bader had given me. ...

The book also showed how to differentiate parameters under the integral sign-It’s a certain operation. It turns out that’s not taught very much in the universities; they don’t emphasize it. But I caught on how to use that method, and I used that one damn tool again and again. So because I was self-taught using that book, I had peculiar methods of doing integrals.

The result was that, when guys at MIT or Princeton had trouble doing a certain integral, it was because they couldn’t do it with the standard methods they had learned in school. If it was contour integration, they would have found it; if it was a simple series expansion, they would have found it. Then I come along and try differentiating under the integral sign, and often it worked. So I got a great reputation for doing integrals, only because my box of tools was different from everybody else’s, and they had tried all their tools on it before giving the problem to me.” (Richard P. Feynman[F])

Introduction

The method of differentiating under the integral sign can be described as follows. Given a function F (x, y) of x and y, one is interested in evaluating

R(x) := ∫_{−∞}^{∞} F (x, y)dy .  (1)

By differentiating R(x) once or more with respect to x, and using integration by parts or change of variable, one gets a differential equation for R(x). Then one evaluates R(x) and its derivatives in an ”easy” point, and then solves the differential equation.

Although termed a ”method”, differentiating under the integral sign could hardly have been considered more than a trick, and the examples given in the few textbooks that treat it ([Ed]) are rather simple. One notable exception is Melzak’s[Me] delightful book, that devotes a whole section to this method.

In this paper we are going to present a unified theory of differentiation under the integral sign for the important and wide class of so-called holonomic functions. We will recall the ”slow” elimination algorithm of [Z1], and present it in a form that will hopefully lead to a fast version using Buchberger’s method of Grobner bases adapted to the Weyl algebra. For the subclass of ”hyperexponential” functions we will give a fast algorithm, that given F (x, y), finds a linear differential equation satisfied by R(x). This ”fast” algorithm is a continuous analog of Zeilberger’s [Z2] [Z3]

* Department of Mathematics, University of Lund, Box 118, 22100 Lund, Sweden. ** Department of Mathematics and Computer Science, Drexel University, Philadelphia, PA 19104.

Supported in part by NSF grant DMS8800663. Address in 1999: Dept. of Mathematics, Temple Univ., Philadelphia, PA 19122, zeilberg@math.temple.edu (see also [W-Z]) algorithm for hypergeometric summation, that uses Gosper’s[Go] algorithm as a subroutine. In Appendix 1 (which is published in [A-Z]) we give a listing of a MAPLE program that implements this algorithm, while Appendix 2 contains sample inputs and outputs.

While the problem of indefinite integration in closed form is completely solved and implemented ([Ri], [No], [Tr], [DST] ), there is currently no well-developed theory of definite integration ( see [No], 3.2). The present paper can be viewed as a first step toward developing a theory of definite integration.

The Theoretical Foundation

In [Z1] Bernstein’s class of holonomic functions is considered, and it is shown how the methodology of holonomic functions and systems forms a natural framework for special function identities. Roughly a function F (x_{1}, ..., x_{n}) is holonomic if it ”satisfies as many linear differential equations with polynomial coefficients as possible” (see appendix 4, footnote 1 for the precise definition).

The algebra of linear differential operators with polynomial coefficients in (x_{1}, ..., x_{n}) is called the Weyl algebra. It is generated by the ”indeterminates” x_{1}, ..., x_{n}, D_{1}, ..., D_{n}, that obey the commutation relations

x_{i}x_{j} = x_{j}x_{i ,}D_{i}D_{j} = D_{j}D_{i ,}x_{i}D_{j} = D_{j}x_{i} (i ≠ j)

and the ”uncertainty principle”

D_{i}x_{i} − x_{i}D_{i} = 1 ,

where 1 is the identity operator.

Excellent expositions on the Weyl algebra and Bernstein’s theory can be found in [Bj] and [Eh]. It was proved in [Z1](section 6.4), using ideas of Bernstein[Be], that if F (x, y) is holonomic and P (x, y, D_{x}, D_{y}) and Q(x, y, D_{x}, D_{y}) are operators that generate a ”holonomic ideal” then one can ”eliminate” y. In other words, it is possible to find operators A(x, y, D_{x}, D_{y}) and B(x, y, D_{x}, D_{y}) such that

S(x, D_{x}, D_{y}) = A(x, y, D_{x}, D_{y})P (x, y, D_{x}, D_{y}) + B(x, y, D_{x}, D_{y})Q(x, y, D_{x}, D_{y})

is independent of y. It follows that F (x, y) is also annihilated by S: S(x, D_{x}, D_{y})F (x, y) ≡ 0. One can write

S(x, D_{x}, D_{y}) = S̄(x, D_{x}) + D_{y}S′(x, D_{x}, D_{y}) .  (2)

It is possible to show that it is always possible to have S̄ ≠ 0. It follows that

S̄(x, D_{x})F (x, y) = −D_{y}S′(x, D_{x}, D_{y})F (x, y) .  (3)

In the applications we have in mind, we will have that lim_{|y|→∞}x^{n}y^{m}D_{x}^{s}D_{y}^{t}F (x, y) = 0 for all integers n, m, s and t. Therefore integrating with respect to y from −∞ to ∞, we get

S̄(x, D_{x})R(x) = 0 .

It was observed in [Z1] that the slow elimination algorithm for the discrete case, that was an adaptation of Sylvester’s dialytic elimination method from commutative algebra, goes verbatim to the continuous case. The hope was raised that Buchberger’s fruitful concept of Grobner bases could be used to get a fast elimination algorithm. We now reiterate this hope and make one more observation that should make the elimination even faster.

Observation: In order for things to work we don’t need that S′ in (2) be independent of y! Our previous insistence probably made S̄ of higher order than necessary.

This leads to the following

Research Problem: Develop a fast algorithm, using Grobner bases in the Weyl algebra ([Ga]), that inputs two operators P (x, y, D_{x}, D_{y}) and Q(x, y, D_{x}, D_{y}) and outputs operators A, B, C, S̄ such that

S̄(x, D_{x}) := AP + BQ + D_{y}C  (4)

only depends on x and D_{x}.

Basically we have to find an element S̄(x, D_{x}) that does not involve y and D_{y}, that belongs to the ”ideal” generated by P,Q and D_{y}. But note that this ideal is neither ”left” nor ”right”, and therefore not an ideal at all, in the usual sense of the word. It is ”left” with respect to P and Q and ”right” with respect to D_{y}, and so one has to develop an ambidextrous Grobner bases algorithm in the Weyl algebra.

A major step towards this goal was recently achieved by Takayama[Ta].

A related notion to holonomicity is that of D-finiteness, which was considered by Lipshitz[Li]. A formal power series f(x_{1}, ..., x_{n}) is D-finite if the vector space generated by all its partial derivatives is finite dimensional over the field of rational functions in (x_{1}, ..., x_{n}). Every holonomic function is ipso facto D-finite, but the converse is not known (see [Ta]). For our present purposes, the elegant elementary approach of Lipshitz[Li] suffices and there is no need to invoke the theory of holonomic functions.

A Detailed Example

In order to make the ideas clear, we will use ad hoc elimination to treat the following integral.

R(x) := ∫_{−∞}^{∞} exp(−\frac{x^{2}}{y^{2}} − y^{2})dy .

With F (x, y) = exp(−x^{2}/y^{2} − y^{2}) we have

\{{D_{x}F = (−2x/y^{2})F \atop D_{y}F = (2x^{2}/y^{3} − 2y)F}

or

\{{(y^{2}D_{x} + 2x)F = 0 \atop (y^{3}D_{y} + 2y^{4} − 2x^{2})F = 0} .

We want to apply ∫_{−∞}^{∞} dy on these equations. Therefore we want to get rid of all terms containing y, except those of the type D_{y}G. This depends on the fact that

∫_{−∞}^{∞} D_{y}Gdy = [G]_{−∞}^{∞} = 0 ,

since we assume that everything in sight tends to zero as y → ±∞. This means that after integration we can ignore all terms starting with D_{y} ( i.e. put D_{y} = 0). Now

y^{3}D_{y} = D_{y}y^{3} − 3y^{2 ,}

hence the second equation becomes

QF = (D_{y}y^{3} − 3y^{2} + 2y^{4} − 2x^{2})F = 0 .

Apply D_{x} on the left and use the first equation

y^{2}D_{x}F = −2xF

several times. We obtain

D_{x}QF = D_{y}D_{x}y^{3}F − 3D_{x}y^{2}F + 2y^{2}D_{x}y^{2}F − 2D_{x}x^{2}F

= D_{y}D_{x}y^{3}F − 3(−2x)F + 2y^{2}(−2x)F + −2D_{x}x^{2}F .

D_{x}^{2}QF = D_{y}D_{x}^{2}y^{3}F + 6D_{x}xF − 4y^{2}(D_{x}x)F − 2D_{x}^{2}x^{2F}

= D_{y}D_{x}^{2}y^{3}F + 6D_{x}xF − 4y^{2}(xD_{x} + 1)F − 2D_{x}^{2}x^{2}F

= D_{y}D_{x}^{2}y^{3}F + 6D_{x}xF + 8x^{2}F − 4y^{2}F − 2D_{x}^{2}x^{2}F .

D_{x}^{3}QF = D_{y}D_{x}^{3}y^{3}F + 6D_{x}^{2}xF + 8D_{x}x^{2}F − 4D_{x}y^{2}F − 2D_{x}^{3}x^{2}F =

= D_{y}D_{x}^{3}y^{3}F + 6D_{x}^{2}xF + 8D_{x}x^{2}F + 8xF − 2D_{x}^{3}x^{2}F .

Differentiating

S(x) := ∫_{−∞}^{∞} F (x, y)dy

with respect to x, we obtain

D_{x}^{r}S(x) = ∫_{−∞}^{∞} D_{x}^{r}F (x, y)dy .

Integrating D_{x}^{3}QF, the first term vanishes, and we get

(6D_{x}^{2}x + 8D_{x}x^{2} + 8x − 2D_{x}^{3}x^{2})S(x) ≡ 0 .

Using the rules of ”multiplication” in the Weyl algebra, this factors to:

(xD_{x} + 3)(D_{x} − 2)(D_{x} + 2)S(x) = 0 .

This differential equation has the general solution

S(x) = C_{1}e^{2x} + C_{2}e^{−2x} + C_{3}y_{3 ,}

where y_{3} satisfies

y′′_{3} − 4y_{3} = A/x^{3} for A ≠ 0 .

One can show that y_{3} cannot be bounded at 0. But

S(0) = ∫_{−∞}^{∞} e^{−y^{2}}dy = \sqrt{π} ,

hence C_{3} = 0. Assume that x > 0. Then S(x) → 0 as x → ∞ implies that also C_{1} = 0. It follows that

S(x) = \sqrt{π}e^{−2x} .

The above example has only pedagogical value, since it is done much easier by performing the substitution t = y − x/y on the ”obvious” integral ∫_{−∞}^{∞} exp(−t^{2}) = π (see Appendix 4, footnote 2) or by an argument that combines differentiation under the integral sign and substitution,that is given in p. 220 of Edwards book [Ed] (reproduced in Appendix 4, footnote 3) . The advantage of the present method is that it is purely mechanical and has been programmed. We really hope that the experts on Grobner bases will develop a fast elimination algorithm.

The Class of Hyperexponential Functions

The discrete analog of the exponential function e^{cx} is the geometric sequence a_{n} = c^{n}, that has the property that a_{n+1}/a_{n} = is constant. An important generalization of the class of geometric sequences is the class of hypergeometric sequences, that consist of sequences a_{n} for which a_{n+1}/a_{n} is a rational function of n. Hypergeometric sequences are important because they are the terms of hypergeometric series, by which most of the classical special functions can be expressed.

Since the exponential function can be defined as a function f(x) for which f′(x)/f(x) is a constant, it is natural to define, as the continuous analog of hypergeometric sequences, an hyperexponential function to be a function f(x) for which the logarithmic derivative f′(x)/f(x) is a rational function of x.

A function of several variables F (x_{1}, ..., x_{n}) is hyperexponential if the logarithmic derivatives D_{x_{i}}F/F, with respect to every variable x_{i}, are rational functions of (x_{1}, ..., x_{n}).

Consider equation (3) where F (x, y) is hyperexponential. Since a hyperexponential function is obviously D-finite (see [Li]), it follows from a result of Lipshitz ([Li], Lemma 3) that there exist operators S̄ and S′ such that (3) holds. Now set

G(x, y) := −S′(x, D_{x}, D_{y})F (x, y) .  (5)

It is easily proved by induction that for hyperexponential F, D_{x}^{i}D_{j}^{j}F/F is a rational function for every nonnegative i and j. It follows that G(x, y) is a multiple of F (x, y) by a certain rational function, and thus hyperexponential on its own right. Now equation (3) can be rewritten as

S̄(x, D_{x})F (x, y) = D_{y}G(x, y) .  (6)

Thus we can forget about S′ and concentrate our attention on finding the winning pair S̄ and G, now that we know that they exist. Suppose you already knew, or guessed, S̄. Then (6) asserts that S̄(x, D_{x})F (x, y) can be integrated in closed form with respect to y, and Risch’s algorithm [Ri], that is implemented on many symbolic packages, can be used to find G(x, y). However, since S̄ is not given, but has to be found, we chose to develop our own algorithm for indefinite integration, that is, in some sense, a special case of Risch’s algorithm, but that we will be able to extend to find both S̄ and G of (6). This algorithm, that is a continuous analog of Gosper’s ([Go],see also [La], [GKP]) decision procedure for indefinite hypergeometric summation, will now be presented.

or: Gosper’s Algorithm Translated To The Continuous”

This section follows closely Gosper’s[Go] classical paper, and the numbering of the equations in this section correspond to that of Gosper’s paper. To distinguish them from the other equations in this paper, they are prefixed with a ”G”.

Input: A hyperexponential function a(y), i.e. a rational function P (y)/Q(y) such that a′(y)/a(y) = P (y)/Q(y). (Throughout this section a′(y) means D_{y}a(y) = \frac{d}{dy} a(y).)

Output:: A hyperexponential function S(y) such that S′(y) = a(y) if such a hyperexponential function S(y) exists, or the statement ”a(y) does not have a hyperexponential indefinite integral”.

We first find polynomials p(y), q(y), r(y) such that

\frac{(a(y)/p(y))′}{(a(y)/p(y))} = \frac{q(y)}{r(y)} ,  (G5)

where q(y) and r(y) satisfy the condition

gcd(r(y), q(y) − jr′(y)) = 1 for every non − negative integer j.  (G6)

It is always possible to do this for if gcd(r(y), q(y) − jr′(y)) = g(y), then this common factor can be eliminated by

p(y) ← p(y)g(y)^{j ,} r(y)← r(y)/g(y) , q(y)← j(r(y)/g(y))′ + (q(y) − jr′(y))/g(y) .  (G6‘)

The values of j for which such g(y)′s exist can be readily detected as the non-negative integer roots of the resultant of r(y) and q(y) − jr′(y) with respect to y.

We now write

S(y) = r(y)f(y)a(y)/p(y) ,  (G7)

where f(y) is to be determined. By using S′(y) = a(y),

f(y) = (S(y)/S′(y))p(y)/r(y) ,

so f(y) is a rational function whenever S(y)/S′(y) is, (i.e. whenever S(y) is hyperexponential.) By substituting Eq. (G7) into S′(y) = a(y), we get

a(y) = r(y)f(y)(a(y)/p(y))′ + r(y)f′(y)(a(y)/p(y)) + r′(y)f(y)(a(y)/p(y))

Multiplying this through by p(y)/a(y), and using Eq. (G5), we have

p(y) = (q(y) + r′(y))f(y) + r(y)f′(y) ,  (G8)

the differential equation for f(y).

Theorem: If S′(y)/S(y) is a rational function of y, then f(y) is a polynomial.

Proof: We already know that f(y) is a rational function when S′(y)/S(y) is, so suppose

f(y) = \frac{c(y)}{d(y)} ,  (G9)

where d(y) is a polynomial of positive degree, and

gcd(c(y), d(y)) = 1 .  (G9a)

Substituting (G9) in (G8) we get

qcd + r′cd + r(c′d − cd′) = d^{2}p .  (G10)

It is always possible to write

d(y) = A(y)^{k}d̄(y) , where gcd(A(y), d̄(y)) = 1, gcd(A(y), A′(y)) = 1 .  (G11)

Substituting (G11) in (G10) it is readily seen that A(y) divides r(y), i.e.. for some polynomial r̄(y)

r(y) = A(y)r̄(y) .  (G12)

Substituting (G12) in (G10) easily yields that A(y) is also a factor of q(y) − (k − 1)r′(y), but this contradicts the assumption (G6), and thus d(y) must be constant, QED.

All that remains is to look for a polynomial f(y) satisfying equation (G8), given p(y), q(y), and r(y). To do this, we must find an upper bound k, for the degree of the polynomial f(y).

Case 1:

deg(r(y)) ≤ deg(q(y) + r′(y)) = l .

Setting f(y) = c_{k}y^{k} + O(y^{k−1}), we get p(y) = Lc_{k}y^{k+l} + O(y^{k+l−1}), where L is a nonzero constant. Since the two sides of the equation must be of equal degree,

k = deg(f(y)) = deg(p(y)) − l .

Case 2:

deg(q(y) + r′(y)) < deg(r(y)) = l .

Estimating

f(y) = c_{k}y^{k} + c_{k−1}y^{k−1} + O(y^{k−2})

and

f′(y) = kc_{k}y^{k−1} + (k − 1)c_{k−1}y^{k−2} + O(y^{k−3)}

in Eq. (G8) leads to

p(y) = L(k)c_{k}y^{k+l−1} + O(y^{k+l−2}) ,

where L(k) is linear in k and free of c_{k}. Let k_{0} be the root of L(k). Then, the largest k for which c_{k} can be nonzero is

k = max{ k_{0} , deg(p(y)) − l + 1 }, if k_{0} is an integer , or

k = deg(p(y)) − l + 1 , if not .

In either case, if k < 0, the indefinite integral S(y) is not hyperexponential; otherwise we can substitute a kth degree polynomial, with k + 1 undetermined coefficients, for f(y) in Eq. (G8). Equating like powers of y, we get a system of linear equations whose consistency is equivalent to the existence of a hyperexponential S(y) such that S′(y) = a(y) and whose solution provides f(y), whereupon the indefinite integral, S(y), is given by Eq. (G7).

A Fast Algorithm For Definite Hyperexponential Integration

Consider (1) where F (x, y) is hyperexponential. By the theory of holonomic systems, we know for sure that there exist S̄(x, D_{x}) and a hyperexponential G(x, y) such that

S̄(x, D_{x})F (x, y) = D_{y}G(x, y) .  (6)

We now consider x as an auxiliary parameter and set

S̄(x, D_{x}) = ∑_{0}^{m} a_{i}(x)D_{x}^{i}

where m is the guessed order of S̄. Letting the a_{i} be unknown, we evaluate the left of (6), with the undetermined a_{i}(x), and apply the algorithm of section 5 to find the G(x, y). In addition to the unknown coefficients of f(y) (that now depend on x) in this algorithm, we also have the unknowns a_{i}(x). Since they appear linearly, we can modify the above algorithm and require that the resulting system of linear equations be consistent, and when it is, we get as the solution, not only the information about G(x, y), but also the successful S̄(x, D_{x}). We first try m = 0 and increase m until we get a solution. Note that there is no ”trial and error” involved here, and we are guaranteed to eventually get a successful m. Furthermore, it is not hard to derive an a priori upper bound for m.

When we use the subroutine of section 5, the left side of (6) already has the form of a polynomial times a hyperexponential function, so the starting value of p(y) is not 1, and we already have a head start.

A listing of a MAPLE program that implements both the algorithms of section 6 and 7 is given in appendix 1, and appendix 2 contains examples of input and output.

The method of differencing under the integral sign The same idea can be used to find a linear recurrence relation satisfied by

a_{n} := ∫_{−∞}^{∞} F_{n}(y)dy ,

where both F′_{n}(y)/F_{n}(y) and F_{n+1}(y)/F_{n}(y) are rational functions of n and y, that possibly depend on other variables and parameters. We try a generic recurrence operator applied to F_{n}(y), and enquire, by the above algorithm, when it has a hyperexponential primitive function with respect to y.

Simple example (Courtesy of L. Euler): If F_{n}(y) = y^{n}e^{−y}, we have

(n + 1 − N)F_{n}(y) = D_{y}(yF_{n}(y)) ,

where Na_{n} := a_{n+1}. Integrating w.r.t to y from −∞ to ∞ yields (n + 1 − N)a_{n} = 0, i.e. a_{n+1} = (n + 1)a_{n}, and so a_{n} = n!.

Appendix 3 gives a computer-generated proof of the recurrence that leads to the evaluation of Euler’s beta integral as well as four other, less trivial, recurrences.

Applications To The Theory Of Special Functions

It is perhaps surprising that while many deep hypergeometric sum identities are known([Ba]), all of which can be handled by the algorithm of [Z2] and [Z3], relatively few hyperexponential integral identities of the form (1) are known, where both F (x, y) and R(x) are hyperexponential. In fact, the only non-trivial cases known to us are those that are obtained by applying the change of variable

t = x − \frac{a}{x}

to the integral ∫_{−∞}^{∞} exp(−t^{2})dt = \sqrt{π}, and iterating (see appendix 4, footnote 2). Furthermore, the differential equation that the algorithm outputs is never, in our experience, first order, so one has to resort to comparing initial values, that is much harder to do then in the discrete case([Z3]). In the discrete case the sums are usually terminating, and checking initial conditions is trivial.

Many definite integral formulas that one can find in integral tables can be handled by the holonomic ”slow” algorithm, but few fall into the hyperexponential class. On the other hand, the present algorithm is very useful for proving theorems in the theory of special functions. All the classical families of orthogonal polynomials may be defined in terms of a generating function

∑_{n=0}^{∞} P_{n}(x)y^{n} = A(x, y) ,  (7)

where A(x, y) turns out to be hyperexponential ( except for the Jacobi polynomials). Then (7) can be rewritten as

P_{n}(x) = \frac{1}{2πi} ∫_{|y|=ρ} \frac{A(x, y)}{y^{n+1}}dy .

Since the integrand is hyperexponential, one can apply our algorithm to find, and prove, the differential equation satisfied by the sequence P_{n}(x) defined by (7). Similarly, one can use our algorithm for differencing under the integral sign, outlined in the previous section, to find and prove the linear recurrence equation satisfied by the family. Examples are given in appendices 2 and 3.

Another application is to proving Rodrigues-type formulas( e.g. [Ho], p.9). Here we use Cauchy’s formula

\frac{d^{n}}{dx^{n}} f(x) = \frac{n!}{(2πi)^{n}} ∫_{|y−x|=ρ} \frac{f(y)}{(y − x)^{n+1}} dy .

By using the algorithms of sections 6 and 7 we can get a differential equation and a recurrence equation respectively.

Another application of the algorithm of section 7 is for finding a pure recurrence relation in n for sequences a_{n} given by

a_{n} = coeff. of x^{b} in P (x)^{n ,}

where P (x) is a Laurent polynomial. This is useful in one-dimensional random walk problems.

Acknowledgement: We would like to thank Viggo Kann for helpful suggestions. We also thank Shalosh B. Ekhad and Sol Tre for executing our programs. Finally we would like to thank the two referees for their numerous helpful comments.

Appendix 1: A Maple Program Implementing The Algorithm Of Section 6

A.1.1 The program itself:

The listing for the program is given in [A-Z]. It is available from Zeilberger via e-mail (Please send your e-mail address via regular mail). Place that program in a file, and name it WHATYOUWILL.

A.1.2: A Typical Input File:

We will describe how to use the program in the UNIX operating system. Readers who use other operating systems will easily make the appropriate modifications.

Place the following input file in a file called INPUT. You should replace WHATYOUWILL by the full path name of the program file.

#Beginning Of Typical Input File read ‘WHATYOUWILL‘: expli:=1: A:=exp(-x2/y2-y2): B:=exp(- 2*x): ORDER:=1: duis(expli,conj,A,B,ORDER,disorcon): quit: #End of Typical Input File

In order to use the program, you type on your terminal

maple < INPUT > OUTPUT

After the execution is completed you can read the output file by entering ”cat OUTPUT”.

The program can find both a differential equation and a recurrence equation. The default is the former case, where disorcon is not set to anything. Then A(x,y) is the integrand, B(x) is the conjectured right side (if known), in which case expli is set to 1, and the program itself finds ”conj”, which is the first order differential operator annihilating B(x). If B(x) is not known explicitly, then we don’t set anything for B, but set expli:=0. If B(x) is conjectured to satisfy a certain differential equation conj(x, D)B(x) = 0, then we set conj := conj(x, D). If we don’t know what differential equation B(x) satisfies, then we set conj:=1. In every case the output is the differential operator gorem(x, D) and the hyperexponential function G(x, y) such that

gorem(x, D)conj(x, D)A(x, y) = D_{y}G(x, y) .

In other words S̄(x, D) = gorem(x, D)conj(x, D). ORDER is the guessed order of gorem(x, D) Note that A(x, y) must always be a function of x and y, and it is always assumed that y is the integration variable (so y is a global variable).

If a recurrence equation is desired, then the integrand, A(n, y), must depend on n and y, and the program finds a recurrence operator, expressed in terms of the shift operator N , Nf(n) := f(n+1). The above input file is the same, only now one sets disorcon:=discrete. In the discrete case set expli:=0 always.

Appendix 2: Some Sample Inputs and Outputs For Finding Differential Equations

1. With the following input the program proves

∫_{−∞}^{∞} exp(−x^{2}/y^{2} − y^{2})dy = \sqrt{π}e^{−2x ,}

after the appropriate initial conditions are checked.

#Begin input read ‘WHATYOUWILL‘: expli:=1: conj:=1: A:=exp(-x2/y2-y2): B:=exp(-2*x): ORDER:=1: duis(expli,conj,A,B,ORDER,disorcon): quit: #End input

The output was:

the conjectured rhs B is, exp(- 2 x) conj, the first order operator annihilating B(x) is, D + 2 Let R(x) be the integral of the following w.r.t y

exp(−\frac{x^{2}}{y^{2}} − y^{2})

The operator annihilating R(x) is the following times conj - 2 + D G(x,y) is the integrand times 2 1/y

2. With the following input, the program finds, and then proves, that

R(x) := ∫_{−∞}^{∞} exp(−x^{2}/y^{2} + ay^{2})dy

satisfies the differential equation (4a + D^{2})R(x) = 0.

#Begin input read ‘WHATYOUWILL‘: expli:=0: conj:=1: A:=exp(-x2/y2+a*y2): ORDER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: #End input The output was: conj is 1 Let R(x) be the integral of the following w.r.t y

exp(−\frac{x^{2}}{y^{2}} + ay^{2})

The operator annihilating R(x) is the following times conj

4a + D^{2}

G(x,y) is the integrand times 2 1/y 3. #Begin of input read ‘WHATYOUWILL‘: expli:=0: conj:=1: A:=exp(-x6/y4-y2): ORDER:=3: duis(expli,conj,A,B,ORDER,disorcon): quit: #End of input The output was conj is 1 Let R(x) be the integral of the following w.r.t y

exp(−\frac{x^{6}}{y^{4}} − y^{2})

The operator annihilating R(x) is the following times conj

216x^{5} + 7D − 12xD^{2} + 4x^{2}D

G(x,y) is the integrand times

−54\frac{(4x^{6} − y^{4} + 2y^{6})x^{5}}{y^{7}}

4. With the following input, the program finds, and then proves, the differential equation satisfied by the Legendre polynomials, when the later are defined by their generating function:

∑_{n=0}^{∞} P_{n}(x)t^{n} = (1 − 2xt + t^{2})^{−1/2} .

#Begin of input #The diff. eq. for Legendre out of the gen. fun. read ‘WHATYOUWILL‘: expli:=0: A:=1/(y(n+1)*(1-2*x*y+y2)(1/2)): ORDER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: #End of input The output was:

#The diff. eq. for Legendre out of the gen. fun. conj is 1 Let R(x) be the integral of the following w.r.t y

\frac{(1 − 2xy + y^{2})^{1/2}}{y^{n+1}}

The operator annihilating R(x) is the following times conj

n (n + 1) − 2 x D − (x − 1)(x + 1) D^{2}

G(x,y) is the integrand times

− \frac{(− n − 1 + 2 n x y + x y − n y^{2}) y}{(− 1 + 2 x y − y^{2}) n}

words used=156843, alloc=88064, time=29.37

5.

With the following input, the program proves the Mehler formula for Hermite polynomials (e.g. [Ho], p. 24, ex. 4).

The Mehler formula came to the limelights when Foata[Fo1] gave a beautiful combinatorial proof that inaugurated the currently very active area of combinatorial special function theory. Unlike the computer-generated proof below, that is a verification, Foata’s proof is a derivation, where every part of the generating function arises naturally. Furthermore, Foata’s proof lead naturally to a proof of a multi-variate generalization of the Mehler formula that was accomplished by Foata and Garsia[FG](see also [Fo2]). The theory of holonomic systems is incapable, even in principle, of proving multi-variate identities with an arbitrary number of variables.

#Begin input #The Mehler formula for Hermite polynomials read ‘WHATYOUWILL‘: expli:=0: A:=exp((2*x*y*z-y2*(x2+z2))/(2*(1-y2)))/((1-y2)(1/2)*y(n+1)): ORDER:=2:

duis(expli,conj,A,B,ORDER,disorcon): quit: #End of input

The output was:

#The Mehler formula for Hermite polynomials conj is 1 Let R(x) be the integral of the following w.r.t y

exp(\frac{2 x y z −y^{2} (x^{2} +z^{2})}{2 − 2 y^{2}})
(1 − y^{2})^{1/2} y^{(n +1)}

The operator annihilating R(x) is the following times conj

− n + x D − D^{2} G(x,y) is the integrand times y words used=138796, alloc=79872, time=26.87

6. With the following input the program proves the Doetsch identity (see [Fo2]). #To Prove the Doetsch formula for Hermite polynomials(see Foata,Adv. #Appl. Math. 2(1981),250-259.) read ‘WHATYOUWILL‘: expli:=0: A:=y(-n-1)*(1+2*y2)(-3/2)*(1+y*x+2*y2)*exp(y2*x2/(1+2*y2)): ORDER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: The output was: conj is , 1 Let R(x) be the integral of the following w.r.t y

\frac{y^{(− n −1)}(1 + y x + 2 y^{2})exp(\frac{y^{2}x^{2}}{1 + 2y^{2}} )}{(1 + 2 y^{2})^{3/2}}

The operator annihilating R(x) is the following times conj

−n + x D − D^{2}

G(x,y) is the integrand times y words used=194380, alloc=112640, time=40.10 7. With the following input, the program proves the Rodrigues formula for Jacobi polynomials ([Ra], p. 257) #Begin input #Proves Rodrigues’ Formula for Jacobi polynomials by finding the #diff. eq read ‘WHATYOUWILL‘: expli:=0: A:=(y2-1)n*(1-y)a*(1+y)b*(1-x)(-a)*(1+x)(-b)/(y-x)(n+1): OR-DER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: #End input The output was: #Proves Rodrigues’ Formula for Jacobi polynomials by finding the #diff. eq conj is the identity 1 Let R(x) be the integral of the following w.r.t y

\frac{(y^{2} − 1)^{n} (1 − y)^{a} (1 + y)^{b} (1 − x)^{(− a)} (1 + x)^{(− b)}}{(y − x)^{(n + 1)}}

The operator annihilating R(x) is the following times conj

−(n + a + 1 + b)n + (a + ax − b + xb + 2x)D + (−1 + x)(1 + x)D^{2}

G(x,y) is the integrand times

\frac{(n + 1)(− 1 + y) (1 + y)}{(− y + x)}

words used=397108, alloc=155648, time=96.80 8. With the following input, the program proves Euler’s integral representation for the hypergeometric function _{2}F_{1}(a, b; c; x) ([Ra], Th. 16, p. 47), if the latter is defined by its differential equation, or alternatively, finds the differential equation if it is defined in terms of the integral representation. #Begin of input #Proves Euler’s integral representation for 2F1 read ‘WHATYOUWILL‘: expli:=0: A:=(1-x*y)(-a)*y(b-1)*(1-y)(c-b-1): ORDER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: #End of input The output was: #Proves Euler’s representation for 2F1 conj is , 1 Let R(x) be the integral of the following w.r.t y

(1 − x y)^{(− a)} y^{(b − 1)} (1 − y)^{(c − b − 1)}

The operator annihilating R(x) is the following times conj

b a + (x + a x + x b − c) D + x (x − 1) D^{2}

G(x,y) is the integrand times

\frac{ay (− 1 + y)}{− 1 + x y}

words used=132760, alloc=77824, time=25.07 9. Rational functions can always be integrated either by contour-integration or by partial fractions, at least in principle. Here is an example how our program finds a differential equation satisfied by such an integral. #Begin input read ‘WHATYOUWILL‘: expli:=0: A:=x2/((x3+y3)*(1+y3)): ORDER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: #End input The program outputted the following differential operator. We omit G(x, y).

3 x^{2} + D + 5 x^{3} D + x^{4} D^{2} − x D^{2}

Appendix 3: Sample Inputs And Outputs For Finding Recurrence Equations In all the following, the same program was used, and the input files have the same form, only now the parameter ”disorcon” is set to discrete , and A should involve n and y. N is the shift operator in the n direction: Nf(n) := f(n + 1).

With the following input the program finds, and proves, the recurrence equation satisfied by the Legendre polynomials, if the latter are defined by the Rodrigues formula. #Beginning of input read ‘WHATYOUWILL‘: expli:=0: A:=(y2-1)n/(y-x)(n+1): ORDER:=2: duis(expli,conj,A,B,ORDER,discrete): quit: #End of input The output was: conj is , 1 Let R(n) be the integral of the following w.r.t y

\frac{(y^{2} − 1)^{n}}{(y − x)^{(n + 1)}}

The operator annihilating R(n) is the following times conj

− 4 n − 4 + 2 x (3 + 2 n) N + (− n − 2) N^{2}

G(n,y) is the integrand times

−\frac{(y − 1) (y + 1) (− y^{2} + 2 y x − 1)}{− y + x}

words used=161099, alloc=94208, time=31.43 2. With the following input, the program finds the recurrence equation satisfied by the Jacobi polynomials, when defined in terms of their Rodrigues’s formula. #Begin of input #Finds, and proves, the recurrence equation for the Jacobi polynomials #when they are defined in terms of the Rodrigues formula. read ‘WHATYOUWILL‘: expli:=0: A:=(y2-1)n*(1-y)a*(1+y)b*(1-x)(-a)*(1+x)(-b)/(y-x)(n+1): ORDER:=2: duis(expli,conj,A,B,ORDER,discrete): quit: #End of input The output is omitted( do it on your computer!). 3. With the following input the program finds, and proves the recurrence satisfied by the sequence R(n) :=constant term of (y + 1 + 1/y)^{n}. #Begin of input #Finds the recurrence for the constant term of (1/y+1+y)n: read ‘WHATY-OUWILL‘: expli:=0: A:=(1/y+1+y)n/y: ORDER:=2: duis(expli,conj,A,B,ORDER,discrete): quit: #End of input The output was: #Finds the recurrence for the constant term of (1/y+1+y)n: conj is , 1 Let R(n) be the integral of the following w.r.t y

\frac{(1/y + 1 + y)^{n}}{y}

The operator annihilating R(n) is the following times conj

− 3 n − 3 + (− 2 n − 3) N + (n + 2) N^{2}

G(n,y) is the integrand times

\frac{(y − 1)(y + 1)(1 + y + y^{2})}{y}

words used=130521, alloc=73728, time=23.02 4. With the following input the program finds, and proves, the recurrence equation satisfied by the Legendre polynomials, when they are given in terms of their generating function([Ra], p. 157, (1)). #The recurrence eq. for Legendre out of the gen. fun. read ‘WHATYOUWILL‘: expli:=0: A:=1/(y(n+1)*(1-2*x*y+y2)(1/2)): ORDER:=2: duis(expli,conj,A,B,ORDER,discrete): quit: The output was: #The recurrence eq. for Legendre out of the gen. fun. conj is , 1 Let R(n) be the integral of the following w.r.t y

\frac{1}{y^{(n + 1)} (1 − 2 x y + y^{2})^{1/2}}

The operator annihilating R(n) is the following times conj

n + 1 − x (3 + 2 n) N + (n + 2) N^{2}

G(n,y) is the integrand times

\frac{− 1 + 2 x y − y^{2}}{y}

words used=117996, alloc=65536, time=21.25

5. Euler’s Beta integral One of the referees suggested that we present the evaluation of the recurrence equation that leads to Euler’s derivation of the beta function. Note that the recurrence, and the proof, hold when m and n are thought of as real arguments, and then by standard uniqueness theorems, factorials can be replaced by Gamma functions.

The input file was: read ‘WHATYOUWILL‘ expli:=0: conj:=1: A:=yn*(1-y)m: ORDER:=1: duis(expli,conj,A,B,ORDER,discrete): quit: The output was: conj is , 1 Let R(n) be the integral of the following w.r.t y n m y (1 - y) The operator annihilating R(n) is the following times conj - n - 1 + (n + m + 2) N G(n,y) is the integrand times y (- 1 + y) words used=74417, alloc=43008, time=12.17 Appendix 4: Footnotes 1.The definition of ”holonomic function” The precise definition of holonomicity is as follows([Bj],p. viii). Recall that the principal symbol of a linear partial differential operator

P (x, D) = ∑_{|α|≤m} a_{α}(x)D^{α}

is the ”highest homogeneous part in the D’s”:

σ(P)(x, ζ) = ∑_{|α|=m} a_{α}(x)ζ^{α}

A function or distribution F is holonomic if the ideal of linear partial differential operators with polynomial coefficients annihilating it

I_{F} := {P (x, D) ; P (x, D)F = 0 }

has the property that the ”variety of common zeroes of the principal symbols of the elements of I_{F} ”

{(x, ζ)ϵC^{2n}; σ(P)(x, ζ) = 0 , P (x, D)ϵI_{F} }

has the minimal dimension it could have, namely n. (That such a variety cannot have dimension less than n follows from the celebrated Bernstein inequality ([Bj]). As an application we will show that a wide class of hyperexponential functions is holonomic. Lemma A.4.1: Let S(x, y) be a rational function: S(x, y) = P (x, y)/Q(x, y) (P, Q polynomials) such that the varieties (x, y)ϵC^{2} ; P (x, y) = 0, Q(x, y) = 0 and (x, y)ϵC^{2} ; Q_{y}(x, y)/G(x, y) = 0 , Q_{x}(x, y)/G(x, y) = 0 (G := gcd(Q_{x}, Q_{y})) and (x, y)ϵC^{2} ; G(x, y) = 0, Q(x, y) = 0 are zerodimensional (i.e. consists of a finite set of points) then F (x, y) := exp(S) is holonomic. Proof: F (x, y) = exp(S(x, y)) is easily seen to satisfy the following differential equations

(Q^{2}D_{x} − (P_{x}Q − P Q_{x}))F (x, y) ≡ 0  (i)

(Q^{2}D_{y} − (P_{y}Q − P Q_{y}))F (x, y) ≡ 0 .  (ii)

(P_{y}Q − P Q_{y})(i) − (P_{x}Q − P Q_{x})(ii) yields that F (x, y) also satisfies ( since F (x, y) is C^{∞} we can ”divide” by Q):

((P_{y}Q − P Q_{y})D_{x} − (P_{x}Q − P Q_{x})D_{y})F (x, y) ≡ 0 .  (iii)

Thus the variety of common zeroes of the principal symbols of I_{F} is a subset of the variety

{(x, y, ζ_{1}, ζ_{2})ϵC^{4} ; Q^{2}ζ_{1} = 0, Q^{2}ζ2 = 0, ((P_{y}Q − P Q_{y})ζ_{1} − (P_{x}Q − P Q_{x})ζ_{2} = 0} ,

which is easily seen to be two dimensional. It follows that the variety of zeroes has dimension ≤ 2, and by Bernstein’s inequality it is actually 2 and F is holonomic.

The same method of proof yields that under the condition of the lemma (P (x, y)/Q(x, y))^{α} is always holonomic, and, since holonomicity is closed under multiplication([Z1]) also any product of such functions.

Let us remark that since holonomic functions are D-finite ([Z1], 4.1), the existence of (3) could have also been established by Lipshitz’s[Li] method.

2. The Secret Behind Many Definite Integrals

The secret behind many definite integrals is a lemma that can be found in [He]. It is probably not the earliest reference. It is likely that the trick was known to Euler or at least to Cauchy (Herr quotes Cauchy, among others , in the preface ).

Lemma:

∫_{−∞}^{∞} f(x)dx = ∫_{−∞}^{∞} f(x − \frac{a}{x})dx , a > 0

Proof: t = x − \frac{a}{x} if and only if

x = \{{t/2 + \sqrt{t^{2}/4 + a} ifx > 0 \atop t/2 − \sqrt{t^{2}/4 + a} ifx < 0}

so,

dx = \{{1/2 + t/ \sqrt{t^{2}/4 + a} ifx > 0 \atop 1/2 − t/ \sqrt{t^{2}/4 + a} ifx < 0}

We thus have

∫_{0}^{∞} f(x − \frac{a}{x})dx = ∫_{−∞}^{∞} f(t)(1/2 + t/ \sqrt{)}dt ,

∫_{−∞}^{0} f(x − \frac{a}{x})dx = ∫_{−∞}^{∞} f(t)(1/2 − t/ \sqrt{)}dt .

The proof of the lemma is completed by adding the above two equations.

3. How Feynman would have evaluated the integral of section 3 Since the derivation of the integral

R(x) := ∫_{−∞}^{∞} exp(−\frac{x^{2}}{y^{2}} − y^{2})dy ,

given in ([Ed], p. 220) is ”not as well known as it should be”, and is ”the argument that Feynman would have used directly”, one of the referees suggested that we reproduce this derivation. It runs as follows.

R(x) = 2 ∫_{0}^{∞} exp(−\frac{x^{2}}{y^{2}} − y^{2})dy ,

so

R′(x) = −4x ∫_{0}^{∞} exp(−\frac{x^{2}}{y^{2}} − y^{2})\frac{dy}{y^{2}} ,

= −4x ∫_{∞}^{0} exp(−t^{2} − \frac{x^{2}}{t^{2}} )(\frac{−xdt}{t^{2}})\frac{t^{2}}{x^{2}} ,

= −2R(x) .

This integrates to R(x) = Ae^{−2x} and x = 0 gives R(0) = \sqrt{π}. References [A-Z] Almkvist, G., and Zeilberger, D., A MAPLE program that finds, and proves, recurrences and differential equations satisfied by hyperexponential definite integrals, SIGSAM Bulletin, to appear.

[Ba] Bailey, W. N. , ”Generalized Hypergeometric Series”, Cambridge Math. Tracts 32, Cambridge University Press, London, 1935. (Reprinted: Hafner, New York, 1964 .)

[Be] Bernstein, I. N. , Modules over a ring of differential operators, study of the fundamental solutions of equations with constant coefficients, Functional Analysis and Its Applications, Vol 5, No. 2(1971), Russian original :1-16, English translation:89-101.

[Bj] Bjo:rk, J. -E., ”Rings Of Differential Operators”, North-Holland, Amsterdam, 1979.

[Bu] Buchberger ,B. Grobner bases - an algorithmic method in polynomial ideal theory, Chapter 6 in N.K.Bose (ed.): ”Multidimensional System Theory”, D.Reidel, Publ. Comp., 1985.

[Co] Coxeter, H. S. M., ”Twelve Geometric Essays”, Southern Illinois University Press, Carbondale, 1968, preface.

[DST] Davenport, J. H., Siret, Y., Tournier, E., ”Computer Algebra”, Academic Press, London, 1988.

[Ed] Edwards, J., ”A treatise on the Integral Calculus”, vol. 2, Macmillan, New York, 1921.

[Eh] Ehlers, F. The Weyl algebra, in: ”Algebraic D-modules”, by A. Borel et. al., Academic Press, Boston, 1987.

[F] Feynman, R. P., ”Surely You are Joking, Mr. Feynman!”, pp. 71-72, paperback, Bantam books 1986. originally published by W. W. Norton, 1985

[Fo1] Foata, D., A combinatorial proof of the Mehler formula, J. Comb. Theory, ser. A, 24(1978), 250-259.

[Fo2] Foata, D., Some Hermite polynomial identities and their combinatorics, Advances in App. Math. 2(1981), 250-258.

[FG] Foata, D., and Garsia, A. M., A combinatorial approach to the Mehler formulas for Hermite polynomials, in: ”Relations between combinatorics and other parts of mathematics” (D. K. Ray-Chaudhury, Ed.), Proceedings of Symposia in Pure Math. No. 34, 163-179, Amer. Math. Soc., Providence, 1979.

[Gal] Galligo A., Some algorithmic questions on ideals of differential operators, in: Lecture Notes in Computer Science 204, ( edited by Bob F.Caviness, EUROCAL ’85, vol.2), Springer, Berlin.

[Go] Gosper, R. W., Jr., Decision Procedure of Indefinite Summation, Proc. Natl. Acad. Sci. USA75, 40-42 , 1978.

[GKP] Graham, R., Patashnik, O., and Knuth, D. E., ”Concrete Mathematics”, Addison Wesley, Reading, 1989.

[He] Herr, J. P., ”Lehrbuch der hoeren Mathematik”, Bd 2, p.356, Wien, 1857.

[Ho]Hochstadt, H., ”Special Functions of Mathematical Physics”, Holt, Rinehart and Winston, New York, 1961.

[La] Lafon, J. C., Summation in finite terms, in: ”Computer Algebra Symbolic and Algebraic Computation, edited by B. Buchberger, G. E. Collins, and R. Loos , in cooperation with R. Albrecht, Springer-Verlag, Vienna, second edition, 1983.

[Li] Lipshitz, L., The diagonal of a D-finite power series is D-finite, J. Algebra 113(1988), 373-378.

[Me] Melzak, Z. A.,”Companion to Concrete Mathematics”, John Wiley, New York, 1973.

[No] Norman, A. C., Integration in finite terms, in: ”Computer Algebra Symbolic and Algebraic Computation, edited by B. Buchberger, G. E. Collins, and R. Loos , in cooperation with R. Albrecht, Springer-Verlag, Vienna, second edition, 1983.

[Ra] Rainville, E. D., ”Special Functions”, Chelsea, Bronx, 1971. Originally published by Macmillan, 1960.

[Ri] Risch, R. H. , The solution of the problem of integrating in finite terms, Bull. AMS 76 (1970).

[Ta] Takayama, N., An approach to the zero recognition problem by the Buchberger algorithm, preprint, Kobe University.

[Tr] Trager, B. M.,”On the Integration of Algebraic Functions”, Ph. D. thesis, MIT, 1985. [W-Z] Wilf, H. S., and Zeilberger, D., Towards computerized proofs of identities, Bulletin of the Amer. Math. Soc., to appear.

[Z1] Zeilberger, D., A holonomic systems approach to special function identities, preprint, Drexel University.

[Z2] Zeilberger, D., A fast algorithm for proving terminating hypergeometric identities, Discrete Math., to appear.

[Z3] Zeilberger, D., The method of creative telescoping, preprint, Drexel University.