The Basic Formula

Recall that the trapezoidal approximation of  involved the weighted sum

( +2*( +...+ ))

where the interval [a,b] was divided into n  subintervals of length  =    and , ,...  are the endpoints of the subintervals. In the trapezoidal rule the function  was approximated by a linear function having the same values at the endpoints of each subinterval. In Simpson's rule we approximate the function  on subintervals by quadratic functions or parabolas. Since it takes three points to determine a unique parabola, we work with pairs of subintervals, and therefore, it is necessary to assume that the number of subintervals, n, is even . Consider the first two subintervals having endpoints ,  and . Let ,  and  denote the values of the function f  at the points ,  and  respectively. Let  denote the unique quadratic function whose graph passes through the points ( , ), ( ) and ( ). The idea behind Simpson's Rule is to approximate  by . It turns out that the last integral can be calculated in terms of h , ,  and . The result is

 = ( + + )

So  is approximated by ( + + ). The next pair of intervals involves the endpoints ,  and . By the result above  is approximated by ( + + ). Adding these results we approximate  by

( )

Adding up the approximations for all pairs of subintervals leads to the following approximation for :

( +...+ )

Again we obtain a weighted sum involving the values of the function f  at the points in the partition of the interval [a,b]. This approximation is known as Simpson's approximation with n  subintervals. Remember that n  is assumed to be even .

Example with 10 Subintervals

We will study the integral of the function  on the interval [0,5] with 10 subintervals.

>	restart:f:=x->x^2*exp(-x);

>

a:=0;

>

b:=5;

>

n:=10;

>	h:=(b-a)/n;

>	x:=i->a+i*h;

We now calculate the Simpson approximation for  with 10 subintervals.

>	SimpApprox:=h/3*(f(x(0))+4*f(x(1))+2*f(x(2))+4*f(x(3))+2*f(x(4))+4*f(x(5))+2*f(x(6))+4*f(x(7))+2*f(x(8))+4*f(x(9))+f(x(10)));

>	evalf(SimpApprox);

Since the integral can be calculated directly with antiderivatives, we compute the actual value of the integral and compute the error.

>	TrueValue:=evalf(int(f(x),x=0..5));

>	SimpError:=evalf(abs(SimpApprox-TrueValue));

We also compute the percentage error.

>	PercentError:=SimpError/TrueValue;

Discuss the result. In an earlier project the same integral was approximated using the trapezoidal rule. Compare the results.

The General Formula

We need a general formula for the Simpson Approximation with n  subintervals. The trick is to group the terms having the same coefficients. If x(i) is defined as , the first and last terms have coefficient 1, the odd terms have coefficient 4 and the even terms not including the first and last terms have coeffieient 2. We let . It follows that the Simpson approximation with n  subintervals can be written as

( + + )

In order to implement this formula in Maple, we note that x(i) must be defined in terms of both i  and n . Also h  is now a function of n.

>	restart:f:=x->x^2*exp(-x);

>

a:=0;

>

b:=5;

>	h:=n->(b-a)/n;

>	x:=(i,n)->a+i*h(n);

>	y:=(i,n)->f(x(i,n));

>	SimpApprox:=n->evalf((h(n)/3)*(y(0,n)+y(n,n)+2*sum(y(2*i,n),i=1..n/2-1)+4*sum(y(2*i-1,n),i=1..n/2)));

>	SimpApprox(10);

>	TrapApprox:=n->evalf((h(n)/2)*(y(0,n)+y(n,n)+2*sum(y(i,n),i=1..n-1)));

>	TrapApprox(10);

Determine the Simpson approximation for the same function on the same interval with 16 subintervals. Also calculate the error and percent error with 16 subintervals.

Project