Simpson's Rule and Beyond

kw: numerical integration, algorithms, derivations

On July 28 I posted a derivation of the (1-4-1)/6 method known as Simpson's Rule. It is based on projecting a parabola through three points to determine the area between the curve defined by those points and the X-axis. Simply put, the area under any portion of a parabola defined by three points (X0,Y0), (X1,Y1), and (X2,Y2), for which X1 is the midpoint of the interval (so X2-X1 = X1-X0), and h=X2-X0, is

As = h(Y0+4Y1+Y2)/6

I have also seen in the literature, again presented without proof, that performing a convergence acceleration on a Trapezoid rule integration using one, then two steps in an interval is equivalent to Simpson's rule. This is quite simple to prove.

The Trapezoid rule is based on using the area between a line segment and the X-axis to approximate the area under a curve that passes through the segment's end points. For a single segment, you multiply the width of the interval (h) by the average of the end-points. Formally, A = ½h•(Y0+Y1). If you add a third point at the center of the interval, making two segments, and call the Y at the midpoint Ym, the area is A = ¼h•(Y0+Y1+2Ym).

Now let us rename the points so we can derive the rule: Y0, Y1, Y2 in order. Let us sum areas for Y = 1/X, from X=1 to X=2, to show the convergence. First, we need the three points with which we'll work:

(X0,Y0) = (1,1)
(X1,Y1) = (1.5,0.6666667)
(X2,Y2) = (2,0.5)

Now, the 1-step area is

A1 = ½h•(Y0+Y2) = 0.5*1.5 = 0.75

and the 2-step area is

A2 = ¼h•(Y0+2Y1+Y2) = 0.25*2.833333 = 0.7083333

The actual area under the curve 1/X within X={1,2} is ln(2)=0.6931472.

The error in A1 is +0.0568528 and that in A2 is +0.0151861. We are gratified to find that there is indeed convergence; A2 is much smaller than A1. Now, A1/A2 = 3.74..., and if we were to pursue further analysis (such as by using four steps, 8 steps, etc.), we find that doubling the number of steps in an interval increases the accuracy by a factor of about four. This is second-order convergence, characteristic of the Trapezoid rule.

The second-order nature of the convergence means that we can combine these two areas to produce a third, better estimate:

Ae = (4•A2-A1)/3 = (4*0.7083333-0.75)/3 = 0.6944444, which has an error of 0.001297.

Firstly, note that the combined error is 1/44th of the error in A1. It would take a Trapezoid rule summation with sixteen steps in the same interval to achieve this level of accuracy. Now, let us use the formula for Ae, substituting back the Y values:

Ae = h(4•A2-A1)/3 = h(4•¼(Y0+2Y1+y2)-½(Y0+Y2))/3 = h(2Y0+4Y1+2Y2-Y0-Y2)/6 = h(Y0+4Y1+Y2)/6. This is Simpson's rule. QED

Clue to a later post: Simpson's rule also converges, at an even higher order...