Newton-Cotes Integration Rules.

Deriving the Basic Rules with Interpolating Polynomials.

This demonstration shows some examples of deriving integration rules with Mathematica. The basic idea for all of these rules is the same: we fit a polynomial to the function f[x] at evenly spaced points and then integrate the polynomial. 

The error term comes from manipulating Taylor series.

FitData[n_,a_, h_,f_] := Table[{a + i*h,f[a + i*h]},{i,0,n}];
Poly[n_,a_,h_,f_] := Collect[Simplify[InterpolatingPolynomial[FitData[n,a,h,f],x]],x];
NewtonCotes[n_,a_,h_,f_] := Integrate[Poly[n,a,h,f],{x,a,a+n*h}];

Now we demonstrate how these work for the simple case of Simpson's rule:

FitData[2,a,h,f]

Poly[2,a,h,f]

NewtonCotes[2,a,h,f]

Simplify[NewtonCotes[2,a,h,f]]

The other rules follow the same lines:

Simplify[NewtonCotes[3,a,h,f]]

This is called Simpson's 3/8 rule.

Simplify[NewtonCotes[4,a,h,f]]

This is called Boole's rule.

Simplify[NewtonCotes[5,a,h,f]]

This is just called the six-point closed Newton-Cotes rule.

Composite versions of the Rules

We now work out the composite versions of the various rules. Since the number of intervals has to be divisible by n for the n+1 point Newton-Cotes rule, it is easier to just let kn be the total number of intervals. 

With this convention, we see that

CompositeNewtonCotes[k_,n_,a_,h_,f_] := Sum[NewtonCotes[n,a+i*n*h,h,f],{i,0,k-1}];

So, for instance, for the composite Simpson's rule with 4 intervals, we get

Simplify[CompositeNewtonCotes[2,2,a,h,f]]

and for the composite 3/8 rule with 6 intervals we get

Simplify[CompositeNewtonCotes[2,3,a,h,f]]

Comparing the numerical performance of the various rules.

In order to compare the actual numerical performance of the various rules directly, we'll make sure to use the composite rules with the same total number of intervals. In order to get the number of intervals divisible by 2, 3, 4, and 5, we'll have to use 60 intervals overall.

NCData = Table[N[CompositeNewtonCotes[60/i,i,0,Pi/60,Sin],20],{i,2,5}]

Column[NCData]

We now compute the number of correct digits in each of these integrations:

Digits[a_,b_] := N[Abs[Log[10,Abs[a -b]]]];

NCDigits = Digits[NCData,2]

ListLinePlot[NCDigits,Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]

Notice that we're still not particularly close to full machine precision. If we resurrect the Romberg code from a few lectures ago, we get (for 64 intervals):

Trap[{f_,a_,b_,n_}] := Module[ {h},

h = (b - a)/n; 
h*((1/2)*(f[b] + f[a]) + Sum[f[a + i*h],{i,1,n-1}])

];

RombergArray[{f_,a_,b_,NN_,digits_}] := Module[{A,RA,m,i,n},

(* Code to generate an nxn array of R(n,m) terms. *)

RA = ConstantArray[0.0,{NN,NN}];

(* Mathematica arrays are indexed by 1--which is annoying from our point of view. We solve this by computing the 0th column of the array containing the R(n,0) terms obtained from the Trapezoid rule separately *)

R0 = Table[N[Trap[{f,a,b,2^(i-1)}],digits],{i,1,NN}];

(* Now we fill in the first column of the array by hand *)
	
RA[[2;;NN,1]] = 
Table[
(4/3)*R0[[n]] - (1/3)*R0[[n-1]],
{n,2,NN}];

(* The funny [[ 2 ;; N, 1]] means "the submatrix from row 2 to row N and column 1" *)

(* Now we can fill the remaining columns with a loop *)

Do[

RA [[m+1;;NN,m]] = 

Table[

(4^m/(4^m -1))*RA[[n,m-1]] - (1/(4^m - 1))*RA[[n-1,m-1]],

{n,m+1,NN}] ,

{m,2,NN}];

(* We now want to shift the columns of the array. After all,    the 0th column is stored in R0 and the Nth column is empty. *)

Do[ RA[[All,m]]= RA[[All,m-1]],{m,NN,2,-1}];
RA[[All,1]] = R0;

RA

];

Romberg[{f_,a_,b_,NN_,digits_}]  := RombergArray[{f,a,b,NN,digits}][[NN,NN]];

Digits[Romberg[{Sin,0,Pi,6,20}],2]

which is considerably better than the best of the Newton-Cotes rules.

Error Formulas for the (non-composite) Newton-Cotes Rules

We now work out the error formulas for the Newton-Cotes rules. Since Mathematica won't have a nice way to simplify the error formula for the composite rules, we'll just derive the error formula for the simple rule.

To start with, we just rederive the error formula for Simpson's rule:

SimpsonSeries = Series[NewtonCotes[2,a,h,f],{h,0,5}]

IntegralSeries = Series[Integrate[f[a+x],{x,0,2*h}],{h,0,5}]

IntegralSeries - SimpsonSeries

And now we repeat the process for the other terms:

NewtonCotesError[n_,a_,h_] := Series[Integrate[f[a+x],{x,0,n*h}],{h,0,n +3}]  - Series[NewtonCotes[n,a,h,f],{h,0,n +3}] ;

NewtonCotesError[2,a,h]

NewtonCotesError[3,a,h]

Notice in particular that this is WORSE than the error formula for Simpson's rule!

NewtonCotesError[4,a,h]

NewtonCotesError[5,a,h]

Likewise, since

N[8/945]

is less than

N[275/12096]

we see that the six-point rule is actually worse than the five point rule as well.

Integrands for which polynomial interpolation is particularly poorly behaved.

Recall that when we were first studying polynomial interpolation, we saw that some functions were unusually difficult to model with an interpolating polynomial based on equally spaced gridpoints.

ff[x_] := 1/(1 + (10*x)^2);

a = -1;
b =    1;
n  =   40;

EqualNodes := Table[{a + (b-a)i/n,ff[a+(b-a)i/n]},{i,0,n}]
EqP := InterpolatingPolynomial[EqualNodes,x];

Here's the equally-spaced nodes interpolating polynomial:

Show[Plot[{ff[x],EqP},{x,a,b}],ListPlot[EqualNodes,PlotStyle->{PointSize[Large]}],PlotRange->{{-1,1},{-2,2}}]

If we reduce the number of nodes to six, we get a better looking polynomial graph, but you can still see there's substantial error.

ff[x_] := 1/(1 + (10*x)^2);

a = 0;
b = 1;
n  =  6;

EqualNodes := Table[{a + (b-a)i/n,ff[a+(b-a)i/n]},{i,0,n}]
EqP := InterpolatingPolynomial[EqualNodes,x];

Show[Plot[{ff[x],EqP},{x,-0.1,1.1},Filling->{2}],ListPlot[EqualNodes,PlotStyle->{PointSize[Large]}]]

NewNCData = Table[N[CompositeNewtonCotes[60/i,i,0,1/60,ff],20],{i,2,5}]

Column[NewNCData]

N[Integrate[ff[x],{x,0,1}],20]

We now compute the number of correct digits in each of these integrations, revealing the terrible truth: for these integration methods, the blowup in derivatives means that the number of correct digits actually goes DOWN with more complicated NC rules.

NewNCDigits = Digits[NewNCData,Integrate[ff[x],{x,0,1}]]
ListLinePlot[NewNCDigits,Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]

Even the much-hyped (by me!) Romberg algorithm has trouble with this integrand, performing worse than the worst of the Newton-Cotes rules. The problem is that this function simply doesn't obey the nice properties that make the clever cancellation of error terms work in Romberg and these other algorithms.

Romberg[{ff,0,1,6,20}]

N[Integrate[ff[x],{x,0,1}],20]

Digits[Romberg[{ff,0,1,6,20}],N[Integrate[ff[x],{x,0,1}],20]]