Chapter 32 — Multiple Integrals

32.1 From Single to Double Integrals

For thirteen chapters of integration you accumulated along a line. The definite integral $\int_a^b f(x)\,dx$ chops the interval $[a,b]$ into tiny pieces of width $\Delta x$, weights each piece by a height $f(x_i)$, sums, and takes a limit. The whole machine lived on a one-dimensional stage.

Now we raise the curtain on a two-dimensional stage. Instead of accumulating over an interval, we accumulate over a region $R$ in the plane. Instead of slicing into little widths $\Delta x$, we tile $R$ with tiny patches of area $\Delta A$. Instead of a graph, we have a surface $z = f(x,y)$ floating above the region. The result is the double integral:

$$\iint_R f(x, y) \, dA = \lim_{\|P\| \to 0} \sum_{i,j} f(x_i^*, y_j^*)\,\Delta A_{ij}.$$

The right-hand side is a double Riemann sum. We partition $R$ into small rectangles, pick a sample point $(x_i^*, y_j^*)$ in each, multiply the height $f$ there by the patch area $\Delta A_{ij} = \Delta x\,\Delta y$, and add up all the resulting little box volumes. As the partition $P$ gets finer (the mesh $\|P\| \to 0$), the sum converges — for any continuous $f$ on a reasonable region — to a single number, the double integral. The symbol $dA$ is "infinitesimal area," the limit of $\Delta A$.

Geometric Intuition. Picture a hilly surface $z = f(x,y)$ floating above a patch of floor $R$. Each term $f(x_i^*, y_j^*)\,\Delta A_{ij}$ is the volume of a thin rectangular column: floor area $\Delta A$ times height $f$. Stack up all the columns and you have a "skyline" approximation to the solid trapped between the floor and the surface. Refine the tiling and the staircase of columns melts into a smooth solid. When $f \ge 0$, the double integral $\iint_R f\,dA$ is exactly the volume of that solid — the three-dimensional sibling of "area under a curve."

This single picture carries three of the book's recurring themes at once: the double integral accumulates change over a region (Theme 1), it fuses an algebraic limit with a geometric volume (Theme 2), and it is built by approximation — boxes refining toward a solid (Theme 6).

Triple integrals $\iiint_E f \, dV$ push the same idea into 3D: tile a solid region $E$ with tiny boxes of volume $\Delta V$, weight each by $f$, sum, and take the limit. We will reach those in §32.7.

The Key Insight. A multiple integral is one number: the limit of a sum of "value times patch-size" over a region. Everything in this chapter — Fubini's theorem, polar coordinates, the Jacobian factor $r$ — is machinery for computing that number. The meaning never changes; only the bookkeeping does.

Linearity and the basic properties

Before computing anything, record the properties that come for free from the limit definition, exactly mirroring the single-variable case:

$$\iint_R \big(\alpha f + \beta g\big)\,dA = \alpha\iint_R f\,dA + \beta\iint_R g\,dA \quad\text{(linearity)},$$

and if $R = R_1 \cup R_2$ with $R_1, R_2$ overlapping only on their boundary,

$$\iint_R f\,dA = \iint_{R_1} f\,dA + \iint_{R_2} f\,dA \quad\text{(additivity over regions)}.$$

The second property is the one we exploit constantly: to integrate over an awkward region, cut it into friendly pieces and add.

32.2 Fubini's Theorem and the Double Integral over a Rectangle

The Riemann-sum definition tells us what the double integral means but gives no practical way to evaluate it — summing over a 2D grid by hand is hopeless. The escape is Fubini's theorem, which says we may integrate one variable at a time, collapsing a double integral into two ordinary single integrals.

Start with the simplest region, a rectangle $R = [a, b] \times [c, d]$.

Fubini's Theorem (rectangular case). If $f$ is continuous on $R = [a,b]\times[c,d]$, then $$\iint_R f(x, y) \, dA = \int_a^b \left(\int_c^d f(x, y) \, dy\right) dx = \int_c^d \left(\int_a^b f(x, y) \, dx\right) dy.$$

The two nested integrals on the right are called iterated integrals, and the theorem makes two promises at once: that an iterated integral equals the double integral, and that the order of integration does not matter. You may sweep horizontally first or vertically first; you get the same number.

The geometric reason is the method of slices from Chapter 18. Fix $x$; then $A(x) = \int_c^d f(x,y)\,dy$ is the cross-sectional area of the solid in the plane at that $x$. Integrating those slice-areas, $\int_a^b A(x)\,dx$, sweeps the slices through and recovers the volume — just as $V = \int A(x)\,dx$ did for solids of revolution. Fubini is the slice method dressed in two-variable clothing.

Procedure. For the inner integral, treat the other variable as a constant (it is, after all, frozen while you integrate). Then feed the result into the outer integral.

Example 1. Evaluate $\iint_R xy \, dA$ over $R = [0, 1] \times [0, 2]$.

Inner first over $y$, holding $x$ constant:

$$\int_0^1 \left(\int_0^2 xy \, dy\right) dx = \int_0^1 x\left[\frac{y^2}{2}\right]_0^2 dx = \int_0^1 x\cdot 2 \, dx = \big[x^2\big]_0^1 = 1.$$

Now confirm the other order, integrating over $x$ first:

$$\int_0^2 \left(\int_0^1 xy \, dx\right) dy = \int_0^2 y\left[\frac{x^2}{2}\right]_0^1 dy = \int_0^2 \frac{y}{2} \, dy = \left[\frac{y^2}{4}\right]_0^2 = 1. \checkmark$$

Both orders give $1$, as Fubini guarantees.

A shortcut worth knowing. When the integrand factors as $f(x,y) = g(x)\,h(y)$ and the region is a rectangle with constant limits, the double integral splits into a product of single integrals:

$$\iint_{[a,b]\times[c,d]} g(x)h(y)\,dA = \left(\int_a^b g(x)\,dx\right)\left(\int_c^d h(y)\,dy\right).$$

For Example 1, $\big(\int_0^1 x\,dx\big)\big(\int_0^2 y\,dy\big) = \tfrac12 \cdot 2 = 1$ in one stroke. This trick fails the moment the limits depend on the other variable — so it is a rectangle-only convenience.

Check Your Understanding. Evaluate $\displaystyle\iint_R (x + 2y)\,dA$ over $R = [0,2]\times[0,1]$.

Answer

By linearity, integrate the inner over $y$ first: $\int_0^1 (x+2y)\,dy = x\cdot 1 + [y^2]_0^1 = x + 1$. Then $\int_0^2 (x+1)\,dx = \big[\tfrac{x^2}{2}+x\big]_0^2 = (2+2) = 4$. (The factoring shortcut does not apply here because $x+2y$ is a sum, not a product.)

Historical Note. The theorem is named for Guido Fubini (1879–1943), the Italian mathematician who proved the general version for Lebesgue integrals in 1907. The idea of iterating is far older — Cauchy and others swapped integration orders routinely in the 1800s — but Fubini pinned down exactly when the swap is legal. Fubini fled Mussolini's Italy in 1938 (he was Jewish) and spent his last years at the Institute for Advanced Study in Princeton.

Warning. Fubini's permission to swap orders relies on $f$ being "nice" — continuous, or at least absolutely integrable ($\iint |f|\,dA < \infty$). For unbounded integrands on infinite regions the two iterated integrals can genuinely disagree. The standard counterexample $f(x,y)=\frac{x^2-y^2}{(x^2+y^2)^2}$ on the unit square integrates to $+\tfrac{\pi}{4}$ one way and $-\tfrac{\pi}{4}$ the other. For every example in this chapter the hypotheses hold, but never swap orders blindly on an improper integral.

32.3 Double Integrals over General Regions

Rectangles are a luxury; real regions are triangles, disks, and curved blobs. The fix is to let the inner limits be functions while the outer limits stay constant. Two standard descriptions cover most regions.

Type I (vertically simple): $R = \{(x, y) : a \leq x \leq b,\ g_1(x) \leq y \leq g_2(x)\}$

Here $x$ ranges over a fixed interval $[a,b]$, and for each $x$ the variable $y$ runs between a bottom curve $g_1(x)$ and a top curve $g_2(x)$:

$$\iint_R f \, dA = \int_a^b \int_{g_1(x)}^{g_2(x)} f(x, y) \, dy \, dx.$$

Geometric Intuition. A Type I region is one you can sweep with vertical strips. Slide a vertical line from $x=a$ to $x=b$; at each position it enters the region at $y=g_1(x)$ and exits at $y=g_2(x)$. The inner integral $\int_{g_1}^{g_2} f\,dy$ accumulates along that strip; the outer integral sweeps the strip across. The inner limits are functions precisely because the top and bottom of the strip move as $x$ changes.

Type II (horizontally simple): $R = \{(x, y) : c \leq y \leq d,\ h_1(y) \leq x \leq h_2(y)\}$

Now the roles swap: $y$ ranges over $[c,d]$, and for each $y$, $x$ runs between a left curve $h_1(y)$ and a right curve $h_2(y)$. Sweep with horizontal strips:

$$\iint_R f \, dA = \int_c^d \int_{h_1(y)}^{h_2(y)} f(x, y) \, dx \, dy.$$

The Key Insight. The outer limits of an iterated integral are always constants; the inner limits may be functions of the outer variable, never of the inner variable. If you ever write an iterated integral whose outer limits contain a variable, you have made an error — the final answer is a number, and the outermost integral must produce one.

Always sketch the region first

The single most reliable way to get the limits right is to draw the region, then read the limits off the picture. Decide whether vertical strips (Type I) or horizontal strips (Type II) are simpler, mark where a strip enters and exits, and the limits write themselves.

Example 2 (a triangle, both ways). Integrate over the triangle $T$ with vertices $(0, 0)$, $(1, 0)$, $(1, 1)$. The slanted edge is the line $y = x$.

Type I (vertical strips): $x$ runs from $0$ to $1$. For a fixed $x$, the strip enters at the bottom edge $y=0$ and exits at the slant $y=x$. So $0 \le y \le x$: $$\iint_T f\,dA = \int_0^1 \int_0^x f(x,y) \, dy \, dx.$$

Type II (horizontal strips): $y$ runs from $0$ to $1$. For a fixed $y$, the strip enters at the slant $x=y$ and exits at the right edge $x=1$. So $y \le x \le 1$: $$\iint_T f\,dA = \int_0^1 \int_y^1 f(x,y) \, dx \, dy.$$

Both describe the same triangle and give the same integral. Choose whichever makes the integrand easier — which is exactly the subject of the next subsection.

Let us test the setup with a concrete integrand, $f(x,y)=x$, by Type I: $$\int_0^1\int_0^x x\,dy\,dx = \int_0^1 x\cdot[y]_0^x\,dx = \int_0^1 x\cdot x\,dx = \int_0^1 x^2\,dx = \tfrac13.$$

Reversing the order of integration

Sometimes one order is merely tedious and the other is impossible — because the inner integrand has no elementary antiderivative. Reversing the order can rescue you.

Example 3 (an "impossible" integral made easy). Evaluate $$\int_0^1 \int_x^1 \sin(y^2) \, dy \, dx.$$

The inner integral asks for $\int \sin(y^2)\,dy$, a Fresnel integral with no elementary antiderivative (Chapter 17 catalogues such integrals). We are stuck — in this order.

So reverse it. First, read the region from the given limits: $0 \le x \le 1$ and $x \le y \le 1$. This is the triangle bounded by $y = x$, $y = 1$, and $x = 0$. Now re-describe it with $y$ on the outside. For each $y \in [0,1]$, $x$ runs from the left edge $x=0$ to the slant $x = y$:

$$\int_0^1 \int_0^y \sin(y^2) \, dx \, dy = \int_0^1 \sin(y^2)\,[x]_0^y\,dy = \int_0^1 y\,\sin(y^2) \, dy.$$

Now the extra factor of $y$ — handed to us for free by the reversal — is exactly the chain-rule factor we need. Substitute $u = y^2$, $du = 2y\,dy$, so $y\,dy = \tfrac12\,du$, with limits $u : 0 \to 1$:

$$\int_0^1 y\sin(y^2)\,dy = \frac{1}{2}\int_0^1 \sin u \, du = \frac{1}{2}\big[-\cos u\big]_0^1 = \frac{1}{2}(1 - \cos 1) \approx 0.2298.$$

An integral that was impossible in one order is elementary in the other. Recognizing this is one of the signature skills of multivariable integration.

Common Pitfall. When you reverse the order, students very often just swap the limits in place — writing $\int_x^1\!\!\int_0^1$, leaving the variable $x$ stranded in a limit it can no longer reach. You must re-derive the limits from the region, not shuffle the old ones. Sketch the region, choose the new strip direction, and read off fresh limits. The new outer limits must be pure numbers, and the new inner limits must involve only the new outer variable.

Computational Note. A reliable algorithm for reversing order: (1) sketch the region from the original limits; (2) describe it as the other type; (3) write the new inner limits as the entry/exit curves of the new strip, and the new outer limits as the overall range of the new outer variable. When the region is a single triangle or a lens between two curves, this is mechanical. When it is more complex, additivity (§32.1) lets you split it into pieces and reverse each piece separately.

32.4 Area and Volume

Two special integrands turn the double integral into a measuring tape.

Area. Set $f \equiv 1$. Each column has height $1$, so its volume is just its floor area, and summing recovers the area of $R$ itself:

$$\text{Area}(R) = \iint_R 1 \, dA.$$

For a Type I region this unwinds to something familiar:

$$\int_a^b \int_{g_1(x)}^{g_2(x)} 1 \, dy \, dx = \int_a^b \big(g_2(x) - g_1(x)\big) \, dx,$$

the area between two curves from Chapter 18. The double integral did not invent a new idea; it contains the single-variable area formula as the case $f=1$. That is the right way to see every formula in this chapter — as a 2D generalization that collapses back to something you already trust.

Volume. When $f \ge 0$, the double integral is the volume under the surface, as §32.1 argued.

Example 4 (volume under a paraboloid). Find the volume under $z = x^2 + y^2$ over the square $[0,1]\times[0,1]$. The region is a rectangle, so:

$$V = \int_0^1 \int_0^1 (x^2 + y^2) \, dy \, dx = \int_0^1 \left[x^2 y + \frac{y^3}{3}\right]_{y=0}^{y=1} dx = \int_0^1 \left(x^2 + \frac{1}{3}\right) dx = \frac{1}{3} + \frac{1}{3} = \frac{2}{3}.$$

The bowl-shaped surface encloses a volume of $\tfrac23$ above the unit square.

Average value. Dividing the integral by the area gives the average value of $f$ over $R$, the exact 2D analog of the §14.9 average:

$$\bar f = \frac{1}{\text{Area}(R)}\iint_R f(x,y)\,dA.$$

For Example 4, the average height of the paraboloid over the unit square is $\frac{2/3}{1} = \tfrac23$.

32.5 Polar Coordinates for Double Integrals

Some regions and integrands are built around a center, not a corner: disks, rings, anything featuring $x^2 + y^2$. For these, Cartesian rectangles are the wrong tiles. Polar coordinates tile the plane with the region's own grain.

Recall from Chapter 26 the polar substitution $x = r\cos\theta$, $y = r\sin\theta$, so that $x^2 + y^2 = r^2$. The crucial new ingredient is the polar area element:

$$dA = r \, dr \, d\theta.$$

The extra factor of $r$ is not optional, and forgetting it is the defining error of polar integration. Here is why it is there.

Geometric Intuition — where the $r$ comes from. Tile the plane not with little Cartesian rectangles but with little polar patches: pick a thin wedge of angular width $d\theta$ and a thin ring of radial width $dr$. Their intersection is a curved patch. Its radial side has length $dr$. Its arc side, at radius $r$, has length $r\,d\theta$ (arc length $=$ radius $\times$ angle). So the patch is approximately a rectangle of sides $dr$ and $r\,d\theta$, giving area $dA \approx r\,dr\,d\theta$. The patches are bigger far from the origin — a one-degree wedge spans more area at radius 10 than at radius 1 — and the factor $r$ is precisely that growth. This is the first instance of a general phenomenon: a coordinate change rescales area by a local factor. Chapter 33 will name that factor the Jacobian determinant and prove that for polar coordinates it equals exactly $r$. For now, accept it on the strength of the picture.

The conversion recipe: replace $x^2+y^2$ by $r^2$, replace the region by its $(r,\theta)$ description, and replace $dA$ by $r\,dr\,d\theta$.

Example 5 (a Gaussian over a disk). Evaluate $\iint_D e^{-(x^2+y^2)} \, dA$ over the unit disk $D : x^2 + y^2 \leq 1$.

In Cartesian coordinates this is hopeless — $\int e^{-x^2}dx$ has no elementary form. In polar coordinates the disk is simply $0 \le r \le 1$, $0 \le \theta \le 2\pi$, and the integrand becomes $e^{-r^2}$:

$$\iint_D e^{-(x^2+y^2)} \, dA = \int_0^{2\pi} \int_0^1 e^{-r^2}\, r \, dr \, d\theta.$$

The factor of $r$ from $dA$ is exactly what the integrand was missing. Substitute $u = r^2$, $du = 2r\,dr$:

$$\int_0^1 e^{-r^2}\,r\,dr = \frac{1}{2}\int_0^1 e^{-u}\,du = \frac{1}{2}\big(1 - e^{-1}\big).$$

The angular integral is just $\int_0^{2\pi}d\theta = 2\pi$, so

$$\iint_D e^{-(x^2+y^2)} \, dA = 2\pi \cdot \frac{1}{2}(1 - e^{-1}) = \pi\left(1 - \frac{1}{e}\right) \approx 1.986.$$

Polar coordinates did not just simplify the algebra — they turned a non-elementary integral into a one-line substitution.

Common Pitfall. The single most frequent error in this entire chapter is writing $dA = dr\,d\theta$ and dropping the $r$. The area element is $r\,dr\,d\theta$, always. A quick sanity check: the area of the unit disk should be $\pi$. With the $r$: $\int_0^{2\pi}\!\int_0^1 r\,dr\,d\theta = 2\pi\cdot\tfrac12 = \pi$. ✓ Without it: $\int_0^{2\pi}\!\int_0^1 dr\,d\theta = 2\pi$, which is wrong. If your disk areas come out as $2\pi$ instead of $\pi$, you forgot the $r$.

The Gaussian integral $\int_{-\infty}^\infty e^{-x^2}\,dx = \sqrt{\pi}$

We now reach one of the most beautiful computations in all of calculus — and one this book has been promising you since Chapter 13. The function $e^{-x^2}$ is the unnormalized bell curve, the heart of the normal distribution (our anchor example since the definite integral first appeared). Its integral over the whole line, $\int_{-\infty}^\infty e^{-x^2}\,dx$, is the normalization constant of all of probability theory. Yet $e^{-x^2}$ has no elementary antiderivative (Liouville's theorem, mentioned in §14.12), so no amount of one-variable cleverness can crack it.

The trick — due to Poisson in the early 1800s — is to refuse to compute the integral and instead compute its square, by promoting it to two dimensions. Call the unknown value $I$:

$$I = \int_{-\infty}^\infty e^{-x^2}\,dx.$$

Square it, and rename the dummy variable in the second copy from $x$ to $y$ (legal — it is just a label):

$$I^2 = \left(\int_{-\infty}^\infty e^{-x^2}\,dx\right)\left(\int_{-\infty}^\infty e^{-y^2}\,dy\right) = \int_{-\infty}^\infty\int_{-\infty}^\infty e^{-x^2}e^{-y^2}\,dx\,dy.$$

The two factors recombine into a single double integral over the whole plane, because $e^{-x^2}e^{-y^2} = e^{-(x^2+y^2)}$:

$$I^2 = \iint_{\mathbb{R}^2} e^{-(x^2+y^2)} \, dA.$$

And this integral — unlike its one-dimensional parent — succumbs to polar coordinates. The whole plane is $0 \le r < \infty$, $0 \le \theta \le 2\pi$:

$$I^2 = \int_0^{2\pi} \int_0^\infty e^{-r^2}\, r \, dr \, d\theta.$$

The radial integral converges (it is an improper integral of the Chapter 17 kind): with $u = r^2$,

$$\int_0^\infty e^{-r^2}\,r\,dr = \frac{1}{2}\int_0^\infty e^{-u}\,du = \frac{1}{2}\big[-e^{-u}\big]_0^\infty = \frac{1}{2}(0 - (-1)) = \frac{1}{2}.$$

Therefore

$$I^2 = \int_0^{2\pi}\frac{1}{2}\,d\theta = 2\pi\cdot\frac{1}{2} = \pi \quad\Longrightarrow\quad \boxed{\;I = \int_{-\infty}^\infty e^{-x^2}\,dx = \sqrt{\pi}\;}.$$

(We take the positive root because $e^{-x^2} > 0$ makes $I > 0$.)

The Key Insight. A one-dimensional integral that cannot be done was solved by lifting it into two dimensions, where rotational symmetry makes the magic factor $r$ appear and the substitution $u = r^2$ finishes the job. This is the payoff of the area-under-the-normal-curve anchor example: Chapter 13 introduced the bell curve, Chapter 17 will use $\sqrt{\pi}$ to declare the normal integral convergent, and Chapter 23 will approximate the partial areas with Taylor series. The total area, the keystone that normalizes the whole distribution, comes from this polar double-integral trick and nowhere else.

Real-World Application — the normal distribution (statistics & data science). The standard normal density is $\phi(x) = \frac{1}{\sqrt{2\pi}}e^{-x^2/2}$. The constant $\frac{1}{\sqrt{2\pi}}$ exists for exactly one reason: it forces $\int_{-\infty}^\infty \phi(x)\,dx = 1$, so that total probability is $1$. A quick substitution $x = \sqrt{2}\,t$ turns $\int e^{-x^2/2}dx$ into $\sqrt2\int e^{-t^2}dt = \sqrt2\cdot\sqrt\pi = \sqrt{2\pi}$, which is precisely the reciprocal of that constant. Every $z$-score, every confidence interval, every hypothesis test in statistics rests on the number $\sqrt{\pi}$ you just computed with a polar double integral. The same Gaussian, in 2D, becomes the point-spread function of optical and image-processing systems — the blur kernel behind a camera's bokeh.

32.6 Triple Integrals

Everything generalizes one dimension upward. For a solid region $E \subset \mathbb{R}^3$,

$$\iiint_E f(x, y, z) \, dV = \lim_{\|P\|\to 0}\sum_{i,j,k} f(x_i^*,y_j^*,z_k^*)\,\Delta V_{ijk},$$

the limit of sums of $f$ times tiny box-volumes $\Delta V = \Delta x\,\Delta y\,\Delta z$. With $f = 1$ this computes the volume of $E$; with $f = \rho$ a density, it computes mass (§32.10). When $f \ge 0$, a triple integral has no "volume under a graph" picture — that would need a fourth dimension — so we lean on the mass interpretation instead: think of $f$ as density and the integral as total mass.

Fubini extends, too: a triple integral becomes three nested single integrals. The standard "$z$-simple" setup integrates $z$ first, between a bottom surface and a top surface, then integrates the shadow region $D$ in the $xy$-plane as a double integral:

$$E = \{(x, y, z) : (x, y) \in D,\ u_1(x, y) \leq z \leq u_2(x, y)\} \implies \iiint_E f \, dV = \iint_D \left(\int_{u_1(x,y)}^{u_2(x,y)} f(x, y, z) \, dz\right) dA.$$

Example 6 (volume between a paraboloid and a plane). Find the volume of the solid $E$ trapped between the paraboloid $z = x^2 + y^2$ below and the plane $z = 4$ above.

The two surfaces meet where $x^2 + y^2 = 4$, a circle of radius $2$; its interior disk $D$ is the shadow. For each $(x,y)\in D$, $z$ runs from the paraboloid up to the plane. Integrate $f=1$ over $z$ first:

$$V = \iiint_E 1\,dV = \iint_D \int_{x^2+y^2}^{4} 1\, dz\, dA = \iint_D \big(4 - (x^2+y^2)\big)\, dA.$$

The shadow $D$ is a disk and the integrand depends only on $x^2+y^2$ — a flashing neon sign for polar coordinates. With $D: 0\le r\le 2$, $0\le\theta\le2\pi$ and $dA = r\,dr\,d\theta$:

$$V = \int_0^{2\pi}\int_0^2 (4 - r^2)\,r\,dr\,d\theta = 2\pi\int_0^2 (4r - r^3)\,dr = 2\pi\left[2r^2 - \frac{r^4}{4}\right]_0^2 = 2\pi(8 - 4) = 8\pi.$$

Notice the structure: the triple integral reduced to a double integral, and the double integral was best handled in polar coordinates. That layering — peel off one variable, then change coordinates for the rest — is the everyday rhythm of multiple integration.

Check Your Understanding. Set up (do not evaluate) the triple integral for the volume of the solid bounded below by the plane $z = 0$, above by the plane $z = 1 - x - y$, with shadow the triangle $0 \le x \le 1$, $0 \le y \le 1 - x$.

Answer

$\displaystyle V = \int_0^1\int_0^{1-x}\int_0^{1-x-y} 1\,dz\,dy\,dx$. The innermost limits use the two $z$-surfaces; the middle limits use the triangle's slanted edge $y = 1-x$; the outer limits are the constant range of $x$. (Evaluating gives $\tfrac16$, the volume of a corner tetrahedron.)

32.7 Cylindrical Coordinates

When a solid has an axis of symmetry — a cylinder, a cone, a paraboloid of revolution — wrap the $xy$-plane in polar coordinates and leave $z$ alone. These are cylindrical coordinates $(r, \theta, z)$:

$$x = r\cos\theta, \qquad y = r\sin\theta, \qquad z = z.$$

Since cylindrical coordinates are "polar in $x,y$, untouched in $z$," the volume element inherits the polar factor:

$$dV = r \, dr \, d\theta \, dz.$$

Example 7 (Example 6 again, in cylindrical coordinates). The paraboloid–plane solid is the natural test case, since it is a surface of revolution. In cylindrical coordinates the paraboloid $z = x^2+y^2$ is simply $z = r^2$, and the solid is $0\le\theta\le2\pi$, $0\le r\le 2$, $r^2 \le z \le 4$:

$$V = \int_0^{2\pi}\int_0^2\int_{r^2}^4 r \, dz\, dr\, d\theta = \int_0^{2\pi}\int_0^2 r(4 - r^2)\,dr\,d\theta = 8\pi. \checkmark$$

Same answer as Example 6, but the setup was almost effortless because cylindrical coordinates match the solid's symmetry.

32.8 Spherical Coordinates

When the symmetry is about a point — balls, shells, anything featuring $x^2+y^2+z^2$ — use spherical coordinates $(\rho, \phi, \theta)$. Here $\rho$ is the distance from the origin, $\phi$ is the polar angle measured down from the positive $z$-axis, and $\theta$ is the same azimuthal angle as before:

$$x = \rho\sin\phi\cos\theta, \qquad y = \rho\sin\phi\sin\theta, \qquad z = \rho\cos\phi,$$

with $\rho \geq 0$, $0 \leq \phi \leq \pi$, $0 \leq \theta < 2\pi$. The crowning fact is the spherical volume element:

$$dV = \rho^2 \sin\phi \, d\rho \, d\phi \, d\theta.$$

Geometric Intuition — the $\rho^2\sin\phi$ factor. Build a tiny spherical box from increments $d\rho$, $d\phi$, $d\theta$. Its radial edge has length $d\rho$. As $\phi$ changes by $d\phi$, the point sweeps an arc on a circle of radius $\rho$, length $\rho\,d\phi$. As $\theta$ changes by $d\theta$, it sweeps an arc on a circle of radius $\rho\sin\phi$ (the distance from the $z$-axis, which shrinks toward the poles), length $\rho\sin\phi\,d\theta$. Multiply the three nearly-perpendicular edges: $dV \approx d\rho\cdot\rho\,d\phi\cdot\rho\sin\phi\,d\theta = \rho^2\sin\phi\,d\rho\,d\phi\,d\theta$. The $\rho^2$ says far-away patches are larger; the $\sin\phi$ says patches shrink to nothing at the poles ($\phi=0,\pi$), where all longitude lines converge. Chapter 33 confirms this is the Jacobian of the spherical map.

Example 8 (the volume of a ball, finally derived). A ball of radius $R$ is $0\le\rho\le R$, $0\le\phi\le\pi$, $0\le\theta\le2\pi$ — constant limits, the dream scenario. Because the integrand $\rho^2\sin\phi$ factors and the limits are constants, the triple integral splits into a product of three single integrals:

$$V = \int_0^{2\pi}\!\!\int_0^\pi\!\!\int_0^R \rho^2\sin\phi \, d\rho\, d\phi\, d\theta = \left(\int_0^{2\pi}d\theta\right)\!\left(\int_0^\pi \sin\phi\,d\phi\right)\!\left(\int_0^R \rho^2\,d\rho\right).$$

Evaluate each factor: $\int_0^{2\pi}d\theta = 2\pi$; $\int_0^\pi\sin\phi\,d\phi = [-\cos\phi]_0^\pi = -(-1)-(-1) = 2$; $\int_0^R \rho^2\,d\rho = R^3/3$. Multiply:

$$V = 2\pi \cdot 2 \cdot \frac{R^3}{3} = \frac{4}{3}\pi R^3.$$

The formula you memorized in grade school, derived from a triple integral. There is real satisfaction in watching $\tfrac43\pi R^3$ emerge from $\rho^2\sin\phi$ — geometry and algebra agreeing once again (Theme 2).

Common Pitfall. Two warnings about spherical coordinates. First, conventions differ: physics texts often swap the names of $\phi$ and $\theta$. Always state which angle is which. In this book $\phi$ is the angle down from the $z$-axis and $\theta$ is the azimuth. Second, the volume element is $\rho^2\sin\phi$, not $\rho^2$ and not $\rho\sin\phi$ — both pieces are essential. A fast check: the unit ball must have volume $\tfrac43\pi$. Get $\tfrac43\pi$ and your element is right; get anything else and re-examine it.

32.9 Choosing a Coordinate System

A short field guide, because picking the right coordinates is half the battle:

The choice can transform days of work into minutes. When in doubt, look at the region's boundary and the integrand together: if both speak the language of circles or spheres, switch coordinates; if both are flat, stay Cartesian.

32.10 Applications: Mass, Center of Mass, and Moment of Inertia

Multiple integrals are how we weigh, balance, and spin continuous objects. The unifying idea is always the same: mass is the integral of density.

Mass

A flat plate occupying region $R$ with variable area-density $\rho(x,y)$ (mass per unit area) has total mass

$$M = \iint_R \rho(x, y) \, dA.$$

For a 3D solid $E$ with volume-density $\rho(x,y,z)$,

$$M = \iiint_E \rho(x, y, z) \, dV.$$

Each tiny patch contributes density-times-size; the integral sums them into the total mass.

Center of mass

The center of mass $(\bar x, \bar y, \bar z)$ is the balance point — the average position weighted by mass. Each coordinate is a moment (mass times lever arm) divided by the total mass:

$$\bar{x} = \frac{1}{M} \iiint_E x\,\rho \, dV, \qquad \bar{y} = \frac{1}{M} \iiint_E y\,\rho \, dV, \qquad \bar{z} = \frac{1}{M} \iiint_E z\,\rho \, dV.$$

When the density is constant, $\rho$ cancels between numerator and denominator and the center of mass reduces to the purely geometric centroid, $\bar x = \frac{1}{\text{Area}}\iint_R x\,dA$ (and likewise for the other coordinates).

Example 9 (centroid of a triangle). Find the centroid of the triangle $T$ with vertices $(0,0)$, $(1,0)$, $(1,1)$, taking uniform density. Its area is $\tfrac12$. Using Type I limits $0\le x\le1$, $0\le y\le x$:

$$\iint_T x\,dA = \int_0^1\int_0^x x\,dy\,dx = \int_0^1 x^2\,dx = \frac{1}{3}, \qquad \iint_T y\,dA = \int_0^1\int_0^x y\,dy\,dx = \int_0^1 \frac{x^2}{2}\,dx = \frac{1}{6}.$$

So $\bar x = \frac{1/3}{1/2} = \frac{2}{3}$ and $\bar y = \frac{1/6}{1/2} = \frac{1}{3}$. The centroid is $\big(\tfrac23, \tfrac13\big)$ — which is exactly the average of the three vertices, $\frac{1}{3}\big((0,0)+(1,0)+(1,1)\big) = \big(\tfrac23,\tfrac13\big)$, as it must be for a triangle. The integral and the vertex-average agree.

Moment of inertia

Where center of mass weights position by mass, moment of inertia weights squared distance from an axis by mass — it measures how hard an object is to spin. For rotation about the $z$-axis, the distance from the axis is $\sqrt{x^2+y^2}$, so:

$$I_z = \iiint_E (x^2 + y^2)\,\rho \, dV.$$

The integrand is "(distance from axis)$^2 \times$ density." Moments of inertia are why a figure skater spins faster when she pulls her arms in (reducing $x^2+y^2$ for her mass) and why flywheels store rotational energy. In cylindrical coordinates, $x^2+y^2 = r^2$, so $I_z$ becomes $\iiint r^2\rho\,(r\,dr\,d\theta\,dz)$ — the integrand acquires a clean $r^3$ when $\rho$ is constant.

Real-World Application — engineering & structures (moments of inertia). A beam's resistance to bending and a shaft's resistance to twisting are governed by area moments of inertia, $I = \iint_R y^2\,dA$, computed by exactly these double integrals. This is why steel I-beams put most of their material far from the center (large $y^2$, large $I$, stiff beam) while using little material near the neutral axis. Every bridge, crane, and aircraft wing is sized by an engineer evaluating moment-of-inertia integrals over the cross-section.

32.11 Application: Probability and Joint Distributions

In one variable, a probability density $p(x)$ satisfies $p \ge 0$ and $\int_{-\infty}^\infty p\,dx = 1$, and probabilities are areas (the Chapter 14 net-change reading of the CDF). In two variables, a joint probability density $p(x,y)$ describes a pair of random quantities $(X, Y)$ at once, and probabilities become volumes under the density surface:

$$P\big((X, Y) \in R\big) = \iint_R p(x, y) \, dA, \qquad \iint_{\mathbb{R}^2} p \, dA = 1.$$

From the joint density you recover the behavior of $X$ alone by integrating out $y$ — the marginal density:

$$p_X(x) = \int_{-\infty}^{\infty} p(x, y) \, dy.$$

This integrates over every possible value of the variable you do not care about, leaving the distribution of the one you do.

The 2D Gaussian. The independent standard bivariate normal density is

$$p(x, y) = \frac{1}{2\pi} \, e^{-(x^2+y^2)/2},$$

the same bell shape spun about the origin into a smooth hill. That it integrates to $1$ is the Gaussian trick of §32.5 in disguise — switch to polar and the $r$ in $dA$ makes the integral elementary. The normalizing constant $\frac{1}{2\pi}$ is, once again, a consequence of $\int e^{-x^2}dx = \sqrt\pi$. This same Gaussian hill is the blur kernel in image processing: a Gaussian blur replaces each pixel by a $p$-weighted average of its neighbors, and the weights must integrate to $1$ so the image neither brightens nor darkens.

Check Your Understanding. A point $(X,Y)$ is uniformly distributed on the unit square $[0,1]^2$, meaning $p(x,y) = 1$ there and $0$ outside. What is $P(X + Y \le 1)$?

Answer

$P(X+Y\le 1) = \iint_R 1\,dA$ where $R$ is the part of the square below the line $x+y=1$ — a right triangle with legs $1$. Its area is $\tfrac12$, so the probability is $\tfrac12$. Because the density is the constant $1$, the probability is literally the area of the region, a clean illustration of $\iint_R 1\,dA = \text{Area}(R)$.

32.12 Computation: Multiple Integration in Python

Hand computation builds the understanding; Python supplies the power (Theme 4). For double and triple integrals, scipy provides dblquad and tplquad. The one quirk to remember: scipy expects the integrand's arguments in inner-to-outer order, so a function for $\iint f\,dy\,dx$ is written f(y, x).

# Double and triple integrals with scipy, checked against hand answers.
from scipy.integrate import dblquad, tplquad

# (1) Double integral of (x*y + x + y) over the rectangle [0,1] x [0,1].
#     Note argument order: dblquad calls f(y, x) -- inner variable first.
f = lambda y, x: x*y + x + y
val2, _ = dblquad(f, 0, 1, lambda x: 0, lambda x: 1)   # x in [0,1], y in [0,1]
print(f"double integral = {val2:.6f}")   # double integral = 1.083333  (= 13/12)

# (2) A NON-rectangular region: integrate f(x,y) = 1 over the triangle
#     0 <= x <= 1, 0 <= y <= x. The answer is the triangle's area, 1/2.
area, _ = dblquad(lambda y, x: 1.0, 0, 1, lambda x: 0, lambda x: x)
print(f"triangle area   = {area:.6f}")   # triangle area   = 0.500000

# (3) Triple integral of (x + y + z) over the unit cube [0,1]^3.
g = lambda z, y, x: x + y + z
val3, _ = tplquad(g, 0, 1, 0, 1, 0, 1)
print(f"triple integral = {val3:.6f}")   # triple integral = 1.500000  (= 3/2)

The non-rectangular case (2) is the important one: the inner limits are functions lambda x: 0 and lambda x: x, exactly the Type I description from §32.3. We can also verify the Gaussian payoff numerically:

# Verify the Gaussian double-integral trick: I^2 = pi, so I = sqrt(pi).
import numpy as np
from scipy.integrate import dblquad

gauss2d = lambda y, x: np.exp(-(x**2 + y**2))
# Integrate over a large box approximating all of R^2 (the tails are negligible).
I2, _ = dblquad(gauss2d, -8, 8, lambda x: -8, lambda x: 8)
print(f"I^2 (numeric) = {I2:.6f},  pi = {np.pi:.6f}")  # I^2 (numeric) = 3.141593, pi = 3.141593
print(f"I = sqrt(pi)? {np.sqrt(I2):.6f} vs {np.sqrt(np.pi):.6f}")  # 1.772454 vs 1.772454

The numerical double integral over a large box reproduces $\pi$ to six digits — machine and polar-coordinate hand computation shaking hands, exactly as our FTC checks did in Chapter 14.

Computational Note — the curse of dimensionality. Deterministic methods like tplquad work beautifully through three or four dimensions but collapse beyond that: a grid of $n$ points per axis needs $n^d$ total points, exploding exponentially with dimension $d$. For high-dimensional integrals — common in physics (path integrals), finance (option pricing), and Bayesian statistics (posterior expectations) — the standard escape is Monte Carlo integration: sample $N$ random points in the region and estimate $\iint_R f \approx \text{vol}(R)\cdot\bar f$, the volume times the sample mean of $f$. Its error shrinks like $1/\sqrt{N}$ regardless of dimension, which is why a 100-dimensional integral is tractable by sampling but utterly hopeless by grids.

32.13 A Worked Triple Integral: Mass of a Hemisphere

To tie the threads together, here is one fully worked applied problem combining spherical coordinates with a variable density.

Problem. A solid hemisphere of radius $R$ (the top half of the ball, $z \ge 0$) has density proportional to the distance from the center, $\rho = k\rho_{\text{dist}} = k\rho$ — that is, $\rho_{\text{density}}(\rho) = k\rho$ where $\rho$ is the spherical radial coordinate and $k$ is a constant. Find its mass.

Setup. The hemisphere in spherical coordinates is $0 \le \rho \le R$, $0 \le \phi \le \tfrac{\pi}{2}$ (the polar angle only reaches the equator, not the south pole, because $z \ge 0$), and $0 \le \theta \le 2\pi$. The mass integral is

$$M = \iiint_E \rho_{\text{density}}\,dV = \int_0^{2\pi}\!\int_0^{\pi/2}\!\int_0^R (k\rho)\,\underbrace{\rho^2\sin\phi}_{dV}\,d\rho\,d\phi\,d\theta.$$

Evaluate. The integrand $k\,\rho^3\sin\phi$ factors and the limits are constant, so split into three single integrals:

$$M = k\left(\int_0^{2\pi}d\theta\right)\!\left(\int_0^{\pi/2}\sin\phi\,d\phi\right)\!\left(\int_0^R \rho^3\,d\rho\right) = k\cdot 2\pi\cdot\big[-\cos\phi\big]_0^{\pi/2}\cdot\frac{R^4}{4}.$$

Now $[-\cos\phi]_0^{\pi/2} = -\cos\tfrac\pi2 - (-\cos 0) = 0 + 1 = 1$, so

$$M = k\cdot 2\pi\cdot 1\cdot\frac{R^4}{4} = \frac{\pi k R^4}{2}.$$

Note the $R^4$ rather than $R^3$: the extra power comes from the density growing linearly with radius, piling more mass into the outer shells. Dimensional sanity: $k$ carries units of density-per-length, so $k\cdot R^4$ has units of (density)$\times$(length$^3$) $=$ mass. ✓

Add to Your Modeling Portfolio. Add a multiple-integral accumulation to your model — a total quantity gathered over a 2D or 3D region. Biology: model the mass of an organism (or total biomass of a habitat) as $\iiint_E \rho\,dV$, integrating a tissue- or population-density over the body or territory; report the center of mass as a balance point. Economics: compute an aggregate over a continuum of agents, e.g. total demand $\iint_R q(p,w)\,dA$ integrated over a 2D distribution of prices $p$ and wealth levels $w$. Physics: compute the moment of inertia $I_z = \iiint_E (x^2+y^2)\rho\,dV$ of a continuous body about its spin axis, and connect it to rotational kinetic energy $\tfrac12 I_z\omega^2$. Data Science: verify that a 2D joint density integrates to $1$ over $\mathbb{R}^2$, then compute a marginal $p_X(x) = \int p(x,y)\,dy$ and a region probability $P((X,Y)\in R) = \iint_R p\,dA$ — the foundation of every joint model you will build.

Looking Ahead

We have leaned on a magic factor at every coordinate change — $r$ for polar and cylindrical, $\rho^2\sin\phi$ for spherical — and justified each by a geometric area-or-volume argument. Chapter 33 makes this rigorous and general. The Jacobian determinant of any change of variables measures exactly how that transformation stretches or shrinks area (in 2D) or volume (in 3D) locally, and the general change-of-variables formula reads $\iint_R f\,dx\,dy = \iint_S f\big(x(u,v),y(u,v)\big)\,|J|\,du\,dv$. You will see that $r$ and $\rho^2\sin\phi$ are simply the Jacobians of the polar and spherical maps — special cases of one master formula, the multivariable analog of $u$-substitution from Chapter 15.

Beyond that, Part VII turns integration loose on vector fields: line integrals (Chapter 35), surface integrals (Chapter 36), and the great theorems of Green, Stokes, and Gauss (Chapters 35, 37) that, as Chapter 14 foretold, are all higher-dimensional Fundamental Theorems of Calculus — each one trading an integral over a region for information on its boundary.

Reflection

You began this book accumulating along a line and have now learned to accumulate over a plane and through a solid. The leap is conceptual as much as technical: a multiple integral is still "sum of value times patch-size, in the limit," but the patches are now areas and volumes, and the right choice of patch — rectangular, polar, cylindrical, spherical — can turn an intractable integral into a one-line computation. That choice is where geometry pays for itself. And in the Gaussian integral $\int_{-\infty}^\infty e^{-x^2}\,dx = \sqrt\pi$, you saw the deepest lesson of the chapter: sometimes the only way to solve a problem is to make it bigger — to lift it into a higher dimension where its hidden symmetry can finally do the work.