How do I evaluate a double or triple integral over a region?

Question

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.

Answer 1

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.

Answer 2

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

Here are the corresponding regions:

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} *)

Answer 3

Using M10+ region functionality:

region = ImplicitRegion[
    1 < x^2+y^2 < 4,
    {{x, -2, 2}, {y, -2, 2}, {z, 1, 5}}
];

Integrate[
    x^4 + y^2 + z,
    {x, y, z} ∈ region
] //AbsoluteTiming

{0.016709, 165 π/2}

Note that one could also define the region without giving x and y explicit limits (the limits on z are necessary):

region = ImplicitRegion[
    1 < x^2+y^2 < 4,
    {x, y, {z, 1, 5}}
];

Integrate[
    x^4 + y^2 + z,
    {x, y, z} ∈ region
] //AbsoluteTiming

{0.017081, 165 π/2}