Maple Integration Error

The following give very different answers.

> evalf( int( ln(2+2*cos(t)^2+cos(t)^4), t = 0 .. Pi));

   -8.710344354+19.73920881*I;

> evalf( Int( ln(2+2*cos(t)^2+cos(t)^4), t = 0 .. Pi));

   3.658972529

Is this part of a larger problem?

 

Doug Meade's picture

int vs Int and the perils of unanalyzed antiderivatives

This type of behavior is rather common, and is NOT a bug. It simply shows the differences between symbolic calculations, without careful consideration of the problem, and numeric calculations.

Let's consider the two results separately. First,

evalf( int( ln(2+2*cos(t)^2+cos(t)^4), t = 0 .. Pi));
                        -8.710344354 + 19.73920881 I

The int command produces an exact expression for the antiderivative:

F := int( ln(2+2*cos(t)^2+cos(t)^4), t);
                                          /
                                          |
4 I ln(2) ln(exp(I t)) - I ln(exp(I t)) ln|
                                          |
                                          \

               4                2                6                 8\
  54 (exp(I t))  + 12 (exp(I t))  + 12 (exp(I t))  + 1 + (exp(I t)) |
  ------------------------------------------------------------------|
                                       4                            |
                             (exp(I t))                             /

                             /                     -----                      
                             |                      \                         
                     2       |                       )                        
   - 2 I ln(exp(I t))  + 8 I |                      /                         
                             |                     -----                      
                             |            /     4        2        6         8\
                             \_R1 = RootOf\54 _Z  + 12 _Z  + 12 _Z  + 1 + _Z /

                                                               \
                                                               |
  /1                /_R1 - exp(I t)\   1      /_R1 - exp(I t)\\|
  |- ln(exp(I t)) ln|--------------| + - dilog|--------------|||
  \8                \     _R1      /   8      \     _R1      //|
                                                               |
                                                               /

This is not a very useful result. But, if we continue, evaluating this antiderivative at each endpoint and subtracting gives:

F1 := simplify(eval( F, t=Pi ));
F2 := simplify(eval( F, t=0  ));
                       /Sum(/       /_R1 + 1\        /_R1 + 1\\,             /
                       |    |I Pi ln|-------| + dilog|-------||  _R1 = RootOf\
                 2     |    \       \  _R1  /        \  _R1  //               
Pi ln(5) + 2 I Pi  + I |                                                      
                       |                                                      
                       |                                                      
                       \                                                      

       4        2        6         8\ .. )\
  54 _Z  + 12 _Z  + 12 _Z  + 1 + _Z /     |
                                          |
                                          |
                                          |
                                          |
                                          /
       /                     -----                                    \
       |                      \                                       |
       |                       )                             /_R1 - 1\|
     I |                      /                         dilog|-------||
       |                     -----                           \  _R1  /|
       |            /     4        2        6         8\              |
       \_R1 = RootOf\54 _Z  + 12 _Z  + 12 _Z  + 1 + _Z /              /
evalf( F1-F2 );
                        -8.710344359 + 19.73920881 I

The complex valued result is not expected. It is obtained because we have not give careful consideration to the branch cuts of the antiderivative.

On the other hand, the second option:

evalf( Int( ln(2+2*cos(t)^2+cos(t)^4), t = 0 .. Pi));
                                 3.658972529

does not try to do any symbolic evaluation, it immediately resorts to a numerical integration routine. Here, the integrand is always positive, so the real-valued and positive result is consistent with our expectations.

It is cases like this that make the distincition between int and Int important. With evalf( int( ) ), Maple will first attempt a symbolic evaluation of the definite integral. If this is not successful, then it will use numerical quadrature. But, if you want to avoid the time spent trying to find an antiderivative (or, worse, finding a useless antiderivative as in this example) use evalf( Int( ) ).

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/
JacquesC's picture

Still a bug though

If Maple insists on using FTOC and discontinuous anti-derivatives for definite integration, then it has to be able to figure out all the points of discontinuity in the domain, as well as properly figuring out all the limits involved (which is very very hard in this particular case).  My guess is that sufficient digging into this particular issue would reveal an actual bug in either discont or limit.  I have fixed many bugs in int in my time that were just like this one, and pretty much all of them reduced to a discont or limit bug.  As I have stated before, it is now my belief that although FTOC is a wonderful theorem of "theoretical calculus", its actually computational usefulness is much more limited [a statement with which you seem to agree in this particular instance].

In case of doubt plot

F:=int( ln(2+2*cos(t)^2+cos(t)^4), t );
plot(Re(F),t = 0 .. Pi);

If you see discontinuities like here do not use blindly the definite integral. If there were a Maple wiki, this should be in the FAQ.

Axel Vogt's picture

it is a bug i would say

  Int( ln(2+2*cos(t)^2+cos(t)^4), t = 0 .. Pi);
  S:=simplify(combine(value(%)));

       / -----
       |  \
  S := |   )
       |  /
       | -----
       \_R1 = %1

                                                             \
        /       _R1 - 1          _R1 + 1            _R1 + 1 \|
        |-dilog(-------) + Pi ln(-------) I + dilog(-------)|| I
        \         _R1              _R1                _R1   /|
                                                             |
                                                             /

                            2
         + Pi ln(5) + 2 I Pi

                    4        2        6         8
  %1 := RootOf(54 _Z  + 12 _Z  + 12 _Z  + 1 + _Z )

  S1:=select(has,S,sum);
  S2:=S-S1;

Now evaluating S1 ( = the sum from above) shows the error, since the integral is real.

that is another side of the bug

I would say. A plot:

F:=int( ln(2+2*cos(t)^2+cos(t)^4), t );
plot([Re(F),Im(F)],t = 0 .. Pi);

shows that this imaginary part is the imaginary part of a piecewise "integration" constant.

Axel Vogt's picture

ok

the jumps where to cut and take the directed limits for the anti-derivate

1/2*Pi-1/2*arccos(5^(1/2)-2^(1/2)) and 1/2*Pi+1/2*arccos(5^(1/2)-2^(1/2))
(i used 'discont', stared at the solution and with it, searching a Real one)

however M11 needs evalf to get the limits

seems MMA does not give the anti-derivative (online version)

PS: the new editor is a mess for line breaks, deleting and back space
only usable if typing first in a usual editor, then do copy+paste

Maple Integration Error

From the Wikipedia: "A software bug (or just "bug") is an error, flaw, mistake, failure, or fault in a computer program that prevents it from behaving as intended (e.g., producing an incorrect result)." When Maple produces an incorrect result, it's a bug. 

As JacquesC says, the error arises because Maple is applying the Fundamental Theorem of Calculus incorrectly. The function that the Int command produces is not an antiderivative of our function on the entire interval [0, Pi] because it is discontinuous on the interval. Therefore, the FTC doesn't apply.

Just because Maple has been making this error for as long as we can all remember, doesn't excuse it. To tell the user not to blindly use the definite integral function is good advice since Maple refuses to acknowledge that there is a problem. To suggest that the user is at fault for not properly considering the branch cuts lays the bug on the wrong doorstep. How hard would it be for Maple to check that the answer from the definite integral is (essentially) the same as the numerical result? It's far better to give no result than an incorrect one. Or at least a warning.

I have no problem with the indefinite integral because Maple makes it clear for the beginning that the it will not bother with the constants of integration.

I guess most of my frustration arises from trying to convince students that learning Maple is worth the effort . And then every term Maple undercuts me by making an elementary error like this.

Axel Vogt's picture

complains and bugs

I understand your complains. But would not agree in all your points.

It is certainly not true that the system ignores the assumptions of
the fundamental theorem. More likely it is a bug (may even be some
coding bug) and not a fundamental fault (i.e. ignoring).

In some cases Maple returns unevaluated, because it can not decide
whether the FT can be properly applied (have no example at hand now).

The other thing: automatical numerical cross checks. From time to
time this also is brought up at sci.math.symbolic - mainly there are
2 reasons why not to do: it will slow down the system (may be needed
in many intermediate steps) and in certain cases it will be difficult
to judge (say numerical almost impossible to achieve results or the
contrary (there is an article by Fateman(?) concerning the latter)).

The integral itself is far from being easy (remind that MMA does not
give an indefinite answer) and may be seen as ln((1+p(x)^2))/sqrt(1-x^2),
p(x)=1+x^2, to be integrated over the unit interval, by x=cos(t).

For other troublesome tasks you may for example search for Zacharias
"Some intriguing definite integrals" (1999).

So your frustration is understandable, but you picked up a peculiar
problem, I would not look at a CAS as a universal solver ... and
hope that encourages you again, at least a bit.

http://www.wias-berlin.de/main/publications/wias-publ/run.cgi?template=abstract&type=Preprint&year=1999&number=493

JacquesC's picture

I disagree

I have read the code as well as the book on which the code is based.  They both definitely ignore some of the side-conditions in FTOC!  This is a well-known flaw of basing definite integration on the differential algebra view of anti-derivatives which has been amply documented in the literature.  There are current research projects ongoing which seek to remedy this flaw in the theory, but have only provided partial answers.  The current integrator tries to back-patch this issue by reverse-engineering the new discontinuities introduced in the integral and computing limits to patch over the jumps, but that is often a much harder problem than the original integral was -- and so predictably one encounters bugs.  Note that the theory underlying the computation of limits has holes in it which as just as large as those in integration [the problems are that the theory relies on a theory of series which is too weak, as well as the fact that the theory requires a decision procedure for zero recognition, which is a well-known undecidable problem].

two bugs in a rug

Axel says that Mathematica does not evaluate the indefinite integral, but in fact it does. It is the mess below. Note it does not include the Sum. However, Mathematica's antiderivative is also discontinuous, and when used to evaluate the definite integral, yields the same incorrect answer

-8.710344354 + 19.73920881 I

f:=x->x*ln(2 + 2*cos(x)^2 + cos(x)^4) + 2*(-(x*ln((3 - 2*I) + cos(2*x)))/2 + ((2*I)*x^2 - 4*arcsin(sqrt(2 - I))* arctanh(((1 + I)*tan(x))/ sqrt(1 - 3*I)) - (2*x + 2*arcsin(sqrt(2 - I)))* ln(1 + ((3 - 2*I) - 2*sqrt(1 - 3*I))* exp((2*I)*x)) - (2*x - 2*arcsin(sqrt(2 - I)))* ln(1 + ((3 - 2*I) + 2*sqrt(1 - 3*I))* exp((2*I)*x)) + 2*x*ln((3 - 2*I) + cos(2*x)) + I*(polylog(2, -(((3 - 2*I) - 2*sqrt(1 - 3*I))* exp((2*I)*x))) + polylog(2, -(((3 - 2*I) + 2*sqrt(1 - 3*I))* exp((2*I)*x)))))/4) + 2*(-(x*ln((3 + 2*I) + cos(2*x)))/2 + ((2*I)*x^2 - (4*I)*arcsin( sqrt(2 + I))*arctan( ((1 + I)*tan(x))/ sqrt(1 + 3*I)) - (2*x + 2*arcsin(sqrt(2 + I)))* ln(1 + ((3 + 2*I) - 2*sqrt(1 + 3*I))* exp((2*I)*x)) - (2*x - 2*arcsin(sqrt(2 + I)))* ln(1 + ((3 + 2*I) + 2*sqrt(1 + 3*I))* exp((2*I)*x)) + 2*x*ln((3 + 2*I) + cos(2*x)) + I*(polylog(2, -(((3 + 2*I) - 2*sqrt(1 + 3*I))* exp((2*I)*x))) + polylog(2, -(((3 + 2*I) + 2*sqrt(1 + 3*I))* exp((2*I)*x)))))/4);

 

Axel Vogt's picture

perils ...

(Alex: I used the online version and got nothing, but thx for the result!)

I have no reason to doubt Jacques' answer ... but played a bit more,
may be the following may help (but does not solve the problem).

Problems already occure for Int(1/(cos(t)-(-1-I)^(1/2))*t*sin(t),t = 0 .. Pi)
or -Int(cos(t)*t*sin(t)/(-cos(t)-(1+I)^(1/2)),t = 0 .. Pi).
Axel Vogt's picture

the exact value

I think the eaxt value is

  Pi*ln(1/16*(4+2*10^(1/2))*(3+10^(1/2))^(1/2)+13/16+1/4*10^(1/2));
  evalf(%);


                          1/2         1/2 1/2          1/2
                 (4 + 2 10   ) (3 + 10   )      13   10
           Pi ln(---------------------------- + -- + -----)
                              16                16     4

                           3.6589725292852

and I should find the time to write it up (with continuous anti-derivatives).
Axel Vogt's picture

my solution

Within Int( ln(2+2*cos(t)^2+cos(t)^4), t = 0 .. Pi) we have a (trigonometric)
polynomial of degree 4, which splits over the Reals:

  r1,r2:=(-1+I)^(1/2), -(-1-I)^(1/2);
  p1,p2:= 'eval((x-r)*(x-conjugate(r)), r= r1)', 'eval((x-r)*(x-conjugate(r)), r= r2)';
  'p1*p2': '%'=evalc(collect(expand(%),x));

                                   1/2           1/2
                 r1, r2 := (-1 + I)   , -(-1 - I)


                              _ |                     _ |
       p1, p2 := (x - r) (x - r)|      , (x - r) (x - r)|
                                |r = r1                 |r = r2


                                     4      2
                        p1 p2 = 2 + x  + 2 x


The parabolas p have real coefficients and no roots on the interval, so log
turns the product to a sum, thus we have 2 integrals of the following form:

  Int(ln((cos(t)-conjugate(r))*(cos(t)-r)),t=0..Pi);
  ``=PDEtools[dchange](t=arccos(x),%,[x]);
  ``=value(rhs(%)): # <------- in terms of hypergeometrics
  ``=simplify(rhs(%),size);
  
  J:=rhs(%):

                   Pi
                  /
                 |                 _
                 |    ln((cos(t) - r) (cos(t) - r)) dt
                 |
                /
                  0


                          1
                         /                  _
                        |   ln((x - r) (x - r))
                     =  |   ------------------- dx
                        |             2 1/2
                       /        (1 - x )
                         -1


         /                           /     1  \1/2
         |                           |1 - ----|
         |                           |      2 |
         |             _             \     r  /
    = Pi |ln(-r) + ln(-r) + ln(1/2 + -------------)
         |                                 2
         \

                                      \
                                      |
                                      |
                        /     1  \1/2 |
         + ln(1/2 + 1/2 |1 - ----|   )|
                        |     _2 |    |
                        \     r  /    /


Hence the value is given by

  'eval(J,r=r1) + eval(J,r=r2)';
  ``=evalc(%):
  combine(%,ln);
  evalf(%);

                         J|       + J|
                          |r = r1    |r = r2


                     1/2          1/2 1/2   1/2          1/2 1/2
              13   10      (3 + 10   )    10      (3 + 10   )
      = Pi ln(-- + ----- + -------------------- + --------------)
              16     4              8                   4


                             = 3.658972529


Numerical cross check:

  Int( ln(2+2*cos(t)^2+cos(t)^4), t = 0 .. Pi); evalf(%);

                    Pi
                   /
                  |                   2         4
                  |    ln(2 + 2 cos(t)  + cos(t) ) dt
                  |
                 /
                   0


                             3.658972529


I was using Maple 11.02 with Digits:=14

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
}