Examples in common notation 

>	Int(1/(1+x^2), x) = int(1/(1+x^2), x);

 

>	Int(cos(x), x) = int(cos(x), x);

 

>	Int(1/sqrt(1-x^2), x) = int(1/sqrt(1-x^2), x);

 

>	Int(arccosh(x), x) = int(arccosh(x), x);

 

>

 

The same examples in exp-log notation 

>	Int(1/(1+x^2), x) = convert(int(1/(1+x^2), x), 'expln');

 

>	Int(1/2*exp(I*x) + 1/2*exp(-I*x), x) = int(1/2*exp(I*x) + 1/2*exp(-I*x), x);

 

>	Int(1/sqrt(1-x^2), x) = convert(int(1/sqrt(1-x^2), x), 'expln');

 

>	Int(convert(arccosh(x),'expln'), x) = int(convert(arccosh(x),'expln'), x);

 

>

 

Example 1.2 

>

restart;

 

>	Int(1/(x^2+2*x+1), x) = int(1/(x^2+2*x+1), x);

 

Remark: No terms required in the above example. 

>	Int(1/(x^3+x), x) = int(1/(x^3+x), x);

 

Remark: No algebraic number extensions required in the above example. 

>	Int(1/(x^2-2), x) = `int/risch`(1/(x^2-2), x);

 

Remark: In the latter example, the constant field needed to be extended by the algebraic number . 

>

 

Minimal algebraic extensions are desired for efficiency. 

For example, from above noting that    would lead to the following representation of the integral: 

 . 

Goal:  Avoid introducing algebraic numbers except when necessary. 

Hermite Reductions for the Rational Part 

We want to reduce the problem as follows: 

where    and where    is monic and square-free. 

Hermite's method is a classical method to achieve this reduction. The reduction is accomplished using basic polynomial operations and without introducing any algebraic number extensions.  The method is summarized here. 

Step 1:  Euclidean division with remainder yields 

where    or   . 

Step 2:  Square-free factorization yields 

where    for all , 

and where    for . 

Step 3:  Partial fraction expansion yields 

where   . 

We now have 

. 

Step 4:  Hermite reductions must be applied for each integrand where the denominator is a power > of a square-free factor. 

After this step, each integral that remains will have a square-free denominator. 

Consider one integral    with > . 

Apply the Extended Euclidean Algorithm to solve:    . 

This yields    . 

Integration by parts yields 

 . 

Note the contribution to the rational part of the result. 

Hermite reductions can be continued until only square-free denominators remain. 

Example 1.3 

>	p := x^6+5*x^5+10*x^4+50*x^3+25*x^2+127*x-1:

 

>	q := x^4+10*x^2+25:

 

The integral we wish to compute is 

>	Int(p/q, x);

 

>	r := rem(p, q, x, 'P');

 

>

P;

 

>	int_P := int(P, x);

 

Euclidean division with remainder has yielded: 

>	Int(p/q, x) = int_P + Int(r/q, x);

 

>	q_sqrfree := convert(q, sqrfree, x);

 

>	q2 := x^2 + 5;

 

In this case, the square-free factorization of the denominator is 

 with    and   . 

The problem has been reduced to 

>	Int(p/q, x) = int_P + Int(r/q2^2, x);

 

Apply Hermite reduction to the remaining integral which is of the form    with . 

The Extended Euclidean Algorithm is invoked as follows. 

>	gcdex(q2, diff(q2,x), r, x, 's', 't');

 

>

s;

 

>

t;

 

>

j := 2:

 

The Hermite reduction formula takes the form 

>	Int(r/q2^j, x) = -t/((j-1)*q2^(j-1)) + Int((s+diff(t,x)/(j-1))/q2^(j-1), x);

 

and therefore the original integral has been reduced to 

>	Int(p/q, x) = int_P + rhs(%);

 

>

 

Rothstein-Trager Theorem for the Logarithmic Part 

The problem is to express    where    and    is square-free. 

Classical approach:  Factor B(x) over its splitting field, and then a partial fraction expansion determines the terms. 

But such a factorization may be much more expensive than necessary. 

Theorem 2  (Rothstein-Trager Theorem) 

where are the distinct roots of   

and where are the polynomials   . 

Moreover, this expresses the integral using the minimal algebraic extension field. 

Remark: The roots of may introduce new algebraic numbers. 

Example 1.4 

From Example 1.3 the remaining integral is 

>	A := -1/10:  B := x^2 + 5:

 

>	Int(A/B, x);

 

>	R := resultant(A - z*diff(B,x), B, x);

 

>	(c1,c2) := solve(R=0, z);

 

>	v1 := gcd(A - c1*diff(B,x), B);

 

>	v2 := gcd(A - c2*diff(B,x), B);

 

Therefore the integral can be expressed as follows. 

>	Int(A/B, x) = c1*log(v1) + c2*log(v2);

 

>

 

In this case, the result is identical to that obtained by the "classical approach" where the denominator is completely factored. 

Example 1.5 

One of the integrals from Example 1.2 is 

>	A := 1:  B := x^3+x:

 

>	Int(A/B, x);

 

>	R := resultant(A - z*diff(B,x), B, x);

 

>	solve(R=0, z);

 

The distinct roots are 

>	(c1,c2) := (1,-1/2);

 

>	v1 := gcd(A - c1*diff(B,x), B);

 

>	v2 := gcd(A - c2*diff(B,x), B);

 

Therefore the integral can be expressed as follows. 

>	Int(A/B, x) = c1*log(v1) + c2*log(v2);

 

In contrast, the "classical" approach would require a complete factorization of the denominator 

>	factored_B := factor(B,I);

 

>	integrand := convert(A/factored_B, parfrac, x);

 

from which we see that the integral is 

>	Int(A/B, x) = map(`int/risch`, integrand, x);

 

>

 

The Rothstein-Trager method avoids introducing the algebraic number extension    in this example. 

Computational cost increases significantly with each new algebraic number extension, so it is important to use the minimal number of extensions required by the problem. 

Example 1.6 

For the integral 

>

restart;

 

>	Int(1/ln(x), x);

 

the integrand lies in the function field where is a logarithmic extension of the rational function field . 

In this field the integrand is 

>

1/theta;

 

>

 

Example 1.7 

For the integral 

>	Int((x^3+2*x)*cos(x^2)+x*sin(x^2), x);

 

by converting the integrand to complex - form, the problem can be expressed in the form 

>	Int((x^3+2*x)*(exp(I*x^2)+exp(-I*x^2))/2 - 1/2*I*x*(exp(I*x^2)-exp(-I*x^2)), x);

 

and we can view the integrand as an element of the function field    where .  (Note: .) 

In this field the integrand is 

>	1/2*(x^3+2*x)*(theta + 1/theta) - 1/2*I*x*(theta - 1/theta);

 

>

 

Example 1.8 

For the integral presented in Example 1.1 

>	Int(x*(x+1)*( (x^2*exp(2*x^2) - ln(x+1)^2)^2 + 2*x*exp(3*x^2)*( x - (2*x^3+2*x^2+x+1)*ln(x+1) )) / ((x+1)*ln(x+1)^2 - (x^3+x^2)*exp(2*x^2) )^2, x);

 

we can view the integrand as an element of the function field    where and . 

Note that the two exponential terms    and    cannot be considered as independent extensions because there is an algebraic relationship between them. Rather, we choose to represent these terms in the form    and    which involves a single transcendental extension   . 

>

 

Hermite Reductions for the Rational Part (log extension) 

This case mimics the Hermite method for rational functions.  Applying square-free factorization of the denominator yields    where this means "square-free as a polynomial in ."  For application in subsequent steps of the algorithm, we need to know that the following theorem holds. 

Theorem 3:  If is a monic polynomial which is square-free as a polynomial in the symbol and if    for some function then 

where the latter differentiation is with respect to . 

Then applying partial fraction expansion, we get 

   where     . 

For each integral with a denominator containing a power   > , apply the Extended Euclidean Algorithm to solve 

                     . 

As in rational function integration, apply integration by parts to obtain the Hermite reduction: 

. 

Continuing such reductions until no remaining integral has a denominator containing a power   > , we reduce the integration problem to one in which the denominator is square-free. 

Rothstein-Trager Theorem: log extension 

The problem is to express    where   ,    is square-free, and where    is a logarithmic extension. 

Theorem 4  (Rothstein-Trager Theorem: log extension) 

(i)    is elementary if and only if all of the roots of 

are constants.  In other words, the primitive part of with respect to the variable must have constant coefficients. 

(ii)  If    is elementary then 

where are the distinct roots of and where are defined by 

          . 

Moreover, this expresses the integral using the minimal algebraic extension field. 

Example 1.9 

Consider the following integral. 

>	Int(1/ln(x), x);

 

The integrand is    in the function field where . 

A resultant computation (Rothstein-Trager Theorem) is applicable. 

>	A := 1:  B := theta:

 

>	diff_B := diff(ln(x), x);

 

>	R := resultant(A - z*diff_B, B, theta);

 

The single root of is which is not a constant. 

Therefore the original integral is not elementary. 

>

 

Example 1.10 

The integral 

>	Int(1/(x*ln(x)), x);

 

has integrand    in the function field where . 

This time, the resultant computation is as follows. 

>	A := 1:  B := x*theta:

 

>	diff_B := subs(ln(x)=theta, diff(subs(theta=ln(x), B), x));

 

>	R := resultant(A - z*diff_B, B, theta);

 

The single root of is and the integral is elementary.  Specifically, 

>	c[1] := 1;  v[1] := gcd(A - c[1]*diff_B, B);

 

and the integral is 

>	c[1]*ln(v[1]);

 

In other words, 

>	Int(1/(x*ln(x)), x)  = subs(theta=ln(x), c[1]*ln(v[1]));

 

>

 

Polynomial Part (log extension) 

For the integration of the "polynomial part", the integrand is a polynomial in the log extension : 

with    in    and we can conclude that 

. 

Note that the last term indicates that some new log extensions may arise. 

Differentiating and equating coefficients of powers of the transcendental symbol , we get: 

                 . . . 

                 . . . 

where   . 

We can solve these equations successively, from the top down.  At each step there is an integration operation which requires a recursive invocation of the Risch algorithm.  This works because the integrand lies in the field    involving one less extension.  At any step, if the integral to be computed is not elementary then the algorithm halts with the conclusion that the original integral is not elementary.  Moreover, at any step except the last step, if the result of the recursive integration introduces one or more new log extensions then the algorithm halts with the conclusion that the original integral is not elementary.  At the last step, new log extensions may appear. 

Example 1.11 

The integral 

has integrand    in the field    where   .  If the integral is elementary then 

where the equations to be satisfied are 

            . 

From the first equation we conclude that   , an arbitrary constant of integration.  Plugging this into the second equation and integrating both sides yields 

or 

where is an arbitrary constant of integration. Equating coefficients of on both sides of the equation, we conclude that and . 

The last equation now becomes 

or rearranging (we always leave the term    on the right hand side since integration of this term is trivial): 

            . 

At this point we use the fact that    (i.e. )  on the left hand side, and integrate both sides to get 

from which we conclude (by equating coefficients of on both sides of the equation) that and , ignoring the arbitrary constant of integration in this final step. Putting it all together we have 

>	theta := ln(x):

 

>	B[2] := 0:  B[1] := x:  B[0] := -x:

 

>	Int(theta, x) = B[2]*theta^2 + B[1]*theta + B[0];

 

>

 

Example 1.12 

The integral 

has integrand    in the field    where and .  If the integral is elementary then 

where the equations to be satisfied are 

            . 

From the first equation we conclude that   , an arbitrary constant of integration.  Plugging this into the second equation and integrating both sides yields 

or 

where is an arbitrary constant of integration. Equating coefficients of on both sides of the equation, we conclude that and . The last equation now becomes 

or rearranging: 

            . 

Since   ,  i.e. , substituting for on the left hand side and integrating yields 

            . 

From Example 1.9 we know that the integral appearing here is not elementary. Hence we conclude that    is not elementary. 

Hermite Reductions for the Rational Part (exp extension) 

This case is very similar to the logarithmic case.  Applying square-free factorization of the denominator yields    where this means "square-free as a polynomial in ."  For application in subsequent steps of the algorithm, we need to know that the following theorem holds. 

Theorem 5:  If is a monic polynomial which is square-free as a polynomial in the symbol and if    for some function then 

where the latter differentiation is with respect to ,  provided that    does not divide . 

Then, exactly as in the logarithmic case, applying partial fraction expansion we get 

   where     . 

For each integral with a denominator containing a power   > , apply the Extended Euclidean Algorithm to solve 

                     . 

As before, apply integration by parts to obtain the Hermite reduction: 

. 

Continuing such reductions until no remaining integral has a denominator containing a power   > , we reduce the integration problem to one in which the denominator is square-free. 

Rothstein-Trager Theorem: exp extension 

The problem is to express    where   ,    is square-free, and where    is an exponential extension. 

Theorem 6  (Rothstein-Trager Theorem: exp extension) 

(i)    is elementary if and only if all of the roots of 

are constants.  In other words, the primitive part of with respect to the variable must have constant coefficients. 

(ii)  If    is elementary then 

where are the distinct roots of ,   are defined by 

and where   . 

Moreover, this expresses the integral using the minimal algebraic extension field. 

Example 1.13 

Consider the following integral. 

>

restart;

 

>	f := (x*exp(2*x^2)+8*sqrt(2)*x*(exp(x^2)+1)+2*x)/     (exp(3*x^2)+exp(2*x^2)+2*exp(x^2)+2):

 

>

Int(f,x);

 

The integrand is    in the function field where . 

A resultant computation (Rothstein-Trager Theorem) is applicable. 

Keep in mind that    where   . 

>

u := x^2:

 

>	A := x*theta^2 + 8*sqrt(2)*x*(theta+1) + 2*x:

 

>	B := theta^3 + theta^2 + 2*theta + 2:

 

>	diff_B := subs(exp(u)=theta, diff(subs(theta=exp(u), B), x));

 

>	R := resultant(A - z*diff_B, B, theta);

 

>	solve(R=0, z);

 

The distinct roots of are    and    and the integral is elementary.  Specifically, 

>	c[1] := -1/2:  c[2] := -sqrt(2):

 

>	v[1] := gcd(A - c[1]*diff_B, B);

 

>	v[2] := gcd(A - c[2]*diff_B, B);

 

The function appearing in Theorem 6 is 

>	g := add(-c[i]*degree(v[i],theta), i=1..2) * u;

 

and the desired integral is 

>	result := g + add(c[i]*ln(v[i]), i=1..2);

 

In other words, 

>	Int(f,x) = subs(theta=exp(u), result);

 

>

 

Polynomial Part (exp extension) 

For the integration of the "polynomial part", the integrand is of the form    with    in   , where    is an exponential extension.  We can conclude that 

. 

Note that the latter summation indicates that some new log extensions may arise. 

Differentiating and equating coefficients of powers of the transcendental symbol , we get: 

           ,   for  

where   . 

Solving the latter equation for    requires a recursive integration:    where the integrand lies in the field   .  For the general case,   , the equation to be solved for    is known as a Risch differential equation.  This may sound like a backwards step to have reduced an integration problem to the solution of a differential equation!  However, the problem can be solved because we must look for a solution in the field    only.  Nonetheless, this is a difficult problem, in general, for which a complete solution was presented in [Bronstein90a].  See also [Bronstein97]. 

Example 1.14 

The integral 

>

restart;

 

>	Int(exp(-x^2), x);

 

has integrand    in the function field where . 

The solution must take the form    where in is a solution of the Risch differential equation   . 

Since must be a rational function in , one first argues that it must be a polynomial in . This follows because if   had a nontrivial denominator then  would have a higher-degree denominator and these denominators could not cancel out to correspond to the right hand side of the differential equation which is the constant 1 .  So assume that    is a nonzero polynomial in with    where   .  But then taking the degree of each side of the differential equation, we get    which is a contradiction.  Hence the given Risch differential equation has no solution in the field    from which we conclude that the original integral is not elementary. 

>

 

Example 2.1 

>

restart;

 

>	f := ln(x):

 

We wish to compute the definite integral 

>	Int(f, x=0..2);

 

We can compute the antiderivative 

>	F := int(f, x);

 

but direct evaluation of the antiderivative at encounters a logarithmic singularity. However, by taking limits we can compute the result as follows. 

>	limit(F, x=2, left) - limit(F, x=0, right);

 

Of course, the integration procedure in Maple will do this computation automatically. 

>	int(f, x=0..2);

 

>

 

Example 2.2 

>	g := sin(x)/x;

 

>	series(g, x=0);

 

From the series expansion it is clear what is the limit at . 

>	Limit(g, x=0) = limit(g, x=0);

 

Another example: 

>

g := x^x;

 

>	series(g, x=0, 3);

 

>	Limit(g, x=0) = limit(g, x=0);

 

And another example: 

>	g := ln(1-cos(x^2))/(ln(x)*arctan(sqrt(1-x^2)));

 

>	series(g, x=0, 4);

 

>	Limit(g, x=0) = limit(g, x=0);

 

>

 

Example 2.3 

The following example is discussed in [Bronstein97]. Suppose that the integral to be computed is 

>	f := (x^4-3*x^2+6)/(x^6-5*x^4+5*x^2+4):

 

>	Int(f, x=1..2);

 

It can be verified that the integrand is continuous and positive on the real line, and hence the value of the definite integral must be a positive real number. Applying the integration algorithm in the form discussed in this presentation, the indefinite integral gets expressed in the following form. 

>	F := I/2*ln(x^3+I*x^2-3*x-2*I) - I/2*ln(x^3-I*x^2-3*x+2*I);

 

That is, 

>	Int(f, x) = F;

 

It can easily be verified that the derivative of    is the integrand    and therefore this is a correct formal antiderivative. However, the differential algebra point of view used in developing the integration algorithm takes no account of the analytic concept of branch cuts of logarithmic functions. In this particular example, if we blindly apply the Fundamental Theorem of Calculus we get the following incorrect result for the definite integral. 

>	Int(f, x=1..2) = limit(F, x=2, left) - limit(F, x=1, right);

 

>	evalf(rhs(%));

 

This has yielded a value for the integral which is real, but negative, whereas we know that the correct answer must be positive. The presentation in [Bronstein97] shows how to ensure that the integration algorithm computes a continuous antiderivative for this type of problem. In this case, a continuous antiderivative is 

>	F_continuous := arctan((x^5 - 3*x^3 + x)/2) + arctan(x^3) + arctan(x);

 

By differentiating, we see that this is a correct antiderivative of the original integrand . 

>	simplify( diff(F_continuous, x) );

 

This time, direct application of the Fundamental Theorem of Calculus is valid and we obtain the following result for the definite integral. 

>	Int(f, x=1..2) = limit(F_continuous, x=2, left) - limit(F_continuous, x=1, right);

 

>	evalf(rhs(%));

 

>

 

Example 2.4 

In Example 2.3 we had the following discontinuous antiderivative. 

>	Int(f, x) = F;

 

In Maple, the following command will determine the points of discontinuity over the whole real line. 

>	discont(F, x);

 

>	allvalues(%);

 

The above notation is expressing the fact that there are two points of discontinuity which are the two roots of a quadratic polynomial. Since we are interested in the definite integral over the interval there is just one relevant point of discontinuity: 

>	c := sqrt(2);

 

and the endpoints of integration are 

>	a := 1;  b := 2;

 

The correct value of the definite integral on can therefore be computed by using the fact that 

and noting that the antiderivative    is continuous on each of the two subintervals. 

This leads to the following limit computations for the definite integral. 

>	Int(f, x=a..b) = limit(F, x=c, left) - limit(F, x=a, right) +                 limit(F, x=b, left) - limit(F, x=c, right);

 

>	evalf(rhs(%));

 

>

 

Example 2.5 

>

restart;

 

>	assume(c::real);  interface(showassumed=0):

 

>

 

>	f1 := x^(s-1)*exp(-p*x^4)*erfi(c*x)*erf(c*x):

 

>	Int(f1, x=0..infinity);

 

>	result1 := value(%);

 

The method used was first to convert the integrand into form. 

>	f1new := convert(f1, MeijerG, x);

 

Now when we apply the command, it uses a known formula to express such a definite integral involving the product of two functions. 

>	int(f1new, x=0..infinity);

 

which is equal to . 

>	simplify(% - result1);

 

Of course, for particular values of the parameters we may compute a numerical value. 

>	eval(result1, {s=4, p=1, c=1});

 

>

evalf(%);

 

>

 

Example 2.6 

>

restart;

 

>	assume(b>0, y>0);  interface(showassumed=0):

 

>

 

>	f2 := x^v/(x^2+y^2)*BesselK(v,b*x):

 

>	Int(f2, x=0..infinity);

 

>	result2 := value(%);

 

Note that the result of converting the integrand into form is as follows. 

>	f2new := convert(f2, MeijerG, x);

 

>	int(f2new, x=0..infinity);

 

which is equal to . 

>	simplify(% - result2);

 

>