The Peano Kernel Error Formula

We have seen that the error for the trapezoid rule is given by

	Trap(f,n) - (∫_a^b)f(x) dx = (1/12)f''(ξ)((((b-a))^3)/(n^2))=((b-a)/12)f''(ξ) (h^2)

where ξ is some point in (a,b).

Of course, this formula is exact, so it can’t be improved by itself. However, since we don’t know where ξ is in the interval (a,b), to use the formula we’re forced to make the worst-case assumption that ξ happens to fall at the point where |f''(x)| is maximized.

This gives us the error bound

	|Trap(f,n) - (∫_a^b)f(x) dx| ≤ ((b-a)/12)M (h^2)

where M is chosen so that M ≥ |f''(x)| for all x ∈ (a,b).

Max second derivative? Or average second derivative?

Is the error in the trapezoid rule really controlled by the maximum absolute value of the second derivative of the integrand? Let’s do an example: if we consider the two functions

f[x_] := 20(((Sin[(x/π)]))^5)

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

and

g[x_] :=(x^2)/2

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

on the interval [0,1], we see that the functions look about the same (when scaled appropriately!), and the maximum absolute values of their 2nd derivatives are approximately equal:

NMaximize[{Abs[f''[x]],0 ≤ x ≤ 1},x]

NMaximize[{Abs[g''[x]],0≤x ≤ 1},x]

However, if we look at the plot of these second derivatives:

Plot[{Abs[f''[x]],Abs[g''[x]]},{x,0,1},PlotLegends->"Expressions",Filling->Axis,ImageSize->Large]

we see that |f''(x)|is usually much less than |g''(x)|, even though the maximum values are about the same. 
We are now going to see how well the standard trapezoid rule error bound predicts the actual error for these two functions.

First, we give code for the trapezoid rule and the actual (approximation) error:

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

ActualTrapezoidError[f_,{a_,b_},lgn_]:= 
Block[{exactf,x,h,$MaxExtraPrecision=20000},
exactf = Integrate[f[x],{x,a,b}];
h = ((b-a)/(Ceiling[(10^lgn)])); 
{lgn,N[Abs[exactf-Trap[f,{N[a,1000],N[b,1000]},Ceiling[(10^lgn)]]]]}
]

Next, we code the standard bound:

TrapezoidErrorBound[f_,{a_,b_},lgn_] := Block[{M,x},
M = First@NMaximize[{Abs[f''[x]],a ≤ x ≤ b},x];
{lgn,((1/12)M ((((b-a))^3)/((Ceiling[(10^lgn)])^2))) }
]

We can measure the tightness of the bound by taking the ratio of the bound to the actual error. We know this will be at least one, because the error bound is larger. But a value closer to one means the bound is close to the best one possible.

TightnessOfTrapezoidErrorBound[f_,{a_,b_},lgn_]:= 
{lgn,((Last@TrapezoidErrorBound[f,{a,b},lgn])/(Last@ActualTrapezoidError[f,{a,b},lgn]))}

If we plot the results:

ListLinePlot[{Table[TightnessOfTrapezoidErrorBound[f,{0,1},lgn],{lgn,0.2,3,0.2}],
Table[TightnessOfTrapezoidErrorBound[g,{0,1},lgn],{lgn,0.2,3,0.2}]},
PlotMarkers->Automatic,Filling->Axis,AxesOrigin->{0,0},AxesLabel->{"\!\(\*SubscriptBox[\(Log\), \(10\)]\)(n)","Ratio"},PlotLegends->{"\!\(\*FractionBox[\((usual error bound for f)\), \((actual error for f)\)]\)","\!\(\*FractionBox[\((usual error bound for g)\), \((actual error for g)\)]\)"},ImageSize->Large]

We see that the traditional error bound works very well for g (it’s almost perfect!), but overestimates the trapezoid error for f by a factor of about 4.

The Peano kernel bound.

The Peano kernel formulation of the error, as we worked out in the notes, is given by defining the function

	K(x) = (1/2)(x - (Subscript[x, i]))(x-(Subscript[x, i+1]))   for (Subscript[x, i])≤ x≤ (Subscript[x, i+1])
	
where (Subscript[x, 0]),..., (Subscript[x, n]) are the endpoints of the subdivisions in the trapezoid rule. (This function is “the Peano kernel”) and then proving that 

	Trap(f,n) - (∫_a^b)f(x) dx = (∫_a^b)K(x) f''(x) dx

Again, we have to replace something by a bound, and we prove that |K(x)|≤ ((h^2)/8) so that 

         |Trap(f,n) - (∫_a^b)f(x) dx| ≤ ((h^2)/8)(∫_a^b)|f''(x)| dx.
         
This formula is the “Peano kernel bound” for the error in the trapezoid rule.

The Peano kernel bound can be somewhat easier to work with than the usual error bound. For instance, if you know that f''(x) ≥ 0 on (a,b) then 

	(∫_a^b)|f''(x)|dx = (∫_a^b)f''(x)dx = f'(b) - f'(a).
	
It can also provide a better bound than the usual one for functions like f where the second derivative is variable:

PeanoErrorBound[f_,{a_,b_},lgn_] := Block[{h},
h = ((b-a)/(Ceiling[(10^lgn)]));
{lgn,N[((h^2)/8)(∫_a^b)Abs[f''[x]]\[DifferentialD]x]}
]

TightnessOfPeanoErrorBound[f_,{a_,b_},lgn_]:= 
{lgn,((Last@PeanoErrorBound[f,{a,b},lgn])/(Last@ActualTrapezoidError[f,{a,b},lgn]))}

Now we can see how well the Peano bound does for both of our test functions:

ListLinePlot[{Table[TightnessOfPeanoErrorBound[f,{0,1},lgn],{lgn,0.2,3,0.2}],
Table[TightnessOfPeanoErrorBound[g,{0,1},lgn],{lgn,0.2,3,0.2}]},
]

We see that it’s about the same for both functions -- significantly better than the traditional bound for the function f, but significantly worse for the function g. This gives us an idea about when the Peano bound might be useful in practice.