The Trapezoid Rule.

This demonstration shows some examples of computing numerical integrals via the trapezoid rule. We see some interesting phenomena.

First, we demonstrate Mathematica's modular programming style. Here h is a local variable. Notice that we can pass a function to Trap as an argument.

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}])

];

f[x_] := Sin[x];
N[Trap[{f,0,Pi,100}],20]

We now estimate the number of points required to compute the integral with given error in terms of a bound on the second derivative.

Npts[M_,eps_,a_,b_] := N[Ceiling[Abs[n /. Solve[(1/12)*M*(b-a)^3/n^2 == eps, n]][[1]]]];

Npts[1,10^(-5),0,Pi]

Ok, let's try it out:

Abs[N[Trap[{f,0,Pi,509}],20] - N[Integrate[Sin[x],{x,0,Pi}],20]]

Now let's try this another way. We can use the error formula to predict the number of correct digits (as a function of n) for a given integrand.

PredictedDigits[{M_,a_,b_,n_}] := N[Abs[Log[10,Abs[(1/12)*(b-a)^3/n^2]]]];

PredictedDigits[{1,0,Pi,509}]

Our next goal is to compare the number of predicted digits of accuracy to the actual number of digits of accuracy for a bunch of n values, to see how far off our bound is.

Digits[x_] := N[Abs[Log[10,Abs[N[x,200]]]]];

N[Pi,16]

PredictedData = Table[{n,PredictedDigits[{1,0,Pi,n}]},{n,100,20000,100}];
ActualData =Table[{n,Digits[Trap[{f,0,Pi,n}] - 2]},{n,100,20000,100}]; 
NActualData = Table[{n,Digits[Trap[{f,0,3.141592653589793,n}] - 2]},{n,100,20000,100}];

Show[ListLinePlot[{PredictedData,ActualData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}],
ListPlot[{PredictedData,ActualData}]]

There are a couple of things worth noting about this plot. First, the number of predicted and actual correct digits is really small (the actual number is a little larger than the predicted number, of course, but by less than one digit in all cases).

We're doing this calculation exactly, so roundoff error is not an issue here. (Here's the plot, if we use the calculation done in machine arithmetic.)

Show[ListLinePlot[{PredictedData, NActualData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}],
ListPlot[{PredictedData,NActualData}]]

We now try another integrand.

g[x_] := Exp[Sin[x]];
Plot[g[x],{x,0,4*Pi},Filling->Axis]

Of course, we need a second derivative bound. So let's compute:

SecondDeriv = D[g[x],x,x];
Plot[SecondDeriv,{x,0,2*Pi},Filling->Axis]
NMaximize[{SecondDeriv,x>0, x < 2*Pi},x]

We now have a function and a bound on its second derivative. So let's recompute the predicted and actual correct digits and see what happens!

Ans = Integrate[g[x],{x,0,2*Pi}]

N[Trap[{g,0,2*Pi,10}],20]

N[Ans,20]

Hang on! We just got a lot more correct digits than we planned to! Let's take this opportunity to make a table and plot things:

ActualGData =Table[{n,Digits[Trap[{g,0,2*Pi,n} ] - Ans]},{n,2,15}]

(Yikes! 40 correct digits with 15 function evaluations!)

PredictedGData = Table[{n,PredictedDigits[{1,0,2*Pi,n}]},{n,2,15,1}]

NActualGData = Table[{n,Digits[Trap[{g,0,2*3.141592653589793,n} ]- Ans]},{n,2,15,1}]

Show[ListLinePlot[{PredictedGData,ActualGData,NActualGData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}],
ListPlot[{PredictedGData,ActualGData,NActualGData}],Plot[16,{x,0,15}]]

This plot is rich in fascinating information. 

The bottom blue curve is the number of predicted digits, which turns out to hover around 1 for this number of function evaluations. 

The pink curve is the number of correct digits if we evaluate with no roundoff error. We can see that this goes to about 40, and that it drops sharply on the even N and essentially doubles on the odd N. 

The brown curve (trapezoid rule with machine arithmetic) follows the pink one when it can, but maxes out at about 16 digits, as expected.

Hmm. The two integrands don't look that different:

Plot[{f[x],g[x]},{x,0,2*Pi},Filling->Axis]

They have comparable second derivatives:

Plot[{-Sin[x],SecondDeriv},{x,0,2*Pi},Filling->Axis]

So what's going on?
Obviously, some deep mathematics is coming into play here...