The Trapezoid Rule (Revisited)

Now that we have some more insight into the error in the trapezoid rule, we’re going to revisit our examples. First, let’s recall the (new) error formula. 

Theorem. For any k ≥ 1, we have

	Trap(f,n) - (∫_a^b)f(x)dx = (∑_(i=1)^k)(Subscript[c, i]) (h^(2i)) + o((h^(2k)))

where (Subscript[c, i ])= -((Subscript[B, 2i])/((2i)!)) ((f^((2i-1)))(b) - (f^((2i-1)))(a)) and (Subscript[B, i]) is the Bernoulli number.

Bernoulli numbers.

I admit that this formula looks particularly mysterious when you haven’t computed any Bernoulli numbers yet! So let’s take a moment and do a few examples. First, we recall that the Bernoulli numbers are supposed to relate to coefficients in the Taylor expansion of (t/((e^t)-1)) by

	(t/((e^t)-1))= (∑_(i=0)^∞)(Subscript[B, i])((t^i)/(i!)) 
	
We can get a table of these easily enough:

Series[(t/((E^t)-1)),{t,0,12}]

Table[i!,{i,0,12}]*CoefficientList[Normal[Series[(t/((E^t)-1)),{t,0,12}]],t]

The other formula we had for the Bernoulli numbers was:

BernoulliNumber[i_] := (∑_(k=0)^i)((∑_(v=0)^k)(((-1))^v)Binomial[k,v] ((v^i)/(k+1)) )

And we can check that it works (at least for i>0):

Table[BernoulliNumber[i],{i,0,12}]

Notice that it would have worked even for i=0 if we’d figured out how to tell Mathematica to use the convention (0^0)=1.

Somewhat periodic functions.

Recovering our previous code for the trapezoid rule:

Trap[f_,{a_,b_},n_] := ((b-a)/(2n))(f[a] + (∑_(i=1)^(n-1))2f[a+i(((b-a)/n))] + f[b])

We are now going to try a new example:

(∫_0^(2π))Sin[(x/4)]\[DifferentialD]x

Since the first derivative (d/dx)sin (x/4) = (1/4)cos x has f'(b) = (1/4)cos(π/2) = 0 and f'(a) = (1/4)cos(0) = (1/4), we expect the leading term in the error to be of order (h^2).  We can check this by plotting Log[h] against Log[error], doing the computation in an absurd number of digits of accuracy to ensure that roundoff error plays no role here:

TrapezoidError[f_,{a_,b_},n_]:= 
Block[{exact,x,h,$MaxExtraPrecision=20000},
exact = Integrate[f[x],{x,a,b}];
h = ((b-a)/n); 
{N@Log10[h],N@Log10[Abs[exact-Trap[f,{N[a,1000],N[b,1000]},n]]]}
]

TrapezoidErrorData1 = Table[TrapezoidError[Sin[((#)/4)]&,{0,2π},Ceiling[10^x]],{x,1,2,0.1}]

lm = LinearModelFit[TrapezoidErrorData1,lgh,lgh][lgh]

We can see from the plot that the data are indeed very linear, and the slope is close to 2.

Show[{Plot[lm,{lgh,-1.5,0},],
ListPlot[TrapezoidErrorData1,]},]

Now according to our theory, 

	Trap(f,n) - (∫_a^b)f(x)dx =  -((Subscript[B, 2])/(2!)) (f'(b) - f'(a)) (h^2)  -((Subscript[B, 4])/(4!)) (f'''(b) - f'''(a)) (h^4) + ... 

So if we integrate a function where the first derivative agrees at the start and end of the interval, but the third derivative disagrees at the start and end of the interval, the trapezoid rule should converge at rate (h^4). For instance:

f[x_] := (x^2)(((2π-x))^3)

Plot[f[x],{x,0,2π},Filling->Axis]

We can check directly whether

{f'[0] == f'[2π],f'''[0] == f'''[2π]}

Now we construct the Log10[h], Log10[error] table:

TrapezoidErrorData2 = Table[TrapezoidError[f,{0,2π},Ceiling[10^x]],{x,1,2,0.1}]

and fit it to a linear model:

lm2 = LinearModelFit[TrapezoidErrorData2,lgh,lgh][lgh]

We can compare the two error rates in a single plot:

Show[{Plot[{lm,lm2},{lgh,-1.5,0},Filling->Axis,AxesOrigin->{0,0}],
ListPlot[{TrapezoidErrorData1,TrapezoidErrorData2},
PlotLegends->{"Error for (∫_0^(2  π))sin (x/4) dx","Error for (∫_0^(2  π))(x^2)(2π-x)^3dx"},
PlotMarkers->{Automatic,Small}]},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)(h)","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(error)"},ImageSize->Large]

Of course, we can continue this game. The function

g[x_] :=(x^4)(((2π-x))^6)

has

{g'[0] == g'[2π],g'''[0] == g'''[2π],g'''''[0] == g'''''[2π]}

and doesn’t look particularly different from the other integrand:

Plot[g[x],{x,0,2π},Filling->Axis]

But using the trapezoid rule to approximate

Integrate[g[x],{x,0,1}]

should have error O((h^6)). Trying out some examples:

TrapezoidErrorData3 = Table[TrapezoidError[g,{0,2π},Ceiling[10^x]],{x,1,2,0.1}]

and fitting them to a linear model:

lm3 = LinearModelFit[TrapezoidErrorData3,lgh,lgh][lgh]

We now compare all three error rates in a single plot:

Show[{Plot[{lm,lm2,lm3},{lgh,-1.6,0},Filling->Axis,AxesOrigin->{0,0}],
ListPlot[{TrapezoidErrorData1,TrapezoidErrorData2,TrapezoidErrorData3},PlotMarkers->{Automatic,Small},
PlotLegends->{"Error for \!\(\*SubsuperscriptBox[\(∫\), \(0\), \(2π\)]\)sin \!\(\*FractionBox[\(x\), \(4\)]\) dx","Error for \!\(\*SubsuperscriptBox[\(∫\), \(0\), \(2π\)]\)\!\(\*SuperscriptBox[\(x\), \(2\)]\)(2π-x\!\(\*SuperscriptBox[\()\), \(3\)]\)dx","Error for \!\(\*SubsuperscriptBox[\(∫\), \(0\), \(2π\)]\) \!\(\*SuperscriptBox[\(x\), \(4\)]\)(2π-x\!\(\*SuperscriptBox[\()\), \(6\)]\) dx"}
]},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)(h)","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(error)"},ImageSize->Large]

We can see that the error really is O((h^6)) for the g(x) integrand.

Periodic functions.

Now if f(x) is a periodic function with period b-a, then 

(f^((k)))(b) - (f^((k)))(a) = 0 

for all derivatives (not just the odd order ones). It’s really tempting to conclude that our error theorem 

Theorem. For any k ≥ 1, we have

	Trap(f,n) - (∫_a^b)f(x)dx = (∑_(i=1)^k)(Subscript[c, i]) (h^(2i)) + o((h^(2k)))

where (Subscript[c, i ])= -((Subscript[B, 2i])/((2i)!)) ((f^((2i-1)))(b) - (f^((2i-1)))(a)) and (Subscript[B, i]) is the Bernoulli number.

allows us to conclude that for a periodic f, there is NO error in the trapezoid rule. Of course, that’s not quite right: what the theorem really tells you is that “the error drops faster than any power of h”. We can check this explicitly by computing

FullSimplify@Trap[Exp[Sin[#]]&,{0,2π},3]

and

FullSimplify@Trap[Exp[Sin[#]]&,{0,2π},4]

and checking

FullSimplify@Trap[Exp[Sin[#]]&,{0,2π},3] == FullSimplify@Trap[Exp[Sin[#]]&,{0,2π},4]

So there is still some error in a trapezoid rule approximation, even if the function is smooth and periodic. We can get data on how fast this error drops in practice:

TrapezoidErrorData4 = ParallelTable[TrapezoidError[Exp[Sin[#]]&,{0,2π},n],{n,10,150,1}]

Plotting this against the other error estimates, we see that the log-log graph of h vs error is simply not a straight line of any slope! That’s what it means to say that the error drops “faster than any power of h”.

Show[{Plot[{lm,lm2,lm3},{lgh,-1.5,0},Filling->Axis,AxesOrigin->{0,0}],
ListPlot[{TrapezoidErrorData1,TrapezoidErrorData2,TrapezoidErrorData3,TrapezoidErrorData4},PlotRange->{{-1.6,0},All},PlotMarkers->{Automatic,Small},
PlotLegends->{"Error for \!\(\*SubsuperscriptBox[\(∫\), \(0\), \(2π\)]\)sin (x/4) dx","Error for \!\(\*SubsuperscriptBox[\(∫\), \(0\), \(2π\)]\)\!\(\*SuperscriptBox[\(x\), \(2\)]\)(2π-x\!\(\*SuperscriptBox[\()\), \(3\)]\)dx","Error for \!\(\*SubsuperscriptBox[\(∫\), \(0\), \(2π\)]\) \!\(\*SuperscriptBox[\(x\), \(4\)]\)(2π-x\!\(\*SuperscriptBox[\()\), \(6\)]\) dx",
"Error for (∫_0^(2  π))\!\(\*SuperscriptBox[\(e\), \(sin\\ x\)]\)\[DifferentialD]x"}
]},PlotRange->{{-1.5,0},All},AxesOrigin->{0,0},
AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)(h)","\!\(\*SubscriptBox[\(Log\), \(10\)]\)(error)"},ImageSize->Large]

The moral of the story.

What we see from this is that in situations where the function is periodic (or even repeats some derivatives!) at the points a and b, the performance of the trapezoid is much better than you might expect. This is true for even for functions where f(a) ≠ f(b) such as

g2[x_] := 50x+(x^2)(((2π - x))^3)

Plot[g2[x],{x,0,2π},Filling->Axis]

Now we construct the Log10[h], Log10[error] table:

TrapezoidErrorData5 = Table[TrapezoidError[g2,{0,2π},Ceiling[10^x]],{x,1,2,0.1}]

and fit it to a linear model:

lm5 = LinearModelFit[TrapezoidErrorData5,lgh,lgh][lgh]

Now we’d like to achieve these kinds of impressive results for functions which don’t have this special property. We’ll learn shortly that one can use Richardson extrapolation to improve the results of the Trapezoid rule for every function, and that this improvement is better than the simple trapezoid rule for basically everything except periodic functions.