Double Riemann Sums 

 

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. 

The Code

 

 

Example:  f(x, y) = VectorCalculus:-`+`(VectorCalculus:-`+`(4, VectorCalculus:-`-`(`*`(`^`(x, 2)))), VectorCalculus:-`-`(`*`(`^`(y, 2)))).  Let us explore  where R is the rectangle Typesetting:-delayCrossProduct([0, 2], [-1, 1]). 

 

> f := (x,y) -> 4-x^2-y^2:
'f(x,y)' = f(x,y);
 

f(x, y) = `+`(4, `-`(`*`(`^`(x, 2))), `-`(`*`(`^`(y, 2)))) (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);
 

Int(Int(`+`(4, `-`(`*`(`^`(x, 2))), `-`(`*`(`^`(y, 2)))), x = 0 .. 2), y = -1 .. 1) = `/`(28, 3) (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);
 

7.060000000 (3)
 

 

The corresponding (double) Riemann sum is: 

 

> doubleRiemann(f,x=0..2,y=-1..1,4,5,[Right,Left],RiemannSum);
 

Sum(Sum(`+`(`/`(4, 5), `-`(`*`(`/`(1, 5), `*`(`^`(`+`(`*`(`/`(1, 2), `*`(i)), `/`(1, 2)), 2)))), `-`(`*`(`/`(1, 5), `*`(`^`(`+`(`-`(1), `*`(`/`(2, 5), `*`(j))), 2))))), i = 0 .. 3), j = 0 .. 4) (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);
 

Plot_2d
 

 

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);
 

 

Sum(Sum(`+`(`/`(4, 5), `-`(`*`(`/`(1, 5), `*`(`^`(`+`(`*`(`/`(1, 2), `*`(i)), `/`(1, 4)), 2)))), `-`(`*`(`/`(1, 5), `*`(`^`(`+`(`-`(`/`(4, 5)), `*`(`/`(2, 5), `*`(j))), 2))))), i = 0 .. 3), j = 0 .. 4...
Plot_2d
 

 

Example:  f(x, y) = VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(`*`(`^`(x, 2)), y), VectorCalculus:-`-`(VectorCalculus:-`*`(2, y))), 3).  Let us explore  where R is the rectangle Typesetting:-delayCrossProduct([-1, 3], [-2, 2]). 

 

> f := (x,y) -> x^2*y-2*y+3:
'f(x,y)' = f(x,y);
 

f(x, y) = `+`(`*`(`^`(x, 2), `*`(y)), `-`(`*`(2, `*`(y))), 3) (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);
 

Int(Int(`+`(`*`(`^`(x, 2), `*`(y)), `-`(`*`(2, `*`(y))), 3), x = -1 .. 3), y = -2 .. 2) = 48 (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);
 

 

Sum(Sum(`+`(`*`(8, `*`(`^`(`+`(`-`(1), `*`(2, `*`(i))), 2), `*`(j))), `-`(`*`(16, `*`(j))), 12), i = 0 .. 1), j = 0 .. 1)
Plot_2d
 

 

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);
 

 

Sum(Sum(`+`(`*`(4, `*`(`^`(`+`(`-`(1), `*`(2, `*`(i))), 2), `*`(`+`(`-`(1), `*`(2, `*`(j)))))), 20, `-`(`*`(16, `*`(j)))), i = 0 .. 1), j = 0 .. 1)
Plot_2d
 

 

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);
 

 

Sum(Sum(`+`(`*`(4, `*`(`^`(`+`(`-`(1), `*`(2, `*`(i))), 2), `*`(`+`(`-`(2), `*`(2, `*`(j)))))), 28, `-`(`*`(16, `*`(j)))), i = 0 .. 1), j = 0 .. 1)
Plot_2d
 

 

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);
 

 

 

 

 

 

 

 

 

 

Sum(Sum(`+`(`*`(16, `*`(`^`(`+`(`-`(1), `*`(4, `*`(i))), 2), `*`(`+`(2, `*`(4, `*`(j)))))), `-`(16), `-`(`*`(128, `*`(j)))), i = 0 .. 0), j = 0 .. 0)
Plot_2d
Sum(Sum(`+`(`*`(8, `*`(`^`(`+`(`-`(1), `*`(2, `*`(i))), 2), `*`(j))), `-`(`*`(16, `*`(j))), 12), i = 0 .. 1), j = 0 .. 1)
Plot_2d
Sum(Sum(`+`(`*`(`^`(`+`(`-`(1), i), 2), `*`(`+`(`-`(1), j))), 5, `-`(`*`(2, `*`(j)))), i = 0 .. 3), j = 0 .. 3)
Plot_2d
Sum(Sum(`+`(`*`(`/`(1, 4), `*`(`^`(`+`(`-`(1), `*`(`/`(1, 2), `*`(i))), 2), `*`(`+`(`-`(`/`(3, 2)), `*`(`/`(1, 2), `*`(j)))))), `/`(3, 2), `-`(`*`(`/`(1, 4), `*`(j)))), i = 0 .. 7), j = 0 .. 7)
Plot_2d
Sum(Sum(`+`(`*`(`/`(1, 16), `*`(`^`(`+`(`-`(1), `*`(`/`(1, 4), `*`(i))), 2), `*`(`+`(`-`(`/`(7, 4)), `*`(`/`(1, 4), `*`(j)))))), `/`(13, 32), `-`(`*`(`/`(1, 32), `*`(j)))), i = 0 .. 15), j = 0 .. 15)
Plot_2d
 

 

Example: Use the "doubleRiemann" function defined in the initialiazation above to approximate the area under f(x, y) = VectorCalculus:-`+`(cos(sqrt(VectorCalculus:-`+`(`*`(`^`(x, 2)), VectorCalculus:-`*`(2, `*`(`^`(y, 2)))))), 1) over the rectangle R = Typesetting:-delayCrossProduct([-4, 4], [0, 3]).  

 

> f := (x,y) -> cos(sqrt(x^2+2*y^2))+1:
'f(x,y)' = f(x,y);
 

f(x, y) = `+`(cos(`*`(`^`(`+`(`*`(`^`(x, 2)), `*`(2, `*`(`^`(y, 2)))), `/`(1, 2)))), 1) (7)
 

 

First, let's create a plot of the surface over R 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);
 

Plot_2d
 

 

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);
 

 

Plot_2d
Sum(Sum(`+`(`*`(`/`(3, 25), `*`(cos(`+`(`*`(`/`(1, 40), `*`(`^`(`+`(20754, `-`(`*`(9216, `*`(i))), `*`(1024, `*`(`^`(i, 2))), `*`(72, `*`(`^`(j, 2))), `*`(72, `*`(j))), `/`(1, 2)))))))), `/`(3, 25)), ... (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);
 

 

Plot_2d
Sum(Sum(`+`(`*`(`/`(3, 25), `*`(cos(`+`(`*`(`/`(1, 20), `*`(`^`(`+`(4096, `-`(`*`(2048, `*`(i))), `*`(256, `*`(`^`(i, 2))), `*`(18, `*`(`^`(j, 2)))), `/`(1, 2)))))))), `/`(3, 25)), i = 0 .. 9), j = 0 ... (9)
 

 

Example: Let's approximate the volume under the part of the graph of  z = f(x, y) (defined below) which lies above the rectangle `and`(`<=`(0, x), `<=`(x, 1)) and `and`(`<=`(0, y), `<=`(y, 2)) 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);
 

 

f(x, y) = `+`(`*`(`^`(x, 3)), `*`(2, `*`(y)))
Plot_2d
 

 

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  z = g(x, y) (defined below) which lies above the rectangle `and`(`<=`(-1, x), `<=`(x, 1)) and `and`(`<=`(-1, y), `<=`(y, 1)) 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);
 

 

g(x, y) = `+`(9, `-`(`*`(`^`(x, 2))), `-`(`*`(`^`(y, 2))))
Plot_2d
 

 

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);
 

 

 

 

 

 

33.01037037
33.01037037
33.01037037
33.01037037
33.05037037
33.05037037 (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);
 

 

 

 

33.45481482
33.45481482
33.45481483
33.45481483 (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.