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

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.

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.

• 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
• @Melab, you misunderstand how Boole[] (i.e. Iverson brackets) works; the key here is that it is zero if its argument is false, and one 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.'s technical difficulties Apr 19 '13 at 2:20

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[], order[], order[]}, 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} *)


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}

• ...or alternatively put the limit into the region condition: region = ImplicitRegion[1 <= x^2 + y^2 <= 4 && 1 <= z <= 5, {x, y, z}] – Alexander Erlich Oct 17 '17 at 21:09