Symbolic definite integration: methods and open issues

Daniel Lichtblau

Wolfram Research, Inc.

Abstract

The computation of definite integrals presents one with a variety of choices. There are various methods such as Newton-Leibniz or Slater's convolution method. There are issues such as whether to split or merge sums, how to search for singularities on the path of integration, when to issue conditional results, how to assess (possibly conditional) convergence, and more. These various considerations moreover interact with one another in a multitude of ways. Herein we discuss these various issues and illustrate with examples.

Introduction

In principle, algorithmic integration is fairly well understood. There are algorithms to handle various sorts of integrands and integration bounds. But, in the words of librettist Ira Gershwin, "It ain't necessarily so" [Gershwin 1935],  and actual practice is an altogether different situation. Troublesome areas include: different situation. Troublesome areas include:

Determining (possibly conditional) convergence

Recognizing and handling cancellation of additive singularities

Intricacies of the algorithmic implementation, e.g. need for recursion or reliance on powerful transformations

Handling of parameters

Working with legacy code.

Basic Structure of Integrate Code in Mathematica

Indefinite Integration

Uses: Risch algorithm

Table lookup for elliptic integrals, antiderivatives that require special functions, and antiderivatives of special function.

Definite Integration

Uses: Contour integration

Newton-Leibniz: Find antiderivative, evaluate at endpoint

Slater convolution: Write integrand as product of two MeijerG functions, integrated from 0 to infinity. Form result as a MeijerG.

Extra considerations of definite integration

Check for convergence

Check for singularities along path of integration

Formulate provisos to be satisfied by parameters in order that result hold

Brief Descriptions and Examples of the Two Methods

Newton-Leibniz

First we find an antiderivative to the integrand. We then look at the integration path, finding candidate singular points. Either issue a proviso on parameters that guarantee the point is not a "bad" one or else split the integration path at that point and use limits as we approach from either side. Sum appropriate signed values to get the resulting definite integral.

Here is a classical example that incorrectly gave zero in a prior version of Mathematica. Note that one path, from -1 + i to -1 - i, must be split where it crosses the (negative) real axis.

Integrate[1/z, {z, 1 + I, -1 + I, -1 - I, 1 - I, 1 + I}]

Slater Convolution via MeijerG Products

The idea in this case is to make suitable transformations so that we have an integrand in, say, the variable z, as a product (z^k)MeijerG[(Subscript[expression, 1])]MeijerG[(Subscript[expression, 2])] where we integrate from 0 to infinity.

Several common functions such as trigs, various exponentials, Bessel functions, and certain common algebraics e.g. (((a + b (x^n)))^m) may be represented as MeijerG functions. Others may be brought to MeijerG form after some suitable transformation. Do a table driven conversion in an effort to obtain the desired form.

Once in the factored form above we go into convolution code. This formulates the integrated result as a new MeijerG expression via the Slater theorem. It then performs some manipulations in an effort to recast in terms of better known functions.

Here is a simple example.

val = 
      Integrate[
        BesselJ[2, x] (x^2) (E^(-x + 2)), {x, 
          0, ∞}, GenerateConditions -> False]
nval = 
      NIntegrate[
        BesselJ[2, x] (x^2) (E^(-x + 2)), {x, 
          0, ∞}]
Chop[val - nval, (1/(10^8))]

This operated by first converting the integrand to an "inert" form representing the integrand product as

(x^2) MeijerG[{{}, {}}, {{0}, {}}, 
          x] MeijerG[{{}, {}}, {{1}, {-1}}, ((x^2)/4)]

Software Engineering Problems Associated with this Body of Code

List of problem spots:

Different modules, written by different people, are often not on the best of speaking terms.

Much is legacy code. It was written over a period of about 15 years, by several different peopleNo one person understands all of it, and some parts are no doubt not understood by anyone.

Some of the steps used are quite fragile. Specifically they may be sensitive to small changes in seemingly unrelated functions such as Together, Apart (which is in fact stitched at the hip to Integrate), Factor, Solve, Simplify, and others. I have come to conclude that this issue really is endemic to symbolic integration and not simply an artifact of implementation details.

There are tradeoffs to be made between speed and power.

Many classes of integration problem can make use of assumptions regarding, say, behavior of integrand at infinity. The code must be sufficiently careful to delimit what inputs enter handlers that may utilize such assumptions.

Symbolic integration is one of the most complicated pieces of machinery to be found in algorithmic mathematics. It is (even m ore than usual) vital that the pieces be adequately documented.

When I began work on this body of code, more than one out of every four open bugs in the Mathematica kernel (with around 2000 functions) was in the category of definite integration Such a mass poses several problems in and of itself: difficult to start from scratch due to widespread usage, difficult to perform adequate triage, the massive overhaul needed implies considerable damage will ensue for some time, and the likelihood that the scope of trouble exceeds the capabilities of any one developer.

Issues in the Implementation of Indefinite Integration

We will discuss in brief some issues that tend to be specific to indefinite integration.

The Curse of Recursing

Often an integral may be broken into two parts. The technical term for the first part is "done" while that for the second is "not done". The second term might be further rewritten in ways that return one to the original problem or a variant thereof, thus leading to recursive splitting.

Another common cause of descent into the infinite is the utilization of pairs of inverse transformations.

Out on a Limb with a Cut Branch

Frequently we require transformations that bring into play multivalued functions. The consequence is that we may arrive at an antiderivative that is only correct up to a piecewise multiplicative constant. Hence we now attempt to restore the proper factor. It is not trivial to recognize in all cases how to correctly reverse effects of such transformations.

Transformations

In many cases it is well understood how one might obtain an antiderivative for a particular class of integrand e.g. rational functions of trigs. But the needed transformations must be applied carefully in order to avoid potential explosion in intermediate complexity or that of the final result.

Issues in the Implementation of Definite Integration

In this section we arrive at the main focus of this report. Within symbolic definite integration one encounters a wide array of issues. In the subsections below I will endeavor to present and illustrate a few.

To Expand or Not Expand, Definitely

When one has certain types of input it makes sense to expand over summands and integrate term-by-term. In other cases this can be a very bad thing to do. We will illustrate with a few examples.

First, a "simple" integrand involving a rational function.

p = 
        7 (x^13) + 10 (x^8) + 4 (x^7) - 7 (x^6) - 4 (x^3) - 4 (x^2) + 3 x + 
          3;
q = 
        (x^14) - 2 (x^8) - 2 (x^7) - 2 (x^4) - 4 (x^3) - (x^2) + 2 x + 
          1;
result = (∫_2^5) (p/q) \[DifferentialD]x

We do a numeric check that this is correct.

N[result] - NIntegrate[p/q, {x, 2, 5}]

Now we expand and integrate over each term, then try to simplify the result. It is a mess.

result2 = 
      Simplify[Map[Integrate[#, {x, 2, 5}] &, Expand[(p/q)]]]

N[result2 - result]

We see that the mess at least is correct. But it is obvious that we are better served by avoiding such expansion in this example.

Here is an example where one must not expand. The reason is that the summands have cancelling singularities.

ii = (∫_0^∞)((E^((I/x))) - 
                1) Sin[x] \[DifferentialD]x

Again we validate this with a numeric check.

N[ii] - 
      NIntegrate[((E^((I/x))) - 1) Sin[
            x], {x,  .001, 1000}, WorkingPrecision -> 25]

Clearly we cannot expand and handle summands separately.

Map[Integrate[#, {x, 0, Infinity}] &, 
      Expand[((E^((I/x))) - 1) Sin[
            x]]]

By contrast, here is a case where preexpansion helps considerably.

integrand = ((a c Cos[t])/(g s)) + ((b q Cos[2 t])/(c f)) + ((c Cos[3 t])/(a d)) + ((d f Cos[4 t])/(a h n)) + ((e q Cos[5 t])/(a g)) + ((f l Cos[6 t])/(m r w)) + ((b g Cos[7 t])/(n o x)) + ((h Sin[t])/(b c)) + ((i Sin[2 t])/(e h r)) + ((j y Sin[3 t])/(l p)) + ((d k Sin[4 t])/(c j)) + ((a l m Sin[5 t])/(b f h s)) + ((m p Sin[6 t])/(j k)) + ((n q Sin[7 t])/(c x));
Timing[(∫_0^(2 π))integrand \[DifferentialD]t]

Earlier versions which did not expand  took about a minute to handle the same integration. Transformations involving trig summands can be costly.

Heuristics Involving Method Choices and Transformations of Integrand

Quite often one encounters an integrand that might be handled in different ways. In particular we might try for an antiderivative evaluation or use a convolution approach.

Here is an example where use of convolution approach gives a result faster.

Timing[(∫_0^∞) (1/(x^8)) Sin[x] Sin[
                (x/3)] Sin[(x/5)] Sin[(x/7)] Sin[(x/9)] Sin[(x/11)] Sin[
                (x/13)] Sin[(x/15)] \[DifferentialD]x]

A related issue of heuristics is that frequently an integrand may appear amenable to any of several different transformations. For example, for a quadratic radical function of the integration variable we might want to transform to a MeijerG function, or do  a linear change of coordinates in an effort to simplify or remove the radical.

Another related issue is that some transformations may not be valid for all possible values of a parameter. For example, one might wish to write (((a + b (x^d)))^n) as (a^n) (((1 + ((b (x^d))/a)))^n) (for integration variable x) in order to convert to a MeijerG form. This transformation is of course not valid for all values of the parameter a. We then must choose between issuing a conditional answer, or abandoning this approach and trying another that might give an unconditional result.

Generations of Generation of Provisos (Results That Depend on Conditions)

Below I indicate ways in which generation of conditions can be problematic. Several are illustrated with older versions of Mathematica because improvements render the specific examples obsolete. The ideas behind them are general and still lurk.

Propagation of Conditions for Multiple Integrals

For multiple integration one may want to propagate conditions from inner to outer integrations. This is problematic if the conditions generated cannot adequately be simplified. On the other hand if we fail to generate and propagate them we run a greater risk of obtaining a bad result.

Here is an example where in version 5.0 we required conditions in order to propagate singularity information.

result = 
      Integrate[Abs[x - y]^n, {x, 0, 1}, {y, 0, 1}, 
        GenerateConditions -> True]

A numeric check validates this.

(result /. n -> 3.2)  -  
      NIntegrate[Abs[x - y]^3.2, {x, 0, 1}, {y, 0, 1}]

Note that if we do not insist on generation of conditions we get a result that is quite obviously incorrect; the correct one must lie between zero and one for positive n.

Integrate[Abs[x - y]^n, {x, 0, 1}, {y, 0, 1}, GenerateConditions->False]

(∫_0^1)(∫_0^1) Abs[
                  x - y]^n \[DifferentialD]y \[DifferentialD]x /. 
      n -> 3.2

Excessive Conditions

Some code may give generation of unneeded conditions. For example, MeijerG convolution will require that certain values lie in "wedges" emanating from the origin, and Newton-Leibniz methods may issue conditions based on the specific form of the antiderivative.

Here is a simple example using version 5.0. Special case code for handling exponentials via convolution wants to insist that a parameter take a negative real part.

Integrate[(E^(-((a (((x - b))^2))/(sigma^2)))), {x, -∞, ∞}, 
      Assumptions -> {a > 0, sigma ∈ Reals}]

The current behavior is nicer.

Integrate[(E^(-((a (((x - b))^2))/(sigma^2)))), {x, -∞, ∞}, 
      Assumptions -> {a > 0, sigma ∈ Reals}]

Necessity of Condition Generation

Often one may wish to ignore conditions. A problem is that results can depend on what path is taken in the code. If the integration path is split into parts and if the parts take routes in the code that make different assumptions about parameters, then results might be entirely incorrect. One common manifestation is incorrect cancellation. We illustrate with the preceding example.

Integrate[(E^(-((a (((x - b))^2))/(sigma^2)))), {x, -∞, ∞}, GenerateConditions -> False, 
      Assumptions -> {a > 0, sigma ∈ Reals}]

The upshot is that partial results cancelled, with no indication that they require conflicting conditions to hold.

Genericity of Generated Conditions

Conditions generated may be only generically correct. Here is an example that exhibits this phenomenon. Specifically, when Im[a] == Im[b] (that is, a vertical integration path) we have a problem because that their difference appears in some denominators.

(∫_a^b) Log[x] \[DifferentialD]x

Assessment of Convergence

Testing for convergence is closely related to generation of conditions insofar as the former can depend on parameter values. But convergence testing, even in the absence of symbolic parameters, is no easy matter (complicated by oscillatory factors, possible cancellation of singular terms, and so forth). Code that tests integrands for convergence ad hoc.

Here are some simple examples that may help to give an idea of the difficulties lurking within convergence assessment. Mathematica will evaluate them correctly but this was not always the case.

(∫_0^∞) Cos[(x^2)]/Log[
              x] \[DifferentialD]x

(∫_0^∞) (E^x) Sech[
            a x] x \[DifferentialD]x

The first diverges, but not due to problems at either endpoint (there is a pole at 1). The second converges conditionally, and version 4 even knows this, but then gets the condition wrong.

Below is an example that gives an idea of what is involved in sorting out conditions for convergence. We need to look at possibilities of oscillatory, exponentially growing, and exponentially damped factors in order to figure out the correct conditions on parameters. In this case the original assumption suffices to guarantee convergence.

Integrate[(Sin[betty*r]*Sin[gramma*r])/
        E^(alf*r^2), {r, -Infinity, Infinity}, 
      Assumptions -> Re[alf] > 0]

Here is a slightly more difficult variant, one that requires nontrivial conditions on parameters beyond what is given in the assumptions option.

Integrate[(Sin[betty*r^2]*Sin[gramma*r^2])/
        E^(alf*r^2), {r, -Infinity, Infinity}, 
      Assumptions -> Re[alf] > 0]

As another example we show an integral that converges. The integrand is a sum of two terms, each of which separately will diverge. If the code splits this then it must recognize that there will be a cancellation of singular parts.

Integrate[(Log[x] - Log[a])/(x - a), {x, 0, a}, 
      Assumptions -> a > 1]

Difficulties Involving Parameters and Detection of Singularities

One issue is in finding singular points on the integration path. Even something as simple as a trig function can be problematic. There are other pits lurking beneath the surface of the swamp.

Parametric Singularities

The presence of elliptic functions in the antiderivative usually removes any hope of correctly detecting parameter dependent singularities.  Below is an example.

(∫_0^((π/2)))(-(m MeijerG[{{1, 1}, {}}, {{(1/2), 
                          (3/2)}, {}}, (1/(1 - m))]) + π EllipticK[
                  m] (1 + (-1 + 
                          m) Sin[y]^2))/(π (Sqrt[1 + (-1 + m) Sin[y]^2])) \[DifferentialD]y

Whether the result is worth the screen real estate it occupies is open to debate.

Transcendentals and Singularity Detection

Often finding singularities depends on finding roots of expressions. These need not be algebraic, and a general purpose transcendental root finder is a nontrivial undertaking (also it might be slow). In the example below, we find the bad point at the origin, and correctly decide the integral is divergent.

(∫_(-(1/2))^((1/2)))((-2 + (x^2) + 
              20 Cos[x])/(((-6 x + (x^3) + 60 Sin[x]))^2)) \[DifferentialD]x

But minor modifications give rise to a singularity that is not algebraic number nor a "convenient" multiple thereof (e.g. π times a rational).

ii = (∫_1^5)((-2 + (x^2) + 20 Cos[x])/(-6 x + 
                (x^3) + 60 Sin[x])) \[DifferentialD]x
N[ii]
NIntegrate[((-2 + (x^2) + 20 Cos[x])/(-6 x + (x^3) + 60 Sin[x])), {x, 1, 5}]

Clearly this diverges. The denominator has a pole in the integration path, and it is a (nonalgebraic) root of a transcendental equation.

FindRoot[-6 x + (x^3) + 60 Sin[x] == 0, {x, 2}]

Parametrized Singularities in Multiple Integrals

Parametrized singularities can cause trouble, particularly in multidimensional integration. Here is an example from [Math Horizons 2002]. We begin with a small variation on the actual problem.

(∫_0^π)(∫_0^π) Abs[
            Sin[x - y]] \[DifferentialD]y \[DifferentialD]x

Now we double the region of integration in both directions.

(∫_0^(2 π))(∫_0^(2 π))Abs[
            Sin[x - y]] \[DifferentialD]y \[DifferentialD]x

The second one is off by a factor of 2. The problem is that one singularity line, x == y, is recognized. But the broken line x == y + π modulo 2  π goes unrecognized. Hence we do not handle ranges correctly, splitting only into two rather than three pieces in the first level of integration. Note that this is in version 5.0 of Mathematica; later versions do get the correct value of 8 π.

Branch Cuts Intersecting the Integration Path

It is important to assess whether an integration path crosses a branch cut of an antiderivative (so that we might split the path into segments). The example below does this in order to get the correct result.

i1 = (∫_(-1 - I)^(-2 + I))Log[x] \[DifferentialD]x
N[i1]
NIntegrate[
      Log[x], {x, -1 - I, -2 + I}]

Version 4 of Mathematica failed to catch the crossing and gave a result with imaginary part off by 3 π. Needless to say, this problem becomes vastly more difficult if the crossing might or might not exist based on parameter values. One might try the example below to get an indication of what might be a reasonable result.

(∫_(-1 - I)^y) Log[
          x] \[DifferentialD]x

Problems with Algebraic Manipulation

Here is an example that was quite wrong prior to version 5.1 (we show the result from version 5.0).

ii = (∫_0^∞) 
            2 y (E^y) BesselK[
              0, (Sqrt[3]) y Sign[y]] \[DifferentialD]y
N[ii]
NIntegrate[
      2 y (E^y) BesselK[0, (Sqrt[3]) y Sign[y]], {y, 
        0, ∞}]

The problem was in the innards of the convolution code. At a key step we made use of an exponential of the form exponential((-2 i) π - 
              log((4/3)) + 2 log(2)) and we subsequently required a square root. Once we replace exponential with Exp we lose track of a factor of (-1) and this gave rise to a bad result. This example is now repaired but the general problem remains of how to find and forestall this phenomenon.

The next examples have bugs in current versions due to branch cut issues not correctly handled in convolution code. They show that not all changes can be viewed as "progress".

Integrate[(E^(-I x t)) ((E^((1/3) (-I) (x^3))) - 
            1) (x^(-(3/2))), {x, 0, ∞}, 
      Assumptions -> t > 0]

A numerical check will reveal that this is simply not correct.

The next example does not even involve a symbolic parameter.

ii = Integrate[ArcSin[v]*Log[v]/v, {v, 0, 1}]
N[ii]
NIntegrate[ArcSin[v]*Log[v]/v, {v, 0, 1}]

General Considerations

Here are some of the questions that require thought in the design and implementation of definite integration. They have emerged from prolonged overhaul of our existing code base.

Given a choice of methods, which should one attempt first? This can have serious repercussions in terms of speed. For example, what might be a fast MeijerG convolution can be very slow to evaluate as an indefinite integral followed by extraction of limiting values. And of course the opposite can happen. Or the methods might be comparable in speed but give results of vastly different complexity of form.

When or how should one presimplify the input?

When or how should one simplify the result?

Special case methods may be very helpful for speed or simplicity of result. But they also proliferate opportunities for bugs. Decisions of what classes to handle in these ways requires careful consideration.

Some methods require intrinsically "slow" technology. For example, refinement of conditions (sometimes essential in order that they not blow up) may require some level of CAD support. Even limit extraction for Newton-Leibniz methods can be slow. When to apply such technology? How to prevent hangs?

We have use time constraints to prevent hangs. New set of problems:
Asynchronous interrupt handling, required by TimeConstrained, is imperfect and in rare cases will cause a kernel crash. 
Results now can be platform dependent.

Indefinite integration is (generally) more powerful in version 5 than in the past. This has come at the price of speed: it simply tries more transformations and the like. This in turn means all Newton-Leibniz code is at risk of getting slower.

Provisional moral, or maybe conundrum: The more things improve, the more opportunity for closely related things to deteriorate.

Summary and Directions for Future Work

Development and overhaul of Integrate has had successes and at least some setbacks. One should hope these latter may prove to be temporary.

Among the successes, virtually all infinite recursion problems from earlier versions have been fixed. Premature bailout has likewise been addressed. Numerous problems in convergence assessment and singularity detection have likewise been fixed. Much of the idiosyncratic ad hoc code has been removed or rewritten. Most open problems fall into the various categories described above, which, while large, is better understood than the integration swamp that we had four years ago.

This has come at a cost of speed. Moreover we now have some level of platform dependent behavior. Work remains to iron out bugs in convolution code and elsewhere. Heuristics to determine use of expansion, simplification, refinement of conditions, and the like are crude and certainly imperfect. Modularity of the code, which affects usability, maintaince and further development, has improved but is likewise still far from perfect. The same may be said for its documentation.

Future work in definite integration within Mathematica might go in a few directions.

Strengthening the MeijerG convolution code:
Range of integrands covered
Correct handling of convergence and branch cut issues.

Better assessment of convergence conditions. This will include splitting the GenerateConditions option into two or three distinct ones so that meanings are more clear. At present the code may use it to test for convergence and look for parameter conditions that guarantee convergence

Figuring out a platform independent way to measure work done, and base an interruption scheme thereon so that results do not depend on processor speed.

Continued debugging based on existing known bugs in that code. At the time of this writing we are engaged in a sweep over all the open bug reports in definite integration.
Approximately forty percent have been closed
It is anticipated that perhaps another thiry percent will be closed by the time the sweep is concluded.
This will bring the quantity to a point where a careful categorization of specific problem areas becomes realistic. Some of these specific areas will simply need point fixes. What remains should be amenable to rewriting.

Continued addition to the test suite so as to avoid large scale breakage during overhaul. In past a problem that led to the massive bugs list was that we did not have an adequate test suite. Hence bug fixes might well cause breakage elsewhere that went long undetected (code decay). At present, virtually all fixed bugs get made into new tests.