Numerical Integration - Simpson's Rule

Review

So far we've approximated the integrand over a subinterval by a constant and by a line. Let's continue the pattern, and use approximation by a parabola. This introduces a bit of complication, however: it takes three points to determine a parabola. Thus, we will approximate over two adjacent subintervals, using the three points , , and . Clearly, this will require that n, the number of subintervals, be even. We then integrate this approximation exactly, from to ; it may be shown that the result is

Then, adding over all of the subintervals, we obtain

which may be written

Application to the Potential Integral

The pseudocode below can be used as a basis for a true computer-code-language program (in Fortran, C, etc.) A well-written code will check for improper values of the input values (e.g., n must be a positive even integer).

        input (positive even) integer n
        input x and y (the coordinates of the location at which 

             Vs is to be computed)
        xsq = x*x
        yl = -1.0
        yu = 1.0
        dy = (yu - yl)/n
        sum = 0.0
       
        for i = 1 thru (n/2)-1 do
            ym = yl + (2*i - 1)*dy
            yp = yl + 2*i*dy
            sum = sum + 2.0*F(ym)/sqrt(xsq + (ym-y)^2) 
                   + F(yp)/sqrt(xsq + (yp-y)^2)
        end do
   
        ans = (2.0*sum + F(yl)/sqrt(xsq + (yl-y)^2)
                       + 4.0*F(yu-dy)/sqrt(xsq + (yu-dy-y)^2)
                       + F(yu)/sqrt(xsq + (yu-y)^2))*dy/3.0
        output ans

In the above, F is the numerator of the integrand of the integral we are approximating.

Trapezoidal Rule vs. Simpson's Rule: Which to Use?

A comparison of the pseudocodes for the trapezoidal and Simpson's rules shows that the Simpson's code is only slightly more complex than the trapezoidal code. Most importantly, each requires the same number of function evaluations (which is important if F is expensive to compute). So is there a good reason for choosing one method over the other?

Recall that the trapezoidal rule was derived by approximating the integrand by a piecewise-linear function. Well, suppose that the integrand is a linear function; then the trapezoidal rule gives an exact result (if we ignore computer roundoff error). Since Simpson's rule was derived by approximating the integrand by a piecewise-quadratic function, we expect an exact result (again, ignoring roundoff error) if the integrand is a quadratic function. But Simpson's does an even better job than that: by a fortunate accident, it turns out that Simpson's rule gives an exact result for integrating a cubic function. So if the integrand is a "well-behaved" function (for example, if it is continuously differentiable throughout the interval of integration), then Simpson's rule is to be preferred. The only other possible drawback is the requirement that n be even; this requirement is rarely of any consequence whatever (and can be circumvented by an application of the so-called "three-eighths rule" -- but that's another story).

Go back to the Electric Potential module