Laboratory 6
Double Integrals and Riemann Sums

Introduction

We generalize the idea of Riemann sum to functions of two variables and give a module which calculates Riemann sums and also gives a graphical illustration. We use  Mathematica to evaluate iterated (double) integrals. Such integrals are interpreted as limits of Riemann sums and are applied to finding volumes bounded by surfaces of the form z=f(x,y).

6.1 Riemann Sums and Their Three Dimensional Representations

Suppose we have a function of two variables z=f(x,y)≥0, defined on a rectangle R=[(Subscript[a, 1]),(Subscript[b, 1])]×[(Subscript[a, 2]),(Subscript[b, 2])], then it is meaningful to talk about the volume bounded by the surface z, the xy-plane and the planes x=(Subscript[a, 1]),y=(Subscript[a, 2]),x=(Subscript[b, 1]) and y=(Subscript[b, 2]), (i.e., the volume beneath the surface and above the xy-plane). Just as we did for functions of one variable, we can approximate such volumes by Riemann sums. Let us subdivide each interval into n subintervals of equal width, Δx=((Subscript[b, 1])-(Subscript[a, 1]))/n  and Δy=((Subscript[b, 2])-(Subscript[a, 2]))/n. These subdivisions partition our rectangle into (n^2) equal subrectangles. Next we choose at a random sample point in each subrectangle, denoting ((Subscript[x, ij]^( * )),(Subscript[y, ij]^( * ))) the point in the (i,j)-th subrectangle. A Riemann sum is a sum of the form
		(∑_(i=1)^n)(∑_(j=1)^n)f((Subscript[x, ij]^( * )),(Subscript[y, ij]^( * )))Δx Δy.
If the function z=f(x,y) is continuous and n is sufficiently large, such a sum will give a good approximation to the volume in questione. This leads to the definition of the double integral of f over the given rectangle (even when f is not non-negative) as
		∫(∫_R) f(x,y)d A=(lim_(n->∞))(∑_(i=1)^n)(∑_(j=1)^n)f((Subscript[x, ij]^( * )),(Subscript[y, ij]^( * )))Δx Δy.

Let us illustrate the notion of Riemann sum with a specific function, say f(x,y)=4-(x^2)-(y^2) which satisfies f(x,y)≥0 on the square [0,1]×[0,1] in the xy-plane. The graph z=f(x,y) is a paraboloid of revolution, though that is hardly obvious from the plot over the given region.

Clear[f]
f[x_,y_]:=4-x^2-y^2
Plot3D[f[x,y],{x,0,1},{y,0,1}]

Below is a module, showRiemannSum, which enables us to calculate and visualize the Riemann sum approximations to a function of two variables, over a rectangular region of integration. There are six input parameters: the function f, two two-dimensional lists {a1,b1} and {a2,b2} representing the rectangle [(Subscript[a, 1]),(Subscript[b, 1])]×[(Subscript[a, 2]),(Subscript[b, 2])], and an integer max representing the maximum value of n that we want to observe. The module first sets n=max, evaluates the Riemann sum with randomly generated sample points and (max^2) subrectangles and displays a solid (explained below) whose volume equals the value of the Riemann sum. The output is dynamic and allows to change the value of n between 1 and max.

showRiemannSum[f_,{a1_,b1_},{a2_,b2_},max_]:=Module[
{g,n,sampleX,sampleY,Δx,Δy,p1,p2,riemannSum},
g[x_,y_]=f[x,y];
Manipulate[
Δx=(b1-a1)/n;Δy=(b2-a2)/n;
sampleX=Table[a1+i*Δx+RandomInteger[{1,7}]*Δx/8,{i,0,n-1},{j,0,n-1}];
sampleY=Table[a2+j*Δy+RandomInteger[{1,7}]*Δy/8,{i,0,n-1},{j,0,n-1}];
p1=Table[
Graphics3D[
{PointSize[.02],
Point[{sampleX[[i]][[j]],sampleY[[i]][[j]],g[sampleX[[i]][[j]],sampleY[[i]][[j]]]}]
}],
{i,1,n},{j,1,n}];
p2=Table[
Graphics3D[
{Boxed->False,
Cuboid[{a1+(i-1)*Δx,a2+(j-1)*Δy,0},{a1+i*Δx,a2+j*Δy,g[sampleX[[i]][[j]],sampleY[[i]][[j]]]}]}],
{i,1,n},{j,1,n}];
riemannSum=N[Sum[g[sampleX[[i]][[j]],sampleY[[i]][[j]]],{i,1,n},{j,1,n}]*Δx*Δy,10];
Show[p1,p2,PlotLabel->Style[riemannSum,Blue]],
{{n,max,"n"},1,max,1,Appearance->"Labeled"}]]

For example, to calculate Riemann sums for f(x,y)=4-(x^2)-(y^2), obtained by subdividing the interval [0,1]×[0,1] into 1,4,...,10×10=100 equal squares, we may use the following command.

showRiemannSum[f,{0,1},{0,1},10]

By increasing the number of rectangles, we should obtain a better approximation:

showRiemannSum[f,{0,1},{0,1},25]

What are we actually finding when we form the Riemann sum for a function of two variables f? Let us take a closer look at the terms f((Subscript[x, ij]^( * )),(Subscript[y, ij]^( * )))Δx Δy of the Riemann sum. The product ΔA =Δx Δy is the common area of all the subrectangles into which we have broken the original rectangular region [(Subscript[a, 1]),(Subscript[b, 1])]×[(Subscript[a, 2]),(Subscript[b, 2])]. The value f((Subscript[x, ij]^( * )),(Subscript[y, ij]^( * ))) is the height of the point on the graph z=f(x,y) at the sample point ((Subscript[x, ij]^( * )),(Subscript[y, ij]^( * ))) belonging to the (i,j)-th rectangle. Hence, the term f((Subscript[x, ij]^( * )),(Subscript[y, ij]^( * )))Δx Δy is the volume of a rectangular box (or what Mathematica calls a cuboid)that has the (i,j)-th subrectangle as its base and has height equal to f((Subscript[x, ij]^( * )),(Subscript[y, ij]^( * ))). The full Riemann sum is the sum of the volumes of (n^2) such cuboids. Thus the Riemann sum represents the volume of the solid obtained by taking all (n^2) (non-intersecting) rectangular boxes and combining them into one polyhedron. It is this polyhedron that the module showRiemannSum displays. (The black points on the top surface of the polyhedron represent the points ((Subscript[x, ij]^( * )),(Subscript[y, ij]^( * )),f((Subscript[x, ij]^( * )),(Subscript[y, ij]^( * )))) on the graph z=f(x,y) that determine the heights of the cuboids.)   

As another example, let us take a look at the function f(x,y)=(Sqrt[2-(x^2)-(y^2)]) bounded by the planes x=±1, y=±1:

Clear[f]
f[x_,y_]:=Sqrt[2-x^2-y^2]

showRiemannSum[f,{-1,1},{-1,1},6]

The function whose Riemann sum approximation we are giving is the top part of the sphere (x^2)+(y^2)+(z^2)=2, over the rectangular region [-1,1]×[-1,1]. The above approximation used 6×6=36 rectangular boxes. By increasing the number of rectangles (say, to 30×30=900 rectangles), we should obtain a better approximation.

showRiemannSum[f,{-1,1},{-1,1},30]

Ultimately, we see that as n increases without limit, the cuboids become "good approximations" to the solid below the graph of the function f(x,y). In particular, the Riemann sums (the volume of all the rectangular boxes added together), should give good approximations to the volume beneath the surface. 

Next, let us try the same thing with the function g(x,y)=ln((x^2)+ (y^2))/x on the rectangle [1/2,2]×[1/2,2]. This function is not strictly positive on this rectangle.

Clear[g]
g[x_,y_]:= Log[x^2+ y^2]/x

showRiemannSum[g,{1/2,2},{1/2,2},20]

In this case, the total volume of the boxes below the xy-plane (they correspond to negative sample values) appears in the Riemann sum with negative sign, so that volume is subtracted from the volume of the boxes above the xy-plane (they correspond to positive sample values).

Exercises

Exercise 1

Use the module showRiemannSum above, to obtain Riemann sum approximations to the volumes determined by the following functions. In each case, start with n=100 and adjust this value up or down until the values of the Riemann sums vary by not more than 2%. Once you have found an acceptable value of n, use the dynamic feature of showRiemannSum to compute 5 different Riemann sums with that value of n and approximate the double integral ∫(∫_R)f(x,y) dA by their average. You are to turn in only one plot, but you should write as your answer the values of all five Riemann sums and of their average. 

(a) f(x,y)=x sin(x y), on the rectangle [0,π]×[0,π]
(b) f(x,y)=(y^2)(e^(-(x+y))), on the rectangle [0,1]×[0,1]. 
(c) f(x,y)=y+sin(x y), on the rectangle [0,π]×[0,π].
(d) f(x,y)=(1/5)sin(x+y) (e^(Sqrt[x+y])), on the rectangle [0,π]×[0,π]. 
(e) f(x,y)=y(x^2)cosx+xsin y, on the rectangle [1,2]×[0,3].
(f) f(x,y)=(x^3)(y^2)+x y, on the rectangle [0,1]×[1,2].

Credits

These laboratories were created by Raouf Boules, Geoffrey Goodson, Ohoe Kim, Alexei Kolesnikov, and Michael O'Leary for use in the Calculus courses of Towson University. Commercial use is prohibited without permission of the authors. Non-commercial use is permitted, provided this credit section is retained.