Double Riemann SumsThis 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.? student,multivariatecalculus,approximateint\342\207\223Click 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.restart;
with(plots):
with(plottools):
with(VectorCalculus):
#
# doubleRiemann computes Riemann sums for double integrals.
#
# f is a function Ex: f might be "(x,y) -> x^2 + y^2;"
# X and Y are variables with range Ex: X might be "x=0..4"
# The net volume under f restricted to the rectangle defined
# by X and Y's ranges is approximated using m * n rectangles
# (where m and n are positive integers).
#
# EvalType is a list describing the evaluation scheme. The
# default setting is midpoint rule. Ex: EvalType could be
# something like "[Left,Right]" in this case rectangles would
# be evaluated in their upper left hand corners.
# OutputType can be "Value" "RiemannSum" or "Plot". "Value"
# causes the function to return the riemann sum's value.
# "RiemannSum" causes the function to return an unevaluated
# double sum. "Plot" returns a plot of the rectangular boxes
# and the surface. The plot title states the exact and approximate
# values of the double integral.
# T and Tp are transparency settings for the boxes and graph
# respectively. 0 means opaque and 1 means totally transparent.
#
doubleRiemann := proc(f,X,Y,m,n,EvalType := [Midpoint,Midpoint],OutputType := Plot, T := 0, Tp := 0.95)
local i,j,a,b,c,d,xV,yV,Dx,Dy,xS,yS,x,y,rSum,approxVal,actualVal,pTitle,blocks,graph,points;
# grab variable names and range data: "x=a..b & y=c..d"
x := lhs(X); y := lhs(Y);
a := lhs(rhs(X)); b := rhs(rhs(X));
c := lhs(rhs(Y)); d := rhs(rhs(Y));
# compute delta x and delta y (lengths and widths of subrectangles).
Dx := (b-a)/m; Dy := (d-c)/n;
# x & y coords of box bases.
xV := [seq(a+i*Dx,i=0..m-1)];
yV := [seq(c+j*Dy,j=0..n-1)];
# Set delta x according to mid, right or left hand rule
if EvalType[1] = Midpoint then
xS := Dx/2;
elif EvalType[1] = Right then
xS := Dx;
elif EvalType[1] = Left then
xS := 0;
else
RETURN(cat("Invalid evaluation type: ",EvalType[1]));
end if:
# Set delta y according to mid, right or left hand rule
if EvalType[2] = Midpoint then
yS := Dy/2;
elif EvalType[2] = Right then
yS := Dy;
elif EvalType[2] = Left then
yS := 0;
else
RETURN(cat("Invalid evaluation type: ",EvalType[2]));
end if:
# return an unevaluated Riemann sum
rSum := Sum(Sum(f(a+i*Dx+xS,c+j*Dy+yS)*Dx*Dy,i=0..nops(xV)-1),j=0..nops(yV)-1);
if OutputType = RiemannSum then RETURN(rSum); end if:
# compute the Riemann sum
approxVal := sum(sum(evalf(f(xV[i]+xS,yV[j]+yS)*Dx*Dy),i=1..nops(xV)),j=1..nops(yV));
if OutputType = Value then RETURN(approxVal); end if;
# compute the actual value
actualVal := evalf(int(int(f(x,y),x=a..b),y=c..d));
# create the plot title
pTitle := cat("Actual value = ",convert(actualVal,string),"\134n Approximate Value = ",convert(approxVal,string));
# create plots of the boxes
blocks := {seq(seq(cuboid([xV[i],yV[j],0],[xV[i]+Dx,yV[j]+Dy,f(xV[i]+xS,yV[j]+yS)],color=gray,transparency=T),
i=1..nops(xV)),j=1..nops(yV))}:
# create plots of the evaluation points
points := {seq(seq(point([xV[i]+xS,yV[j]+yS,f(xV[i]+xS,yV[j]+yS)],symbol=diamond,symbolsize=10,color=red,transparency=T),
i=1..nops(xV)),j=1..nops(yV))}:
# create a plot of the graph itself
graph := plot3d(f(x,y),x=a..b,y=c..d,color=blue,transparency=Tp):
display(blocks,graph,points,title=pTitle);
end proc:Example: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYsLUkjbWlHRiQ2JVEiZkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYsNiVRInhGJ0YvRjItSSNtb0dGJDYtUSIsRicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYsNiVRInlGJ0YvRjJGQUZBLUY+Ni1RIj1GJ0ZBRkMvRkdGRUZIRkpGTEZORlAvRlNRLDAuMjc3Nzc3OGVtRicvRlZGam4tSSNtbkdGJDYkUSI0RidGQS1GPjYtUSgmbWludXM7RidGQUZDRmhuRkhGSkZMRk5GUC9GU1EsMC4yMjIyMjIyZW1GJy9GVkZkby1JJW1zdXBHRiQ2JUY6LUYjNiUtRl1vNiRRIjJGJ0ZBRi9GMi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGYG8tRmdvNiVGWEZpb0ZecC1GLDYjUSFGJ0ZB. Let us explore LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzY2LUkjbW9HRiQ2LVErJkludGVncmFsO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNC8lKXN0cmV0Y2h5R1EldHJ1ZUYnLyUqc3ltbWV0cmljR0Y0LyUobGFyZ2VvcEdGOS8lLm1vdmFibGVsaW1pdHNHRjQvJSdhY2NlbnRHRjQvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZELUklbXN1YkdGJDYlRistRiM2JS1JI21pR0YkNiVRIlJGJy8lJ2l0YWxpY0dGOS9GMFEnaXRhbGljRidGUEZSLyUvc3Vic2NyaXB0c2hpZnRHUSIwRictRk02JVEiZkYnRlBGUi1JKG1mZW5jZWRHRiQ2JC1GIzYmLUZNNiVRInhGJ0ZQRlItRiw2LVEiLEYnRi9GMi9GNkY5L0Y4RjRGOi9GPUY0Rj5GQEZCL0ZGUSwwLjMzMzMzMzNlbUYnLUZNNiVRInlGJ0ZQRlJGL0YvLUYsNi1RIn5GJ0YvRjJGNUZgb0Y6RmFvRj5GQEZCRkUtRk02JVEjZHhGJ0ZQRlJGZ28tRk02JVEjZHlGJ0ZQRlJGZ28tRiw2LVEiPUYnRi9GMkY1RmBvRjpGYW9GPkZAL0ZDUSwwLjI3Nzc3NzhlbUYnL0ZGRmRwRmdvLUkobXN1YnN1cEdGJDYoRistRiM2Ji1GLDYtUSomdW1pbnVzMDtGJ0YvRjJGNUZgb0Y6RmFvRj5GQC9GQ1EsMC4yMjIyMjIyZW1GJy9GRkZfcS1JI21uR0YkNiRRIjFGJ0YvRlBGUi1GIzYlRmFxRlBGUi8lMXN1cGVyc2NyaXB0c2hpZnRHRlZGVC9JK21zZW1hbnRpY3NHRiRRLFtub25lLG5vbmVdRictRmdwNihGKy1GIzYlLUZicTYkRlZGL0ZQRlItRiM2JS1GYnE2JFEiMkYnRi9GUEZSRmdxRlRGaXFGV0ZaRmdvRmpvRmdvRl1wRi8= where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiUkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= is the rectangle LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkobWZlbmNlZEdGJDYmLUYjNiYtSSNtbkdGJDYkUSIwRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEiLEYnRjQvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHUSV0cnVlRicvJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRjE2JFEiMkYnRjRGNEY0LyUlb3BlbkdRIltGJy8lJmNsb3NlR1EiXUYnLUY4Ni1RKCZ0aW1lcztGJ0Y0RjsvRj9GPUZBRkNGRUZHRkkvRkxRLDAuMjIyMjIyMmVtRicvRk9GaW4tRiw2Ji1GIzYnLUY4Ni1RKiZ1bWludXMwO0YnRjRGO0ZnbkZBRkNGRUZHRklGaG5Gam4tRjE2JFEiMUYnRjRGN0Zib0Y0RjRGVEZXRjQ=.f := (x,y) -> 4-x^2-y^2:
'f(x,y)' = f(x,y);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);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);The corresponding (double) Riemann sum is:doubleRiemann(f,x=0..2,y=-1..1,4,5,[Right,Left],RiemannSum);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: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYuLUkjbWlHRiQ2JVEiZkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYsNiVRInhGJ0YvRjItSSNtb0dGJDYtUSIsRicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYsNiVRInlGJ0YvRjJGQUZBLUY+Ni1RIj1GJ0ZBRkMvRkdGRUZIRkpGTEZORlAvRlNRLDAuMjc3Nzc3OGVtRicvRlZGam4tSSVtc3VwR0YkNiVGOi1GIzYlLUkjbW5HRiQ2JFEiMkYnRkFGL0YyLyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJ0ZYLUY+Ni1RKCZtaW51cztGJ0ZBRkNGaG5GSEZKRkxGTkZQL0ZTUSwwLjIyMjIyMjJlbUYnL0ZWRlxwRmFvLUY+Ni1RIn5GJ0ZBRkNGaG5GSEZKRkxGTkZQRlIvRlZGVEZYLUY+Ni1RIitGJ0ZBRkNGaG5GSEZKRkxGTkZQRltwRl1wLUZibzYkUSIzRidGQUZB. Let us explore LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzY2LUkjbW9HRiQ2LVErJkludGVncmFsO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNC8lKXN0cmV0Y2h5R1EldHJ1ZUYnLyUqc3ltbWV0cmljR0Y0LyUobGFyZ2VvcEdGOS8lLm1vdmFibGVsaW1pdHNHRjQvJSdhY2NlbnRHRjQvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZELUklbXN1YkdGJDYlRistRiM2JS1JI21pR0YkNiVRIlJGJy8lJ2l0YWxpY0dGOS9GMFEnaXRhbGljRidGUEZSLyUvc3Vic2NyaXB0c2hpZnRHUSIwRictRk02JVEiZkYnRlBGUi1JKG1mZW5jZWRHRiQ2JC1GIzYmLUZNNiVRInhGJ0ZQRlItRiw2LVEiLEYnRi9GMi9GNkY5L0Y4RjRGOi9GPUY0Rj5GQEZCL0ZGUSwwLjMzMzMzMzNlbUYnLUZNNiVRInlGJ0ZQRlJGL0YvLUYsNi1RIn5GJ0YvRjJGNUZgb0Y6RmFvRj5GQEZCRkUtRk02JVEjZHhGJ0ZQRlJGZ28tRk02JVEjZHlGJ0ZQRlJGZ28tRiw2LVEiPUYnRi9GMkY1RmBvRjpGYW9GPkZAL0ZDUSwwLjI3Nzc3NzhlbUYnL0ZGRmRwRmdvLUkobXN1YnN1cEdGJDYoRistRiM2Ji1GLDYtUSomdW1pbnVzMDtGJ0YvRjJGNUZgb0Y6RmFvRj5GQC9GQ1EsMC4yMjIyMjIyZW1GJy9GRkZfcS1JI21uR0YkNiRRIjJGJ0YvRlBGUi1GIzYlRmFxRlBGUi8lMXN1cGVyc2NyaXB0c2hpZnRHRlZGVC9JK21zZW1hbnRpY3NHRiRRLFtub25lLG5vbmVdRictRmdwNihGKy1GIzYmRltxLUZicTYkUSIxRidGL0ZQRlItRiM2JS1GYnE2JFEiM0YnRi9GUEZSRmdxRlRGaXFGV0ZaRmdvRmpvRmdvRl1wRi8= where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiUkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= is the rectangle LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkobWZlbmNlZEdGJDYmLUYjNictSSNtb0dGJDYtUSomdW1pbnVzMDtGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRjkvJSlzdHJldGNoeUdGOS8lKnN5bW1ldHJpY0dGOS8lKGxhcmdlb3BHRjkvJS5tb3ZhYmxlbGltaXRzR0Y5LyUnYWNjZW50R0Y5LyUnbHNwYWNlR1EsMC4yMjIyMjIyZW1GJy8lJ3JzcGFjZUdGSC1JI21uR0YkNiRRIjFGJ0Y0LUYxNi1RIixGJ0Y0RjcvRjtRJXRydWVGJ0Y8Rj5GQEZCRkQvRkdRJjAuMGVtRicvRkpRLDAuMzMzMzMzM2VtRictRkw2JFEiM0YnRjRGNEY0LyUlb3BlbkdRIltGJy8lJmNsb3NlR1EiXUYnLUYxNi1RKCZ0aW1lcztGJ0Y0RjdGOkY8Rj5GQEZCRkRGRkZJLUYsNiYtRiM2J0YwLUZMNiRRIjJGJ0Y0Rk9GYm9GNEY0RmVuRmhuRjQ=.f := (x,y) -> x^2*y-2*y+3:
'f(x,y)' = f(x,y);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);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 LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYqLUkjbWlHRiQ2JVEiZkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYsNiVRInhGJ0YvRjItSSNtb0dGJDYtUSIsRicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYsNiVRInlGJ0YvRjJGQUZBLUY+Ni1RIj1GJ0ZBRkMvRkdGRUZIRkpGTEZORlAvRlNRLDAuMjc3Nzc3OGVtRicvRlZGam4tRiw2JVEkY29zRicvRjBGRUZBLUY2NiQtRiM2Ji1GLDYjUSFGJy1GIzYlRmRvLUkmbXNxcnRHRiQ2Iy1GIzYoLUklbXN1cEdGJDYlRjotRiM2JS1JI21uR0YkNiRRIjJGJ0ZBRi9GMi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRictRj42LVEiK0YnRkFGQ0ZobkZIRkpGTEZORlAvRlNRLDAuMjIyMjIyMmVtRicvRlZGXnFGY3AtRj42LVEifkYnRkFGQ0ZobkZIRkpGTEZORlBGUi9GVkZULUZfcDYlRlhGYXBGZ3BGQUZBRmRvRkFGQUZqcC1GZHA2JFEiMUYnRkFGQQ== over the rectangle LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYqLUkjbWlHRiQ2JVEiUkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1GNjYtUSI9RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjI3Nzc3NzhlbUYnL0ZORlNGNS1JKG1mZW5jZWRHRiQ2Ji1GIzYnLUY2Ni1RKiZ1bWludXMwO0YnRjlGO0Y+RkBGQkZERkZGSC9GS1EsMC4yMjIyMjIyZW1GJy9GTkZobi1JI21uR0YkNiRRIjRGJ0Y5LUY2Ni1RIixGJ0Y5RjsvRj9GMUZARkJGREZGRkhGSi9GTlEsMC4zMzMzMzMzZW1GJ0ZqbkY5RjkvJSVvcGVuR1EiW0YnLyUmY2xvc2VHUSJdRictRjY2LVEoJnRpbWVzO0YnRjlGO0Y+RkBGQkZERkZGSEZnbkZpbi1GVjYmLUYjNiYtRltvNiRRIjBGJ0Y5Rl5vLUZbbzYkUSIzRidGOUY5RjlGZG9GZ29GOQ==. f := (x,y) -> cos(sqrt(x^2+2*y^2))+1:
'f(x,y)' = f(x,y);First, let's create a plot of the surface over LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiUkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= 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);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);Example: Let's approximate the volume under the part of the graph of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2JVEiekYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1GLDYlUSJmRidGL0YyLUkobWZlbmNlZEdGJDYkLUYjNiYtRiw2JVEieEYnRi9GMi1GNjYtUSIsRidGOUY7L0Y/RjFGQEZCRkRGRkZIL0ZLUSYwLjBlbUYnL0ZOUSwwLjMzMzMzMzNlbUYnLUYsNiVRInlGJ0YvRjJGOUY5Rjk= (defined below) which lies above the rectangle LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbW5HRiQ2JFEiMEYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1JI21vR0YkNi1RJiZsZXE7RidGLy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGOC8lKXN0cmV0Y2h5R0Y4LyUqc3ltbWV0cmljR0Y4LyUobGFyZ2VvcEdGOC8lLm1vdmFibGVsaW1pdHNHRjgvJSdhY2NlbnRHRjgvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZHLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnL0YwUSdpdGFsaWNGJ0YyLUYsNiRRIjFGJ0YvRi8= and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbW5HRiQ2JFEiMEYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1JI21vR0YkNi1RJiZsZXE7RidGLy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGOC8lKXN0cmV0Y2h5R0Y4LyUqc3ltbWV0cmljR0Y4LyUobGFyZ2VvcEdGOC8lLm1vdmFibGVsaW1pdHNHRjgvJSdhY2NlbnRHRjgvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZHLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnL0YwUSdpdGFsaWNGJ0YyLUYsNiRRIjJGJ0YvRi8= 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 LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2JVEiekYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1GLDYlUSJnRidGL0YyLUkobWZlbmNlZEdGJDYkLUYjNiYtRiw2JVEieEYnRi9GMi1GNjYtUSIsRidGOUY7L0Y/RjFGQEZCRkRGRkZIL0ZLUSYwLjBlbUYnL0ZOUSwwLjMzMzMzMzNlbUYnLUYsNiVRInlGJ0YvRjJGOUY5Rjk= (defined below) which lies above the rectangle LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbW9HRiQ2LVEqJnVtaW51czA7RicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y0LyUpc3RyZXRjaHlHRjQvJSpzeW1tZXRyaWNHRjQvJShsYXJnZW9wR0Y0LyUubW92YWJsZWxpbWl0c0dGNC8lJ2FjY2VudEdGNC8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRkMtSSNtbkdGJDYkUSIxRidGLy1GLDYtUSYmbGVxO0YnRi9GMkY1RjdGOUY7Rj1GPy9GQlEsMC4yNzc3Nzc4ZW1GJy9GRUZOLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnL0YwUSdpdGFsaWNGJ0ZKRkZGLw== and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbW9HRiQ2LVEqJnVtaW51czA7RicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y0LyUpc3RyZXRjaHlHRjQvJSpzeW1tZXRyaWNHRjQvJShsYXJnZW9wR0Y0LyUubW92YWJsZWxpbWl0c0dGNC8lJ2FjY2VudEdGNC8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRkMtSSNtbkdGJDYkUSIxRidGLy1GLDYtUSYmbGVxO0YnRi9GMkY1RjdGOUY7Rj1GPy9GQlEsMC4yNzc3Nzc4ZW1GJy9GRUZOLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnL0YwUSdpdGFsaWNGJ0ZKRkZGLw== 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);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);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.