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.

Say I need to evaluate the integral $\iiint_W f(x,y,z) dx dy dz$ and $W$ is a region given to me like $W = \{ (x,y,z) : 1 \leq x^2 + y^2 \leq 4, 1 \leq z \leq 5\}$. I don't how to do this with a triple integral in Mathematica code.

share|improve this question
5  
Use Boole. Look here for more –  Spawn1701D Apr 18 '13 at 19:45
add comment

2 Answers

As it's been rightly pointed out in the comments, you can use Boole. Here's a simple example:

f[x_, y_, z_] = x^4 + y^2 + z;
Integrate[f[x, y, z]*Boole[1 < x^2 + y^2 < 4],
  {x, -2, 2}, {y, -2, 2}, {z, 1, 5}] // Timing

(* Out: {14.240965, (165 Pi)/2} *)

It certainly is well worth understanding the underlying transformations, though. In this example, cylindrical coordinates are very natural.

Integrate[f[r*Cos[t], r*Sin[t], z] r,
  {r, 1, 2}, {z, 1, 5}, {t, 0, 2 Pi}] // Timing

(* {0.498571, (165 Pi)/2} *)

Note that we got the same answer in much less time.

share|improve this answer
    
How are the maximums and minimums for $x$, $y$, and $z$ generated? –  Melab Apr 18 '13 at 20:48
    
@Melab Well, I'm not sure which input you're talking about. In the first input, using Boole, the bounds simply need to be large enough to contain the region; in fact, $-\inf$ to $\inf$ will work. In the second input, again, you need to understand the coordinate system. –  Mark McClure Apr 18 '13 at 21:24
    
Is it unable to determine the bounds dynamically since they do, after all, vary across the region? –  Melab Apr 18 '13 at 22:29
1  
@Melab, you misunderstand how Boole[] (i.e. Iverson brackets) works; the key here is that it is zero if its argument is false, and zero if its argument is true (i.e. within your domain of interest); one can then allow the integration limits to be infinite, since Boole[] zeroes out the parts outside the domain. –  J. M. Apr 19 '13 at 2:20
add comment

Clearly the integral over the region $W$ in the question is most easily set up in polar coordinates. But in case one wants to do it in rectangular coordinates or for other regions, Reduce can help. Suppose we want the integral in a particular order say {x, y, z}. Then Reduce will yield inequalities corresponding to the limits of the integral.

tmp = Reduce[{1 <= x^2 + y^2 <= 4, 1 <= z <= 5}, {x, y, z}, Reals]

(* 1 <= z <= 5 &&
    ((y == -2 && x == 0) ||
     (-2 < y < -1 && -Sqrt[4 - y^2] <= x <= Sqrt[4 - y^2]) ||
     (-1 <= y <= 1 && (-Sqrt[4 - y^2] <= x <= -Sqrt[1 - y^2] ||
                        Sqrt[1 - y^2] <= x <=  Sqrt[4 - y^2])) ||
     (1 < y < 2 && -Sqrt[4 - y^2] <= x <= Sqrt[4 - y^2]) ||
     (y == 2 && x == 0)) *)

Of course one sees that the order of the inequalities is complicated. One can eliminate sets of measure zero (in which == appears), and get a normal form of sorts with LogicalExpand:

Select[List @@ ((tmp /.
     {Inequality -> \[FormalI], Less -> \[FormalL], LessEqual -> \[FormalE]}) // 
   LogicalExpand) /.
     {\[FormalI] -> Inequality, \[FormalL] -> Less, \[FormalE] -> LessEqual},
 FreeQ[#, Equal] &]

(* {-2 < y < -1 && 1 <= z <= 5 && -Sqrt[4 - y^2] <= x <= Sqrt[4 - y^2],
    -1 <= y <= 1 && 1 <= z <= 5 && Sqrt[1 - y^2] <= x <= Sqrt[4 - y^2],
    -1 <= y <= 1 && 1 <= z <= 5 && -Sqrt[4 - y^2] <= x <= -Sqrt[1 - y^2], 
     1 <= z <= 5 && 1 < y < 2 && -Sqrt[4 - y^2] <= x <= Sqrt[4 - y^2]}  *)

I had to temporarily disable Inequality etc. since LogicalExpand breaks them apart. One can take advantage of the resulting inequalities to set up integrals automatically.

intLimits[eqns_List, order_List] := intLimits[And @@ eqns, order];
intLimits[eqns_, order_List] := Module[{redEq}, 
   redEq = Reduce[{eqns}, {order[[3]], order[[2]], order[[1]]}, Reals];
   Function[ineq,
     SortBy[(List @@ #)[[{3, 1, 5}]] & /@ List @@ ineq, 
      Position[Reverse@order, First@#] &]] /@ (If[
      FreeQ[redEq, Or], {redEq}, 
      Select[List @@ ((redEq /.
           {Inequality -> \[FormalI], Less -> \[FormalL], LessEqual -> \[FormalE]}) // 
         LogicalExpand) /.
           {\[FormalI] -> Inequality, \[FormalL] -> Less, \[FormalE] -> LessEqual},
        FreeQ[#, Equal] &]])];

We can get the limits in the z, y, x order:

limm = intLimits[1 <= x^2 + y^2 <= 4 && 1 <= z <= 5, {z, y, x}]

(* {{{x, -2, -1}, {y, -Sqrt[4 - x^2], Sqrt[4 - x^2]}, {z, 1, 5}},
    {{x, -1, 1}, {y, Sqrt[1 - x^2], Sqrt[4 - x^2]}, {z, 1, 5}},
    {{x, -1, 1}, {y, -Sqrt[4 - x^2], -Sqrt[1 - x^2]}, {z, 1, 5}},
    {{x, 1, 2}, {y, -Sqrt[4 - x^2], Sqrt[4 - x^2]}, {z, 1, 5}}}  *)

One can set up the integrals as follows:

int = HoldForm@Integrate[f[x, y, z], ##] & @@@ limm // Total

Integrals over hollow cylinder

Here are the corresponding regions:

Regions of integration

Here we test it on @MarkMcClure's test function:

int /. f[x, y, z] -> x^4 + y^2 + z // ReleaseHold // Timing

(* {4.645733, (165 π)/2} *)

And compare it to using Boole:

ff[x_, y_, z_] = x^4 + y^2 + z;
Integrate[ff[x, y, z]*Boole[1 < x^2 + y^2 < 4],
  {x, -2, 2}, {y, -2, 2}, {z, 1, 5}] // Timing

(* {11.202080, (165 π)/2} *)
share|improve this answer
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.