Trapezoidal Rule

Consider the result of Theorem 10 with the linear polynomial fit to two data points:

We can take the integral of both sides of this expression. First, we get

> P1:=x->((x-x[1])/(x[0]-x[1]))*f(x[0])+((x-x[0])/(x[1]-x[0]))*f(x[1]);

> x[1]:=x[0]+h:x[0]:=x0:

When we use the notation , for equidistant points, the integral becomes:

> iP1:=Int(P1(x),x=x0..x[1]);

> iP1:=factor(value(iP1));

The error term becomes:

> err:=(D@@2)(f)(xi(x))*(x-x[0])*(x-x[1])/2;

> Ierr:=Int(err,x=x[0]..x[1]);

In this integrand, the term does not change sign over the interval so that we can apply the Weighted Mean Value Theorem for Integrals:

> Ierr=(D@@2)(f)(xi)/2*Int((x-x0)*(x-x0-h),x=x0..x0+h);

with in the interval . This becomes

> Ierr2:=rhs(%);

> Ierr3:=value(Ierr2);

>

The formula then becomes:

> Int(f(x),x=x[0]..x[1])=iP1+Ierr3;

>

which is known as the Trapezoidal Rule:

> pl1:=plot(3+sin(2*x+.5)/2,x=0..3,y=0..4):y0:=3.+sin(2.5)/2.:y1:=3.+sin(4.5)/2.:

> pl2:=plot({[1,t,t=0..y0],[2,t,t=0..y1],[t,(t-1)*(y1-y0)+y0,t=1..2]},x=0..3,y=0..4,color=black):

> with(plots):display({pl1,pl2});

>

>

>

If the second derivative of is bounded over the interval , i.e. , then the maximum error is given by .

For example,

> f:=x->3+sin(2*x+.5)/2;ddf:=(D@@2)(f);

> for i from 0 to 5 do h:=2.^(-i): int1:=int(f(x),x=1..1.+h):int2:=h*(f(1.+h)+f(1.))/2.:print(h, int2,abs(int2-int1),h^3*2./12.); od:

>

> f:=x->3*x-5;

> for i from 0 to 5 do h:=2.^(-i): int1:=int(f(x),x=1..1.+h):int2:=h*(f(1.+h)+f(1.))/2.:print(h, int2,abs((int2-int1)/int1),0); od:

>