Error Estimates for Numerical Integration

Helmut Knaust, Department of Mathematical Sciences, UTEP, El Paso TX 79968
hknaust@utep.edu
Created:8/21/2000; last edits:4/11/2012

The command will output the integral approximations and the errors in the approximation for the various methods first in the case of n subdivisions and then in the case of 10n subdivisions. Afterwards it will output the factor by which the error decreases when the number of subdivisions increases from n to 10n.

From the theory we expect this factor to be about 10 for the Left and Right Riemann Sums, about (10^2)=100 for the Midpoint and the Trapezoidal Rule, and about (10^4)=10,000 for Simpson's Rule.
If you compute these factors for a few examples, you will realize that the “ideal” numbers above are very rough estimates for the factor. When the number of subdivisions is "too small", the points are sampled too much at random to give results close to the expaected ones; on the other hand, when the number of subdivisions is "too large", other considerations such as the round-off errors committed may distort the results.

The following command will only work for functions, for which Mathematica can symbolically compute the definite integral!

Error[func_,{var_,a_,b_},k_Integer:10]:=Module[{f,n,int,listend,listmid,leftriemann,rightriemann,err1,err2,midpoint,err3,trapez,err4,simpson,err5,err10,err20,err30,err40,err50},f[x_]=func/.{var->x};n=k;int=Integrate[f[x],{x,a,b}];listend=Table[f[x],{x,a,b,(b-a)/n}];listmid=Table[f[x],{x,a+(b-a)/(2n),b,(b-a)/n}]
Print[TableForm[Table[{{StringForm["[`` Subdivisions]",n],"Integral Approximation","Error"},{"","",""},{"Left Riemann Sum:",leftriemann=N[Apply[Plus,Drop[listend,-1]] (b-a)/n],err1=Abs[leftriemann-int]},
			{"Right Riemann Sum:",rightriemann=N[Apply[Plus,Drop[listend,1]] (b-a)/n],err2=Abs[rightriemann-int]},{"Midpoint Rule:",midpoint=N[Apply[Plus,listmid](b-a)/n],err3=Abs[midpoint-int]},
			{"Trapezoidal Rule:",trapez=1/2 (rightriemann+leftriemann),err4=Abs[trapez-int]},{"Simpson's Rule:",simpson=2/3 midpoint +1/3 trapez,err5=Abs[simpson-int]}}]]];
n=10k;listend=Table[f[x],{x,a,b,(b-a)/n}];listmid=Table[f[x],{x,a+(b-a)/(2n),b,(b-a)/n}];
Print[TableForm[Table[{{StringForm["[`` Subdivisions]",n],"Integral Approximation","Error"},{"","",""},"Left Riemann Sum:",leftriemann=N[Apply[Plus,Drop[listend,-1]] (b-a)/n],err10=Abs[leftriemann-int]},{"Right Riemann Sum:",rightriemann=N[Apply[Plus,Drop[listend,1]] (b-a)/n],err20=Abs[rightriemann-int]},
			{"Midpoint Rule:",midpoint=N[Apply[Plus,listmid](b-a)/n],err30=Abs[midpoint-int]},{"Trapezoidal Rule:",trapez=1/2 (rightriemann+leftriemann),err40=Abs[trapez-int]},
			{"Simpson's Rule:",simpson=2/3 midpoint +1/3 trapez,err50=Abs[simpson-int]}}]]];
		Print[TableForm[Table[{{"","Error decreases approximately by a factor:" },{"",""},{"Left Riemann Sum:",Round[err1/err10]},
{"Right Riemann Sum:",Round[err2/err20]},{"",""},{"Midpoint Rule:",Round[err3/err30]},{"Trapezoidal Rule:",Round[err4/err40]},{"",""},{"Simpson's Rule:",Round[err5/err50]}}]]];]

Error[Sin[5x],{x,0,2},50]