Simpson's Rule & More

Many problems in calculus textbooks are contrived. In fact, most of them are. In reality, the sorts of integrals people like engineers and applied mathematicians see are not (practically) solvable analytically. You can't use the tricks from Calculus II to antidifferentiate the integrand, because it's simply too complicated. So what do we do?

The Quest for Quadrature

When integrals can't be antidifferentiated, or when it's too time-consuming to do so, the area under the curve can be approximated by a computer. These numerical, rather than analytical, methods are called "quadrature." (There is an algorithm to either symbolically integrate a function, or else show it cannot be expressed in elementary functions, but this is far beyond this post's scope!)

Perhaps the simplest quadrature algorithm is to draw a line between the function's value at the endpoints and calculate the area under it, using simple geometry. This would be both A) "interpolation to a linear function," and B) very inaccurate. But, as always, we can make smaller slices.

Trapezoidal Quadrature

Suppose we take the interval \([a,b]\) and slice it into many smaller pieces, drawing a line between the function values at each. This would give us a series of trapezoids under the curve. It seems like this would nicely approximate the curve, but it turns out that the trapezoidal method is nearly always worse than the technique we originally learned early in calculus for estimating area with sums of rectangles!

Simpson's Rule

So, if we have a gently curved integrand which is too complicated to find an antiderivative for, like this, \(\int_a^b f(x),dx\), and we want to integrate it, then let's try interpolating it to something we can calculate. We'll evaluate the function three times, once at \(a\), then at \(b\), and finally at the midpoint.

This curve can be seen as three points on a piece-wise quadratic. Instead of interpolating to a linear function, we've interpolated to a quadratic, using one extra function evaluation. And the area under this quadratic is approximately: \[\frac{1}{3}(f(a)+4f(c)+f(b))\] where \(c\) is the midpoint of \(a\) and \(b\). Then, to gain accuracy, we can split the interval into \(n\) small sections, and interpolate for each one of those, so that:

\[\int_a^b f(x),dx \approx \sum_{1\leq i\leq n} \frac{x_{i+1} - x_i }{6} (f(x_i) + 4f(\frac{x_i + x_{i+1}}{2}) + f(x_{i+1})) \]

All we need is a function that takes a number of intervals to evaluate (the precision) and the bounds. Since I want to pass in any mathematical operation, I'll also take a closure, which is a way of passing a function as an argument to another function. (Yes, functions are data types too!) Here's my Rust implementation:

fn int_simpson<F>(a: f64, b: f64, n: i32, integrand: F) -> f64
where
    F: Fn(f64) -> f64,
{
    // We want to sum from a to b, in n subintervals, where the width delta_x is (a - b) / n.
    // Each subinterval is [a, a+delta_x] up to b.
    let delta_x = (b - a) / n as f64;
    let mut result = 0.0_f64;
    for i in 1..=n {
        result += delta_x
            * (integrand(a + (i as f64 - 1_f64) * delta_x)
                + 4_f64 * integrand(a - delta_x / 2_f64 + i as f64 * delta_x)
                + integrand(a + i as f64 * delta_x))
            / 6_f64;
    }
    result
}

The complete program can be found here, along with an example of using it, and even asserting that the floating-point result is accurate to a certain precision!

Other Methods

I think it's interesting that we can see the rectangle rule, trapezoid rule, and Simpson's rule all have something in common. They're interpolating to polynomials (zeroeth, first, and second degrees, respectively) and estimating. There are certainly better methods available, especially adaptive quadratures, for functions with large high-order derivatives, but this is a quick-and-dirty result that achieves enough accuracy for most use cases. I hope you enjoyed!

Footnote: I've added LaTeX formatting to this post using the excellent KaTeX JS package. Simply add it to your site header, and it will typeset any LaTeX markup it sees. Instructions are here.