Double Riemann Sums Maple Code

Initialization

> 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),"\n 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: