Repeated integration is justified by product-measure theory, not by notation alone. Tonelli’s theorem and Fubini’s theorem answer different questions and impose different hypotheses.

Product measures

Let $(X,\mathcal A,\mu)$ and $(Y,\mathcal B,\nu)$ be $\sigma$-finite measure spaces. Their product $\sigma$-algebra $\mathcal A\otimes\mathcal B$ is generated by measurable rectangles $A\times B$. The product measure is characterized by

\[(\mu\otimes\nu)(A\times B)=\mu(A)\nu(B).\]

The $\sigma$-finiteness assumption ensures the standard existence and uniqueness theorem and supports the section results used below. Lebesgue–Stieltjes measures that are finite on bounded intervals are $\sigma$-finite on $\mathbb R$.

For a measurable set $E\subset X\times Y$, define its sections

\[E_x=\{y\in Y:(x,y)\in E\},\qquad E^y=\{x\in X:(x,y)\in E\}.\]

Then the section measures are measurable and

\[(\mu\otimes\nu)(E) =\int_X\nu(E_x)\,d\mu(x) =\int_Y\mu(E^y)\,d\nu(y).\]

The equality is first verified for rectangles, extended to an appropriate monotone class of sets, and then lifted from indicator functions to general functions.

Tonelli’s theorem: nonnegative functions

Theorem. If $f:X\times Y\to[0,+\infty]$ is measurable, then \(\int_{X\times Y}f\,d(\mu\otimes\nu) = \int_X\left(\int_Y f(x,y)\,d\nu(y)\right)d\mu(x)\) \(= \int_Y\left(\int_X f(x,y)\,d\mu(x)\right)d\nu(y),\) with the common value allowed to be $+\infty$.

Tonelli requires no prior integrability assumption because nonnegative contributions cannot cancel. It is the correct theorem for proving integrability: compute either repeated integral of $ f $; if the result is finite, then Fubini applies to $f$.

Fubini’s theorem: integrable functions

Theorem. If $f\in L^1(\mu\otimes\nu)$, then for almost every $x$, $f(x,\cdot)\in L^1(\nu)$; for almost every $y$, $f(\cdot,y)\in L^1(\mu)$; and both repeated integrals equal the product-space integral.

The decisive hypothesis is

\[\int_{X\times Y}|f(x,y)|\,d(\mu\otimes\nu)<\infty.\]

Finite values of two formal iterated integrals do not by themselves establish this condition.

Why order can fail

Define on $(0,1)^2$

\[f(x,y)=\frac{x^2-y^2}{(x^2+y^2)^2}.\]

For fixed $x>0$,

\[f(x,y)=\frac{\partial}{\partial y} \left(\frac{y}{x^2+y^2}\right),\]

so

\[\int_0^1 f(x,y)\,dy=\frac{1}{x^2+1}, \qquad \int_0^1\!\int_0^1 f(x,y)\,dy\,dx=\frac{\pi}{4}.\]

For fixed $y>0$,

\[f(x,y)=-\frac{\partial}{\partial x} \left(\frac{x}{x^2+y^2}\right),\]
which gives the reverse order as $-\pi/4$. The singularity at the origin makes $ f $ nonintegrable. Fubini’s hypothesis fails, so the disagreement is not a contradiction.

Mixed continuous and discrete factors

Let

\[\mu=dx+\sum_{i=1}^{m}a_i\delta_{x_i}, \qquad \nu=dy+\sum_{j=1}^{n}b_j\delta_{y_j},\]

with nonnegative weights. For a nonnegative or integrable $f$,

\[\begin{aligned} \int f\,d(\mu\otimes\nu) ={}&\iint f(x,y)\,dx\,dy\\ &+\sum_i a_i\int f(x_i,y)\,dy\\ &+\sum_j b_j\int f(x,y_j)\,dx\\ &+\sum_{i,j}a_i b_j f(x_i,y_j). \end{aligned}\]

The product measure automatically creates area, line, and point contributions. This is useful when a two-parameter model contains distributed effects together with feature-aligned or nodal effects.

Tensor-product approximation

Suppose $f(x,y)=u(x)v(y)$. Then, whenever the factors are integrable,

\[\int_{X\times Y}u(x)v(y)\,d(\mu\otimes\nu) =\left(\int_Xu\,d\mu\right) \left(\int_Yv\,d\nu\right).\]

Finite sums of separable functions,

\[f_r(x,y)=\sum_{k=1}^{r}u_k(x)v_k(y),\]

lead to low-rank quadrature and tensor-product discretizations. Their efficiency comes from algebraic separation; their validity still depends on measurable and integrable approximations.

CAD and simulation interpretation

Surface parameter domains, tensor-product splines, mass-property integrals, and variational forms all produce product-space integrals. A reliable implementation distinguishes three issues:

  • mathematical interchange: Tonelli or Fubini supplies the condition;
  • parameterization: a surface map contributes the area Jacobian $|\partial_uS\times\partial_vS|$;
  • execution order: loop ordering affects cache behavior and rounding, although it must not alter the exact integral when Fubini applies.

For a regular parameterized surface $S:D\to\mathbb R^3$,

\[\int_{S(D)}g\,dA =\int_D g(S(u,v)) \left\|\frac{\partial S}{\partial u} \times \frac{\partial S}{\partial v}\right\|du\,dv.\]

Degenerate parameter points and trimmed regions require separate geometric analysis; a rectangular parameter box alone does not describe the physical surface measure.