DEFINING  THE  INTEGRAND  f[x],  THE  LIMITS  OF  INTEGRATION  a  AND  b,  
AND  COMPUTING  THE  EXACT  VALUE  OF  THE  INTEGRAL

f[x_]=Sqrt[x]

a = 4

b = 9

exactValue=Integrate[f[x],{x,a,b}]

LEFT  RIEMANN  SUM

numPoints=10;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  leftRiemannSum=Sum[f[x[i]]*DeltaX,{i,0,numPoints-1}];  N[leftRiemannSum - exactValue,8]

numPoints=100;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  leftRiemannSum=Sum[f[x[i]],{i,0,numPoints-1}]*DeltaX;  N[leftRiemannSum - exactValue,8]

numPoints=1000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ ; x[i]=a+i*DeltaX];  leftRiemannSum=Sum[f[x[i]],{i,0,numPoints-1}]*DeltaX;  N[leftRiemannSum - exactValue,8]

numPoints=10000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ ; x[i]=a+i*DeltaX];  leftRiemannSum=Sum[f[x[i]],{i,0,numPoints-1}]*DeltaX;  N[leftRiemannSum - exactValue,8]

numPoints=100000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0]; 
For[ i = 0 ,i≤ numPoints, i++ ; x[i]=a+i*DeltaX];  leftRiemannSum=Sum[f[x[i]],{i,0,numPoints-1}]*DeltaX;  N[leftRiemannSum - exactValue,8]

RIGHT  RIEMANN  SUM

numPoints=10;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  leftRiemannSum=Sum[f[x[i+1]],{i,0,numPoints-1}]*DeltaX;  N[leftRiemannSum - exactValue,8]

numPoints=100;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  leftRiemannSum=Sum[f[x[i+1]],{i,0,numPoints-1}]*DeltaX;  N[leftRiemannSum - exactValue,8]

numPoints=1000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  leftRiemannSum=Sum[f[x[i+1]],{i,0,numPoints-1}]*DeltaX;  N[leftRiemannSum - exactValue,8]

numPoints=10000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  leftRiemannSum=Sum[f[x[i+1]],{i,0,numPoints-1}]*DeltaX;  N[leftRiemannSum - exactValue,8]

numPoints=100000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  leftRiemannSum=Sum[f[x[i+1]],{i,0,numPoints-1}]*DeltaX;  N[leftRiemannSum - exactValue,8]

MIDPOINT  RULE

numPoints=10;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  midpointRule=Sum[f[(x[i]+x[i+1])/2]*DeltaX,{i,0,numPoints-1}];  N[midpointRule - exactValue,8]

numPoints=100;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  midpointRule=Sum[f[(x[i]+x[i+1])/2]*DeltaX,{i,0,numPoints-1}];  N[midpointRule - exactValue,8]

numPoints=1000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  midpointRule=Sum[f[(x[i]+x[i+1])/2]*DeltaX,{i,0,numPoints-1}];  N[midpointRule - exactValue,8]

numPoints=10000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  midpointRule=Sum[f[(x[i]+x[i+1])/2]*DeltaX,{i,0,numPoints-1}];  N[midpointRule - exactValue,8]

numPoints=100000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  midpointRule=Sum[f[(x[i]+x[i+1])/2]*DeltaX,{i,0,numPoints-1}];  N[midpointRule - exactValue,8]

TRAPEZOIDAL  RULE

numPoints=10;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  midpointRule=(f[x[0]]+2*Sum[f[x[i]],{i,1,numPoints-1}]+f[x[numPoints]])*DeltaX/2;  N[midpointRule - exactValue,8]

numPoints=100;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  midpointRule=(f[x[0]]+2*Sum[f[x[i]],{i,1,numPoints-1}]+f[x[numPoints]])*DeltaX/2;  N[midpointRule - exactValue,8]

numPoints=1000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  midpointRule=(f[x[0]]+2*Sum[f[x[i]],{i,1,numPoints-1}]+f[x[numPoints]])*DeltaX/2;  N[midpointRule - exactValue,8]

numPoints=10000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  midpointRule=(f[x[0]]+2*Sum[f[x[i]],{i,1,numPoints-1}]+f[x[numPoints]])*DeltaX/2;  N[midpointRule - exactValue,8]

numPoints=100000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  midpointRule=(f[x[0]]+2*Sum[f[x[i]],{i,1,numPoints-1}]+f[x[numPoints]])*DeltaX/2;  N[midpointRule - exactValue,8]

SIMPSON'S  RULE

numPoints=10;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  simpsonRule=(f[x[0]]+2*Sum[f[x[2*i]],{i,1,numPoints/2-1}]+4*Sum[f[x[2*i+1]],{i,0,numPoints/2-1}]+f[x[numPoints]])*DeltaX/3;  N[simpsonRule - exactValue,8]

numPoints=100;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  simpsonRule=(f[x[0]]+2*Sum[f[x[2*i]],{i,1,numPoints/2-1}]+4*Sum[f[x[2*i+1]],{i,0,numPoints/2-1}]+f[x[numPoints]])*DeltaX/3;  N[simpsonRule - exactValue,8]

numPoints=1000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  simpsonRule=(f[x[0]]+2*Sum[f[x[2*i]],{i,1,numPoints/2-1}]+4*Sum[f[x[2*i+1]],{i,0,numPoints/2-1}]+f[x[numPoints]])*DeltaX/3;  N[simpsonRule - exactValue,8]

numPoints=10000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  simpsonRule=(f[x[0]]+2*Sum[f[x[2*i]],{i,1,numPoints/2-1}]+4*Sum[f[x[2*i+1]],{i,0,numPoints/2-1}]+f[x[numPoints]])*DeltaX/3;  N[simpsonRule - exactValue,8]

numPoints=100000;  DeltaX=(b-a)/numPoints;  
Array[x,numPoints + 1,0];
For[ i = 0 ,i≤ numPoints, i++ , x[i]=a+i*DeltaX];  simpsonRule=(f[x[0]]+2*Sum[f[x[2*i]],{i,1,numPoints/2-1}]+4*Sum[f[x[2*i+1]],{i,0,numPoints/2-1}]+f[x[numPoints]])*DeltaX/3;  N[simpsonRule - exactValue,8]