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$ ,

$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)$ .

One thought on “Thoughts on Numerical Integration (Part 5): Derivation of Simpson’s Rule”

Pingback: Thoughts on Numerical Integration: Index – Mean Green Math

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

Published by John Quintanilla

One thought on “Thoughts on Numerical Integration (Part 5): Derivation of Simpson’s Rule”

Leave a Reply Cancel reply

Share this:

Like this:

Related

Published by John Quintanilla

One thought on “Thoughts on Numerical Integration (Part 5): Derivation of Simpson’s Rule”

Leave a Reply Cancel reply