Volume under a surface by stacking prisms, slicing it either way round, and integrating in polar coordinates
A single integral adds up thin strips to find an area. A double integral adds up thin prisms to find a volume, a mass or an average over a region of the plane. The definition is a limit of sums, but nobody computes it that way: you slice the solid and integrate twice, in whichever order is easier, or change to coordinates that fit the region.
Raise n and watch the prisms close in on the surface. Switch between corner and midpoint sampling and compare how fast the error falls in the log–log strip.
The exact value is computed by 64 × 64-point Gauss–Legendre quadrature, which is accurate to many digits for smooth integrands. Prisms below the floor (rose) count negatively: the double integral is a signed volume, and for sin(x + y) on [0, π]² the two parts cancel exactly. Heights in 3D are scaled so the tallest prism fits; white dots mark the sample points (hidden above n = 30).
Cut the rectangle R into cells of area ΔA, choose a sample point (xᵢⱼ*, yᵢⱼ*) in each, and add up the volumes of the prisms:
∬_R f(x, y) dA = lim (n → ∞) ∑ᵢ ∑ⱼ f(xᵢⱼ*, yᵢⱼ*) ΔA
For a continuous f the limit exists and does not depend on how the sample points are chosen. Where f is negative the prisms hang below the floor and count negatively, so the double integral is a signed volume. Dividing by the area of R gives the average value of f over R.
How fast the sums converge depends on the sample point, just as for single integrals. On a smooth surface a corner sample is wrong by an amount proportional to the slope times the cell size in every cell, and these errors share a sign, so the total error is of order 1/n. Sampling at the centre cancels the first-order term by symmetry and leaves an error of order 1/n². The log–log plot shows these as lines of slope −1 and −2. A random sample point gives an error that fluctuates, which is the starting point of Monte Carlo integration.
Sweep the cutting plane through the solid. Each slice has an area A(x), plotted below; the area under that graph is the volume. Switch the slicing direction and watch the limits change while the volume stays put.
Slice areas and both iterated integrals are computed by adaptive Simpson quadrature. In 3D the heights are scaled so the tallest point of f on the unit square stands 2.2 units high; the indigo floor is R, and the grey walls rise from its boundary to the surface.
Slice the solid with the plane at a fixed x. The slice is a region under the curve y ↦ f(x, y), with area A(x) = ∫ f(x, y) dy taken across R at that x. Adding up the slices, each of thickness dx, gives the volume ∫ A(x) dx. This is Cavalieri’s principle (Bonaventura Cavalieri, 1635) turned into a computation:
∬_R f dA = ∫ₐᵇ ( ∫ f(x, y) dy ) dx = ∫_c^d ( ∫ f(x, y) dx ) dy
Guido Fubini proved in 1907 that, for integrable functions, the double integral equals the iterated integrals in either order. For a rectangle the limits are constants. For other regions the inner limits depend on the outer variable, and reversing the order means redrawing the region: the triangle 0 ≤ y ≤ x ≤ 1 is ∫₀¹ ∫₀ˣ … dy dx one way and ∫₀¹ ∫ᵧ¹ … dx dy the other.
Reversing the order can turn an impossible integral into an easy one. ∫₀¹ ∫ᵧ¹ e^(x²) dx dy has no elementary inner integral, but over the same triangle in the other order it is ∫₀¹ x e^(x²) dx = (e − 1)/2. Triple integrals work the same way with one more layer: ∭_E f dV is a limit of sums of boxes, computed as three nested single integrals, in any of six orders.
Cut a region into polar cells and see how their size grows with distance from the origin. Click a cell to magnify it, and compare the polar sum with a sum over Cartesian squares.
The polar sum evaluates f at the centre (r̄, θ̄) of each cell; wedges of a region with a curved outer edge use the radius at the middle angle of the wedge. The Cartesian sum uses the squares of an n × n grid over the bounding box whose centres lie inside R. The exact value is nested adaptive quadrature in polar coordinates.
A polar rectangle r₀ ≤ r ≤ r₁, θ₀ ≤ θ ≤ θ₁ is the difference of two circular sectors, so its area is
½r₁²Δθ − ½r₀²Δθ = ½(r₁ + r₀)(r₁ − r₀)Δθ = r̄ Δr Δθ
exactly, with r̄ the middle radius. Its sides are Δr and an arc of length about r̄Δθ, so for small cells it is nearly a rectangle with those sides. Summing f times these areas and letting the cells shrink gives
∬_R f dA = ∫∫ f(r cos θ, r sin θ) r dr dθ
The extra r is the stretching factor of the map (r, θ) ↦ (x, y), its Jacobian determinant. It is why polar coordinates tame integrals like ∬ e^(−(x²+y²)) dA over the whole plane: with r dr dθ the inner integral ∫ r e^(−r²) dr has an elementary antiderivative, the answer is π, and taking a square root gives the famous ∫ e^(−x²) dx = √π over the whole line.
Spot a mistake or have an idea?