. Let us explore 
where 
is the rectangle ![Typesetting:-delayCrossProduct([0, 2], [-1, 1])](images/double_Riemann_4.gif){kind=small}
.
This worksheet explores double Riemann sums in Maple. Although Maple does have a built-in command for approximating double integrals (in the "Student[MultivariateCalculus]" package there is the function "ApproximateInt") but I like my function better.
Maple worksheet: double_Riemann.mw
| > | ? student,multivariatecalculus,approximateint |
⇓Click on the triangle below to expand the initialization code. Then execute the code to load it into memory.
Initialization
Don't forget to evaluate these commands or the rest of the sheet will not work properly.
Example: . Let us explore where is the rectangle .
| > | f := (x,y) -> 4-x^2-y^2:
'f(x,y)' = f(x,y); |
 |
(1) |
The actual integral:
| > | Int(Int(f(x,y),x=0..2),y=-1..1) = int(int(f(x,y),x=0..2),y=-1..1); |
 |
(2) |
Here is the Riemann sum using a 4 by 5 grid of rectangles. We evaluate our function in the lower righthand side of each sub-rectangle. This amounts to "[Right,Left]" In general given "[A,B]", A = Left corresponds to lefthand side, A = Right to righthand side, B = Left corresponds to bottom, and B = Right corresponds to top.
| > | doubleRiemann(f,x=0..2,y=-1..1,4,5,[Right,Left],Value); |
 |
(3) |
The corresponding (double) Riemann sum is:
| > | doubleRiemann(f,x=0..2,y=-1..1,4,5,[Right,Left],RiemannSum); |
 |
(4) |
Finally, we have the following plot:
| > | display(doubleRiemann(f,x=0..2,y=-1..1,4,5,[Right,Left],Plot,0,0.75),viewpoint=circleleft); |
 |
Same set up using midpoint rule and a 3 by 3 grid.
| > | doubleRiemann(f,x=0..2,y=-1..1,4,5,[Midpoint,Midpoint],RiemannSum);
display(doubleRiemann(f,x=0..2,y=-1..1,4,5,[Midpoint,Midpoint],Plot,0,0.75),viewpoint=circleleft); |
 |
|
 |
Example: 
. Let us explore 
where 
is the rectangle ![Typesetting:-delayCrossProduct([-1, 3], [-2, 2])](images/double_Riemann_15.gif){kind=small}
.
| > | f := (x,y) -> x^2*y-2*y+3:
'f(x,y)' = f(x,y); |
 |
(5) |
The actual integral:
| > | Int(Int(f(x,y),x=-1..3),y=-2..2) = int(int(f(x,y),x=-1..3),y=-2..2); |
 |
(6) |
Here are some Riemann sum using a 2 by 2 grid of rectangles. First, left-top evaluation...
| > | doubleRiemann(f,x=-1..3,y=-2..2,2,2,[Left,Right],RiemannSum);
display(doubleRiemann(f,x=-1..3,y=-2..2,2,2,[Left,Right],Plot,0,0.75),viewpoint=circleleft); |
 |
|
 |
Next, left-midpoint...
| > | doubleRiemann(f,x=-1..3,y=-2..2,2,2,[Left,Midpoint],RiemannSum);
display(doubleRiemann(f,x=-1..3,y=-2..2,2,2,[Left,Midpoint],Plot,0,0.75),viewpoint=circleleft); |
 |
|
 |
Now, left-bottom...
| > | doubleRiemann(f,x=-1..3,y=-2..2,2,2,[Left,Left],RiemannSum);
display(doubleRiemann(f,x=-1..3,y=-2..2,2,2,[Left,Left],Plot,0,0.75),viewpoint=circleleft); |
 |
|
 |
Now let's see how increasing the number of boxes increases accuracy... [I'll continue to use left-top evaluation.]
| > | doubleRiemann(f,x=-1..3,y=-2..2,1,1,[Left,Right],RiemannSum);
display(doubleRiemann(f,x=-1..3,y=-2..2,1,1,[Left,Right],Plot,0,0.75),viewpoint=circleleft); doubleRiemann(f,x=-1..3,y=-2..2,2,2,[Left,Right],RiemannSum); display(doubleRiemann(f,x=-1..3,y=-2..2,2,2,[Left,Right],Plot,0,0.75),viewpoint=circleleft); doubleRiemann(f,x=-1..3,y=-2..2,4,4,[Left,Right],RiemannSum); display(doubleRiemann(f,x=-1..3,y=-2..2,4,4,[Left,Right],Plot,0,0.75),viewpoint=circleleft); doubleRiemann(f,x=-1..3,y=-2..2,8,8,[Left,Right],RiemannSum); display(doubleRiemann(f,x=-1..3,y=-2..2,8,8,[Left,Right],Plot,0,0.75),viewpoint=circleleft); doubleRiemann(f,x=-1..3,y=-2..2,16,16,[Left,Right],RiemannSum); display(doubleRiemann(f,x=-1..3,y=-2..2,16,16,[Left,Right],Plot,0,0.75),viewpoint=circleleft); |
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
Example: Use the "doubleRiemann" function defined in the initialiazation above to approximate the area under 
over the rectangle ![R = Typesetting:-delayCrossProduct([-4, 4], [0, 3])](images/double_Riemann_35.gif){kind=small}
.
| > | f := (x,y) -> cos(sqrt(x^2+2*y^2))+1:
'f(x,y)' = f(x,y); |
 |
(7) |
First, let's create a plot of the surface over 
using plot3d options "axes=boxed" and "color=x^2+y^2".
| > | plot3d(f(x,y),x=-4..4,y=0..3,axes=boxed,color=x^2+y^2,viewpoint=circleleft); |
 |
Now let's approximate the volume under the surface using midpoint rule and 10 x 20 rectangles using the function "doubleRiemann". Then we'll make the function output the approximation with corresponding plot and then make it output the corresponding double summation.
| > | display(doubleRiemann(f,x=-4..4,y=0..3,10,20,[Midpoint,Midpoint],Plot,0.5,0),viewpoint=circleleft); doubleRiemann(f,x=-4..4,y=0..3,10,20,[Midpoint,Midpoint],RiemannSum); |
 |
|
 |
(8) |
Now let's repeat the last example except we'll use a lower-righthand rule (in each subrectangle evaluate at the lower-right corner). Note: This means right sides of x-intervals and left sides of y-intervals.
| > | display(doubleRiemann(f,x=-4..4,y=0..3,10,20,[Right,Left],Plot,0.5,0),viewpoint=circleleft); doubleRiemann(f,x=-4..4,y=0..3,10,20,[Right,Left],RiemannSum); |
 |
|
 |
(9) |
Example: Let's approximate the volume under the part of the graph of 
(defined below) which lies above the rectangle 
and 
using a 4 by 3 grid of rectangles and evaluating subrectangles in the lower lefthand corner.
| > | f := (x,y) -> x^3+2*y:
'f(x,y)' = f(x,y); display(doubleRiemann(f,x=0..1,y=0..2,4,3,[Left,Left]),orientation=[-125,80,-170],viewpoint=circleleft); |
 |
|
 |
For the above function, double Riemann sums evaluated in lower lefthand corners always give underestimate. This happens because (in one dimension) the left hand rule gives underestimates when a function is increasing. This function is increasing in both the x and y-axis directions so we get an underestimate.
Example: Let's approximate the volume under the part of the graph of 
(defined below) which lies above the rectangle 
and 
using a 5 by 7 grid of rectangles and evaluating the function according to the midpoint rule.
| > | g := (x,y) -> 9-x^2-y^2:
'g(x,y)' = g(x,y); display(doubleRiemann(g,x=-1..1,y=-1..1,5,7,[Midpoint,Midpoint],Plot,0.75,0),viewpoint=circleleft); |
 |
|
 |
Some underestimates...
| > | doubleRiemann(g,x=-1..1,y=-1..1,3,10,[Left,Left],Value);
doubleRiemann(g,x=-1..1,y=-1..1,3,10,[Left,Right],Value); doubleRiemann(g,x=-1..1,y=-1..1,3,10,[Right,Left],Value); doubleRiemann(g,x=-1..1,y=-1..1,3,10,[Right,Right],Value); doubleRiemann(g,x=-1..1,y=-1..1,3,10,[Right,Midpoint],Value); doubleRiemann(g,x=-1..1,y=-1..1,3,10,[Left,Midpoint],Value); |
 |
|
 |
|
 |
|
 |
|
 |
|
 |
(10) |
A few overestimates...
| > | doubleRiemann(g,x=-1..1,y=-1..1,3,10,[Midpoint,Right],Value);
doubleRiemann(g,x=-1..1,y=-1..1,3,10,[Midpoint,Left],Value); doubleRiemann(g,x=-1..1,y=-1..1,10,3,[Right,Midpoint],Value); doubleRiemann(g,x=-1..1,y=-1..1,10,3,[Left,Midpoint],Value); |
 |
|
 |
|
 |
|
 |
(11) |
For the above function, double Riemann sums evaluated using midpoint rule always give overestimates. Why? Recall that (in one dimension) the midpoint rule gives overestimates for concave down functions. This function is concave down in both the x and y-axis directions so we get an overestimate. Now if we switch to left and right hand evaluations, there is no reason to believe this will still work. Left and right hand rule approximations are over or underestimates when functions are only increasing or decreasing. If a function goes up and down, there is no guarantee that we would have an over or underestimate when using right and left hand evaluations. However, because of symmetry it happens that the evaluations built from mixtures of left and right hand rules will always give underestimates. Mixtures with midpoints can go either way depending on the chosen grid size.