Assume the region of integration is bounded by the x axis below, f(x) above, and the vertical lines x = a and x = b on the left and right. Let h(x,y) be the function to integrate. Assume h is integrable over the region.
Consider a nested integral, i.e. the integral of an integral. Integrate g(x), as x runs from a to b, where g(x) is the integral of h(x,y) as y runs from 0 to f(x). Assume the inner integral always exists, so that each g(x) is well defined. This is the case when h is continuous.
We will prove that this nested integral is equal to the double integral. This is rather technical, so hang on. Remember that the order of integration is arbitrary in this proof. Integrating with respect to y, and then x, is the same as integrating with respect to x, and then y, since both are equal to the double integral. So here we go, back to Riemann nets and lower and upper sums.
Select horizontal and vertical nets, thus building a 2 dimensional net over the region. Focus on one of the x coordinates in the horizontal net. We know that g(x), the integral of h(x,y) at x, is trapped between the lower sum and upper sum, according to the vertical net. Call these sums l and u. Let r be the gap between x and the next point in the net. If h had the same profile, from x to x+r, the lower sum over the squares in this strip, within the two dimensional lower sum, would be l×r, and the upper sum would be u×r. Actually the lower sum is probably less and the upper sum is probably greater. Thus the lower and upper sums of the strip of squares from x to x+r bracket r×g(x). This is the term used to compute the Riemann sum of the integral of g(x). The Riemann sum of g(x), defined by the horizontal net, is trapped between the lower and upper sums, as defined by the two dimensional net. Remember that the lower and upper sums of all 2 dimensional nets converge to the double integral. Therefore the Riemann sums of g(x), for finer and finer nets, are squashed between the converging lower and upper sums of the double integral. The nested integral equals the double integral.
Note that a nested integral does not imply a double integral. Let h = 1 when x is rational and -1 when x is irrational, over the unit square. As y runs from 1 to 2, let h(x,y) = - h(x,y-1). The second square, pasted on top of the first, is the opposite of the first. For every x, the integral of g(x) is 0, hence the nested integral is 0. However, the double integral does not exist.
Assuming h is integrable over the region, we can redefine or undefine entire vertical or horizontal lines. This does not change the double integral or the nested integral.
If the region of integration is bounded by one function above and another function below, the nested integral still matches the double integral. Let g(x) be the integral of h(x,y) as y runs between the two functions, and apply the same proof.
When computing a double integral, we sometimes break the region up into manageable subregions, apply a nested integral to each, and add the results together. Remember that the nesting order is arbitrary, and it is sometimes easier to integrate g(y), rather than g(x).
Return to the region under the parabola y = x2 from 0 to 1. Let h(x) = x+y, which is certainly integrable over the region. In fact it looks like a slanted roof rising up over the first quadrant. As we compute the integral of h, we are really finding the volume of a shape with a floor bounded by the parabola, having area 1/3, and a ceiling that slants up from the origin. You'll probably never see this shape in real life, but it's a nice little example.
Fix x, and integrate x+y as y runs from 0 to x2. Think of x as a constant, so the integral of x is xy. The integral of y is y2/2 as usual. Evaluate at 0 and x2 and get:
g(x) = x3 + x4/2
Integrate this from 0 to 1 to get 1/4 + 1/10 = 0.35.
Given the same shape, reverse the order of integration. Let y be a constant and integrate x+y as x runs from sqrt(y) to 1. This gives:
g(y) = y + 1/2 - y/2 - y×sqrt(y)
Integrate this from 0 to 1 to get 1/2+1/2-1/4-2/5 = 0.35. We got the same answer, but as you can see, the problem can be easy or hard, depending on the order of integration.
Generalize this theorem to higher dimensions by induction. If we're dealing with a 3 dimensional region, let the x coordinate determine a cross section whose double integral is given by g(x). Show that the integral of g(x) is the same as the triple integral of h over the entire region. As before, g(x)×r is trapped between the lower and upper sums for that layer. The Riemann sums that determine the integral of g(x) are held arbitrarily close to the triple integral, and our nested integral equals the triple integral. Of course we will probably compute each double integral by a nested integral, giving a triply nested integral.