16. Double integrals over rectangles and general regions
§4.2 built the definite integral as a limit of rectangle sums approximating the area under a curve. This module runs the identical construction one dimension up: instead of rectangles under a curve, use thin boxes under a surface z=f(x,y) — and the sum's limit computes a volume, the direct two-variable sibling of everything Module 4 built.
The double integral over a rectangle
For f(x,y)\ge0 defined on a rectangle R=[a,b]\times[c,d], chop R into a grid of small rectangles, and over each one build a thin box of height f(x_i,y_j). Summing the volumes of every box and taking the limit as the grid gets finer:
\iint_Rf(x,y)\,dA=\lim_{m,n\to\infty}\sum_{i=1}^m\sum_{j=1}^nf(x_i,y_j)\,\Delta x\,\Delta y
— exactly §4.1's Riemann sum, with a second index and a second dimension of chopping added. When f\ge0, this computes the volume between the surface z=f(x,y) and the xy-plane over R, the direct generalization of §4.2's "definite integral equals area."
Fubini's theorem: reducing to two ordinary integrals
Computing a double sum's limit directly is impractical. Fubini's theorem says the double integral over a rectangle can always be computed as an iterated integral — two ordinary single-variable integrals, done one after another:
\iint_Rf(x,y)\,dA=\int_a^b\left(\int_c^df(x,y)\,dy\right)dx=\int_c^d\left(\int_a^bf(x,y)\,dx\right)dy
and the order never matters, for any continuous f. The inner integral treats the outer variable as a constant — exactly §10.2's partial derivative logic, run in reverse: integrate with respect to y while x is temporarily frozen, producing an ordinary function of x alone, then integrate that result over x.
This is the entire practical content of the lesson: every double integral this module computes reduces to two applications of Module 4's techniques (substitution, integration by parts, whatever the integrand demands), never anything conceptually new.
Double integrals over general regions
Most regions of interest aren't rectangles. A Type I region is bounded above and below by functions of x:
R=\{(x,y):a\le x\le b,\ g_1(x)\le y\le g_2(x)\}
\iint_Rf(x,y)\,dA=\int_a^b\int_{g_1(x)}^{g_2(x)}f(x,y)\,dy\,dx
This is exactly §5.0's "top minus bottom" boundary description, repurposed: the inner integral's limits are the same g_1(x) and g_2(x) that described a region's vertical extent back when only its area was being computed — now the height being integrated is f(x,y) instead of the constant function 1.
A Type II region is the mirror image, bounded left and right by functions of y:
R=\{(x,y):c\le y\le d,\ h_1(y)\le x\le h_2(y)\}\qquad \iint_Rf(x,y)\,dA=\int_c^d\int_{h_1(y)}^{h_2(y)}f(x,y)\,dx\,dy
— the same "integrate against whichever variable makes the boundaries single-valued" choice §5.0 already made for area alone.
Setting f(x,y)=1 recovers ordinary area: \iint_R1\,dA=\text{Area}(R) — exactly §5.0's formula, since integrating a constant height of 1 over a region just accumulates the region's own area, with nothing extra contributed by a nontrivial integrand.
Doing it in Python
Confirming Fubini's theorem directly — both orders of integration over a rectangle giving the same answer:
import sympy as sp
x, y = sp.symbols('x y')
f = x + y
order1 = sp.integrate(sp.integrate(f, (y, 0, 3)), (x, 0, 2))
order2 = sp.integrate(sp.integrate(f, (x, 0, 2)), (y, 0, 3))
print(f"integrate y first, then x: {order1}")
print(f"integrate x first, then y: {order2}")
print(f"match: {order1 == order2}")
A double integral over a Type I region — between y=x^2 and y=x on [0,1]:
import sympy as sp
x, y = sp.symbols('x y')
f = x * y
result = sp.integrate(sp.integrate(f, (y, x**2, x)), (x, 0, 1))
print(f"integral of x*y over the region between y=x^2 and y=x: {result}")
Recovering ordinary area by integrating the constant function 1 over the same region, and confirming it matches §5.0's direct area-between-curves formula:
import sympy as sp
x, y = sp.symbols('x y')
area_via_double_integral = sp.integrate(sp.integrate(1, (y, x**2, x)), (x, 0, 1))
area_via_section_5 = sp.integrate(x - x**2, (x, 0, 1)) # top minus bottom, section 5.0's formula
print(f"area via double integral: {area_via_double_integral}")
print(f"area via section 5.0: {area_via_section_5}")
print(f"match: {area_via_double_integral == area_via_section_5}")
Worked example
Evaluate \displaystyle\iint_Rf(x,y)\,dA for f(x,y)=x+y over the rectangle R=[0,2]\times[0,3], using both orders of integration.
Integrating y first:
\int_0^3(x+y)\,dy=\left[xy+\frac{y^2}2\right]_0^3=3x+\frac92
\int_0^2\left(3x+\frac92\right)dx=\left[\frac{3x^2}2+\frac{9x}2\right]_0^2=6+9=15
Integrating x first:
\int_0^2(x+y)\,dx=\left[\frac{x^2}2+xy\right]_0^2=2+2y
\int_0^3(2+2y)\,dy=\left[2y+y^2\right]_0^3=6+9=15
\boxed{\iint_Rf(x,y)\,dA=15}
Sanity check. Both orders give exactly 15, confirming Fubini's theorem directly rather than just citing it. As a rough magnitude check: f(x,y)=x+y ranges from 0 (at the corner (0,0)) to 5 (at the corner (2,3)) across a rectangle of area 2\times3=6 — an average height somewhere in the middle of that range, times an area of 6, should land in a plausible ballpark; 15/6=2.5, squarely inside [0,5] and close to what a rough average of the four corner values would suggest. ✓
Your turn
1. Evaluate \displaystyle\iint_R xy\,dA over R=[0,1]\times[0,2].
2. Evaluate \displaystyle\int_0^1\int_0^{x}(2x+y)\,dy\,dx (a Type I region, bounded above by y=x and below by y=0).
3. True or false: Fubini's theorem guarantees the two orders of integration always agree, for any function f(x,y), regardless of continuity.
Solutions
1.
\int_0^2xy\,dy=\left[\frac{xy^2}2\right]_0^2=2x
\int_0^12x\,dx=\Big[x^2\Big]_0^1=\boxed1
2.
\int_0^x(2x+y)\,dy=\left[2xy+\frac{y^2}2\right]_0^x=2x^2+\frac{x^2}2=\frac{5x^2}2
\int_0^1\frac{5x^2}2\,dx=\frac52\left[\frac{x^3}3\right]_0^1=\boxed{\frac56}
3. False. Fubini's theorem requires f to be continuous on the region (or, more generally, integrable in a suitable sense) — a genuine hypothesis. There exist specially constructed discontinuous functions where the two orders of integration produce different answers, or where one order fails to exist as an ordinary integral at all. For every function this course encounters — polynomials, exponentials, trig functions, and their combinations — continuity holds automatically, but the hypothesis is doing real work, not just decoration.
Check yourself in code
Evaluate \displaystyle\iint_Rf(x,y)\,dA for f(x,y)=x+y over R=[0,2]\times[0,3], using both orders of integration.
Print exactly this:
order 1 (y then x) = 15
order 2 (x then y) = 15
import sympy as sp
x, y = sp.symbols('x y')
f = x + y
order1 = sp.integrate(sp.integrate(f, (y, 0, 3)), (x, 0, 2))
print("order 1 (y then x) = ...")
order2 = sp.integrate(sp.integrate(f, (x, 0, 2)), (y, 0, 3))
print("order 2 (x then y) = ...")
import sympy as sp
x, y = sp.symbols('x y')
f = x + y
order1 = sp.integrate(sp.integrate(f, (y, 0, 3)), (x, 0, 2))
print(f"order 1 (y then x) = {order1}")
order2 = sp.integrate(sp.integrate(f, (x, 0, 2)), (y, 0, 3))
print(f"order 2 (x then y) = {order2}")
A double integral \iint_Rf(x,y)\,dA is a limit of box-volume sums, computing the volume between a surface and the plane — and Fubini's theorem reduces that limit to two ordinary single-variable integrals, done one after another in either order, never anything beyond Module 4's existing toolkit. Type I and Type II regions extend §5.0's "top-minus-bottom" boundary description from area alone to a full integrand, and setting f=1 collapses the whole machinery back down to §5.0's original area formula.
Next: double integrals over regions that are naturally circular — where switching to polar coordinates turns an awkward boundary into a clean one.