Numerical Integration Techniques

Helmut Knaust, Department of Mathematical Sciences, UTEP, El Paso TX 79968
hknaust@utep.edu
Created:3/16/1996; last edits:3/22/2012

Upper and Lower Darboux Sums

DarbouxSum[f[x], {x, a, b}, n, s] visualizes and computes the lower and upper Darboux Sums for the 
function f(x) on the interval [a, b], using n equal - length subintervals.
The maximum and the minimum on each subinterval is computed by sampling f(x) at (s + 1) equally - spaced points in each 
subinterval. (The values for n and s are optional; the defaults are 20, and 25, respectively.)

DarbouxSum[f_,{x_,a_,b_},n_Integer:20,split_Integer:10,options___Rule]:=
Module[{h,g,fct,min,max},h=(b-a)/n;g[s_]=f/.{x->s};Compile{x,g[x]};
fct=Plot[g[x],{x,a,b},PlotStyle->{Black,Thick}];
min[j_]:=Min[Table[g[j+k*h/split],{k,0,split}]];
max[j_]:=Max[Table[g[j+k*h/split],{k,0,split}]];rect=Flatten[Table[If[N[max[j]]<0,
{{Blue,EdgeForm[Darker[Blue]],Opacity[0.6],Rectangle[{j,0},{j+h,min[j]}]},
{Red,EdgeForm[Darker[Red]],Opacity[0.6],Rectangle[{j,0},{j+h,max[j]}]}},
{{Red,EdgeForm[Darker[Red]],Opacity[0.6],Rectangle[{j,0},{j+h,max[j]}]},
{Blue,EdgeForm[Darker[Blue]],Opacity[0.6],Rectangle[{j,0},{j+h,min[j]}]}}],
{j,a,b-h,h}],1];Print[Show[{fct,Graphics[rect],fct},Axes->False,Frame->True,PlotRange->{{a,b},All},ImageSize->600,options]];Print["The Lower Darboux Sum (in blue) with ",n," intervals is  ",
N[Sum[min[j],{j,a,b-h,h}]*h]];Print["The Upper Darboux Sum (in red) with ",n," intervals is  ",
N[Sum[max[j],{j,a,b-h,h}]*h]]]

DarbouxSum[Sin[t^2],{t,0,3},30]

Left and Right Riemann Sums

LeftRSum[f[x], {x, a, b}, n] visualizes and computes the left Riemann Sum for the function f(x) on the interval [a, b], 
using n equal - length subintervals. (The value for n is optional; the default is 20.)

LeftRSum[f_,{x_,a_,b_},n_Integer:20,options___Rule]:=
Module[{h,g,fct},h=(b-a)/n;g[s_]=f/.{x->s};Compile{x,g[x]};
fct=Plot[g[x],{x,a,b},PlotStyle->{Black,Thick}];
Print[Show[{fct,Graphics[{Blue,Opacity[0.6],EdgeForm[Darker[Blue]],
Table[Rectangle[{j,0},{j+h,g[j]}],{j,a,b-h,h}]}],fct},
Axes->False,Frame->True,PlotRange->{{a,b},All},ImageSize->600,options]];
Print["The Left Riemann Sum with ",n," intervals is  ",
N[Sum[g[j],{j,a,b-h,h}]*h]]]

LeftRSum[Sin[t^2],{t,0,3},50]

RightRSum[f[x], {x, a, b}, n] visualizes and computes the right Riemann Sum for the function f(x) on the interval [a, b], 
using n equal - length subintervals.  (The value for n is optional; the default is 20.)

RightRSum[f_,{x_,a_,b_},n_Integer:20,options___Rule]:=
Module[{h,g,fct},h=(b-a)/n;g[s_]=f/.{x->s};Compile{x,g[x]};
fct=Plot[g[x],{x,a,b},PlotStyle->{Black,Thick}];
Print[Show[{fct,Graphics[
{Red,Opacity[0.6],EdgeForm[Darker[Red]],
Table[Rectangle[{j,0},{j+h,g[j+h]}],{j,a,b-h,h}]}],fct},
Axes->False,Frame->True,PlotRange->{{a,b},All},ImageSize->600,options]];
Print["The Right Riemann Sum with ",n," intervals is  ",
N[Sum[g[j],{j,a+h,b,h}]*h]]]

RightRSum[Sin[t^2],{t,0,3},50]

The Midpoint Rule

MidPtSum[f[x], {x, a, b}, n] visualizes and computes the Midpoint Riemann Sum for the function f(x) on the 
interval [a, b], using n equal - length subintervals.  (The value for n is optional; the default is 20.)

MidPtSum[f_,{x_,a_,b_},n_Integer:20,options___Rule]:=
Module[{h,g,fct},h=(b-a)/n;g[s_]=f/.{x->s};Compile{x,g[x]};
fct=Plot[g[x],{x,a,b},PlotStyle->{Black,Thick}];
Print[Show[{fct,Graphics[{Orange,EdgeForm[Darker[Orange]],Opacity[0.6],
Table[Rectangle[{j,0},{j+h,g[j+h/2]}],{j,a,b-h,h}]}],fct},
Axes->False,Frame->True,PlotRange->{{a,b},All},ImageSize->600,options]];
Print["The Midpoint Riemann Sum with ",n, " intervals is  ",N[Sum[g[j+h/2],{j,a,b-h,h}]*h]]]

MidPtSum[Sin[τ^2],{τ,0,3},25]

The Trapezoidal Rule

TrapRule[f[x],{x,a,b},n] visualizes and computes an approximation for the integral of the function f(x) on the interval [a,b],
 using the Trapezoidal Rule with n equal-length subintervals. (The value for n is optional; the default is 20.)

TrapRule[f_,{x_,a_,b_},n_Integer:20,options___Rule]:=
Module[{h,g,fct},h=(b-a)/n;g[s_]=f/.{x->s};Compile{x,g[x]};
fct=Plot[g[x],{x,a,b},PlotStyle->{Black,Thick}];
Print[Show[{fct,Graphics[{Brown,Opacity[0.6],EdgeForm[Darker[Brown]],Table[Polygon[{{j,0},{j+h,0},{j+h,g[j+h]},{j,g[j]}}],{j,a,b-h,h}]}],fct},
Axes->False,Frame->True,PlotRange->{{a,b},All},ImageSize->600,options]];
Print["The Trapezoidal Rule with ",n," intervals yields the result   ",(2*N[Sum[g[j],{j,a+h,b-h,h}]]+N[g[a]]+N[g[b]])*h/2
]]

TrapRule[Sin[t^2],{t,0,3},17]

Simpson' s Rule

SimpRule[f[x],{x,a,b},n] visualizes and computes an approximation for the integral of the function f[x] on the interval [a,b],
 using Simpson's Rule with n equal-length subintervals. (The (even) value for n is optional; the default is 10.)

SimpRule[f_,{x_,a_,b_},n_Integer:10,options___Rule]:=
Module[ {h,g,fct,m},If[EvenQ[n],m=n,(m=n+1;Print["The number of intervals must be even."];
Print["Number of intervals increased to ",m,"."])];h=(b-a)/m;g[s_]=f/.{x->s};Compile{x,g[x]};fct=Plot[g[x],{x,a,b},PlotStyle->{Black,Thick},options];Print[Show[{fct,Table[Plot[
Evaluate[Fit[{{j,g[j]},{j+h,g[j+h]},{j+2h,g[j+2h]}},{1,x,x^2},x]],
{x,j,j+2h},Filling->0,FillingStyle->Green],{j,a,b-2h,2h}],Graphics[{Black,Table[Line[{{j,0},{j,g[j]}}],{j,a,b-2h,2h}]}],
Graphics[{Black,Dashed,Table[Line[{{j,0},{j,g[j]}}],{j,a+h,b-h,2h}]}],fct},
Axes->False,Frame->True,PlotRange->{{a,b},All},ImageSize->600]];
Print["Simpson's Rule with ", m ," intervals yields the result ",(N[g[a]]+4N[g[b-h]]+N[g[b]]+N[Sum[4*g[j]+2*g[j+h],
{j,a+h,b-2h,2h}]])*h/3]]

SimpRule[Sin[t^2],{t,0,3},8]