Simpson's Rule

Consider again the result of Theorem 10 with the quadratic polynomial fit to three data points:

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

> restart:

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

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

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

> iP2:=Int(P2(x),x=x[0]..x[2]);

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

The error term becomes:

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

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

Since now the term does change sign over the interval we can not apply the Weighted Mean Value Theorem for Integrals. We will derive a more usable error term in the next section.

>

The formula then becomes:

> Int(f(x),x=x[0]..x[2])=iP2+Ierr;

>

which is known as Simpson's Rule. Graphically:

> f:=x->2.+sin(4*x+.5)/2:lagr:=interp([1,1.5,2.],[f(1.),f(1.5),f(2.)],x,quadratic):

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

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

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

For example,

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

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

>

> f:=x->3*x^2+x-5;

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

> f:=x->x^3-4*x^2+5*x-7;

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

> f:=x->-2*x^4+x^3-4*x^2+5*x-7;

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

These results suggest that the error term might be proportional to .

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

which suggest a likely dependence of the step size of the form !

>