The definite integral of a function of a single variable is defined to be the limit over all applicable partitions of a Riemann sum. In imitation of this, the double integral of a function of two variables is defined to be the (bivariate) limit of an appropriate Riemann sum.

In two dimensions, where the domain of integration is rectangular in shape and oriented so the sides are parallel to the coordinate axes, the horizontal and vertical edges are partitioned so that a grid of sub-rectangles covers the domain. The function defined over the domain is evaluated at one point in each sub-rectangle, and this value is multiplied by the area of the corresponding sub-rectangle. The Riemann sum is the sum of all the volumes of the rectangular "blocks" so determined. The limit as the number of blocks becomes infinite is, under the appropriate conditions, taken as the meaning of the double integral.

Figure 5.1.1(a) shows the surface $z\=f\left(x\,y\right)\equiv x^2\+y^2$ graphed over the square $0\le x\,y\le 1$. It is the volume beneath this surface that is to be computed by the double integral.

Figure 5.1.1(b) shows a $5\times 5$ grid for which a Riemann sum is constructed. The red dot at the center of each sub-rectangle is the point at which the height of the surface is sampled. Figure 5.1.1(c) shows the "blocks" whose volumes are analogous to the rectangles that approximate area in a Riemann sum for a single integral.

Computationally, the content of Figures 5.1.1(a - c) can be implemented once the coordinates of the "red dots" in Figure 5.1.1(b) are specified as $\left(.1\+.2 i\,.1\+.2 j\right)\,i\=0\,\dots \,4\,j\=0\,\dots \,4$. The sum of the volumes of the 25 blocks so determined is

$\sum _{i\=0}^4 \sum _{j\=0}^{4}f\left(.1\+.2 i\,.1\+.2 j\right)\cdot \left(.04\right)\=0.660$ 

Passage to the limit as the density of the grid increases requires discretizing so that $\mathrm{\Δ} x\=1\/m$ and $y\=1\/n$, in which case the red-dot midpoints become

 $\left(\frac{1}{2 m}\+\frac{i}{m}\,\frac{1}{2 n}\+\frac{j}{n}\right)\,i\=0\,\dots \,m-1\,j\=0\,\dots \,n-1$ 

and the approximating sum becomes

$\sum _{i\=0}^{m-1}\sum _{j\=0}^{n-1}f\left(\frac{1}{2 m}\+\frac{i}{m}\,\frac{1}{2 n}\+\frac{j}{n}\right)\cdot \left(\frac{1}{m n}\right)\=\frac{2}{3}-\frac{1}{12}\left(\frac{1}{m^2}\+\frac{1}{n^2}\right)$ 

which, in the limit as $m$ and $n$ become infinite, approaches $2\/3$ as a limiting value.

A simpler Riemann sum can be constructed if the evaluation point for each sub-rectangle is taken as the "lower-left" corner where the generic coordinate would be $\left(i\/m\,j\/n\right)$ and the approximating sum would be

$\sum _{i\=0}^{m-1} \sum _{j\=0}^{n-1}f\left(\frac{i}{m}\,\frac{j}{n}\right)\cdot \left(\frac{1}{m n}\right)\=\frac{2}{3}\+\frac{1}{6}\left(\frac{1}{m^2}\+\frac{1}{n^2}\right)-\frac{1}{2}\left(\frac{1}{m}\+\frac{1}{n}\right)$

Again, the limiting value as $m$ and $n$ become infinite is $2\/3$. Finally, if the evaluation point in each sub-rectangle is taken as the "lower-right" corner where the generic coordinate would still be $\left(i\/m\,j\/n\right)$, the approximating sum would be

$\sum _{i\=1}^m \sum _{j\=1}^{n}f\left(\frac{i}{m}\,\frac{j}{n}\right)\cdot \left(\frac{1}{m n}\right)\=\frac{2}{3}\+\frac{1}{6}\left(\frac{1}{m^2}\+\frac{1}{n^2}\right)\+\frac{1}{2}\left(\frac{1}{m}\+\frac{1}{n}\right)$ 

with limiting value $2\/3$.

With the following syntax, the ApproximateInt command in the Student MultivariateCalculus package will draw the graph seen in Figure 5.1.1(e).

With the following syntax, the ApproximateInt command in the Student MultivariateCalculus package will return the approximate signed volume between the surface and the $\mathrm{xy}$-plane.