# Integral with HeavisideTheta takes too long to evaluate

I tried to compute $$\int_{-1}^1 d x_1 \int_{-1}^1 d x_2 \int_{-1}^1 d y_1 \int_{-1}^1 d y_2 \theta(x_1 x_2 + y_1 y_2)\,$$ where $\theta$ is Heaviside's step function, by using

Integrate[HeavisideTheta[x1 x2+y1 y2],{x1,-1,1},{x2,-1,1},{y1,-1,1},{y2,-1,1}]


but it took forever to evaluate. I didn't have the patience to wait until the end. However,

Integrate[Boole[x1 x2+y1 y2 >= 0],{x1,-1,1},{x2,-1,1},{y1,-1,1},{y2,-1,1}]


was very quick. Can anybody explain this strange behavior?

-
What about the case when x1 x2 + y1 y2 is zero? Perhaps you want UnitStep? – Mr.Wizard Nov 11 '12 at 19:59
You might want to try UnitStep. – JohnD Nov 11 '12 at 20:25
@Mr.Wizard it's an integral so that shouldn't matter, and that's the whole point of it not evaluating at 0 I think, right? – Rojo Nov 11 '12 at 20:33
@texasAUtiger, did you get better results with UnitStep? I got bored of waiting in both cases. But I have to admit I got bored quite fast – Rojo Nov 11 '12 at 20:34
This is not a problem with multiple integration: even Integrate[HeavisideTheta[x1 x2 + y1 y2], {x1, -1, 1}] will take extremely long (if it finishes at all). For another clue, giving this one about 10 seconds shows that MMA creates an impossible integrand (for it) after the first integration, suggesting the complications it's running into: Integrate[ HeavisideTheta[x1 x2 + y1 y2], {x1, -1, 0}, {x2, -1, 0}, {y1, -1, 0}, {y2, 0, 1}] (notice the limits of integration). Using UnitStep helps, but is extremely slow compared to Boole (which is used exclusively on Integrate's help page). – whuber Nov 12 '12 at 17:59

As outlined in one of the comments, I think Mathematica is running into boundary value problems with the HeavisideTheta. More specifically, that Heaviside remains unevaluated at 0 itself, but you can integrate a function against the Heaviside in an arbitrarily small, open region including zero. If you fix this

Timing[
x = 1 + \$MachineEpsilon;
Integrate[HeavisideTheta[x1 x2 + y1 y2], {x1, -x, x}, {x2, -x, x}, {y1, -x, x}, {y2, -x, x}]]


the integral computes

Out[7]= {13.1659, 8.}


You figured out that using Boole is faster and the better choice. So I guess the answer to your question is: issues with HeavisideTheta[0] remaining unevaluated.

-