Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Sept 2016
Lecture topic: Composite Quadrature

Deriving Some Formulas for Integration

Clear[a,n,k,x,xp,i,Lag,f]

Here is the wiggly polynomial defined. Notice I am using Mathematica's integral abilities to calculate the integral of the wiggly polynomial, since it is just a polynomial and conceptual not hard to integrate (keeping track of the details might be a pain, of course!).

Lag[n_,k_,x_]:=((∏_(i=0)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))
IntLag[n_,k_,x_,a_,b_]:=Integrate[Lag[n,k,x],{x,a,b}]

Here is the quadrature formula that approximates the integral and the error:

intf[n_,a_,b_]:=(∑_(k=0)^n)f[xp[k]]*IntLag[n,k,x,a,b]
intferror[n_,a_,b_]:=(derf/((n+1)!))Integrate[(∏_(i=0)^n)(x-xp[i]),{x,a,b}]

Note that I am assuming we can apply the mean value theorem for integrals to the error term, so I have removed the (f^((n)))(c(x)) from the integral for the error. This is just to make the error computation easier for Mathematica whereas a hand calculuation would be more correct. The expressions we obtain for the integration formulas are correct (and so are the error estimates, as it turns out, although not rigorously defined).

n=0 (Left Hand Rule)

Define our equally spaced nodes:

xp[0]=a
xp[1]=a+h

Here is the simplification when n=0:

intf[0,xp[0],xp[1]]
intferror[0,xp[0],xp[1]]

n=1 (Trapezoidal Rule Eqn 5.21)

Define our equally spaced nodes:

xp[0]=a
xp[1]=a+h

Here is the simplification when n=1:

intf[1,xp[0],xp[1]]
intferror[1,xp[0],xp[1]]

n=0 (MidPoint Rule Eqn. 5.26)

Define our equally spaced nodes:

xp[0]=w=(a+b)/2

Here is the simplification when n=0 (note we needed to go to higher derivatives again to get the error):

intf[0,a,b]
intferror[0,a,b]
intferror[1,a,b]

n=2 (Simpson's Rule Eqn. 5.22)

Define our equally spaced nodes:

xp[0]=a
xp[1]=a+h
xp[2]=a+2h

Here is the simplification when n=2:

intf[2,xp[0],xp[2]]
intferror[2,xp[0],xp[2]]

What happened to our error estimate? We need to go to higher derivatives to get the error estimate.

intferror[3,xp[0],xp[2]]

n=3 (Simpson's Three-Eigths Rule Example 5.7)

Define our equally spaced nodes:

xp[0]=a
xp[1]=a+h
xp[2]=a+2h
xp[3]=a+3h

Here is the simplification when n=3 (Equation 4.25):

intf[3,xp[0],xp[3]]
intferror[3,xp[0],xp[3]]

n=4 (Boole’s Rule)

Define our equally spaced nodes:

xp[0]=a
xp[1]=a+h
xp[2]=a+2h
xp[3]=a+3h
xp[4]=a+4h

Here is the simplification when n=4 (Problem 5.2.15):

intf[4,xp[0],xp[4]]
intferror[4,xp[0],xp[4]]
intferror[5,xp[0],xp[4]]

Approximate Integration Formulas: New Modules Using Composite Idea

Clear[f,a,b,n,h]

Here we define new modules that act as integrators, using the composite idea to reuse a formula that works over one subinterval.

n=2 (Simpson's Rule)

Here is the simplification when n=2:

xp[0]=a;
xp[1]=a+h;
xp[2]=a+2h;
intf[2,xp[0],xp[2]]

Let's use this to construct a new integration module that uses this formula!

Here is the formula for the single interval between a and b:

simp[f_,a_,b_]:=Module[{h,f0,f1,f2,f3},
h=((b-a)/2);
f0=f/.x->a;
f1=f/.x->a+h;
f2=f/.x->a+2h;
Return[SetPrecision[(1/3) h f0+(4/3) h f1+(1/3) h f2,12]];
]

Here is the composite formula that uses the above and works with n subintervals.

simpcomposite[f_,a_,b_,n_]:=Module[{h,sum,x1,x2},
h=N[((b-a)/n)];
sum=0.0;
x1=a-h;
x2=a;
Do[
x1=x1+h;
x2=x2+h;
sum=sum+simp[f,x1,x2]
,{i,1,n}];
Return[SetPrecision[sum,12]];
]

TableForm[Table[{n,simpcomposite[(x^7),0,1,n]},{n,5,100,5}],TableHeadings->{None,{"Subintervals","Simpson's Approximation"}}]

n=3 (Simpson's Three-Eigths Rule)

Here is the simplification when n=3:

xp[0]=a;
xp[1]=a+h;
xp[2]=a+2h;
xp[3]=a+3h;
intf[3,xp[0],xp[3]]

Let's use this to construct a new integration module that uses this formula!

Here is the formula for the single interval between a and b:

simp38[f_,a_,b_]:=Module[{h,f0,f1,f2,f3},
h=((b-a)/3);
f0=f/.x->a;
f1=f/.x->a+h;
f2=f/.x->a+2h;
f3=f/.x->a+3h;
Return[SetPrecision[(3/8) h f0+(9/8) h f1+(9/8) h f2+(3/8) h f3,12]];
]

Here is the composite formula that uses the above and works with n subintervals. Notice that the only change is what function I call!

simp38composite[f_,a_,b_,n_]:=Module[{h,sum,x1,x2},
h=N[((b-a)/n)];
sum=0.0;
x1=a-h;
x2=a;
Do[
x1=x1+h;
x2=x2+h;
sum=sum+simp38[f,x1,x2]
,{i,1,n}];
Return[SetPrecision[sum,12]];
]

TableForm[Table[{n,simpcomposite[(x^7),0,1,n],simp38composite[(x^7),0,1,n]},{n,5,100,5}],TableHeadings->{None,{"Subintervals","Simpson's Approximation","Simpson's 3/8 Approximation"}}]

n=4 (Boole’s Rule)

Here is the simplification when n=4:

xp[0]=a;
xp[1]=a+h;
xp[2]=a+2h;
xp[3]=a+3h;
xp[4]=a+4h;
intf[4,xp[0],xp[4]]

Let's use this to construct a new integration module that uses this formula!

Here is the formula for the single interval between a and b:

n4[f_,a_,b_]:=Module[{h,f0,f1,f2,f3,f4},
h=((b-a)/4);
f0=f/.x->a;
f1=f/.x->a+h;
f2=f/.x->a+2h;
f3=f/.x->a+3h;
f4=f/.x->a+4h;
Return[SetPrecision[(14/45) h f0+(64/45) h f1+(8/15) h f2+(64/45) h f3+(14/45) h f4,12]];
]

Here is the composite formula that uses the above and works with n subintervals.

n4composite[f_,a_,b_,n_]:=Module[{h,sum,x1,x2},
h=N[((b-a)/n)];
sum=0.0;
x1=a-h;
x2=a;
Do[
x1=x1+h;
x2=x2+h;
sum=sum+n4[f,x1,x2]
,{i,1,n}];
Return[SetPrecision[sum,12]];
]

Comparison

TableForm[Table[{n,simpcomposite[(x^7),0,1,n],simp38composite[(x^7),0,1,n],n4composite[(x^7),0,1,n]},{n,5,100,5}],TableHeadings->{None,{"Subintervals","Simpson's","Simpson's 3/8","Boole's Rule"}}]

Composite Trapezoidal Module that Stops at Given Tolerance

I’ll include all the modules needed here. Each interaction simply halves h across the entire region (halving the step size doubles the number of functional evaluations) until the tolerance is met.

trap[f_,a_,b_]:=Module[{h,f0,f1,f2,f3,f4},
h=b-a;
f0=f/.x->a;
f1=f/.x->a+h;
Return[SetPrecision[(1/2) h f0+(1/2) h f1,12]];
]

trapcomposite[f_,a_,b_,n_]:=Module[{h,sum,x1,x2},
h=SetPrecision[((b-a)/n),16];
sum=0.0;
x1=a-h;
x2=a;
Do[
x1=x1+h;
x2=x2+h;
sum=sum+trap[f,x1,x2]
,{i,1,n}];
Return[SetPrecision[sum,16]];
]

traptolerance[f_,a_,b_,tol_]:=Module[{r1,r2,max=(2^20)},
r1=trapcomposite[f,a,b,(2^1)];
r2=trapcomposite[f,a,b,(2^2)];

i=3;
While[(2^i)≤max&&Abs[r1-r2]>tol,
r1=r2;
r2=trapcomposite[f,a,b,2^i];
i=i+1];
If[i≥max,Print["Did Not Converge to Required Tolerance"]];
Return[r2];
]

f[x_]:=Sqrt[1-Cos[x]^2]/Sqrt[1+x^3];
a=0.0;
b=2.0;
Table[{(10.0^(-k)),traptolerance[f[x],a,b,(10.0^(-k))]},{k,2,10}]//TableForm

traptolerance[f[x],a,b,0.0000001]

Composite Simpson’s Three-Eigths Module that Stops at Given Tolerance

I’ll include all the modules needed here. This destroys the Trapezoidal module! Note this is set to work as a function of x, if you want a function of t or some other variable change the x in the first simp38 module to whatever variable you are using.

simp38[f_,a_,b_]:=Module[{h,f0,f1,f2,f3},
h=((b-a)/3);
f0=f/.x->a;
f1=f/.x->a+h;
f2=f/.x->a+2h;
f3=f/.x->a+3h;
Return[SetPrecision[(3/8) h f0+(9/8) h f1+(9/8) h f2+(3/8) h f3,12]];
]

simp38composite[f_,a_,b_,n_]:=Module[{h,sum,x1,x2},
h=SetPrecision[((b-a)/n),16];
sum=0.0;
x1=a-h;
x2=a;
Do[
x1=x1+h;
x2=x2+h;
sum=sum+simp38[f,x1,x2]
,{i,1,n}];
Return[SetPrecision[sum,16]];
]

simp38tolerance[f_,a_,b_,tol_]:=Module[{r1,r2,max=(2^20)},
r1=simp38composite[f,a,b,(2^1)];
r2=simp38composite[f,a,b,(2^2)];

i=3;
While[(2^i)≤max&&Abs[r1-r2]>tol,
r1=r2;
r2=simp38composite[f,a,b,2^i];
i=i+1];
If[i≥max,Print["Did Not Converge to Required Tolerance"]];
Return[r2];
]

f[x_]:=Sqrt[1-Cos[x]^2]/Sqrt[1+x^3];
a=0.0;
b=2.0;
Table[{(10.0^(-k)),simp38tolerance[f[x],a,b,(10.0^(-k))]},{k,2,10}]//TableForm

simp38tolerance[f[x],a,b,(10^(-11))]
simp38tolerance[f[x],a,b,(10^(-14))]

An Integral

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)]
Plot[f[x],{x,0,1},PlotStyle->Blue,ImageSize->{400,Automatic}]
Integrate[f[x],{x,0,1}]

Here, we use the n=4 rule to approximate the integral:

n4composite[f[x],0,1,100]

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

TableForm[Table[{n,simpcomposite[f[x],0,1,n],simp38composite[f[x],0,1,n],n4composite[f[x],0,1,n]},{n,10,40,2}],TableHeadings->{None,{"Subintervals","Simpson's","Simpson's 3/8","n=4"}}]

The Error

Let's use the difference between two successive approximations as an estimate of the error in the method (useful, since we don't know the exact value of the intergral in this case).

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

error1=Table[{n,Abs[simpcomposite[f[x],0,1,n]-simpcomposite[f[x],0,1,n-1]]},{n,10,60,2}];
plot1=ListLogPlot[error1,PlotStyle->{Black,PointSize[0.02]},ImageSize->{500,Automatic}];
error2=Table[{n,Abs[simp38composite[f[x],0,1,n]-simp38composite[f[x],0,1,n-1]]},{n,10,60,2}];
plot2=ListLogPlot[error2,PlotStyle->{Red,PointSize[0.02]},ImageSize->{500,Automatic}];
error3=Table[{n,Abs[n4composite[f[x],0,1,n]-n4composite[f[x],0,1,n-1]]},{n,10,60,2}];
plot3=ListLogPlot[error3,PlotStyle->{Blue,PointSize[0.02]},ImageSize->{500,Automatic}];

errors=Show[plot1,plot2,plot3,FrameLabel->{"n","Log Error"},PlotRange->{{10,60},Automatic},Frame->True,Axes->False,BaseStyle->20]

After n=50 we start to get some round-off error creeping into the n=4 Boole’s Rule approximation. It's incredibly accurate though! Blows the doors off the other methods.

Appendix: Approximate Integration Formulas: Calculus Definitions

Clear[f,a,b,n]

These are the definitions as you probably remember them from calculus.

Right Hand Rule

R[n_,a_,b_] = ((b-a)/n) (∑_(i=1)^n)f[a+((b-a)/n)*i]

MidPoint Rule

M[n_,a_,b_] =  ((b-a)/n) (∑_(i=1)^n)f[a+((b-a)/(2n))*(2i-1)]

Left Hand Rule

L[n_,a_,b_] = ((b-a)/n) (∑_(i=1)^n)f[a+((b-a)/n)*(i-1)]

Trapezoidal Rule

T[n_,a_,b_] =  ((b-a)/(2n))(f[a]+f[b]+2 (∑_(i=1)^(n-1))f[a+((b-a)/n)*i])

Simpson's Rule

S[n_,a_,b_]=((T[n/2,a,b])/3)+(2/3)M[n/2,a,b]

Trapezoidal Rule as a Module

Trapezoid[fcn_,n_,a_,b_] :=Module[{i,appr}, 
appr=(b-a)/(2n)( (fcn/.x->a)+(fcn/.x->b)+2Sum[(fcn/.x->a+(b-a)/n*i),{i,1,n-1}]);
SetPrecision[appr,16]];