An Elaborate Riemann Sum Demo
This Mathematica notebook is an extension of the work in the notebook RiemannSums.nb posted in the MouseOver of abstractmath.org.  It is produced by and is licensed under a MouseOver , which makes it freely available for others to use, provided they refer to this file. It is the source of the demo in the Gyre&Gimble post An Elaborate Riemann Sum Demo.

Riemann Package (opened for modification)

This package provides PlotRiemann, which shows Riemann sum of a given integral with various choices of parameters, and CloudList, which plots the values of a number of random Riemann sums of the given integral.  It is a work in progress. You are encouraged to improve this work and to use it in your own publications and coursework.

SetOptions[Plot,AspectRatio->Automatic];

PlotBounded[expression_,range_,options___]:=Module[{g,f,left,right,newleft,newright,var,leftheight,rightheight,width,newrange},var=range[[1]];
f=Function[Release[var],expression];
left=range[[2]];
right=range[[3]];
width=right-left;
leftheight=f[left];
rightheight=f[right];
newleft=left-0.25*width;
newright=right+0.25*width;
newrange={var,newleft,newright};
g:=Plot[expression,Release[newrange],DisplayFunction->Identity,PlotStyle->{{RGBColor[0,0,1]}},options];
Show[g,Graphics[{Line[{{left,0},{left,leftheight}}],Line[{{right,0},{right,rightheight}}]}],DisplayFunction->$DisplayFunction]]

RandomPartition[range_,mesh_]:=Module[{u=range[[2]],top=range[[3]],list={range[[2]]},v,newlist,actualmesh=0},While[u+mesh<top,v=u+Random[Real,mesh];
actualmesh=Max[actualmesh,v-u];
newlist:=Append[list,v];
u=v;
list=newlist];
partitionstring=StringJoin["Random partition with maximum norm ",ToString[mesh]];
actualmesh=Max[actualmesh,top-Last[list]];
Return[{Append[list,top],actualmesh}]]

UniformPartition[range_,number_]:=Module[{bottom=range[[2]]//N,top=range[[3]]//N,actualmesh},actualmesh=(top-bottom)/number;
partitionstring=StringJoin["Uniform partition into ",ToString[number]," subintervals"];
Return[{N[Range[bottom,top,actualmesh]],actualmesh}]]

Slices[expression_,variable_,list_,options___Rule]:=Module[{i=1,startlist=list,leftlist,rightlist,widthlist,selectlist,arealist,heightlist,valuelist,opt=Height/.{options}/.Options[RiemannData],f=Function[variable,expression]},leftlist=Drop[list,-1];
rightlist=Drop[list,1];
widthlist=rightlist-leftlist;
selectlist=Switch[opt,LeftSide,(heightstring="Left edge of subinterval used for height";leftlist),RightSide,(heightstring="Right edge of subinterval used for height";rightlist),Midpoint,(heightstring="Midpoint of subinterval used for height";
leftlist+0.5*widthlist),Random,(heightstring="Random point in subinterval used for height";
leftlist+(Random[Real,#]&/@widthlist))];
heightlist=f/@selectlist;
arealist=heightlist*widthlist;
{heightlist,arealist}]

(*Boxx::usage="Boxx[bottom, top, left, right]
produces a series of Line statements which in a
Graphics object will draw a hollow rectangle with
the corners given.  (Note that the built in
command Rectangle draws a solid rectangle.)"*)

Boxx[bottom_,top_,left_,right_]:=Line[{{left,bottom},{right,bottom},{right,top},{left,top},{left,bottom}}]

ShowSlices[partition_,heightlist_]:=Table[Boxx[0,heightlist[[i]],partition[[i]],partition[[i+1]]],{i,1,Length[partition]-1}]

Options[RiemannData]:={Subintervals->12,Partition->Uniform,Mesh->intervallength/6.0,Height->Midpoint}

RiemannData[expression_, range_, options___Rule] :=
	Block[{M=Mesh /. {options} 
			/. Options[RiemannData],
		S=Subintervals /. {options} 
			/. Options[RiemannData],
		P=Partition /. {options} 
			/. Options[RiemannData],
		variable=range[[1]],
		intervallength=range[[3]]-range[[2]]},
		list=Switch[P, Uniform, 
			UniformPartition[range, S],
			Random,
			RandomPartition[range, M]
           ]; 
	{list[[1]]} ~Join~ 
		Slices[expression, variable, list[[1]], options]
		~Join~ {list[[2]]}]

RiemannSum[expression_,range_,options___Rule]:=Apply[Plus,RiemannData[expression,range,options][[3]]]

PlotRiemann[expression_,range_,{plotoptions___},{riemannoptions___}]:=
Block[
{g,u,partition,heightlist,sum,mesh,integralvalue,outstring},(u=RiemannData[expression,range,riemannoptions
];
partition=u[[1]];
heightlist=u[[2]];
sum=Apply[Plus,u[[3]]];
mesh=u[[4]];
integralvalue=NIntegrate[expression,range];
outstring=StringJoin[
"\n",
partitionstring,"\n",
heightstring,"\n",
"Riemann Sum = ",ToString[Chop[sum]],"    ",
"Definite Integral = ",ToString[integralvalue],"\n",
"Norm = ",ToString[mesh],"    ","Percent error =",ToString[
NumberForm[
Chop[Abs[100*(integralvalue-sum)/integralvalue]],
3]
]
];
g:=Plot[expression,range,DisplayFunction->Identity,PlotLabel->outstring,PlotStyle->{{RGBColor[0,0,1]}},plotoptions];
h:=Graphics[{RGBColor[.8,.3,0],ShowSlices[partition,heightlist]}];
Show[g,h,DisplayFunction->$DisplayFunction])]

PlotRiemann[expression_,range_]:=PlotRiemann[expression,range,{},{}];

The Display

SetOptions[Plot,AspectRatio->1];

Manipulate[
PlotRiemann[
fn,
{x,0,1},
{ImageSize->{350,350}},
{Partition -> pt,Height -> ht,Mesh->N[1/n],Subintervals ->n}
],
{{n,4},1,100,1},
{{pt,Uniform, "Partition Type"},{Uniform, Random}},
{{ht,Midpoint, "Height"},{Midpoint,Random, LeftSide, RightSide}},
{{fn,Sqrt[1-x^2],"Function"},{ x, x^2, x^5,Sqrt[1-x^2],Sin[Pi x]}},
SaveDefinitions->True]