Calculus II: Approximate Integration

Approximate Integration Formulas

Clear[f,a,b,n]

Here are some rules for integration. Below each general rule, the rule is evaluated for a=0, b=1 and five intervals. The function f is left unspecified. This shows how the formulas we looked at in class relate to the formulas in the text.

Left Hand Rule

L[n_,a_,b_] := (b-a)/n  Sum[f[a+(b-a)/n*(i-1)],{i,1,n}]
L[n,a,b]
L[6,0,1]

Right Hand Rule

R[n_,a_,b_] := (b-a)/n  Sum[f[a+(b-a)/n*i],{i,1,n}]
R[n,a,b]
R[6,0,1]

MidPoint Rule

M[n_,a_,b_] := (b-a)/n  Sum[f[a+(b-a)/(2n)*(2i-1)],{i,1,n}]
M[n,a,b]
M[6,0,1]

Trapezoidal Rule

T[n_,a_,b_] := (b-a)/(2n) ( f[a]+f[b]+2Sum[f[a+(b-a)/n*i],{i,1,n-1}])
T[n,a,b]
T[6,0,1]

Simpson's Rule

S[n_,a_,b_]:=T[n/2,a,b]/3+2M[n/2,a,b]/3
S[n,a,b]
S[6,0,1]

When f(x) is Given Explicitly: 
You Use a Large Value of n, and can Vary n to See How Good Your Approximation Is

First, define the formulas you want to use (note the use of := (set delayed)):

M[n_,a_,b_] := SetPrecision[(b-a)/n  Sum[f[a+(b-a)/(2n)*(2i-1)],{i,1,n}],10]

We have used the SetPrecision to ensure our answer will always be a decimal (10 decimals in this case).

Now, you can use it to create tables using the following bit of code:

M[n_,a_,b_] := SetPrecision[(b-a)/n  Sum[f[a+(b-a)/(2n)*(2i-1)],{i,1,n}],10]

f[x_]=x^2
a=0
b=1
Table[{n,M[n,a,b]},{n,10,40,2}];
TableForm[%]

If you have a different function, or limits, you just have to change them:

M[n_,a_,b_] := SetPrecision[(b-a)/n  Sum[f[a+(b-a)/(2n)*(2i-1)],{i,1,n}],10]

f[x_]=Cos[x]^2*Exp[-x]
a=0
b=2Pi
Table[{n,M[n,a,b]},{n,10,200,10}];
TableForm[%]

An Application

We can use these rules to numerically approximate integrals. For example, here is an integral Mathematica cannot do:

f[x_]=Exp[x^4+x^3]
a=0
b=1
Integrate[f[x],{x,a,b}]

Here, we use the midpoint rule to approximate the integral:

M[n_,a_,b_] := SetPrecision[(b-a)/n  Sum[f[a+(b-a)/(2n)*(2i-1)],{i,1,n}],10]

f[x_]=Exp[x^4+x^3]
a=0
b=1
M[1000,a,b]

We can construct a table of values and help us see how the approximation changes with increasing n:

M[n_,a_,b_] := SetPrecision[(b-a)/n  Sum[f[a+(b-a)/(2n)*(2i-1)],{i,1,n}],10]

f[x_]=Exp[x^4+x^3]
a=0
b=1
Table[{n,M[n,a,b]},{n,10,40}];
TableForm[%]

T[n_,a_,b_] := SetPrecision[(b-a)/(2n) ( f[a]+f[b]+2Sum[f[a+(b-a)/n*i],{i,1,n-1}]),10]

f[x_]=Exp[x^4+x^3]
a=0
b=1
Table[{n,T[n,a,b]},{n,10,40}];
TableForm[%]

M[n_,a_,b_] := SetPrecision[(b-a)/n  Sum[f[a+(b-a)/(2n)*(2i-1)],{i,1,n}],10]
T[n_,a_,b_] := SetPrecision[(b-a)/(2n) ( f[a]+f[b]+2Sum[f[a+(b-a)/n*i],{i,1,n-1}]),10]
S[n_,a_,b_]:=T[n/2,a,b]/3+2M[n/2,a,b]/3

f[x_]=Exp[x^4+x^3]
a=0
b=1
Table[{n,M[n,a,b],T[n,a,b],S[n,a,b]},{n,40,70,2}];
TableForm[%,TableHeadings->{None,{"n","Midpoint","Trapezoidal","Simpson's"}}]

Mathematica has a built-in command to perform numerical integration:

NIntegrate[f[x],{x,0,1}]

The Error

In the above example, we saw that the integral was about 1.88303 using Simpson's rule. Let's assume that is accurate, and sketch the errors associated with the Midpoint Rule, Trapeziod Rule, and Simpson's Rule.

We can construct a table of values and help us see how the approximation changes with increasing n:

midpointerror=Table[{n,Abs[M[n,0,1.]-1.88303]},{n,10,40}];
TableForm[midpointerror]
meplot=ListPlot[midpointerror,PlotStyle->{Green,PointSize[0.02]}]

traperror=Table[{n,Abs[T[n,0,1.]-1.88303]},{n,10,40}];
TableForm[traperror]
teplot=ListPlot[traperror,PlotStyle->{Red,PointSize[0.02]}]

simperror=Table[{n,Abs[S[n,0,1.]-1.88303]},{n,10,40,2}];
TableForm[simperror]
seplot=ListPlot[simperror,PlotStyle->{Blue,PointSize[0.02]}]

Show[meplot,seplot,teplot,AxesLabel->{"n","Error"}]

Simpson's rule has a smaller error than the other two.

The text has exact error bounds for the different rules. Where there error bounds come from is something that you would learn in a Numerical Methods course, where you would also learn a way of calulating these integration formulas that easily allows you to calculate more impressive formulas (impressive since they are more accurate). These techniques are called quadrature methods.

When f(x) is Given in a Data Table:
You Are Stuck With Whatever n You Get From The Table

Here, there is no explicit y=f(x) to refer to, so we cannot use the formulas that constructed to work with different values of n. You are stuck with whatever value of n you have in the table of values, and you get an approximation and that's that (it might be a good approximation, or a crappy approximation, but since you cannot vary n its the best you can do). Here, it is most useful to use the formulas from the textbook (you typically do not use (Subscript[L, n]), (Subscript[R, n]), or (Subscript[M, n]) in this situation, since (Subscript[T, n]) and (Subscript[S, n]) will give better results):

(Subscript[T, n])=((Δx )/2) ((Subscript[y, 0])+2 (Subscript[y, 1])+2(Subscript[y, 2])+2(Subscript[y, 3])+2(Subscript[y, 4])+2(Subscript[y, 5])+...+2(Subscript[y, n-2])+2(Subscript[y, n-1])+(Subscript[y, n]))

(Subscript[S, n])=((Δx )/3) ((Subscript[y, 0])+4 (Subscript[y, 1])+2(Subscript[y, 2])+4(Subscript[y, 3])+2(Subscript[y, 4])+4(Subscript[y, 5])+...+2(Subscript[y, n-2])+4(Subscript[y, n-1])+(Subscript[y, n]))

Note that n must be even for Simpson’s rule.

I switched to f((Subscript[x, i]))=(Subscript[y, i]) since we usually think of the data as given in ((Subscript[x, i]),(Subscript[y, i])) pairs. Notice that we have n+1 points since the counting begins at zero (that's quite common in numerical analysis).

Here is a sample data set, and let's say we want to figure out (∫_1^8)f(x)dx for it, where f(x) represents the underlying function that gave us the data set, but we don't know what f(x) is explicitly.

data={{1,3},{2,4},{3,7},{4,12},{5,8},{6,2},{7,1},{8,0.5},{9,0}}

Here's a sketch:

ListPlot[data,PlotStyle->PointSize[0.02],AxesOrigin->{0,0}]

Since we have 9 points, n=8. It's really easiest just to use Mathematica as a calculator:

(Subscript[T, 8])=((Δx )/2) ((Subscript[y, 0])+2 (Subscript[y, 1])+2(Subscript[y, 2])+2(Subscript[y, 3])+2(Subscript[y, 4])+2(Subscript[y, 5])+2(Subscript[y, 6])+2(Subscript[y, 7])+(Subscript[y, 8]))
=(1/2) (3+2 (4)+2(7)+2(12)+2(8)+2(2)+2(1)+2(0.5)+0)

(1/2) (3+2 (4)+2(7)+2(12)+2(8)+2(2)+2(1)+2(0.5)+0)

We can do something similar for Simpson's Rule, and get another approximation to (∫_1^8)f(x)dx (notice we have n even, which is required for Simpson's Rule)

(Subscript[S, 8])=((Δx )/3) ((Subscript[y, 0])+4 (Subscript[y, 1])+2(Subscript[y, 2])+4(Subscript[y, 3])+2(Subscript[y, 4])+4(Subscript[y, 5])+2(Subscript[y, 6])+4(Subscript[y, 7])+(Subscript[y, 8]))
=(1/3) (3+4 (4)+2(7)+4(12)+2(8)+4(2)+2(1)+4(0.5)+0)

(1/3) (3+4 (4)+2(7)+4(12)+2(8)+4(2)+2(1)+4(0.5)+0)