Thursday, July 28, 2005

Deriving Simpson's Rule

kw: mathematics, numerical integration, derivations

I've been teaching basic calculus to my son over the summer (can't believe it is almost over). We've had a look at the anti-derivative, AKA the indefinite Integral, and its use for summing a continuous function. Of course, if the function you want to integrate isn't a simple algebraic function, and you can't find it in a "Table of Integrals" anywhere, what do you do? Use numerical integration. So, we did a couple of sessions with the Trapezoid Rule, and how it can be accelerated by repeated use. Then, I told him we'd look at Simpson's rule, which has better accuracy and converges more quickly when accelerated.

I went looking on the web for a derivation, and found dozens of pages that say, "Now, for more efficient summations, we'll introduce Simpson's Rule. We present it without derivation..." or "...derivation of Simpson's Rule is beyond the scope..." So, I thought it over and decided, it can't be that hard...and it isn't! A simple picture illustrates the method. (Click for a larger image)



The figure shows how we analyze the total area in three sections:

  • At the top, an area between the parabolic segment we are fitting through the three points (X0,Y0), (X1,Y1), and (X2,Y2) and the base of the parallelogram shown, with height Hp and width W = X2-X0.
  • Just beneath, a Triangle with Altitude Ht and the same width.
  • At the bottom, the Rectangle that fits between the Triangle and the X-axis, having the same width.

In Analytical Geometry (part of Pre-Calculus at most High Schools), we learn that the sagitta of any parabolic section (Hp in this diagram) of a certain width is the same, anywhere on the parabola. Further analysis shows that the area between the parabolic section and its enclosing parallelogram is constant. Thus the area between the parabolic section shown and the base of the parallelogram—the sloping line through the midpoint (X1,Y1)—is W·Hp/3.

Now we determine the heights of the three sections, Hp, Ht, and Hr:

  • Hp = ½(Y0+Y2) - Y1 = ½Y0 - Y1 + ½Y2
  • Ht = Y2 - Y0
  • Hr = Y0 - Hp = ½Y0 + Y1 - ½Y2

This allows us to find the three areas:

  • Ap = Hp·W/3 = W·(Y0/6 - Y1/3 + Y2/6)
  • At = Ht·W/2 = W·(Y2/2 - Y0/2)
  • Ar = Hr·W = W·(Y0/2 + Y1 - Y2/2)

Adding these, it is best to use W/6 for the common prefix, so

  • A = W·[(Y0 - 2Y1 + Y2) + (3Y2 - 3Y0) + (3Y0 + 6Y1 - 3Y2)]/6

From which we find:

  • A = W·(Y0 + 4Y1 + Y2)/6

And this is Simpson's Rule.

No comments: