19.3 Why Do these Rules Work?

Suppose your integrand is f(x) and f has a power series expansion about the point q whose coefficients do not go wild. We can then write, for x near q

Suppose now we want the area under f between q - d and q + d, which is an interval of width 2d.

We can figure out the answer: exactly, by integrating each and every term, or approximately, by stopping at some point and summing the results of integrating the terms up to that point.

If we do either, you will notice that the terms with an odd power in (x - q) do not contribute at all. (they just tilt the contributions in a way that averages to nothing in any interval that is symmetric about its center.) The area under f will then consist of 2d*f(q) from the first term and then from the third term and then another contribution proportional to from the fifth term, and so on.

Now suppose instead we approximate by the left hand rule: This approximates the area by 2d*f(q - d) which has a contribution from the second term, and that contribution is proportional to d2.

This d2 term is wrong, and the relative error is the ratio of it to the first term. When the contribution from the first term is not 0, this relative error is proportional to d.

On the other hand, if we approximate the area by something symmetric between the endpoints of the interval, the odd terms will cancel out, and the error, coming from the third and later terms will be proportional to d3, and the relative error (since the first term is proportional to d) will be proportional to d2.

There are two obvious choices for symmetric "rules", either of which has error that decreases as d2. These are: take half of the value at q - d and half at q + d; or take the value at q itself. The first is the trapezoid rule and the second is the no-name rule. (In practice, the no-name rule is slightly better than the trapezoid rule in accuracy.)

A sneaky trick is to take a combination of the trapezoid and no-name rules chosen to get the contribution to the area of the interval from the quadratic term in the power series expansion exactly right. Then the first term in the error will come from the quartic term in the power series expansion, and it decline relative to the main contribution as d4 does.

The result is called Simpson's rule. We can deduce what it is by noticing what the quadratic term contributes to the true area between q - d and q + d: integrating between these two points gives you (the combines with the fact that you get contributions at both endpoints to give 1). The no-name rule gets no contribution from the quadratic term while the trapezoid rule gets a contribution lacking the that comes from integrating , from the true answer. In other words, it gets a contribution, from that term, of , which is three times too much.

Simpson's rule therefore says: take of the trapezoid rule, and of the no-name rule. Since the trapezoid rule gives weight of d to each end of the interval and the no-name rule gives weight 2d to the center of it, taking the combination gives weight to each end and to the middle point. If you sum this up over a sequence of intervals you get weights of multiplied by 1 4 2 4 2 4 ... 4 1, where the 4's go in the middle of intervals and the 2's at the ends. The separation between adjacent evaluation points here is d.

How can you do better?

Most of the time you won't need to do better, because the error here goes as d4 which, for d on the order 10-4 is on the order of 10-16 and that should make the error quite negligible. If it does not, your function may have some horrible behavior in your region of integration, and no improvement in method will help it much.

However you can do better.

How? Why do you make me repeat myself?

Calm down.

What we can do is to put in a combination of contributions from evaluating f at q, and q plus and minus d, and also at q plus and minus 2d, so that the quartic term comes out exactly right, and the quadratic term does so as well.

And how, how, how?

Well, if you take 16 parts of Simpson's rule with separation d as above, and 1 part of Simpson's rule with separation 2d, the contributions to the error from the quartic term will cancel. So, if we call the result of applying Simpson's rule with separation j, the combination will give the correct answer in up to the terms of order

Why did you divide by 15 here?

I did so because the combination gives the correct answer 16 -1 or 15 times, as well as any error it has. This is 15 times too much.

You mean to say, you can improve the accuracy of Simpson's rule with evaluation point separation d by combining it with a less accurate (separation 2d) Simpson's rule answer?

Exactly so. And moreover, Simpson's rule is the same kind of improvement of the trapezoid rule. In fact if we define to be the result of using the trapezoid rule with evaluation separation d, then . This comes from the fact that the error in , to a first approximation, is of the error in .

So you are saying, that starting with the trapezoid rule for evaluation separation d, you can combine it with less accurate trapezoid rule results at separations 2d and 4d, and end up with something whose error is the order d6? And the trapezoid rule answer is just the average of the left and right hand rule answers?

Right you are. And this means it is quite easy to do all these things.

OK, how?