Symmetry and Simpson's Rule

Mathematicians appreciate that exploiting symmetry is a powerful tool. We shall show that in Simpson's rule for estimating integrals, better error bounds and a simpler derivation result from symmetry considerations. These ideas may be applied to all Newton-Cotes formulae.

Suppose that h > 0. Recall that Simpson's rule says that if f is continuous on [a,a+2h] and four times differentiable on (a,a+2h) then  

$$\left\vert\int_{a}^{a+2h} f(x) dx - \frac{h}{3}(f(a) + 4f(a+... ... \leq \frac{h^5}{90}\sup_{a < z < a+2h}\vert f^{(4)}(z)\vert.$$ (1)

Generally Simpson's rule is applied to integrals of the form

$$\int_A^B f(x) \;dx$$

by letting N be a positive even integer, letting h= (B-A)/N and writing

$$\int_A^B f(x) \;dx = \int_{A}^{A+2h} f(x)\;dx + \int_{A+2h}^{A+4h} f(x)\;dx + \cdots + \int_{A+2(N-1)h}^B f(x)\;dx$$

to obtain the familiar formula  

$$\left\vert\int_A^B f(x) dx - \frac{h}{3}\left(f(A) + 4\left... ...q \frac{(B-A)^5}{180N^4}\sup_{A < z < B}\vert f^{(4)}(z)\vert.$$ (2)

This bound on the error is clearly not optimum because it assumes that the maximum value of |f⁽⁴⁾| on all of (A,B) occurs on every subinterval, but it is quick and easy (relatively speaking). However, one cannot find the maximum of |f⁽⁴⁾| on every subinterval before one knows how many subintervals are to be used. Hence the cruder estimate in (2) is used to get an upper bound on the number of subintervals needed.

Put

$$S_N(A,B,f) = \frac{h}{3}\left(f(A) + 4\left(\sum_{k=2}^{N-1}f(A+kh)\right) + f(B)\right)$$

We shall show now that that symmetry gives a painless improvement of (2) to  

$$\left\vert\int_A^B f(x) dx - S_N(A,B,f)\right\vert \leq \fr... ... < B}\left\vert\frac{f^{(4)}(z)+f^{(4)}(A+B-x)}{2}\right\vert.$$ (3)

It is easy to check that  

$$\int_A^B f(x) \;dx = \int_A^B f(A+B-x)\;dx = \int_A^B \frac{f(x) + f(A+B-x)}{2}\;dx$$ (4)

and that the Simpson's rule estimate of each of these integrals is given by exactly the same expression:

$$\frac{h}{3}\left(f(A) + 4\left(\sum_{k=2}^{N-1}f(A+kh)\right) + f(B)\right)$$

Applying Simpson's rule to the integral in the right-hand side of (4) then gives (3).

By way of examples, consider two standard calculations:

$$\int_1^2 \frac{1}{x} dx = \ln 2$$

and

$$\int_0^1 \frac{4}{1+x^2} dx = \pi.$$

In the latter, using symmetry guarentees that the difference between $\pi$ and the Simpson's rule estimate of the integral is no more than (42/180)N^-4 compared with the usual bound of (96/180)N^-4. Thus to obtain an error tolerance of not more than 10^-8 we need only take N = 70 versus 86 with the conventional estimate, or more generally, an improvement by a factor of $\sqrt[4]{42/96} \approx 0.8133$. In the former example, the difference between $\ln 2$ and the Simpson's rule estimate of the integral is no more than (24/180)N^-4 by the usual estimate compared with $(99/(8\cdot 180))N^{-4}$ when symmetry is taken into account, allowing us to use a value of N about $85\%$ smaller to be guarenteed the same accuracy of the estimate. Of course, one can produce examples where there is no improvement and where the improvement is overwhelming.

Considering symmetry also simplifies the derivation and proof of (1). We will show that for an even function $g:[-h, h]\rightarrow(-\infty,\infty)$which is sufficiently smooth,  

$$\left\vert\int_{-h}^h g(x) dx - \frac{2h}{3}(2g(0) + g(h))\r... ...ert \leq \frac{h^5}{90}\sup_{0 < u < h}\vert g^{(4)}(u)\vert.$$ (5)

It is clear from the discussion above that (5) leads quickly to an improved version of (1). As we shall now demonstrate, (5) is easier to derive than (1).

The key to deriving and proving (5) is that even polynomials should be used to interpolate even functions, and that it is easy to find such polynomials if only two interpolation points are needed. In fact, all one has to do is look at the form for linear interpolation as in the trapezoid rule and one sees what to do.

If the trapezoid rule were to be employed to estimate

$$\int_{0}^{h}g(u)du$$

we would replace the integral of g with the integral of the affine function L(u) = g(h)(h-u)/h + g(0)u/h since y=L(x) is an equation of the line which passes through (0,g(0)) and (h,g(h)). It is then clear that since g is even, we should try to change the variable u in L(u) to u² to find a quadratic polymomial with the same property. It is not hard to see that what we want is P(u) = g(0)(h²-u²)/h² + g(h)u²/h², and that y=P(x) passes through all three points (0,g(0)) and $(\pm h,g(\pm h))$. Note that we did not have to solve any systems of equations to find P, nor did we need higher order interpolation formulae, we simply proceeded by analogy with the linear case. It is then easy to check that  

$$\int_{-h}^{h} P(u) du = \frac{1}{h^2}\int_{0}^{h}g(0)(u^2-h^2) + g(h)u^2 \; du = \frac{2h}{3}(2g(0) + g(h))$$ (6)

We see that (6) tells us that to complete the proof of (5) we can estimate |g(u)-P(u)| on [0,h]. We proceed in a well-known fashion. Since g -p is an even function, we try to identify the constant C such that

g(u)-P(u) = Cu²(h²-u²),

where u is a fixed element of (0,h). Our goal is to express C in terms of the some derivative of g. For $x\in[-h,h]$ put E(x) = g(x)-P(x)-Cx²(h²-x²). Note that E continuous on [-h,h], four times differentiable on (-h,h) and even. Therefore, not only do we have E(u) = E(0) = E(h) = 0, but also E'(0) = E⁽³⁾(0) = 0. Therefore we may apply Rolle's Theorem successively to E, E', E'' and E⁽³⁾ and we find that for some $z\in(0,h)$ we have E⁽⁴⁾(z) = 0. Of course, E⁽⁴⁾(z) = g⁽⁴⁾(z) + 24C, so $\vert C\vert \leq \sup_{0<u<h}g^{(4)}(u)\vert/24$. Take C=0 if u=0 or u=h. This tells us that

which is (5).