Simpson’s Rule
We have seen how the midpoint rule approximates curves with rectangles and the trapezoidal rule uses straight lines. Simpson’s rule takes this progression further by approximating curves with parabolas (quadratic functions).
This method uses piecewise quadratic functions to follow the curve more accurately than either rectangles or straight-line segments.
How Simpson’s Rule Works
The setup requires partitioning the interval into an even number of subintervals, each with equal width. This even number requirement is essential because we work with pairs of subintervals.
For each pair of subintervals: Over the first pair of subintervals, we approximate [latex]\int_{x_0}^{x_2} f(x)dx[/latex] with [latex]\int_{x_0}^{x_2} p(x)dx[/latex], where [latex]p(x) = Ax^2 + Bx + C[/latex] is the unique quadratic function passing through the three points [latex](x_0, f(x_0))[/latex], [latex](x_1, f(x_1))[/latex], and [latex](x_2, f(x_2))[/latex].
Continuing the pattern: Over the next pair of subintervals, we approximate [latex]\int_{x_2}^{x_4} f(x)dx[/latex] with the integral of another quadratic function passing through [latex](x_2, f(x_2))[/latex], [latex](x_3, f(x_3))[/latex], and [latex](x_4, f(x_4))[/latex]. This process continues with each successive pair of subintervals.

Deriving Simpson’s Rule Formula
To understand the formula for Simpson’s rule, we need to derive the approximation over the first two subintervals. This derivation requires keeping track of several key relationships.
The quadratic function [latex]p(x) = Ax^2 + Bx + C[/latex] must pass through our three points:
[latex]\begin{array}{c}f(x_0) = p(x_0) = Ax_0^2 + Bx_0 + C\\ f(x_1) = p(x_1) = Ax_1^2 + Bx_1 + C\\ f(x_2) = p(x_2) = Ax_2^2 + Bx_2 + C\end{array}[/latex]
We also need these geometric facts:
- [latex]x_2 - x_0 = 2\Delta x[/latex], where [latex]\Delta x[/latex] is the length of a subinterval
- [latex]x_2 + x_0 = 2x_1[/latex], since [latex]x_1 = \frac{(x_2 + x_0)}{2}[/latex]
We can now integrate the quadratic approximation:
[latex]\begin{array}{ccccc}\hfill {\displaystyle\int _{{x}_{0}}^{{x}_{2}}f(x)dx}& \approx {\displaystyle\int _{{x}_{0}}^{{x}_{2}}p(x)dx}\hfill & & & \\ & ={\displaystyle\int _{{x}_{0}}^{{x}_{2}}(A{x}^{2}+Bx+C)dx}\hfill & & & \\ & =\frac{A}{3}{x}^{3}+\frac{B}{2}{x}^{2}+Cx|\begin{array}{c}{}^{{x}_{2}} \\ {}_{{x}_{0}}\end{array}\hfill & & & \text{Find the antiderivative.}\hfill \\ & =\frac{A}{3}({x}_{2}{}^{3}-{x}_{0}{}^{3})+\frac{B}{2}({x}_{2}{}^{2}-{x}_{0}{}^{2})+C({x}_{2}-{x}_{0})\hfill & & & \text{Evaluate the antiderivative.}\hfill\end{array}[/latex]
After factoring and using algebraic identities:
[latex]\begin{array}{ccccc} & =\frac{A}{3}({x}_{2}-{x}_{0})({x}_{2}{}^{2}+{x}_{2}{x}_{0}+{x}_{0}{}^{2}) +\frac{B}{2}({x}_{2}-{x}_{0})({x}_{2}+{x}_{0})+C({x}_{2}-{x}_{0})\hfill & & & \\ & =\frac{{x}_{2}-{x}_{0}}{6}(2A({x}_{2}{}^{2}+{x}_{2}{x}_{0}+{x}_{0}{}^{2})+3B({x}_{2}+{x}_{0})+6C)\hfill & & & \text{Factor out}\frac{{x}_{2}-{x}_{0}}{6}.\hfill\end{array}[/latex]
Through careful algebraic manipulation and substitution of our relationships, this eventually simplifies to:
[latex]\int_{x_0}^{x_2} f(x)dx = \frac{\Delta x}{3}(f(x_0) + 4f(x_1) + f(x_2))[/latex]
When we apply the same method to approximate [latex]\int_{x_2}^{x_4} f(x)dx[/latex], we get:
[latex]\int_{x_2}^{x_4} f(x)dx \approx \frac{\Delta x}{3}(f(x_2) + 4f(x_3) + f(x_4))[/latex]
Combining these approximations over the interval [latex][x_0, x_4][/latex]:
[latex]\int_{x_0}^{x_4} f(x)dx = \frac{\Delta x}{3}(f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + f(x_4))[/latex]
This pattern continues as we add pairs of subintervals, leading us to the general Simpson’s rule formula.
Simpson’s Rule
Assume that [latex]f\left(x\right)[/latex] is continuous over [latex]\left[a,b\right][/latex]. Let n be a positive even integer and [latex]\Delta x=\frac{b-a}{n}[/latex]. Let [latex]\left[a,b\right][/latex] be divided into [latex]n[/latex] subintervals, each of length [latex]\Delta x[/latex], with endpoints at [latex]P=\left\{{x}_{0},{x}_{1},{x}_{2},\ldots ,{x}_{n}\right\}[/latex].
Set
[latex]{S}_{n}=\frac{\Delta x}{3}\left(f\left({x}_{0}\right)+4f\left({x}_{1}\right)+2f\left({x}_{2}\right)+4f\left({x}_{3}\right)+2f\left({x}_{4}\right)+\cdots +2f\left({x}_{n - 2}\right)+4f\left({x}_{n - 1}\right)+f\left({x}_{n}\right)\right)[/latex].
Then,
[latex]\underset{n\to \text{+}\infty }{\text{lim}}{S}_{n}={\displaystyle\int }_{a}^{b}f\left(x\right)dx[/latex].