The Double Riemann Sum on a 16 by 16 Grid

Let us repeat the steps of last section on a 16 by 16 grid. We split the x and y intervals from 0 to 1 into 16 equal intervals. Each interval will have a length of 1/16. We will use the midpoint rule so that x* and y* will be in the middle of the intervals. The sixteen midpoints are given as follows:

> midpts:=seq(i/16-1/32,i=1..16);

Let's make a vector "xvec" out of these points to represent x*, and a vector "yvec" to represent y*.

> xvec:=vector([midpts]);

> yvec:=vector([midpts]);

We now define a function "fxy", which takes two integers "m" and "n" as input, and evaluates our function f(x,y) at the value x=xvec[m] and y=yvec[n]:

> fxy:=(m,n)->evalf(f(xvec[m],yvec[n]));

This function is used to set up a 16 by 16 matrix containing all the values f(x*,y*):

> fmat:=matrix(16,16,fxy):

Note: I've used a colon to suppress the output. Nobody wants to see a 16 by 16 matrix of floating-point numbers.

We can plot this matrix as a 3-D histogram using the following command.

> matrixplot(fmat,heights=histogram,axes=frame,orientation=[70,75]);

If you compare this histogram to the plot of the function two sections back, you'll see that the approximation now much better.

Now let's evaluate the Double Riemann sum, by summing all the entries in the matrix and multipling by , which in this case is 1/(16^2).

> sum(sum(fmat[i,j],i=1..16),j=1..16)/16^2;

Let's compare this to the integral obtained in two sections back:

> intf-%;

The approximation is excellent!

>