Integral with HeavisideTheta takes too long to evaluate

Bug introduced in 9.0 or earlier, persisting through 13.2 or later

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? Commented Nov 11, 2012 at 19:59
• You might want to try UnitStep. Commented Nov 11, 2012 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
Commented Nov 11, 2012 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
Commented Nov 11, 2012 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). Commented Nov 12, 2012 at 17:59

1 Answer

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.

• I really don't understand why Mathematica has decided to keep HeavisideTheta[0] as unevaluated. Lots of problems for users would be avoided. Commented Nov 20, 2016 at 17:25
• In 13.3.1 on Windows 10 this results in "Integrate::idiv: Integral of ((x1 x2+y1) (<<1>>))/y1+<<12>>+<<4>> does not converge on {-1,1}.". Commented Dec 5, 2023 at 6:32