IntBug.mws

Integrals in Maple: the so-called branch-cut problem

Sometimes Maple gives integration results that have to be taken carefully. In particular, anti-derivatives sometimes can be treacherous in that they are valid only in a limited parameter range.

We illustrate the problem on a simple example which may help users to identify and hopefully prevent difficult situations. We can control with infolevel[int]  how much informations is displayed about the progress of the integral evaluations.

>    restart; infolevel[int]:=0;

infolevel[int] := 0

Our integrand contains two parameters: eta  is real, and 0 < e < 1 . We were interested in the answer of the definite integral from 0  to 2*Pi . Maple7 and Maple8 gave correct answers for the definite integral, but was giving a wrong answer of 0  when the antiderivative F  was obtained first, and evaluated for F(2*Pi)-F(0) . Maple9 still gives a wrong answer for the latter, but returns unevaluated for the definite integral.

>    f:=sin(u+eta)*cos(u+eta)/(1+e*cos(u))^4;

f := sin(u+eta)*cos(u+eta)/(1+e*cos(u))^4

>    int(f,u=0..2*Pi) assuming(eta>0,e>0,e<1);

int(sin(u+eta)*cos(u+eta)/(1+e*cos(u))^4,u = 0 .. 2*Pi)

Let us illustrate the problem for a numerical choice of parameters:

>    e:=1/2;

e := 1/2

>    eta:=1;

eta := 1

>    f;

sin(u+1)*cos(u+1)/(1+1/2*cos(u))^4

Definite integration does work now:

>    DI:=int(f,u=0..2*Pi);

DI := 160/81*sin(1)*cos(1)*3^(1/2)*Pi

>    evalf(DI);

4.886764727

We can verify the answer numerically:

>    evalf(Int(f,u=0..2*Pi));

4.886764725

For the general parameter case Maple7 or Maple8 gets:

-5*cos(eta)*sin(eta)*e^2*Pi/((e+1)^3)/((-1+e)^3)/((-e^2+1)^(1/2))

Maple9.03 has no answer in this case.


Try the antiderivative:

>    F:=unapply(int(f,u),u);

F := proc (u) options operator, arrow; 8*cos(1)^2/(2+cos(u))^2-32/3*cos(1)^2/(2+cos(u))^3-704/27*cos(1)*sin(1)*tan(1/2*u)/(tan(1/2*u)^2+3)+160/81*cos(1)*sin(1)*3^(1/2)*arctan(1/3*tan(1/2*u)*3^(1/2))-89...
F := proc (u) options operator, arrow; 8*cos(1)^2/(2+cos(u))^2-32/3*cos(1)^2/(2+cos(u))^3-704/27*cos(1)*sin(1)*tan(1/2*u)/(tan(1/2*u)^2+3)+160/81*cos(1)*sin(1)*3^(1/2)*arctan(1/3*tan(1/2*u)*3^(1/2))-89...

>    B:=F(2*Pi);

B := 40/81*cos(1)^2-40/81*sin(1)^2

>    A:=F(0);

A := 40/81*cos(1)^2-40/81*sin(1)^2

>    B-A;

0

This answer is definitely wrong. The error is most likely caused by the arctan(tan)  part. We verify that there is a problem by graphing the anti-derivative:

>    plot(F(x),x=0..2*Pi,discont=true);

[Maple Plot]

It is a good idea to check the antiderivative for discontinuities in the range defined by the limits of integration. Maple's answer for the antiderivative implicitly made a branch cut choice during the calculation that led to a discontinuous anti-derivative.

The correct answer can only be obtained from an anti-derivative that is continuous over the entire interval.

We might think that the following test would give us a pointer: compare the derivative of the anti-derivative with the original function:

>    Fp:=diff(F(u),u);

Fp := 16*cos(1)^2/(2+cos(u))^3*sin(u)-32*cos(1)^2/(2+cos(u))^4*sin(u)-704/27*cos(1)*sin(1)*(1/2+1/2*tan(1/2*u)^2)/(tan(1/2*u)^2+3)+1408/27*cos(1)*sin(1)*tan(1/2*u)^2/(tan(1/2*u)^2+3)^2*(1/2+1/2*tan(1/2...
Fp := 16*cos(1)^2/(2+cos(u))^3*sin(u)-32*cos(1)^2/(2+cos(u))^4*sin(u)-704/27*cos(1)*sin(1)*(1/2+1/2*tan(1/2*u)^2)/(tan(1/2*u)^2+3)+1408/27*cos(1)*sin(1)*tan(1/2*u)^2/(tan(1/2*u)^2+3)^2*(1/2+1/2*tan(1/2...
Fp := 16*cos(1)^2/(2+cos(u))^3*sin(u)-32*cos(1)^2/(2+cos(u))^4*sin(u)-704/27*cos(1)*sin(1)*(1/2+1/2*tan(1/2*u)^2)/(tan(1/2*u)^2+3)+1408/27*cos(1)*sin(1)*tan(1/2*u)^2/(tan(1/2*u)^2+3)^2*(1/2+1/2*tan(1/2...

This looks different from f . However, with some massaging we find that F  is a legitimate anti-derivative to f :

>    simplify(expand(f-Fp));

0

The fact that our anti-derivative as calculated by Maple has a jump discontinuity is not picked up by the test to compared the derivative of the indefinite integral with the integrand!

We can calculate the correct answer by carefully tip-toeing around the discontinuity in F  at x=Pi :

>    Ra:=limit(F(x),x=Pi,left)-F(0);

Ra := 80/81*cos(1)*sin(1)*3^(1/2)*Pi-472/81*cos(1)^2+8/3+40/81*sin(1)^2

>    Rb:=-limit(F(x),x=Pi,right)+F(2*Pi);

Rb := 80/81*cos(1)*sin(1)*3^(1/2)*Pi-8/3+472/81*cos(1)^2-40/81*sin(1)^2

>    Ra+Rb,DI;

160/81*cos(1)*sin(1)*3^(1/2)*Pi, 160/81*cos(1)*sin(1)*3^(1/2)*Pi

This verifies that once one takes care of the jump discontinuity everything is fine.

It raises the question: how does Maple's definite integration package avoid the problem? It is probably true that Mathematica is better in these matters by making use of other representations of the anti-derivative. When encountering inverse trig functions in answers from int : beware! These problems with Maple's answers have been around for many years. Sometimes Maple yields complex-valued answers when integrating real functions over the real axis. Checking symbolic answers by numerical evaluation for fixed parameter choices provides some assurance that results are reasonable. It would be useful if the Maple development team took these issues more seriously. It seems that they act like mathematicians who are satisfied when it is shown that an answer can be obtained. The least we should expect is a warning that the calculated anti-derivative has potential problems with discontinuities in particular locations.

This particular integration problem was encountered by Michael Horbatsch.

>