Double integrals

A surface $z = f (x, y)$ standing over a region of the plane encloses a volume, and nothing in the calculus of one variable says what that volume is or how to compute it.

The previous lesson finished the differential half of the subject: partial derivatives, gradients, tangent planes and extrema, with and without constraints. The other half of calculus is summing. In one variable the integral $\int_{a}^{b} f (x) d x$ was defined as the limit of sums of thin rectangles, height times width, as the widths shrink to zero, and the Fundamental Theorem then evaluated it with an antiderivative. This lesson builds the two variable integral the same way, from sums of thin boxes, and then shows that it can be evaluated with nothing more than two single integrals, one inside the other. It assumes the single integral with its substitution rule, and the partial derivative's habit of holding one variable fixed while the other moves.

Volume as a limit of sums

Take a function $f (x, y)$ on the rectangle $R$ where $a \leq x \leq b$ and $c \leq y \leq d$ , written $R = [a, b] \times [c, d]$ , and suppose for the moment that $f$ is positive. Cut $[a, b]$ into $m$ equal pieces of width $Δ x$ and $[c, d]$ into $n$ pieces of width $Δ y$ . The grid lines divide $R$ into $m n$ small rectangles, each of area $Δ A = Δ x Δ y$ . Over the small rectangle $R_{i j}$ pick a sample point $(x_{i j}^{*}, y_{i j}^{*})$ and erect a box of that base and height $f (x_{i j}^{*}, y_{i j}^{*})$ . The boxes together approximate the solid under the surface, and their total volume is the double Riemann sum

S_{m n} = \sum_{i = 1}^{m} \sum_{j = 1}^{n} f (x_{i j}^{*}, y_{i j}^{*}) Δ A

The double integral of $f$ over $R$ is the limit of these sums as the grid is refined in both directions, provided the limit exists and is the same for every choice of sample points:

\iint_{R} f d A = lim_{m, n \to \infty} \sum_{i = 1}^{m} \sum_{j = 1}^{n} f (x_{i j}^{*}, y_{i j}^{*}) Δ A

When $f \geq 0$ this is the volume under the surface and above $R$ . When $f$ takes both signs, the boxes below the plane count negatively, exactly as area below the axis did in one variable. As there, every continuous function on a closed rectangle has a double integral; the proof uses the same uniform continuity argument and is not repeated here.

Take the paraboloid $f (x, y) = 16 - x^{2} - 2 y^{2}$ over the square $[0, 2] \times [0, 2]$ , and use an $n$ by $n$ grid. Two sample points are natural. The upper right corner of each small square gives the lowest value of $f$ on it, because $f$ falls as $x$ and $y$ grow, so that sum underestimates. The centre of each square gives the midpoint sum. A few lines of code produce these:

Grid	Rectangles	Corner sum	Midpoint sum
$2 \times 2$	4	34.000	49.000
$4 \times 4$	16	41.500	48.250
$8 \times 8$	64	44.875	48.0625
$16 \times 16$	256	46.469	48.0156
$100 \times 100$	10,000	47.759	48.0004
$1000 \times 1000$	1,000,000	47.976	48.000004

Both columns close on $48$ . The corner sums creep up slowly, their error falling roughly like $\frac{1}{n}$ , and even a million boxes leave them wrong in the second decimal place. The midpoint sums miss by exactly $\frac{4}{n^{2}}$ , so doubling the grid quarters the error. What the table cannot do is prove the answer is $48$ , and a million evaluations for three figures is poor value. The next section finds the exact value with two ordinary antiderivatives.

Slicing: the iterated integral

Freeze $x$ at some value between $0$ and $2$ and cut the solid with the vertical plane at that $x$ . The cut face is the region under the curve $y \mapsto 16 - x^{2} - 2 y^{2}$ for $0 \leq y \leq 2$ , an ordinary area under a curve, with $x$ acting as a constant. Its area is

A (x) = \int_{0}^{2} (16 - x^{2} - 2 y^{2}) d y = {[16 y - x^{2} y - \frac{2 y^{3}}{3}]}_{0}^{2} = 32 - 2 x^{2} - \frac{16}{3}

This is partial integration, the reverse of the partial derivative: integrate in $y$ as if $x$ were a number. The result is a function of $x$ alone.

Now stack the slices. A slab of thickness $Δ x$ at position $x$ has volume close to $A (x) Δ x$ , and adding the slabs and letting $Δ x$ shrink is a single integral:

\int_{0}^{2} A (x) d x = \int_{0}^{2} (\frac{80}{3} - 2 x^{2}) d x = \frac{160}{3} - \frac{16}{3} = 48

The table's $48$ is exact. The integral just computed is written with both steps visible,

\int_{0}^{2} \int_{0}^{2} (16 - x^{2} - 2 y^{2}) d y d x

and called an iterated integral. It is read from the inside out: the inner integral and the inner differential $d y$ belong together, so $y$ runs from $0$ to $2$ first, and the outer pair, $d x$ with its limits, is done last. Slicing the other way, at fixed $y$ , gives $\int_{0}^{2} (16 - x^{2} - 2 y^{2}) d x = \frac{88}{3} - 4 y^{2}$ , and then $\int_{0}^{2} (\frac{88}{3} - 4 y^{2}) d y = \frac{176}{3} - \frac{32}{3} = 48$ again.

The slicing argument is the method of Bonaventura Cavalieri, who in 1635 compared solids by comparing their cross sections. It is persuasive, but it quietly assumes that the volume defined by boxes and the volume assembled from slices are the same number. That needs a theorem.

Fubini's theorem

The theorem is named after Guido Fubini, who proved a general form of it in 1907. For continuous functions on a rectangle it says:

If $f$ is continuous on $R = [a, b] \times [c, d]$ , then $\iint_{R} f d A = \int_{a}^{b} \int_{c}^{d} f (x, y) d y d x = \int_{c}^{d} \int_{a}^{b} f (x, y) d x d y$

The reason is visible in the Riemann sum. Choose the sample points on a grid, $(x_{i}, y_{j})$ , and add the boxes column by column rather than all at once. The $i$ -th column contributes $(\sum_{j} f (x_{i}, y_{j}) Δ y) Δ x$ , and the bracket is a Riemann sum for the slice area $A (x_{i}) = \int_{c}^{d} f (x_{i}, y) d y$ . Adding the columns then gives a Riemann sum for $\int_{a}^{b} A (x) d x$ . Adding row by row instead gives the other order. A finite sum can be added in any order, so the only work in the proof is showing that both approximations become good at the same time as the grid is refined, and for a continuous function uniform continuity supplies exactly that.

The hypothesis is not decoration. On the square $[0, 1] \times [0, 1]$ take $f (x, y) = \frac{x^{2} - y^{2}}{(x^{2} + y^{2})^{2}}$ , which is continuous everywhere except the corner $(0, 0)$ , where it blows up. Because $\frac{\partial}{\partial y} \frac{y}{x^{2} + y^{2}} = \frac{x^{2} - y^{2}}{(x^{2} + y^{2})^{2}}$ , the inner integral in $y$ is $\frac{1}{1 + x^{2}}$ , and then $\int_{0}^{1} \frac{1}{1 + x^{2}} d x = \arctan 1 = \frac{π}{4}$ . Swapping the roles of $x$ and $y$ changes the sign of $f$ , so the other order gives $- \frac{π}{4}$ . Two orders, two answers, and neither is the double integral, which does not exist: the function is unbounded, and the Riemann sums near the corner can be made as large as desired of either sign.

For continuous functions, the practical content of the theorem is freedom. The order of integration can be chosen to make the antiderivatives easy.

Example. Evaluate $\iint_{R} x e^{x y} d A$ over $R = [0, 1] \times [0, 2]$ .

Integrating first in $x$ needs integration by parts and leaves a function of $y$ with $\frac{1}{y^{2}}$ terms that is awkward to integrate next. Integrating first in $y$ is immediate, because $x$ is exactly the factor the chain rule produces: $\frac{\partial}{\partial y} e^{x y} = x e^{x y}$ . So

\int_{0}^{2} x e^{x y} d y = {[e^{x y}]}_{y = 0}^{y = 2} = e^{2 x} - 1

and then

\int_{0}^{1} (e^{2 x} - 1) d x = \frac{e^{2} - 1}{2} - 1 = \frac{e^{2} - 3}{2} \approx 2.19453

A midpoint sum on a $1000 \times 1000$ grid gives $2.194527$ , agreeing to six figures.

Now you. Evaluate $\iint_{R} y e^{x y} d A$ over $R = [0, 3] \times [0, 1]$ , to three decimal places. Choose the order first.

Answer

Integrate in $x$ first, since $y$ is the factor $\frac{\partial}{\partial x} e^{x y}$ produces: $\int_{0}^{3} y e^{x y} d x = e^{3 y} - 1$ . Then $\int_{0}^{1} (e^{3 y} - 1) d y = \frac{e^{3} - 1}{3} - 1 = \frac{e^{3} - 4}{3} \approx 5.362$ . A midpoint sum on a $1000 \times 1000$ grid gives $5.36184$ .

Regions between two curves

Few regions of interest are rectangles. For a bounded region $D$ , enclose it in a rectangle $R$ and define $F$ to equal $f$ on $D$ and $0$ outside it. Then $\iint_{D} f d A$ is defined as $\iint_{R} F d A$ : the boxes outside $D$ contribute nothing. $F$ jumps at the edge of $D$ , so it is not continuous, but when the edge is made of graphs of continuous functions the jump does no harm. The double integral still exists, and a stronger form of Fubini's theorem still turns it into an iterated integral.

Two shapes of region cover most cases. A region of the first type lies between two curves in $x$ : it is described by $a \leq x \leq b$ and $g_{1} (x) \leq y \leq g_{2} (x)$ . On the vertical line at a fixed $x$ , $F$ is zero except between $g_{1} (x)$ and $g_{2} (x)$ , so the inner integral only runs over that stretch:

\iint_{D} f d A = \int_{a}^{b} \int_{g_{1} (x)}^{g_{2} (x)} f (x, y) d y d x

A region of the second type lies between two curves in $y$ , described by $c \leq y \leq d$ and $h_{1} (y) \leq x \leq h_{2} (y)$ , and its integral is $\int_{c}^{d} \int_{h_{1} (y)}^{h_{2} (y)} f d x d y$ . In both, the outer limits are constants and the inner limits may depend on the outer variable, never the other way round, because the outer variable is the one held fixed while the inner one sweeps.

The limits come from a description of the region, which is usually given in words. Take "the region enclosed by the parabola $y = x^{2}$ and the line $y = 2 x$ ". The curves meet where $x^{2} = 2 x$ , at $x = 0$ and $x = 2$ , so $x$ runs from $0$ to $2$ . Between those values one curve lies above the other throughout, and testing $x = 1$ shows which: the line gives $2$ , the parabola $1$ . So the region is $0 \leq x \leq 2$ , $x^{2} \leq y \leq 2 x$ . It is also of the second type: solving each curve for $x$ gives $x = \frac{y}{2}$ on the line and $x = \sqrt{y}$ on the parabola, and at a height $y$ between $0$ and $4$ the line is on the left, so $\frac{y}{2} \leq x \leq \sqrt{y}$ .

Example. Evaluate $\iint_{D} x y d A$ over the region enclosed by $y = x^{2}$ and $y = 2 x$ .

As a region of the first type, the inner integral holds $x$ fixed:

\int_{x^{2}}^{2 x} x y d y = x {[\frac{y^{2}}{2}]}_{x^{2}}^{2 x} = 2 x^{3} - \frac{x^{5}}{2}

and then $\int_{0}^{2} (2 x^{3} - \frac{x^{5}}{2}) d x = 8 - \frac{64}{12} = \frac{8}{3}$ . As a check, the second type gives $\int_{0}^{4} y \cdot \frac{1}{2} (y - \frac{y^{2}}{4}) d y = \frac{64}{6} - \frac{256}{32} = \frac{8}{3}$ as well. A midpoint sum with $1000$ columns, each cut into $1000$ pieces between the curves, gives $2.666669$ .

Now you. Evaluate $\iint_{D} y d A$ over the region enclosed by the parabola $x = y^{2}$ and the line $x = y + 2$ . Which type makes it one integral rather than two?

Answer

The curves meet where $y^{2} = y + 2$ , at $y = - 1$ and $y = 2$ , and at $y = 0$ the line ( $x = 2$ ) is to the right of the parabola ( $x = 0$ ). So the region is $- 1 \leq y \leq 2$ , $y^{2} \leq x \leq y + 2$ , of the second type, and $\int_{- 1}^{2} y (y + 2 - y^{2}) d y = {[\frac{y^{3}}{3} + y^{2} - \frac{y^{4}}{4}]}_{- 1}^{2} = \frac{8}{3} - \frac{5}{12} = \frac{9}{4}$ . As the first type the lower edge changes formula at $x = 1$ , from $y = - \sqrt{x}$ to $y = x - 2$ , forcing two integrals.

Reversing the order

Sometimes the order a problem arrives in cannot be done at all. Consider

\int_{0}^{1} \int_{y}^{1} e^{x^{2}} d x d y

The inner integral asks for an antiderivative of $e^{x^{2}}$ , and there is none among the elementary functions: Joseph Liouville's theory of integration in finite terms, developed in the 1830s, proves that no combination of powers, exponentials, logarithms and trigonometric functions differentiates to it.

Reversing the order means reading the region off the limits and describing it the other way. The limits say $0 \leq y \leq 1$ and $y \leq x \leq 1$ : at each height $y$ , $x$ runs from the diagonal $x = y$ to the vertical line $x = 1$ . That is the triangle with corners $(0, 0)$ , $(1, 0)$ and $(1, 1)$ . Described as a region of the first type, $x$ runs from $0$ to $1$ , and on the vertical line at $x$ the triangle stretches from the bottom edge $y = 0$ up to the diagonal $y = x$ . So

\int_{0}^{1} \int_{y}^{1} e^{x^{2}} d x d y = \int_{0}^{1} \int_{0}^{x} e^{x^{2}} d y d x = \int_{0}^{1} x e^{x^{2}} d x

because $e^{x^{2}}$ is constant in $y$ and the inner integral just multiplies it by the length $x$ . That factor $x$ is exactly what the substitution $u = x^{2}$ needs, and

\int_{0}^{1} x e^{x^{2}} d x = {[\frac{e^{x^{2}}}{2}]}_{0}^{1} = \frac{e - 1}{2} \approx 0.859141

The original order can still be computed numerically, since only the antiderivative is missing, not the integral. A midpoint sum in the order given, $2000$ values of $y$ and for each one $2000$ points from $x = y$ to $x = 1$ , gives $0.8591409$ , against the exact $0.8591409$ . Fubini's theorem turned an integral that has no formula in one order into a substitution in the other.

Area, average value and mass

Integrating the constant $1$ over $D$ builds boxes of height $1$ , whose volume equals the area of their base, so

area of D = \iint_{D} 1 d A

For the region between $x = y^{2}$ and $x = y + 2$ this is $\int_{- 1}^{2} (y + 2 - y^{2}) d y = \frac{9}{2}$ , the familiar area between two curves recovered as a special case.

The average value of $f$ over $D$ is defined as in one variable, the total divided by the size: $\overline{f} = \frac{1}{area of D} \iint_{D} f d A$ . It is the height of the flat-topped solid on the same base with the same volume. For the paraboloid over $[0, 2] \times [0, 2]$ it is $\frac{48}{4} = 12$ , between the value $16$ at the corner $(0, 0)$ and $4$ at the corner $(2, 2)$ .

A thin plate whose density $ρ (x, y)$ , in kilograms per square metre, varies from point to point has mass by the same summing: a small rectangle near $(x, y)$ has mass about $ρ (x, y) Δ A$ , and in the limit $m = \iint_{D} ρ d A$ . Its average density is the mass divided by the area.

Example. A plate occupies the region between the parabola $y = x^{2}$ and the line $y = 1$ , with lengths in metres and density $ρ = 1 + y$ kilograms per square metre. Find its area, mass and average density.

The curves meet at $x = \pm 1$ , and the line is on top, so $- 1 \leq x \leq 1$ , $x^{2} \leq y \leq 1$ . The area is $\int_{- 1}^{1} (1 - x^{2}) d x = \frac{4}{3}$ square metres. For the mass, $\int_{x^{2}}^{1} (1 + y) d y = \frac{3}{2} - x^{2} - \frac{x^{4}}{2}$ , and $\int_{- 1}^{1} (\frac{3}{2} - x^{2} - \frac{x^{4}}{2}) d x = 3 - \frac{2}{3} - \frac{1}{5} = \frac{32}{15}$ , about $2.1333$ kg. The average density is $\frac{32}{15} \div \frac{4}{3} = 1.6$ kilograms per square metre, above the plain average of $1 + y$ over $0 \leq y \leq 1$ (which is $1.5$ ) because the plate is wider near the top, where it is denser. A midpoint sum on $1000$ columns of $1000$ pieces gives $2.13333$ .

Now you. A plate occupies the region under $y = \sin x$ and above the $x$ -axis for $0 \leq x \leq π$ , with density $ρ = y$ kilograms per square metre. Find the mass to four decimal places, and the average density.

Answer

The inner integral is $\int_{0}^{\sin x} y d y = \frac{\sin^{2} x}{2}$ , and $\int_{0}^{π} \frac{\sin^{2} x}{2} d x = \frac{π}{4} \approx 0.7854$ kg, using $\sin^{2} x = \frac{1 - \cos 2 x}{2}$ . The area is $\int_{0}^{π} \sin x d x = 2$ , so the average density is $\frac{π}{8} \approx 0.3927$ kilograms per square metre.

A disc in the wrong coordinates

Every region so far has had edges that were easy to write as $y$ in terms of $x$ or $x$ in terms of $y$ . The unit disc, $x^{2} + y^{2} \leq 1$ , is of both types, and in both it is awkward: $- 1 \leq x \leq 1$ and $- \sqrt{1 - x^{2}} \leq y \leq \sqrt{1 - x^{2}}$ . Even its area, $\int_{- 1}^{1} 2 \sqrt{1 - x^{2}} d x$ , needs a trigonometric substitution to reach $π$ . For $\iint e^{- (x^{2} + y^{2})} d A$ over the disc, the inner integral is $\int e^{- y^{2}} d y$ between square roots, which has no elementary antiderivative, and here reversing the order does not help, since the disc looks the same either way.

The problem is not the function but the grid. Rectangles in $x$ and $y$ are the wrong pieces for a round region; a disc is naturally described by a distance from the centre and an angle. The next lesson cuts the plane along circles and rays instead, works out the area of the small pieces those cuts produce, and with them evaluates this integral and one that single variable calculus could not reach.