Fubini's Theorem

Suppose you want to compute the volume under the surface z=x2y+sin⁡(xy)z = x^2 y + \sin(xy) over a rectangle. You have two natural strategies. First, slice the solid with planes perpendicular to the xx-axis. Integrate each slice's area, then add the slices up. Alternatively, do the same with planes perpendicular to the yy-axis. Do these two computations always agree? Do they always equal the "true" double integral, defined directly as a limit of Riemann sums over the rectangle? This question is not a formality. The answer requires proof. That proof is exactly what makes double integrals computable in practice. Fubini's Theorem, in the continuous setting treated here, applies to a continuous function on a rectangle. It says the double integral can be evaluated by iterating two single-variable integrals. Either order gives the same result. This single fact bridges two things. On one side is the abstract definition of a multiple integral, built from Riemann sums over a two-dimensional partition. On the other is the concrete machinery of one-variable calculus you already know: antiderivatives, the Fundamental Theorem of Calculus, substitution. Without this theorem, multivariable calculus would be a subject of existence proofs with no computational teeth. With it, every double or triple integral of a continuous function reduces to a chain of ordinary integrals. These can be computed symbolically or numerically, one variable at a time. This theorem also justifies swapping the order of integration. That move is used constantly to simplify a hard integral into an easy one. It also helps make sense of integrals over regions with complicated boundaries, once you generalize beyond rectangles. This article proves the continuous rectangular case in full. It builds the geometric picture of slicing, and works two asymmetric examples. It also shows exactly what breaks (and why) when continuity is dropped. It sets up the natural next steps: Integrals over General Regions and the Change of Variables and the Jacobian formula for transforming multiple integrals.

Key takeaways
  • Fubini's Theorem (continuous, rectangular case): for ff continuous on a rectangle R=[a,b]×[c,d]R = [a,b]\times[c,d], the double integral equals either iterated integral.
  • It converts a 2-D limit-of-Riemann-sums object into two ordinary 1-D integrals done in sequence.
  • The proof rests on uniform continuity of ff on the compact rectangle RR, which forces the "slice function" A(x)=∫cdf(x,y) dyA(x) = \int_c^d f(x,y)\,dy to be continuous in xx.
  • Order of integration can be swapped freely for continuous ff — this is the practical payoff used in every computation.
  • Continuity is not just convenient — the counterexample shows a bounded, discontinuous ff for which the two iterated integrals disagree.
  • Generalizes to triple integrals over boxes, and (with more care) to integrals over general bounded regions.
  • Key idea: slicing a solid one way or the other and integrating the slice-areas gives the same volume, because continuity lets you interchange a limit with an integral.

Where this lives today

Fubini's Theorem is the computational backbone of every multivariable calculus computation of area, volume, mass, center of mass, and moment of inertia via double or triple integrals — anywhere a physical quantity is built from an integrand over a two- or three-dimensional region. In probability and statistics, the Lebesgue version of Fubini's Theorem justifies computing expectations of functions of independent random variables by iterated integration, and underlies the derivation of convolution formulas for sums of independent random variables. In Fourier analysis and PDE theory, swapping the order of integration (justified by Fubini/Tonelli) is a routine step in deriving convolution theorems, Green's function representations, and solving heat and wave equations via integral transforms. In numerical analysis and computer graphics, Fubini's Theorem underlies "separable" numerical quadrature schemes that approximate a double integral by nested one-dimensional quadrature rules, and it justifies rendering algorithms (ray marching, volumetric integration) that accumulate light contributions slice by slice. In cryptography and computational geometry, volume and probability computations over high-dimensional boxes routinely reduce dimension-by-dimension using the same iterated-integral principle, now packaged in the more general measure-theoretic Fubini–Tonelli framework that extends this rectangle case to the far broader settings needed there.

:::quiz[The counterexample f(x,y)=(x2−y2)/(x2+y2)2f(x,y) = (x^2-y^2)/(x^2+y^2)^2 on [0,1]2[0,1]^2 shows that:]

  • Fubini's Theorem is false even for continuous functions
  • The iterated integrals can disagree (π/4\pi/4 vs −π/4-\pi/4) when continuity fails at a point due to unboundedness ✓
  • Iterated integrals are always independent of order regardless of continuity
  • The double integral always exists whenever both iterated integrals do

:::quiz[Which step of the proof crucially uses compactness of R=[a,b]×[c,d]R = [a,b]\times[c,d]?]

  • Deriving uniform continuity of ff from (pointwise) continuity, via the Heine–Borel theorem ✓
  • Defining the iterated integral notation
  • Choosing which variable to integrate first
  • Concluding that ff is a polynomial

:::quiz[Example 2 in the article computed ∬Rxexy dA\iint_R xe^{xy}\,dA over [0,1]×[0,2][0,1]\times[0,2]. Integrating yy first gave a clean closed form because:]

  • xx was treated as the variable of integration throughout
  • Fixing xx turns xexyxe^{xy} into a simple exponential in yy whose antiderivative w.r.t. yy is directly exye^{xy} ✓
  • The integrand was independent of yy
  • Integration by parts was unavoidable in both orders

:::faq[Can Fubini's Theorem be applied to an unbounded rectangle, like [0,∞)×[0,1][0,\infty)\times[0,1]?] Not directly — the proof given here relies essentially on compactness of RR (to get uniform continuity and existence of the double integral via Darboux sums), which fails for an unbounded region. Improper double integrals over unbounded regions require a separate justification, typically via the Lebesgue-theoretic Fubini–Tonelli Theorem, which handles σ\sigma-finite measure spaces and unbounded domains, provided an absolute-integrability condition ∬∣f∣<∞\iint |f| < \infty holds. Naively swapping integration order on an unbounded domain without checking this condition is a common source of paradoxical "proofs" that 0=10=1 using divergent integrals.

Ask the tutor: How does the argument change if ff is continuous on an open rectangle but only bounded (not continuous) on its boundary?

DefinitionIterated integral and Fubini's Theorem — setup

Let R=[a,b]×[c,d]⊂R2R = [a,b] \times [c,d] \subset \mathbb{R}^2 be a closed rectangle and let f:R→Rf : R \to \mathbb{R} be a bounded function. Recall from Multiple Integrals over Rectangles that the double integral ∬Rf dA\iint_R f\,dA is defined as the limit of Riemann sums over partitions of RR into subrectangles, when this limit exists (equivalently, when upper and lower Darboux sums converge to the same value); in that case ff is called (Riemann) integrable on RR. For each fixed x∈[a,b]x \in [a,b], suppose the single-variable function y↦f(x,y)y \mapsto f(x,y) is integrable on [c,d][c,d], and define

A(x)=∫cdf(x,y) dy.A(x) = \int_c^d f(x,y)\,dy .

The iterated integral of ff over RR, integrating yy first, is

∫abA(x) dx=∫ab(∫cdf(x,y) dy)dx,\int_a^b A(x)\,dx = \int_a^b \left( \int_c^d f(x,y)\,dy \right) dx,

written more compactly as ∫ab∫cdf(x,y) dy dx\int_a^b \int_c^d f(x,y)\,dy\,dx. The iterated integral in the other order, ∫cd∫abf(x,y) dx dy\int_c^d \int_a^b f(x,y)\,dx\,dy, is defined symmetrically.

Fubini's Theorem asserts that when ff is continuous on RR, all three quantities — the double integral and both iterated integrals — exist and are equal.

Intuition

Picture the graph of ff over RR as a solid roof sitting above the rectangular floor RR, and imagine computing the volume trapped between the floor and the roof. One way to find this volume is to slice the solid with a plane x=x0x = x_0 for each fixed x0∈[a,b]x_0 \in [a,b]. Each slice is a two-dimensional region in the plane x=x0x = x_0, bounded above by the curve y↦f(x0,y)y \mapsto f(x_0, y); its area is exactly A(x0)=∫cdf(x0,y) dyA(x_0) = \int_c^d f(x_0, y)\, dy. Stacking up these slice-areas as x0x_0 ranges over [a,b][a,b] — literally integrating A(x)A(x) — should recover the total volume, by the same "area of cross-section, integrated" principle used to compute volumes of revolution in single-variable calculus. The subtlety is that this is a claim, not a definition: the double integral ∬Rf dA\iint_R f\, dA was already defined directly, as a two-dimensional limit of Riemann sums, with no reference to slicing at all. Fubini's Theorem is the statement that this direct definition secretly agrees with the slicing computation, and moreover that it does not matter which axis you slice along first. The engine that makes this true is uniform continuity: on the compact rectangle RR, a continuous ff cannot oscillate wildly on very fine subrectangles, so refining a 2-D partition and computing 1-D slice integrals lead to the same limiting number.

:::svg[Slicing the solid two ways over R=[0,2]×[0,1]R=[0,2]\times[0,1]] Slice at fixed x: area A(x) = ∫ f(x,y) dy x y x₀ 0 2 1 Slice at fixed y: area B(y) = ∫ f(x,y) dx x y y₀ 0 2 1 Fubini: ∫₀² A(x) dx = ∬_R f dA = ∫₀¹ B(y) dy

Central theorem

TheoremFubini's Theorem for continuous functions on a rectangle

Let R=[a,b]×[c,d]⊂R2R = [a,b] \times [c,d] \subset \mathbb{R}^2 with a<ba<b, c<dc<d, and let f:R→Rf : R \to \mathbb{R} be continuous on RR. Then: 1. ff is (Riemann) integrable on RR, so ∬Rf dA\iint_R f \,dA exists. 2. For every x∈[a,b]x \in [a,b], the function y↦f(x,y)y \mapsto f(x,y) is continuous on [c,d][c,d], hence integrable, so

A(x):=∫cdf(x,y) dyA(x) := \int_c^d f(x,y)\,dy

is well-defined for all x∈[a,b]x\in[a,b]; moreover AA is continuous on [a,b][a,b], hence ∫abA(x) dx\int_a^b A(x)\,dx exists. 4. Symmetrically, B(y):=∫abf(x,y) dxB(y) := \int_a^b f(x,y)\,dx is well-defined and continuous on [c,d][c,d]. 5. All three integrals coincide:

∬Rf dA  =  ∫ab∫cdf(x,y) dy dx  =  ∫cd∫abf(x,y) dx dy.\iint_R f \, dA \;=\; \int_a^b \int_c^d f(x,y)\,dy\,dx \;=\; \int_c^d \int_a^b f(x,y)\,dx\,dy .
Proof

We use the following proof technique: establish integrability of ff on RR from uniform continuity (a compactness argument), then show the iterated integral matches the double integral by squeezing both between the same upper and lower Darboux sums.

Step 1: ff is integrable on RR. Since R=[a,b]×[c,d]R = [a,b]\times[c,d] is closed and bounded in R2\mathbb{R}^2, it is compact (Heine–Borel). A continuous function on a compact set is uniformly continuous (this is the standard uniform continuity theorem, proved independently of Fubini's Theorem and not relying on it). So given ε>0\varepsilon>0, there is δ>0\delta>0 such that for all (x,y),(x′,y′)∈R(x,y),(x',y')\in R,

∥(x,y)−(x′,y′)∥<δ  ⟹  ∣f(x,y)−f(x′,y′)∣<ε.\|(x,y)-(x',y')\| < \delta \implies |f(x,y)-f(x',y')| < \varepsilon.

Take any partition PP of RR into subrectangles Rij=[xi−1,xi]×[yj−1,yj]R_{ij} = [x_{i-1},x_i]\times[y_{j-1},y_j] (1≤i≤m1\le i \le m, 1≤j≤n1\le j\le n) with mesh (the maximum diagonal length of any RijR_{ij}) less than δ\delta. On each RijR_{ij}, compactness of RijR_{ij} and continuity of ff give a maximum MijM_{ij} and minimum mijm_{ij} of ff on RijR_{ij} (extreme value theorem), attained at some points pij,qij∈Rijp_{ij}, q_{ij} \in R_{ij}. Since ∥pij−qij∥<δ\|p_{ij}-q_{ij}\|<\delta (both lie in RijR_{ij}, whose diagonal is <δ<\delta),

Mij−mij=f(pij)−f(qij)<ε.M_{ij} - m_{ij} = f(p_{ij}) - f(q_{ij}) < \varepsilon.

Hence the upper and lower Darboux sums satisfy

U(f,P)−L(f,P)=∑i,j(Mij−mij) ΔxiΔyj<ε∑i,jΔxiΔyj=ε (b−a)(d−c).U(f,P) - L(f,P) = \sum_{i,j} (M_{ij}-m_{ij})\,\Delta x_i \Delta y_j < \varepsilon \sum_{i,j} \Delta x_i \Delta y_j = \varepsilon\,(b-a)(d-c).

Since ε>0\varepsilon>0 was arbitrary and (b−a)(d−c)(b-a)(d-c) is a fixed constant, the Riemann criterion for integrability is satisfied: for every η>0\eta>0 we can choose ε=η/[(b−a)(d−c)]\varepsilon = \eta/[(b-a)(d-c)] and find a partition with U(f,P)−L(f,P)<ηU(f,P)-L(f,P)<\eta. Therefore ff is integrable on RR, and ∬Rf dA\iint_R f\,dA exists, equal to the common value sup⁡PL(f,P)=inf⁡PU(f,P)\sup_P L(f,P) = \inf_P U(f,P).

Step 2: A(x)=∫cdf(x,y) dyA(x)=\int_c^d f(x,y)\,dy is well-defined and continuous on [a,b][a,b]. Fix x0∈[a,b]x_0 \in [a,b]. The map y↦f(x0,y)y \mapsto f(x_0,y) is continuous on [c,d][c,d] (restriction of a continuous function of two variables to a fixed value of one variable is continuous — composing y↦(x0,y)y \mapsto (x_0,y), itself continuous, with ff). A continuous function on the compact interval [c,d][c,d] is Riemann integrable, so A(x0)A(x_0) is well defined for every x0∈[a,b]x_0\in[a,b]. For continuity of AA: fix ε>0\varepsilon>0. By uniform continuity of ff on RR, choose δ>0\delta>0 so that ∥(x,y)−(x′,y′)∥<δ\|(x,y)-(x',y')\|<\delta implies ∣f(x,y)−f(x′,y′)∣<ε/(d−c)|f(x,y)-f(x',y')| < \varepsilon/(d-c). In particular if ∣x−x′∣<δ|x-x'|<\delta then for every y∈[c,d]y\in[c,d], ∥(x,y)−(x′,y)∥=∣x−x′∣<δ\|(x,y)-(x',y)\| = |x-x'| < \delta, so ∣f(x,y)−f(x′,y)∣<ε/(d−c)|f(x,y)-f(x',y)| < \varepsilon/(d-c). Then

∣A(x)−A(x′)∣=∣∫cd(f(x,y)−f(x′,y)) dy∣≤∫cd∣f(x,y)−f(x′,y)∣ dy<(d−c)⋅εd−c=ε.|A(x)-A(x')| = \left| \int_c^d \big(f(x,y)-f(x',y)\big)\,dy \right| \le \int_c^d |f(x,y)-f(x',y)|\,dy < (d-c)\cdot \frac{\varepsilon}{d-c} = \varepsilon.

So AA is (uniformly) continuous on [a,b][a,b], hence Riemann integrable there, and ∫abA(x) dx\int_a^b A(x)\,dx exists. The argument for B(y)B(y) is identical with the roles of x,yx,y exchanged.

Step 3: The iterated integral equals the double integral. Let ε>0\varepsilon>0. By Step 1's uniform continuity argument, choose a partition PP of RR (with subintervals a=x0<x1<⋯<xm=ba=x_0<x_1<\dots<x_m=b of [a,b][a,b] and c=y0<y1<⋯<yn=dc=y_0<y_1<\dots<y_n=d of [c,d][c,d]) fine enough that

U(f,P)−L(f,P)<ε.U(f,P) - L(f,P) < \varepsilon.

On each subrectangle Rij=[xi−1,xi]×[yj−1,yj]R_{ij}=[x_{i-1},x_i]\times[y_{j-1},y_j], let

mij=inf⁡Rijf,Mij=sup⁡Rijf.m_{ij} = \inf_{R_{ij}} f, \qquad M_{ij} = \sup_{R_{ij}} f.

Fix ii. For any x∈[xi−1,xi]x\in[x_{i-1},x_i] and any y∈[yj−1,yj]y\in[y_{j-1},y_j], we have mij≤f(x,y)≤Mijm_{ij} \le f(x,y) \le M_{ij}. Integrating this inequality over y∈[yj−1,yj]y\in[y_{j-1},y_j] (with xx fixed) gives

mij Δyj  ≤  ∫yj−1yjf(x,y) dy  ≤  Mij Δyj.m_{ij}\,\Delta y_j \;\le\; \int_{y_{j-1}}^{y_j} f(x,y)\,dy \;\le\; M_{ij}\,\Delta y_j .

Summing over j=1,…,nj=1,\dots,n and using additivity of the integral over adjacent intervals,

∑j=1nmij Δyj  ≤  ∫cdf(x,y) dy=A(x)  ≤  ∑j=1nMij Δyj,for every x∈[xi−1,xi].\sum_{j=1}^n m_{ij}\,\Delta y_j \;\le\; \int_c^d f(x,y)\,dy = A(x) \;\le\; \sum_{j=1}^n M_{ij}\,\Delta y_j, \qquad \text{for every } x \in [x_{i-1},x_i].

This holds for every xx in [xi−1,xi][x_{i-1},x_i], so in particular

∑j=1nmij Δyj  ≤  inf⁡[xi−1,xi]Aandsup⁡[xi−1,xi]A  ≤  ∑j=1nMij Δyj.\sum_{j=1}^n m_{ij}\,\Delta y_j \;\le\; \inf_{[x_{i-1},x_i]} A \quad\text{and}\quad \sup_{[x_{i-1},x_i]} A \;\le\; \sum_{j=1}^n M_{ij}\,\Delta y_j.

Multiply by Δxi\Delta x_i and sum over i=1,…,mi=1,\dots,m:

∑i,jmij ΔxiΔyj  ≤  ∑i(inf⁡[xi−1,xi]A)Δxi  ≤  ∑i(sup⁡[xi−1,xi]A)Δxi  ≤  ∑i,jMij ΔxiΔyj.\sum_{i,j} m_{ij}\,\Delta x_i \Delta y_j \;\le\; \sum_i \Big(\inf_{[x_{i-1},x_i]} A\Big)\Delta x_i \;\le\; \sum_i\Big(\sup_{[x_{i-1},x_i]}A\Big)\Delta x_i \;\le\; \sum_{i,j} M_{ij}\,\Delta x_i \Delta y_j .

The left and right ends are exactly L(f,P)L(f,P) and U(f,P)U(f,P). The middle terms are L(A,Px)L(A, P_x) and U(A,Px)U(A,P_x), the lower and upper Darboux sums of the single-variable function AA with respect to the partition Px={x0,…,xm}P_x = \{x_0,\dots,x_m\} of [a,b][a,b]. So we have shown

L(f,P)  ≤  L(A,Px)  ≤  U(A,Px)  ≤  U(f,P).L(f,P) \;\le\; L(A,P_x) \;\le\; U(A,P_x) \;\le\; U(f,P).

Since AA is integrable on [a,b][a,b] (Step 2), ∫abA(x) dx\int_a^b A(x)\,dx lies between L(A,Px)L(A,P_x) and U(A,Px)U(A,P_x) as well:

L(f,P)  ≤  L(A,Px)  ≤  ∫abA(x) dx  ≤  U(A,Px)  ≤  U(f,P).L(f,P) \;\le\; L(A,P_x) \;\le\; \int_a^b A(x)\,dx \;\le\; U(A,P_x) \;\le\; U(f,P).

Also ∬Rf dA\iint_R f\,dA lies between L(f,P)L(f,P) and U(f,P)U(f,P) by definition of the double integral. Therefore both ∫abA(x) dx\int_a^b A(x)\,dx and ∬Rf dA\iint_R f\,dA lie in the interval [L(f,P),U(f,P)][L(f,P), U(f,P)], whose length is U(f,P)−L(f,P)<εU(f,P)-L(f,P) < \varepsilon. Hence

∣∫abA(x) dx  −  ∬Rf dA∣<ε.\left| \int_a^b A(x)\,dx \;-\; \iint_R f\,dA \right| < \varepsilon.

As ε>0\varepsilon>0 was arbitrary, the two quantities are equal:

∫ab∫cdf(x,y) dy dx  =  ∬Rf dA.\int_a^b \int_c^d f(x,y)\,dy\,dx \;=\; \iint_R f\,dA.

The identical argument with the roles of xx and yy interchanged (partitioning yy first, bounding B(y)B(y) between analogous Darboux sums) gives

∫cd∫abf(x,y) dx dy  =  ∬Rf dA.\int_c^d \int_a^b f(x,y)\,dx\,dy \;=\; \iint_R f\,dA.

Combining the two displayed equalities proves all three quantities in the theorem are equal. ■\blacksquare

∎

Worked examples

Example 1 — a genuinely asymmetric integrand. Compute ∬R(3x2y+cos⁡(πy)) dA\displaystyle\iint_R (3x^2y + \cos(\pi y))\,dA where R=[1,2]×[0,1/2]R = [1,2] \times [0, 1/2], using both orders of integration and checking they agree. The integrand f(x,y)=3x2y+cos⁡(πy)f(x,y) = 3x^2 y + \cos(\pi y) is a polynomial-plus-cosine, hence continuous on R2\mathbb{R}^2, so continuous on the compact rectangle RR. Fubini's Theorem applies. Order 1: integrate yy first, then xx.

A(x)=∫01/2(3x2y+cos⁡(πy)) dy=[3x2y22+sin⁡(πy)π]y=0y=1/2=3x28+1π.A(x) = \int_0^{1/2} \big(3x^2 y + \cos(\pi y)\big)\,dy = \left[ \frac{3x^2 y^2}{2} + \frac{\sin(\pi y)}{\pi} \right]_{y=0}^{y=1/2} = \frac{3x^2}{8} + \frac{1}{\pi}.

(Using sin⁡(π/2)=1\sin(\pi/2) = 1, sin⁡0=0\sin 0 = 0.) Now integrate over x∈[1,2]x\in[1,2]:

∫12(3x28+1π)dx=[x38+xπ]12=(88+2π)−(18+1π)=1−18+1π=78+1π.\int_1^2 \left( \frac{3x^2}{8} + \frac{1}{\pi} \right) dx = \left[ \frac{x^3}{8} + \frac{x}{\pi} \right]_1^2 = \left(\frac{8}{8}+\frac{2}{\pi}\right) - \left(\frac{1}{8}+\frac{1}{\pi}\right) = 1 - \frac18 + \frac{1}{\pi} = \frac{7}{8} + \frac{1}{\pi}.

Order 2: integrate xx first, then yy.

B(y)=∫12(3x2y+cos⁡(πy)) dx=[x3y+xcos⁡(πy)]x=1x=2=(8y+2cos⁡(πy))−(y+cos⁡(πy))=7y+cos⁡(πy).B(y) = \int_1^2 \big(3x^2 y + \cos(\pi y)\big)\,dx = \left[ x^3 y + x\cos(\pi y)\right]_{x=1}^{x=2} = (8y + 2\cos(\pi y)) - (y + \cos(\pi y)) = 7y + \cos(\pi y).

Now integrate over y∈[0,1/2]y \in [0, 1/2]:

∫01/2(7y+cos⁡(πy)) dy=[7y22+sin⁡(πy)π]01/2=72⋅14+1π=78+1π.\int_0^{1/2} \big(7y + \cos(\pi y)\big)\,dy = \left[ \frac{7y^2}{2} + \frac{\sin(\pi y)}{\pi}\right]_0^{1/2} = \frac{7}{2}\cdot\frac14 + \frac{1}{\pi} = \frac{7}{8} + \frac{1}{\pi}.

Both orders give 78+1π≈0.875+0.3183=1.1933\dfrac{7}{8} + \dfrac{1}{\pi} \approx 0.875 + 0.3183 = 1.1933, confirming Fubini's Theorem and, incidentally, illustrating that neither order was obviously "easier" here — both required exactly one polynomial and one trigonometric antiderivative.

Example 2 — using order-swap to simplify a hard integral. Compute ∬Rxexy dA\displaystyle\iint_R x e^{xy}\,dA where R=[0,1]×[0,2]R = [0,1]\times[0,2]. Here f(x,y)=xexyf(x,y) = xe^{xy} is continuous (product and composition of continuous functions), so Fubini applies and both orders must agree — but they are not equally easy to execute, which is the real lesson. Attempt order 1: integrate xx first. We'd need ∫01xexy dx\int_0^1 xe^{xy}\,dx for fixed yy. This requires integration by parts (treating yy as a constant): with u=x,dv=exydxu = x, dv = e^{xy}dx, so du=dxdu = dx, v=exy/yv = e^{xy}/y (for y≠0y\neq 0):

∫01xexy dx=[xexyy]01−∫01exyy dx=eyy−1y[exyy]01=eyy−ey−1y2.\int_0^1 x e^{xy}\,dx = \left[\frac{x e^{xy}}{y}\right]_0^1 - \int_0^1 \frac{e^{xy}}{y}\,dx = \frac{e^{y}}{y} - \frac{1}{y}\left[\frac{e^{xy}}{y}\right]_0^1 = \frac{e^y}{y} - \frac{e^y - 1}{y^2}.

This is correct but messy, and we'd still need to integrate this expression over y∈[0,2]y\in[0,2] (with a removable singularity to handle carefully at y=0y=0) — doable, but painful. Order 2: integrate yy first. For fixed xx, treat xx as constant:

B(x)=∫02xexy dy.B(x) = \int_0^2 x e^{xy}\,dy.

If x≠0x \ne 0, substitute u=xyu = xy, du=x dydu = x\,dy:

B(x)=∫02xexy dy=[exy]y=0y=2=e2x−1.B(x) = \int_0^2 x e^{xy}\,dy = \left[ e^{xy} \right]_{y=0}^{y=2} = e^{2x} - 1.

(Check via direct antidifferentiation in yy: ddyexy=xexy\frac{d}{dy}e^{xy} = xe^{xy}, so ∫02xexydy=exy∣02=e2x−1\int_0^2 xe^{xy}dy = e^{xy}\big|_0^2 = e^{2x}-1 — no division by xx ever needed, and the formula manifestly also gives the correct value 00 at x=0x=0 since B(0)=∫020 dy=0=e0−1B(0) = \int_0^2 0\,dy = 0 = e^0 - 1.) So

∫01B(x) dx=∫01(e2x−1) dx=[e2x2−x]01=(e22−1)−12=e22−32.\int_0^1 B(x)\,dx = \int_0^1 \big(e^{2x}-1\big)\,dx = \left[\frac{e^{2x}}{2} - x\right]_0^1 = \left(\frac{e^2}{2} - 1\right) - \frac12 = \frac{e^2}{2} - \frac32.

Numerically, e22−32≈7.3892−1.5=3.6945−1.5=2.1945\frac{e^2}{2} - \frac32 \approx \frac{7.389}{2} - 1.5 = 3.6945 - 1.5 = 2.1945. This is dramatically cleaner than order 1. Fubini's Theorem is what licenses choosing the easier order 2 — without it, there would be no guarantee that computing "yy first" answers the same question as the double integral originally posed. (One can verify the messier order-1 computation gives the same numeric answer 2.1945…2.1945\ldots after integrating eyy−ey−1y2\frac{e^y}{y}-\frac{e^y-1}{y^2} from 00 to 22, but this is exactly the computation Fubini lets us avoid.)

CounterexampleWhy continuity cannot be dropped: a bounded, discontinuous ff where the iterated integrals disagree

Consider R=[0,1]×[0,1]R = [0,1]\times[0,1] and the bounded function

f(x,y)={x2−y2(x2+y2)2(x,y)≠(0,0)0(x,y)=(0,0).f(x,y) = \begin{cases} \dfrac{x^2 - y^2}{(x^2+y^2)^2} & (x,y) \ne (0,0) \\[2mm] 0 & (x,y) = (0,0). \end{cases}

This ff is continuous at every point of RR except the origin, where it is discontinuous (indeed unbounded near the origin along most paths, though we can still ask about the iterated integrals as improper/pointwise limits — the point of this example is the classical asymmetric-kernel phenomenon). To keep things elementary and fully rigorous with an actually bounded function, we use the standard sharper example built from a difference of bumps. Let gng_n for n=1,2,3,…n = 1,2,3,\dots be a continuous "tent" function on [0,1][0,1] supported on [1/(n+1),1/n][1/(n+1), 1/n], with gn≥0g_n \ge 0, gn(1/(n+1))=gn(1/n)=0g_n(1/(n+1)) = g_n(1/n) = 0, and ∫01gn(t) dt=1\int_0^1 g_n(t)\,dt = 1 (a triangular spike of height 2n(n+1)2n(n+1) at the midpoint, base 1/n−1/(n+1)=1/(n(n+1))1/n - 1/(n+1) = 1/(n(n+1)), so that its area is 11). Define on R=[0,1]2R=[0,1]^2:

f(x,y)=∑n=1∞(gn(x)gn(y)−gn+1(x)gn+1(y))∣restricted so only one term is nonzero at each (x,y)f(x,y) = \sum_{n=1}^\infty \big(g_n(x)g_n(y) - g_{n+1}(x)g_{n+1}(y)\big)\Big|_{\text{restricted so only one term is nonzero at each }(x,y)}

More concretely and standardly, define ff on the square [0,1]2[0,1]^2 by: partition (0,1](0,1] into intervals In=(1/(n+1),1/n]I_n = (1/(n+1), 1/n] for n≥1n\ge1. On the square In×InI_n \times I_n, let f(x,y)=an⋅hn(x)hn(y)f(x,y) = a_n \cdot h_n(x) h_n(y) alternating in sign so that ∫In∫Inf dy dx=1\int_{I_n}\int_{I_n} f\,dy\,dx = 1 while ∫In∫Inf dx dy=−1\int_{I_n}\int_{I_n} f\,dx\,dy = -1 for each nn, and set f=0f = 0 off the union of these squares (and at the origin). Such ff is bounded (each block is bounded, and blocks shrink to 00 so no single value blows up), but ff is discontinuous at the origin: along the sequence of block-centers ff takes values bouncing between roughly +cn+c_n and −cn-c_n that do not converge to f(0,0)=0f(0,0)=0, so lim⁡(x,y)→(0,0)f(x,y)\lim_{(x,y)\to(0,0)} f(x,y) does not exist. For this ff, direct computation block-by-block gives

∫01 ⁣ ⁣∫01f(x,y) dy dx=∑n=1∞1⋅(weight decaying appropriately, arranged so the series telescopes to)  1≠−1=∫01 ⁣ ⁣∫01f(x,y) dx dy.\int_0^1\!\!\int_0^1 f(x,y)\,dy\,dx = \sum_{n=1}^\infty 1 \cdot (\text{weight decaying appropriately, arranged so the series telescopes to}) \; 1 \ne -1 = \int_0^1\!\!\int_0^1 f(x,y)\,dx\,dy .

The clean, standard instance of this phenomenon (found in Rudin's Principles of Mathematical Analysis and many other texts) uses instead the simple rational function on [0,1]×[0,1]∖{(0,0)}[0,1]\times[0,1]\setminus\{(0,0)\}:

f(x,y)=x2−y2(x2+y2)2,f(x,y) = \frac{x^2-y^2}{(x^2+y^2)^2},

extended by f(0,0)=0f(0,0)=0. One computes directly, for fixed x∈(0,1]x \in (0,1],

∫01f(x,y) dy=∫01x2−y2(x2+y2)2 dy=[yx2+y2]y=0y=1=1x2+1,\int_0^1 f(x,y)\,dy = \int_0^1 \frac{x^2-y^2}{(x^2+y^2)^2}\,dy = \left[ \frac{y}{x^2+y^2} \right]_{y=0}^{y=1} = \frac{1}{x^2+1},

using the antiderivative identity ddy(yx2+y2)=x2−y2(x2+y2)2\frac{d}{dy}\left(\frac{y}{x^2+y^2}\right) = \frac{x^2-y^2}{(x^2+y^2)^2} (verify by the quotient rule: (x2+y2)−y⋅2y(x2+y2)2=x2−y2(x2+y2)2\frac{(x^2+y^2) - y\cdot 2y}{(x^2+y^2)^2} = \frac{x^2-y^2}{(x^2+y^2)^2}, confirmed). Then

∫01 ⁣ ⁣∫01f(x,y) dy dx=∫01dx1+x2=arctan⁡(1)−arctan⁡(0)=π4.\int_0^1\!\!\int_0^1 f(x,y)\,dy\,dx = \int_0^1 \frac{dx}{1+x^2} = \arctan(1) - \arctan(0) = \frac{\pi}{4}.

By the antisymmetry f(x,y)=−f(y,x)f(x,y) = -f(y,x) (visible directly from the formula: swapping x,yx,y negates the numerator, leaves the denominator fixed), the same computation with roles swapped gives

∫01 ⁣ ⁣∫01f(x,y) dx dy=−π4.\int_0^1\!\!\int_0^1 f(x,y)\,dx\,dy = -\frac{\pi}{4}.

So the two iterated integrals are π/4\pi/4 and −π/4-\pi/4 — they disagree, and in fact the double integral ∬Rf dA\iint_R f\,dA does not even exist in this case (the double limit fails to converge, since ff is unbounded near the origin — every neighborhood of (0,0)(0,0) contains points where ∣f∣|f| is arbitrarily large, e.g. along y=0,x→0+y=0,x\to0^+, f(x,0)=1/x2→∞f(x,0) = 1/x^2 \to \infty). This is precisely the failure of hypothesis (continuity, which here fails at the single point (0,0)(0,0) due to unboundedness) that Fubini's Theorem requires: drop continuity, and not only can the iterated integrals fail to match each other, the very double integral they were supposed to compute may cease to exist at all.

Scope, edge cases & common traps

SituationWhat happens
ff continuous on closed rectangle RRFubini applies exactly as stated: double integral = either iterated integral.
ff bounded, integrable, but discontinuous on a "small" (measure-zero) setA more general Fubini–Tonelli-type theorem (Lebesgue framework) still guarantees agreement of the iterated integrals with the double integral, for almost every slice; this is a genuinely stronger theorem beyond the Riemann/continuous setting treated here.
ff bounded but discontinuous on a "large" set, or unbounded near a pointIterated integrals may exist separately yet disagree (see counterexample), or the double integral may fail to exist altogether.
Region is not a rectangle (e.g. a disk, triangle)Iterated integrals still work but the limits of integration become functions of the other variable; this is the subject of Integrals over General Regions.
Triple (or higher) integrals over a box [a1,b1]×[a2,b2]×[a3,b3][a_1,b_1]\times[a_2,b_2]\times[a_3,b_3]The same theorem holds verbatim for continuous f:box→Rf:\text{box}\to\mathbb R, with any of the 3!=63! = 6 orders of iterated single integrals agreeing; proof is a direct induction on the two-variable case.

The mathematician behind it

Origins · 1 minGuido Fubini and how this was discovered
GF
Guido Fubini
1879–1943 · Italian analysis and differential geometry, Turin and later New York

Guido Fubini was an Italian mathematician working chiefly at the University of Turin, contributing broadly to real analysis, differential geometry (including the "Fubini–Study metric" used in complex projective geometry and later in quantum information theory), and mathematical physics. His 1907 theorem on interchange of order of integration, proved in the more general Lebesgue-measure setting, became one of the most quoted results in all of analysis because of its sheer everyday utility. Fubini was forced from his position in Italy in 1938 by the Fascist racial laws targeting Jewish academics, and he emigrated to the United States, joining the Institute for Advanced Study and later teaching at New York University until his death in 1943.

History
- 1907 — Guido Fubini proves the general measure-theoretic order-of-integration theorem for Lebesgue integrals.
- 1900s — Continuous-function, Riemann-integral special cases (essentially the theorem proved above) were already implicit in 19th-century treatments of multiple integrals via iterated slicing, predating Fubini's general measure-theoretic statement.
- 1910 — Leonida Tonelli extends the interchange-of-order result to non-negative measurable functions without requiring absolute integrability, giving the complementary "Tonelli's Theorem" used alongside Fubini's in modern measure theory.
- 1938 — Fubini, targeted by Italy's Fascist racial laws, is dismissed from his Italian academic post and emigrates.
- 1943 — Fubini dies in New York; his name remains permanently attached to the order-of-integration theorem taught in every real analysis and calculus curriculum.
Does Fubini's Theorem require the function to be positive?

No. The version proved here requires only that ff be continuous (hence bounded, since RR is compact) on the rectangle RR; it may take negative values freely. Positivity becomes relevant in the more general Lebesgue-theoretic setting, where Tonelli's Theorem handles non-negative measurable functions (allowing infinite integrals) without requiring absolute integrability, while the full Fubini's Theorem for signed or complex measurable functions requires ∬R∣f∣ dA<∞\iint_R |f|\,dA < \infty. In the continuous-function-on-a-compact-rectangle setting of this article, boundedness and hence absolute integrability come for free from continuity plus compactness, so no extra positivity or integrability hypothesis is needed.

Why does the proof bound the iterated integral between L(f,P)L(f,P) and U(f,P)U(f,P) instead of directly comparing Riemann sums?

Working with the Darboux upper and lower sums (built from suprema and infima of ff on each subrectangle) avoids having to track the behavior of arbitrary sample points inside each subrectangle, which is exactly the extra freedom that ordinary Riemann sums allow and that would complicate the argument. By sandwiching both the iterated integral's value and the double integral's value inside the same shrinking interval [L(f,P),U(f,P)][L(f,P), U(f,P)], the proof shows their difference is smaller than any ε(b−a)(d−c)\varepsilon(b-a)(d-c), forcing equality. This squeeze technique is standard whenever one needs to show two limiting quantities coincide without directly manipulating limits.

Is it true that if the two iterated integrals of some function agree, the function is automatically well-behaved enough for Fubini to apply?

No — this converse is false and is exactly the trap discussed in the warning box. It is possible to construct bounded, badly discontinuous functions where both iterated integrals exist and even happen to be numerically equal, yet the underlying double integral fails to exist because the two-dimensional Riemann sums do not converge. Agreement of iterated integrals is a necessary consequence of continuity (or suitable integrability), not a sufficient replacement for it. Always verify a genuine hypothesis on ff over the whole rectangle — continuity, or an appropriate integrability/measurability condition — before invoking any interchange-of-order argument.

Does Fubini's Theorem extend to triple integrals and higher dimensions?

Yes, and the extension is essentially by induction on dimension using the two-variable case proved here as the base step. For ff continuous on a box [a1,b1]×[a2,b2]×[a3,b3][a_1,b_1]\times[a_2,b_2]\times[a_3,b_3], one first treats the box as [a1,b1]×Q[a_1,b_1]\times Q where Q=[a2,b2]×[a3,b3]Q=[a_2,b_2]\times[a_3,b_3], applies the two-variable theorem to peel off the x1x_1-integral, and then applies it again inside QQ to separate x2x_2 and x3x_3; continuity of the relevant "slice" functions at each stage follows from the same uniform continuity argument used above. The conclusion is that any of the 3!=63! = 6 possible orders of iterating three single-variable integrals over a continuous integrand on a box give the same value, matching the triple integral. This is exactly the tool needed once you study volumes, mass, and moments over three-dimensional regions.

How does Fubini's Theorem relate to swapping order of integration when the region is not a rectangle, like a triangle or disk?

On a non-rectangular region, the limits of the inner integral become functions of the outer variable rather than fixed constants — for instance integrating over a triangle 0≤y≤x≤10\le y\le x\le 1 gives ∫01∫0xf(x,y) dy dx\int_0^1\int_0^x f(x,y)\,dy\,dx, and swapping order requires rewriting the same region as 0≤y≤1, y≤x≤10 \le y \le 1,\, y \le x \le 1 to get ∫01∫y1f(x,y) dx dy\int_0^1\int_y^1 f(x,y)\,dx\,dy. The underlying justification is still a Fubini-type theorem, now stated for general bounded regions with continuous or piecewise-continuous boundary functions; this is developed fully under Integrals over General Regions, building directly on the rectangular case proved here.

Why is uniform continuity, rather than plain pointwise continuity, the key ingredient in the proof?

Plain continuity at each point only guarantees that ff doesn't jump near that specific point, with a δ\delta that may depend on which point you're at. The proof needs a single δ\delta that works uniformly across the entire rectangle RR, so that a partition fine enough (mesh <δ<\delta) controls the oscillation of ff on every subrectangle simultaneously — this is what makes U(f,P)−L(f,P)U(f,P)-L(f,P) small in one shot. The Heine–Borel compactness of RR is precisely what upgrades pointwise continuity to uniform continuity, which is why the hypothesis "ff continuous on the closed, bounded rectangle RR" (not an open or unbounded region) is essential to the argument.

Before this

Where this leads