Simpson's Rule

The second method that we shall develop is called Simpson's Rule. In this method we approximate f(x) with a collection of arcs from quadratic functions and integrate across each of these.

Once again let P be a partition of [a,b] into n subintervals of equal width, , where for . Here we require that n be even. Over each interval , for , we approximate f(x) with a quadratic curve that interpolates the points , , and .

 
Figure 4:   Approximating the graph of y = f(x) with parabolic arcs across successive pairs of intervals to obtain Simpson's Rule.

Since only one quadratic function can interpolate any three (non-colinear) points, we see that the approximating function must be unique for each interval . Note that the following quadratic function interpolates the three points , , and :

Since this function is unique, this must be the quadratic function with which we approximate f(x) on . Also, if the three interpolating points all lie on the same line, then this function reduces to a linear function. Therefore, since for each i,

By evaluating the integral on the right, we obtain

Summing the definite integrals over each interval , for , provides the approximation

By simplifying this sum we obtain the approximation scheme

 

This method of approximation is known as Simpson's Rule.

Example 9 As we did for the Trapezoidal Rule, let us now apply Simpson's Rule with n = 6 to approximate the value of

Once again we have and

for . Simpson's Rule then provides

to nine decimal places. So here we have a better approximation of the value .

Example 10 Let us apply Simpson's Rule with n = 8 to approximate the value of

As in Example 3, our function , but here we have with

for . Rounded to nine decimal places, this gives

Applying Simpson's Rule then provides

Example 11 Recall Example 5, where we found a sequence of approximations of the value of the definite integral

by applying the Trapezoidal Rule with . Here we shall find a sequence of approximations of the same definite integral by applying Simpson's Rule with n = 2, 4, 6, 8, and 10. Once again, we shall build a table of values of for these n.

We see here that the values are converging to the value at a must faster rate than the values in Example 5.

Example 12 The gamma function, , is defined by the integral

It is not difficult to see that . Furthermore, this function satisfies the property

Note that this implies that if n is a positive integer, then

Thus is a function that interpolates all of the points (n,(n-1)!) for positive integers n. With the substitution

the integral can be rewritten so that

Thus we have a definite integral over the finite interval [0,1].

Now consider the function

Note that for x > 1 this function is defined for all s in [0,1] except at s = 1. Since

we can assign f(1) = 0 as part of our definition of f(s), making this function defined and continuous for all s in [0,1].

Let us apply Simpson's Rule with n = 4 to approximate the value of . Then , and for i = 0, 1, 2, 3, 4. Since x = 3/2, our function f(s) is of the form

Rounded to nine decimal places, this gives

When we apply Simpson's Rule we obtain

The actual value of is (to nine places). Thus -.064088633 is the resulting error.

 
Figure 5:   The curve .

A problem that arises in this example is that as , we have

while

In the evaluation of the function f(s), we find the product of the values of each of these. A considerable amount of error may be introduced at this point.

Just as we have the doubling scheme for the Trapezoidal Rule, in which we incorporate (4) to derive from , we can derive from the information gathered in obtaining . First, let us denote

Then from (6),

Note that

so that

Thus to compute , one must beforehand preserve the values of and obtained from the derivation of . The only additional work then required is to find the sum .

Using this same notation, the sum in (4) can be written

We now wish to examine how accurately Simpson's Rule approximates the definite integral of f(x) on [a,b]. In this case the following has been proven:

Theorem 1.2 Suppose that exists on [a,b]. Then for n an even positive integer,

where

and the error is given by

for some point c in [a,b].

Once again, the number c is not specified, and so we are not able to use this result to determine the exact value of for most functions f(x). However, we can see that the magnitude of the error has the bounds

Example 13 We shall now find bounds on the error for the problem of using Simpson's Rule with n = 6 to approximate the value of

For our function f(x) = 1/x, we have , so that

Now, is strictly decreasing for all x > 0, so that

Example 14 In Example 10 we incorporated Simpson's Rule to approximate the value of

with n = 8. For our function , we have

With some work it can be shown that for . Therefore is increasing on [1,e], and so

The following Mathematica program applies Simpson's rule to approximate the area under a given curve y = f(x) on an interval [a,b]. In its present form, it is written for the purpose of finding the area under the curve on the interval [1,e] starting with n = 6 divisions of the interval. The program will continue to double the number of divisions until the magnitude of the difference between consecutive approximations of the definite integral is less than a given value (in this case .0005).

       f[x_]:=Exp[x]*Log[x]
       a:=1.
       b:=Exp[1.]
       n:=6
       If[Mod[n,2]==0,
          Bn:=0;
          Cn:=Sum[f[a+2*i*(b-a)/n],{i,1,(n-1)/2}];
          sum1:=1.0;
          sum2:=0.0;
          While[Abs[sum2-sum1]>=.0005,
             sum1=sum2;
             delta=(b-a)/n;
             Cn=Bn+Cn;
             Bn=Sum[f[a+(2*i-1)*delta],{i,1,n/2}];
             sum2=delta*(f[a]+f[b]+4*Bn+2*Cn)/3;
             Print[n,"  ",N[sum2,20]];
             n=2*n],
          Print["The value of n, n = ",n,", is not even, and we cannot proceed."]]