Double Riemann Sums

Code

generatePoints[p_,deltaR_,deltaTheta_,rPoint_,thetaPoint_]:=
(*Generates a list of points to use with ListPlot3D for one Riemann summand*)
(
innerR=p[[1]]-deltaR*rPoint;
outerR=p[[1]]+deltaR*(1-rPoint);
innerTheta=p[[2]]-deltaTheta*thetaPoint;
outerTheta=p[[2]]+deltaTheta*(1-thetaPoint);
innerPoints=Append[Table[{innerR Cos[theta],innerR Sin[theta],p[[3]]},{theta,innerTheta,outerTheta,0.1}],{innerR Cos[outerTheta],innerR Sin[outerTheta],p[[3]]}];
outerPoints=Append[Table[{outerR Cos[theta],outerR Sin[theta],p[[3]]},{theta,innerTheta,outerTheta,0.1}],{outerR Cos[outerTheta],outerR Sin[outerTheta],p[[3]]}];
Join[innerPoints,outerPoints]
)

graphicsInitial := (
doubleIntegral=NIntegrate[f[r,theta]*r,{r,rLeft,rRight},{theta,thetaLeft,thetaRight}];
surface=If[showSurface,ParametricPlot3D[{r Cos[theta],r Sin[theta],f[r,theta]},{r,rLeft,rRight},{theta,thetaLeft,thetaRight},Mesh->None,PlotStyle->{Red,Opacity[0.5]}],{}];
)

graphicsFinal[rSubdivisions_,thetaSubdivisions_]:=(
deltaR=(rRight-rLeft)/rSubdivisions;
deltaTheta=(thetaRight-thetaLeft)/thetaSubdivisions;
cylindricalPoints=Join@@Table[{r,theta,f[r,theta]},{r,rLeft+deltaR*rPoint,rRight,deltaR},{theta,thetaLeft+deltaTheta*thetaPoint,thetaRight,deltaTheta}];
points={#[[1]]Cos[#[[2]]],#[[1]]Sin[#[[2]]],#[[3]]}&/@cylindricalPoints;
riemannSum=Total[#[[3]]*#[[1]]&/@cylindricalPoints]*deltaR*deltaTheta;

riemannPlot=If[showBoxes,ListPlot3D[generatePoints[#,deltaR,deltaTheta,rPoint,thetaPoint],Filling->Axis,FillingStyle->{Hue[2/3,1,.5,fillingOpacity]},Mesh->None,PlotStyle->{LightBlue,Opacity[0.5]}]&/@cylindricalPoints,{}];
pointsPlot=If[showPoints,Graphics3D[{Green,PointSize[0.02],Point[points]}],{}];

label=If[showNumericalResults,Style[Row[{"Riemann Sum: ",NumberForm[riemannSum,{5,3}],"\nDouble Integral: ",NumberForm[doubleIntegral,{5,3}]}],FontSize->16],""];

Show[riemannPlot,surface,pointsPlot,Axes->Automatic,AxesLabel->{x,y,z},BoxRatios->{2,2,3},ViewPoint->{20,4,8},PlotRange->plotRange,PlotLabel->label,Ticks->axisTicks])

graphics[rSubdivisions_,thetaSubdivisions_]:=(
graphicsInitial;
graphicsFinal[rSubdivisions,thetaSubdivisions]
)

Graphics

Modify the function and options below to change the Riemann sum.
rPoint specifies where in each r interval to take the sample point, and similarly for thetaPoint and θ.  Taking “midpoint” sample points corresponds to using rPoint = 0.5 and thetaPoint = 0.5.
showPoints specifies whether to draw the sample points on the surface.

f[r_,theta_]:= r^2-3r Cos[theta];
rSubdivisions=6;
thetaSubdivisions=6;
rPoint=0.5;
thetaPoint=0.5;
rLeft=0;
rRight=1;
thetaLeft=Pi/2;
thetaRight=Pi;
showSurface=True;
showBoxes=True;
showPoints=True;

(*Display parameters*)
showNumericalResults=True;
fillingOpacity=0.2;
axisTicks=Automatic;
plotRange=Automatic;


graphics[rSubdivisions,thetaSubdivisions]

Animation

This animation shows the polar Riemann sum as the number of subdivisions in both r and θ increases from startingSubdivions to endingSubdivisions.
All options except rSubdivisions and thetaSubdivions are taken from the Graphics section above.
It may take some time to run this command, as all frames of the animation are precomputed so that it can run smoothly.

startingSubdivisions=1;
endingSubdivisions=16;

frames=Table[Image[graphics[t,t]],{t,startingSubdivisions,endingSubdivisions,1}];
ListAnimate[frames]