The Trapezoid Rule

This demonstration shows some examples of the numerical properties of the trapezoid rule.

First, we code the formula

Trap(n) = ((b-a)/(2n))(f((Subscript[x, 0])) + 2 f((Subscript[x, 1])) + ... + 2 f((Subscript[x, n-1])) + f((Subscript[x, n]))).

This is an interesting place to demonstrate Mathematica coding, and you can do it several ways. My favorite is probably to sum the middle values separately from the first and last, like

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

because this looks most like the math.

This requires either the use of a palette or some keyboard shortcuts. Because this is an interactive online class now, let’s go to the “palettes” menu and open up “Basic Math Assistant”. The symbols there can be used in Text boxes to help you typeset math. For instance, we could typeset the identity

 1 + 2+... + n = (∑_(i=1)^n)i = ((n(n+1))/2)

writing the sum by picking the “sum” object from the typesetting menu and filling in the boxes. However, if you put that same “sum” object in the Code box, and shift-enter:

(∑_(i=1)^n)i

it actually evaluates the sum!

A first example.

We now consider an easy example:

	(∫_0^(π/2))sin x dx = (-cos π/2) - (-cos 0) = 1.

Mathematica is happy to do this integral for us:

ExactAnswer = Integrate[Sin[x],{x,0,π/2}]

We can approximate this with

Trap100 = Trap[Sin,{0.0,N[π/2,16]},100]

and compute the number of correct digits as usual with

CorrectDigits[x_,xhat_] :=1.0 Abs[Log10[Abs[(x-xhat)/xhat]]];

In this case, we get:

CorrectDigits[ExactAnswer,Trap100]

As we’d expect from a method with error O((h^2)), the number of correct digits (log error) is a linear function of log(datapoints), with slope very close to 2:

ErrorDataExact = Table[{n,CorrectDigits[1,Trap[Sin,{0.0,1.0N[π/2,50]},10^n]]},{n,1,12}];

ErrorDataExact

lm = LinearModelFit[ErrorData,x,x][x]

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

Show[{Plot[{lm,16},{x,0.5,8.5},],ListPlot[ErrorDataExact]},]

Note: The  symbols just conceal some lengthy plotting options. You do this by selecting “iconize” under “Edit”. If you click, or “uniconize”, you can see all of them written out, but I think it’s just kind of distracting.

There are several things to notice about this plot. The first is that you can actually get very high accuracy doing this calculation in machine precision arithmetic: numerical differentiation is numerically unstable (you lose digits of accuracy to roundoff error when h is small), but numerical integration is numerically stable (you only lose accuracy as you get very close to the 16 digit threshold).

The second thing to notice is that the trapezoid rule requires a fantastic amount of computation to get close to full precision-- at (10^8) gridpoints, we’re doing several hundred million arithmetic operations to compute the sum.

A theoretical error bound.

We saw in the notes that we can derive an error bound for the trapezoid rule in the form:

	||Trap(f,n) - (∫_a^b)f(x) dx|| ≤ (1/12)M ((((b-a))^3)/(n^2))
	
where ||f''(x)|| ≤ M on (a,b).

We can use this in two ways: to solve for the n required to achieve a given error bound ϵ. The algebra is easy enough to do by hand, but it doesn’t hurt to practice how to do it with Mathematica.

Solve[(1/12)M ((((b-a))^3)/(n^2))== ϵ,n]

gives you a weird looking solution! The problem is that Mathematica doesn’t know that (Subscript[b,  ])> a. 
We can make it look a little better (and make it clear which solution to pick!) by replacing b -a with |b-a|.

Solve[(1/12)M ((((Abs[b-a]))^3)/(n^2))== ϵ,n]

We’re now ready to write this down as a function call-- just remember the last point, which is that n must be an integer. So we’ll take the Ceiling of the result to make sure we’re rounding up:

RequiredTrapezoidIntervals[M_,a_,b_,ϵ_] := N@Ceiling[(((Sqrt[M]) ((Abs[-a+b])^(3/2)))/(2 (Sqrt[3]) (Sqrt[ϵ])))]

For our test integration (∫_0^(π/2))sin x dx we know that ((d^2)/(dx^2))sin x = -sin x, so M=1 is the smallest bound that will work over (0,(π/2)), so we can try

RequiredTrapezoidIntervals[1,0,π/2,0.001]

and indeed when we measure the error

Abs[(∫_0^(π/2))Sin[x]\[DifferentialD]x -  Trap[Sin,{0.0,π/2},18]]

We can similarly take it the other way and derive an error bound from the number of intervals:

TrapezoidErrorBound[M_,a_,b_,n_] := (1/12)M ((((Abs[b-a]))^3)/(n^2))

TrapezoidErrorBound[1.0,0.0,π/2.0,18]

Note: The predicted error is less than 0.001 because 18 is the ceiling of the number of intervals required.

How good is the theoretical error bound?

We’re now going to get a sense of how tight the error bound is for the trapezoid rule by trying some experiments. We’ll do this making a log plot of the actual error and the error bound for various n for our test integral. Because there isn’t a “LogLogListLinePlot”, we’ll generate tables of “log of # intervals” vs “log of error” and plot them with “ListLinePlot”.

ErrorBoundData = Table[{n,Log10[TrapezoidErrorBound[1.0,0.0,π/2,Ceiling[(10^n)]]]},{n,1,5,0.2}];

TrapezoidData = Table[{n,Log10[Abs[1 - Trap[Sin,{0.0,π/2.0},Ceiling[(10^n)]]]]},
{n,1,5,0.2}];

We’ll now take a look at our data.

ListLinePlot[{ErrorBoundData,TrapezoidData},]

A (mysterious) example.

I encourage you to retry the code above with different test integrals. You’ll generally see that the error bound is pretty close to the actual error. Here’s an example where that isn’t true. This integral has an exact answer in terms of special functions:

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

Mathematica is happy to get as many digits as we’d like here:

N[(∫_0^(4π))(E^(Sin[x]))\[DifferentialD]x,20]

We can compute the second derivative of the function:

FullSimplify[D[Exp[Sin[x]],x,x]]

and bound the absolute value of the second derivative:

Maximize[{Abs[Exp[Sin[x]](((Cos[x])^2) - Sin[x])],0≤ x ≤ 4π},x]

We can now gather data in a table as before:

SecondErrorBoundData = Table[{n,Log10[TrapezoidErrorBound[E,0.0,4π,Ceiling[(10^n)]]]},{n,0.2,5,0.2}];

ErrorFloor =  Table[{n,-16},{n,0.2,5,0.2}];

SecondTrapezoidData = Select[Table[{n,Log10[Abs[4 π BesselI[0,1] - Trap[Exp[Sin[#]]&,{0.0,4.0π},Ceiling[(10^n)]]]]},
{n,0.2,5,0.2}],NumberQ[#[[2]]]&];

and plot it!

ListLinePlot[{SecondErrorBoundData,SecondTrapezoidData,ErrorFloor},]

This is extremely weird! We can see that the error is always less than the error bound (good, math still works). But it’s a lot less; in fact with only 20 intervals or some it’s already close to the minimum possible error and it simply stays near there. If we repeat the error computation, but have Mathematica do everything with absurd precision, we can see the trend more clearly:

SecondExactErrorBoundData = Table[{n,Log10[TrapezoidErrorBound[E,0.0,4π,n]]},{n,10,150,1}];

ExactSecondTrapezoidData = Block[{$MaxExtraPrecision = 20000},
Table[{n,Log10[N[Abs[4 π BesselI[0,1] - Trap[Exp[Sin[#]]&,{0,4π},n]],2500]]},
{n,10,150,1}]
];

ListLinePlot[{SecondExactErrorBoundData,ExactSecondTrapezoidData},]

You’re reading those gridlines right-- with 49 intervals, the trapezoid rule gets about 180 correct digits!

N[4 π BesselI[0,1],185]

N[Trap[Exp[Sin[#]]&,{0,4π},47],185]

So what’s going on? There should be a theory which explains this weird phenomenon...