# Bug in Integrate for Mathematica

Consider the following:

regFunc[x_,y_]:=Boole[-5<4x+3y<5 && -2<3x+2y<2];
Integrate[regFunc[x,y]*((4x+3y)(3x+2y))^4,{x,-100,100},{y,-100,100}]//N
Integrate[regFunc[x,y]*(12x^2+17x*y+6y^2)^4,{x,-100,100},{y,-100,100}]//N


In the first example, Mathematica seems to figure out the correct substitution, and arrives at 16000, which is the correct answer. In the second integration, (integrating the same expression, but expanded), Mathematica gives 5885078144/382725, which is 15376.8. This is clearly a bug, is this well-known?

I am using Mathematica 8.04, 64bit Linux.

EDIT: I got an email from the tech support, and I took the answer as a confirmed bug.

-
did you really mean xy in the second example? (guess not) –  acl Mar 27 '12 at 12:39
Yes, I know that NIntegrate manages to get the answer correct. What I am curious about, is why the two identical integrals give different answers, depending on if the integrand is factorized or not. –  Per Alexandersson Mar 27 '12 at 12:53
Since this is clearly a bug, the proper place to send it is support@wolfram.com. –  Mark McClure Mar 27 '12 at 13:22
Blah, I assumed it meant confirmed by other users @sjoerd –  nixeagle Mar 27 '12 at 15:49
In Mathematica 9.0.0, I get the same result 16000 from both forms. –  murray Dec 13 '12 at 16:38

## 2 Answers

I think it is indeed a bug specific to version 8 of Mathematica. The same integrals in version 7 give the correct result. Compare this issue with this answer. In the both cases one works with assumptions which make Integrate behaving improperly.

Edit 1

It seems that definite integrals are calculated correctly and if we subtract the limits of integration in the way that the boolean formula is slightly neutralized, then the result is correct, e.g. :

Integrate[ regFunc[x, y]*((4 x + 3 y) (3 x + 2 y))^4, {x, -10, 10}, {y, -10, 10}] // N

Integrate[ regFunc[x, y]*(12 x^2 + 17 x*y + 6 y^2)^4, {x, -10, 10}, {y, -10, 10}] // N

7836.43
7836.43

RegionPlot[ {-5 < 4 x + 3 y && 4 x + 3 y < 5 && -2 < 3 x + 2 y && 3 x + 2 y < 2,
-10 < x < 10 && -10 < y < 10 },
{x, -25, 25}, {y, -25, 25}, PlotPoints -> 150, MaxRecursion -> 4]


It should be emphasized that Integrate doesn't work either when we use insted of Boole for example UnitStep :

regFuncUS[x_, y_] := UnitStep[ 5 + 4 x + 3 y, 5 - 4 x - 3 y, 2 + 3 x + 2 y, 2 - 3 x - 2 y]


Edit 2

In Mathematica 9 this bug has been fixed :

Integrate[ regFunc[x,y] (( 4 x + 3 y )( 3 x + 2 y ))^4,{x, -100, 100},{y, -100, 100}] //N
Integrate[ regFunc[x,y] ( 12x^2 + 17 x y + 6 y^2 )^4,{x, -100, 100},{y, -100, 100}] //N

16000.
16000.

-
Ah, but what are the assumptions in this case? I guess the inequalities forces the variables to be real, but the region one integrates over is also a real domain... –  Per Alexandersson Mar 27 '12 at 15:53
Assumptions are in Boole, or in UnitStep if you want. Did I mention anything about Complexes ??? –  Artes Mar 27 '12 at 16:11
Well, yes, but the same assumptions are made in both integrals. What surprises me is that factorization of the integrand changes the answer. –  Per Alexandersson Mar 27 '12 at 17:42
I mean the problem is with the limits of integration, when you impose additionally boolean expression (assumptions or UnitStep) Integrate appears to work improperly. Note for example when you imopse smaller limits e.g. Integrate[...{x, -10, 10}, {y, -10, 10}] it seems to work properly. –  Artes Mar 27 '12 at 17:51
As an added note, this gives 16000: Integrate[regFunc[x,y] #, {x, -100, 100}, {y, -100, 100}]& /@ Expand[(12 x^2 + 17 x*y + 6 y^2)^4]. So, the problem is not entirely with the limits. –  rcollyer Mar 28 '12 at 1:30

It seems to me it's a roundoff error problem; indeed

Map[(Integrate[regFunc[x, y]*#, {x, -100, 100}, {y, -100, 100}] // N[#, 20] &) &,
((4 x + 3 y) (3 x + 2 y))^4 // Expand]


gives $160000$ correctly. So does the analytic approach:

Integrate[regFunc[x, y] #, {x, -100, 100}, {y, -100, 100}] & /@
Expand[(12 x^2 + 17 x*y + 6 y^2)^4]//N


whereas

Map[(Integrate[regFunc[x, y]*#, {x, -100, 100}, {y, -100, 100}] // N) &,
((4 x + 3 y) (3 x + 2 y))^4 // Expand]


gives the wrong answer.

Beware of numerical cancellations!

-