# NIntegrate using Boole doesn't work

This doesn't want to work for me - I don't know why not. It usually gives me 0, but sometimes also a different number. I know it should be a small, non-zero number.

myO[e1_, e2_] := (e1)/((e1) + (e2));
myU[e1_, e2_] := (1 - e2)/((1 - e1) + (1 - e2));
myV[p_, e1_, e2_] := p (e1)/(p (e1) + (1 - p) (1 - e2));
myW[p_, e1_, e2_] := p (1 - e1)/(p (1 - e1) + (1 - p) e2);

NIntegrate[ Boole[0 < e1 < e2 < 1/2 && 0 < e3 < 1/2 && e1 < p < 1 - e1 &&
e3 < myV[p, e1, e2] < 1 - e3 < myW[p, e1, e2] &&
myO[e1, e2] < e3 < myU[e1, e2]],
{p, 0, 1}, {e1, 0, 1/2}, {e2, 0, 1/2}, {e3, 0, 1/2}]

-
Do you mean that sometimes the same command gives different outputs with everything else the same ? –  b.gatessucks Dec 6 '13 at 15:13
Yes, it does. But more importantly, it doesn't give the correct value. It mostly gives zero, which cannot be true. –  LBogaardt Dec 6 '13 at 15:30
In version 9 I get the same result every time. If I change the method to MonteCarlo, I get a slightly different number every time, as expected. Raising MaxPoints makes it converge to the same result the default method gives: 0.00177529. In version 8 I get the same. What version are you using? –  Szabolcs Dec 7 '13 at 1:06
I am using Version 9. –  LBogaardt Dec 9 '13 at 15:59

One way to proceed is to try out portions and see where the problem lies. For instance:

NIntegrate[Boole[0 < e1 < e2 < 1/2 && 0 < e3 < 1/2 && e1 < p < 1 - e1],
{p, 0, 1}, {e1, 0, 1/2}, {e2, 0, 1/2}, {e3, 0, 1/2}]
0.0416667


is fine. But the inequalities

NIntegrate[Boole[e3 < myV[p, e1, e2] < 1 - e3 < myW[p, e1, e2]],
{p, 0, 1}, {e1, 0, 1/2}, {e2, 0, 1/2}, {e3, 0, 1/2}]
0.0239936 - 3.06317*10^-29 I


are giving (small) complex values. Most likely, this is because there are points where your function is trying to divide by zero. One culprit is the function myW, which, at 0,0,0 gives

myW[0, 0, 0]
Power::infy: "Infinite expression 1/0 encountered."

-

You can use Reduce or CylindricalDecomposition to break down the inequalities into forms that can be used as limits of integration. (Note: One might try BooleanMinimize instead of distr below to distribute And over Or, but it changes the order of the inequalities. Consequently, the intervals of integration end up out of order and NIntegrate does not work. One could sort them, but I opted for writing distr, which distributes from the bottom up.)

 (* distribute And over Or *)
distr[e_And] := Distribute[distr /@ e, Or];
distr[e_Or] := distr /@ e;
distr[e_] := e;
(* convert ineqs to limits of integration *)
sub = {(Less | LessEqual)[a_, b_, c_] :> {b, a, c},
HoldPattern[Inequality[a_, _, b_, _, c_]] :> {b, a, c}};

simp = Reduce[
0 < e1 < e2 < 1/2 && 0 < e3 < 1/2 && e1 < p < 1 - e1 &&
e3 < myV[p, e1, e2] < 1 - e3 < myW[p, e1, e2] &&
myO[e1, e2] < e3 < myU[e1, e2],
{e1, e2, e3, p}];
NIntegrate[1, ##] & @@@ (distr[simp] /. sub /. {And | Or -> List}) // Total

(* 0.00177529 *)

-
Thanks! That's great. –  LBogaardt Dec 9 '13 at 16:01