Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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}]
share|improve this question
    
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
1  
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
add comment

2 Answers 2

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."
share|improve this answer
add comment

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 *)
share|improve this answer
    
Thanks! That's great. –  LBogaardt Dec 9 '13 at 16:01
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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