Thoughts on Numerical Integration (Part 5): Derivation of Simpson’s Rule

Numerical integration is a standard topic in first-semester calculus. From time to time, I have received questions from students on various aspects of this topic, including:

  • Why is numerical integration necessary in the first place?
  • Where do these formulas come from (especially Simpson’s Rule)?
  • How can I do all of these formulas quickly?
  • Is there a reason why the Midpoint Rule is better than the Trapezoid Rule?
  • Is there a reason why both the Midpoint Rule and the Trapezoid Rule converge quadratically?
  • Is there a reason why Simpson’s Rule converges like the fourth power of the number of subintervals?

In this series, I hope to answer these questions. While these are standard questions in a introductory college course in numerical analysis, and full and rigorous proofs can be found on Wikipedia and Mathworld, I will approach these questions from the point of view of a bright student who is currently enrolled in calculus and hasn’t yet taken real analysis or numerical analysis.

In the previous post in this series, I discussed three different ways of numerically approximating the definite integral \displaystyle \int_a^b f(x) \, dx, the area under a curve f(x) between x=a and x=b.

In this series, we’ll choose equal-sized subintervals of the interval [a,b]. If h = (b-a)/n is the width of each subinterval so that x_k = x_0 + kh, then the integral may be approximated as

\int_a^b f(x) \, dx \approx h \left[f(x_0) + f(x_1) + \dots + f(x_{n-1}) \right] \equiv L_n

using left endpoints,

\int_a^b f(x) \, dx \approx h \left[f(x_1) + f(x_2) + \dots + f(x_n) \right] \equiv R_n

using right endpoints, and

\int_a^b f(x) \, dx \approx h \left[f(c_1) + f(c_2) + \dots + f(c_n) \right] \equiv M_n

using the midpoints of the subintervals. We have also derived the Trapezoid Rule:

\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{2} [f(x_0) + 2f(x_1) + \dots + 2f(x_{n-1}) + f(x_n)] \equiv T_n

This last approximation was obtained by connecting adjacent points on the curve by line segments, creating trapezoids:

In this post, we will derive Simpson’s Rule. Instead of connecting two adjacent points with line segments, we will connect three adjacent points with a parabola. In the picture below, the points (x_0, f(x_0)), (x_1, f(x_1)) and (x_2,f(x_2)) are connected with one parabola, while the points (x_2, f(x_2)), (x_3, f(x_3)) and (x_4,f(x_4)) are connected with a different second parabola.

Clearly, for this to work, there has to be an even number of subintervals. (By contrast, for the Trapezoid Rule, the Midpoint Rule, or the endpoint rules, the number of subintervals could be even or odd.)

The derivation of Simpson’s Rule is more complicated than the derivation of the Trapezoid Rule because we need to use calculus to find the area under these parabolas. To begin, we make the simplifying assumption that x_1 = 0. Since each subinterval has width h, this means that x_0 = -h and x_2 = h.

To find the area under this parabola, we first need to find the equation of the parabola y = ax^2 + bx + c connecting the three points (-h,y_0), (0,y_1), and (h,y_2). This entails solving a system of three equations in three unknowns:

a(-h)^2 + b(-h) + c = y_0

a(0)^2+b(0) + c = y_1

ah^2 + bh + c = y_2,

or

ah^2 - bh + c = y_0

c = y_1

ah^2 + bh + c = y_2.

While most 3×3 systems are cumbersome to solve, this system is straightforward. Clearly, c = y_1. Also, subtracting the first equation from the third equation yields

2bh = y_2 - y_0, or b = \displaystyle \frac{y_2 - y_0}{2h}

Finally, we solve for a by substituting into the third equation:

ah^2 + \displaystyle \frac{y_2 - y_0}{2h} h + y_1 = y_2

ah^2 + \displaystyle \frac{y_2 - y_0}{2} + y_1 = y_2

ah^2 = \displaystyle \frac{y_0 - y_2}{2} - \frac{2y_1}{2} + \frac{2y_2}{2}

ah^2 = \displaystyle \frac{y_0 - 2y_1 + y_2}{2}

a = \displaystyle \frac{y_0 - 2y_1 + y_2}{2h^2}

Next, we find the integral of y = ax^2 + bx + c between x = -h and x = h:

\displaystyle \int_{-h}^h (ax^2 + bx + c) \, dx = \left[ \frac{ax^3}{3} + \frac{bx^2}{2} + cx \right]^h_{-h}

= \displaystyle \left[ \frac{ah^3}{3} + \frac{bh^2}{2} + ch \right] - \left[ -\frac{ah^3}{3} + \frac{bh^2}{2} - ch \right]

= \displaystyle \frac{2ah^3}{3} + 2ch

= \displaystyle \frac{(y_0 - 2y_1 + y_2)h}{3} + 2y_1h

= \displaystyle \frac{h(y_0 + 4y_1 + y_2)}{3}.

We now turn to the more general case of finding the area under the parabola passing through (x_0,y_0), (x_1,y_1), and (x_2,y_2), where x_1 = x_0 +h and x_2 = x_1 + 2h. Geometrically, it should be clear that this parabola can be obtained from the above parabola by a horizontal translation. Since the area under the curve is not changed by a horizontal translation, the area (and the formula) will be the same.

More formally, if y = ax^2 + bx + c passes through the points (-h,y_0), (0,y_1), and (h,y_2), then y = a(x-x_1)^2 + b(x-x_1) + c will pass through the points (x_0,y_0), (x_1,y_1), and (x_2,y_2). The area under this curve is

\displaystyle \int_{x_0}^{x_2} \left[ a(x-x_1)^2 + b(x-x_1) + c \right] \, dx.

After using the substitution u = x-x_1, this becomes

\displaystyle \int_{-h}^h (au^2 + bu + c) \, du,

which is the same integral that we saw earlier. Therefore,

\displaystyle \int_{x_0}^{x_2} \left[ a(x-x_1)^2 + b(x-x_1) + c \right] \, dx = \displaystyle \frac{h(y_0 + 4y_1 + y_2)}{3}.

Finally, we need to find the sum of the areas under all of these parabolas. Similarly, the area under the parabola passing through (x_2,y_2), (x_3,y_3), and (x_4,y_4) will be \displaystyle \frac{h(y_2 + 4y_3 + y_4)}{3}. So, for the particular example shown above, the total area under the parabolas will be

\displaystyle \frac{h(y_0 + 4y_1 + y_2)}{3} + \frac{h(y_2 + 4y_3 + y_4)}{3} = \frac{h}{3} (y_0 + 4 y_1 + 2 y_2 + 4 y_3 + y_4).

The coefficients of 4 arose from the above integrals, while the coefficient of 2 came from combining the two areas. In general, if there are n subintervals and n is even, then Simpson’s Rule gives the approximation

S_n = \displaystyle \frac{h}{3} \left(y_0 + 4 y_1 + 2 y_2 + 4 y_3 + \dots + 2y_{n-2} + 4 y_{n-1} +  y_{n} \right).

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.