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?

  • $\begingroup$ What about the case when x1 x2 + y1 y2 is zero? Perhaps you want UnitStep? $\endgroup$
    – Mr.Wizard
    Commented Nov 11, 2012 at 19:59
  • $\begingroup$ You might want to try UnitStep. $\endgroup$
    – JohnD
    Commented Nov 11, 2012 at 20:25
  • $\begingroup$ @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? $\endgroup$
    – Rojo
    Commented Nov 11, 2012 at 20:33
  • $\begingroup$ @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 $\endgroup$
    – Rojo
    Commented Nov 11, 2012 at 20:34
  • 1
    $\begingroup$ 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). $\endgroup$
    – whuber
    Commented Nov 12, 2012 at 17:59

1 Answer 1


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

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.

  • 1
    $\begingroup$ I really don't understand why Mathematica has decided to keep HeavisideTheta[0] as unevaluated. Lots of problems for users would be avoided. $\endgroup$ Commented Nov 20, 2016 at 17:25
  • $\begingroup$ 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}.". $\endgroup$
    – user64494
    Commented Dec 5, 2023 at 6:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.