How To Numerically Integrate Like A Mathematician


I have a guest post that I mean to put up shortly which is a spinoff of the talk last month about calculating logarithms. There are several ways to define a logarithm but one of the most popular is to define it as an integral. That has the advantages of allowing the logarithm to be studied using a lot of powerful analytic tools already built up for calculus, and allow it to be calculated numerically because there are a lot of methods for calculating logarithms out there. I wanted to precede that post with a discussion of a couple of the ways to do these numerical integrations.

A great way to interpret integrating a function is to imagine drawing a plot of function; the integral is the net area between the x-axis and the plot of that function. That may be pretty hard to do, though, so we fall back on a standard mathematician’s trick that they never tell you about in grade school, probably for good reason: don’t bother doing the problem you actually have, and instead do a problem that looks kind of like it but that you are able to do.

Normally, for what’s called a definite integral, we’re interested in the area underneath a curve and across an “interval”, that is, between some minimum and some maximum value on the x-axis. Definite integrals are the kind we can approximate numerically. An indefinite integral gives a function that would tell us what the definite integral on any interval would be, but that takes symbolic mathematics to work out and that’s way beyond this article’s scope.

A squiggly function between two vertical green lines that show the interval to be integrated over.
The red curve here represents some function. Any function will do, really. This is just one of them.

While we may have no idea what the area underneath a complicated squiggle on some interval is, we do know what the area inside a rectangle is. So if we pretend we’re interested in the area of the rectangle instead of the original area, good. Take my little drawing of a generic function here, the wavey red curve. The integral of it from wherever that left vertical green line is to the right is the area between the x-axis, the horizontal black line, and the red curve.

A generic function, with rectangles approximating the function's value along the y-axis based on the function's value at the left and at the right of the area we're integrating over.
For the Rectangle Rule, find the area of a rectangle that approximates the function whose integral you’re actually interested in. Shown here are the left-endpoint (yellow line) and the right-endpoint (blue line) approximations, but any rule can be used to find the approximation.

If we use the “Rectangle Rule”, we draw a horizontal line based on the value of the function somewhere from the left line to the right. The yellow line up top is based on the value at the left endpoint. The blue line is based on the value the function has at the right endpoint. We can use any point, although the most popular ones are the left endpoint, the right endpoint, and the midpoint, because those are nice, easy picks to make. (And if we’re trying to integrate a function whose definition we don’t know, for example because it’s the data we got from an experiment, these will often be the only data points we have.) The area under the curve is going to be something like the area of the rectangle bounded by the green lines, the horizontal black line, and the blue horizontal line or the yellow horizontal line.

Drawn this way you might complain this approximation is rubbish: the area of the blue-topped rectangle is obviously way too low, and that of the yellow-topped rectangle is way too high. The mathematician’s answer to this is: oh, hush. We were looking for easy, not good. The area is the width of the interval times the value of the function at the chosen point; how much easier can you get?

(It also happens that the blue rectangle obviously gives too low an area, while the yellow gives too high an area. This is a coincidence, caused by my not thinking to make my function wiggle up and down quite enough. Generally speaking neither the left- nor the right-endpoints are maximum or minimum values for the function. It can be useful analytically to select the points that are “where the function is its highest” and “where the function is its lowest” — this lets you find the upper and lower bounds for the area — but that’s generally too hard to use computationally.)

A generic function divided along the x-axis, with rectangles approximating the function's value along the y-axis.
For the Composite Rectangle Rule, slice the area you want to integrate into a bunch of small strips — not necessarily all the same width — and find the areas of the rectangles that approximate the function in each of those smaller strips.

But we can turn into a good approximation. What makes the blue or the yellow lines lousy approximations is that the function changes a lot in the distance between the green lines. If we were to chop up the strip into a bunch of smaller ones, and use the rectangle rule on each of those pieces, the function would change less in each of those smaller pieces, and so we’d get an area total that’s closer to the actual area. We find the distance between a pair of adjacent vertical green lines, multiply that by the height of the function at the chosen point, and add that to the running total. This is properly called the “Composite Rectangle Rule”, although it’s really only textbooks introducing the idea that make a fuss about including the “composite”. It just makes so much sense to break the interval up that we do that all the time and forget to explicitly say that except in the class where we introduce this.

(And, notice, in my drawings that in some of the regions behind vertical green lines the left-endpoint and the right-endpoint are not where the function gets its highest, or lowest, value. They can just be points.)

There’s nothing special about the Rectangle Rule that makes it uniquely suited for composition. It’s just easier to draw that way. Any numerical integration rule lets you do the same trick. Also, it’s very common to make all the smaller rectangles — called the subintervals — the same width, but that’s not because the method needs that to work. It’s easier to calculate if all the subintervals are the same width, because then you don’t have to remember how wide each different subinterval is.

A generic function, with a trapezoid approximating the function's value. The trapezoid matches the original function at the left and the right ends of the interval we integrate over.
For the Trapezoid Rule, approximate the area under the function by finding the area of a right trapezoid (or trapezium) which is as wide as the interval you want to integrate over and which matches the original function at the left- and the right-endpoints.

Rectangles are probably the easiest shape of all to deal with, but they’re not the only easy shapes. Trapezoids, or trapeziums if you prefer, are hardly a challenge to find the area for. This gives me the next really popular integration rule, the “Trapezoid Rule” or “Trapezium Rule” as your dialect favors. We take the function and approximate its area by working out the area of the trapezoid formed by the left green edge, the bottom black edge, the right green edge, and the sloping blue line that goes from where the red function touches the left end to where the function touches the right end. This is a little harder than the Rectangle Rule: we have to multiply the width of the interval between the green lines by the arithmetic mean of the function’s value at the left and at the right endpoints. That means, evaluate the function at the left endpoint and at the right endpoint, add those two values together, and divide by two. Not much harder and it’s pleasantly more accurate than the Rectangle Rule.

If that’s not good enough for you, you can break the interval up into a bunch of subintervals, just as with the Composite Rectangle Rule, and find the areas of all the trapezoids created there. This is properly called the “Composite Trapezoid Rule”, but again, after your Numerical Methods I class you won’t see the word “composite” prefixed to the name again.

And yet we can do better still. We’ll remember this when we pause a moment and think about what we’re really trying to do. When we do a numerical integration like this we want to find, instead of the area underneath our original curve, the area underneath a curve that looks like it but that’s easier to deal with. (Yes, we’re calling the straight lines of the Rectangle and Trapezoid Rules “curves”. Hush.) We can use any curve that we know how to deal with. Parabolas — the curving arc that you see if, say, you shoot the water from a garden hose into the air — may not seem terribly easy to deal with, but it turns out it’s not hard to figure out the area underneath a slice of one of them. This gives us the integration technique called “Simpson’s Rule”.

The Simpson here is Thomas Simpson, 1710 – 1761, who in accord with Mathematics Law did not discover or invent the rule named for him. Johannes Kepler knew the rule a century before Simpson got into the game, at minimum, and both Galileo’s student Bonaventura Cavalieri (who introduced logarithms to Italy, and was one of those people creeping up on infinitesimal calculus ahead of Newton) and the English mathematician/physicist James Gregory (who discovered diffraction grating, and seems to be the first person to have published a proof of the Fundamental Theorem of Calculus) were in on it too. But Simpson wrote a number of long-lived textbooks about calculus, which probably is why his name got tied to this method.

A generic function, with a parabola approximating the function's value. The parabola matches the original function at the left, middle, and right of the region we're integrating over.
For Simpson’s Rule, approximate the function you’re interested in with a parabola that exactly matches the original function at, at minimum, the left endpoint, the midpoint, and the right endpoint.

In Simpson’s Rule, you need the value of the function at the left endpoint, the midpoint, and the right endpoint of the interval. You can draw the parabola which connects those points — it’s the blue curve in my drawing — and find the area underneath that parabola. The formula may sound a little funny but it isn’t hard: the area underneath the parabola is one-third the width of the interval times the sum of the value at the left endpoint, the value at the right endpoint, and four times the value at the midpoint. It’s a bit more work but it’s a lot more accurate than the Trapezoid Rule.

There are literally infinitely many more rules you could use, with such names as “Simpson’s Three-Eighths Rule” (also called “Simpson’s Second Rule”) or “Boole’s Rule”[1], but they’re based on similar tricks of making a function that looks like the one you’re interested in but whose area you know how to calculate exactly. For the Simpson’s Three-Eighth Rule, for example, you make a cubic polynomial instead of a parabola. If you’re good at finding the areas underneath wedges of circles or underneath hyperbolas or underneath sinusoidal functions, go ahead and use those. You can find the balance of ease of use and accuracy of result that you like.


[1]: Boole’s Rule is also known as Bode’s Rule, because of a typo in the 1972 edition of the Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, or as everyone ever has always referred to this definitive mathematical reference, Abromowitz and Stegun. (Milton Abromowitz and Irene Stegun were the reference’s editors.)

Advertisements

Author: Joseph Nebus

I was born 198 years to the day after Johnny Appleseed. The differences between us do not end there.

8 thoughts on “How To Numerically Integrate Like A Mathematician”

  1. Great post! I do appreciate such neat little tricks even though software does all those calculations for us today.

    BTW your series on logarithms and this post reminded me of the way Feynman has introduced numerical calculations in his physics lectures – he starts solving the equation of motion for a planet revolving around the sun numerically … actually in parallel with explaining what a differential equation is.
    And in the chapter about numbers Feynman introduces e “naturally” by playing with logarithms numerically first: http://www.feynmanlectures.caltech.edu/I_22.html#Ch22-S4

    Like

Please Write Something Good

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 )

Twitter picture

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

Facebook photo

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

Google+ photo

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

Connecting to %s