Theory

An alternative route to Sinpson's Rule uses Taylor series expansions about the middle point :

> restart;

> tayl:=convert(taylor(f(xx),xx=x1,4),polynom);

The error term is given by

> err:=(D@@4)(f)(xi(xx))*(xx-x1)^4/4!;

so that,

> f(xx)=tayl+err;

We can then integrate this expansion over the interval :

> expr1:=Int(f(xx),xx=x[0]..x[2])=Int(tayl+err,xx=x[0]..x[2]);

For equidistant points:

> x[1]:=x[0]+h:x[2]:=x[0]+2*h:x1:=x[1]:expr1;

The first four terms can be integrated to yield:

> expr2:=int(tayl,xx=x[0]..x[2]);

If we use a centered three point formula to obtain the second derivative,

> expr3:=(D@@2)(f)(x[0]+h)=(f(x[0])-2*f(x[1])+f(x[2]))/h^2-h^2/12*(D@@4)(f)(xi[2]);

we obtain:

> expr4:=2*h*f(x[0]+h)+h^3/3*rhs(expr3);

> expr5:=expand(expr4);

>

Meanwhile, the error term

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

can be simplified since is always positive over the interval , so that

> expr6:=expr5=(D@@4)(f)(xi[1])/24*Int((xx-x[1])^4,xx=x[0]..x[2]);

which evaluates to

> expr7:=(D@@4)(f)(xi[1])/24*int((xx-x[1])^4,xx=x[0]..x[2]);

>

The total error term then becomes

> ierr1:=expr7-1/36*h^5*`@@`(D,4)(f)(xi[2]);

>

This term can be rewritten using only one unknown value in the range as (see Exercises for its derivation):

> ierr2:=-1/90*h^5*(D@@4)(f)(xi);

So the formula for Simpson's Rule becomes:

>

As estimated in the previous section, the error term is indeed proportional to !

>