4.5. Definite Integrals, Part 3: The (Composite) Simpson’s Rule and Richardson Extrapolation

References:

• Section 5.2, Newton-Cotes Formulas for Numerical Integration, of [Sau22], Sub-sections 5.2.2 and 5.2.3.

• Chapter 4, Numerical Differentiation and Integration, of [BFB16], Sections 4.3 and 4.4.

• Section 5.3, Simpson’s Rule and Newton-Cotes Rules in [CK13].

4.5.1. Introduction

The Composite Simpson’s Rule can be be derived in several ways. The traditional approach is to devise Simpson’s Rule by approximating the integrand function with a colocating quadratic (using three equally spaced nodes) and then “compounding”, as seen with the Trapezoid and Midpoint Rules.

We have already seen another approach: using a 2:1 weighted average of the Trapezoid and Midpoint Rules with th goal of cancelling their \(O(h^2)\) error terms.

This section will show a third approach, based on Richardson extrapolation: this will set us up for Romberg Integration.

4.5.2. The Basic Simpson’s Rule by Richardson Extrapolation

Equation (4.17) in the section on The Composite Trapezoid and Midpoint Rules says that the error in the composite trapezoid rule is approximately

\[\int_a^b f(x)\ dx - T_n \approx \frac{f'(b) - f'(a)}{12} h^2\]

In fact it can be shown (but not here) that this error can be expanded in even powers of \(h\), giving

\[T_n = \int_a^b f(x)\ dx + \frac{f'(b) - f'(a)}{12} h^2 + c_4 h^4 + \cdots c_{2p} h^{2p} + o(h^{2p})\]

so long as \(f\) is sufficiently diferentiable.

For now, we just use the first part of this:

\[T_n = I + c_2 h^2 + O(h^4)\]

where \(I = \int_a^b f(x)\ dx\) is the integral to be approximated (the “Q” in the section on Richardson Extrapolation) and \(c_2 = (f'(b) - f'(a))/12\).

Thus, using the “n form” of Richardson Extrapolation in Equation (4.8) with \(p=2\) gives a new approximation

(4.19)\[S_{2n} = \frac{4T_{2n} - T_n}{4-1}\]

The name is the same as given in (4.18), because it can be shown that these give the same value.

To start, look at the simplest case of this:

\[ S_{2} = \frac{4 T_{2} - T_1}{3} \]

Definfing \(h = (b-a)/2\), the ingredients are

\[ T_1 = \frac{f(a) + f(b)}{2} (b-a) = \frac{f(a) + f(b)}{2} 2 h = (f(a) + f(b)) h \]

and

\[ T_2 = \left[ \frac{f(a)}{2} + f(a+h) + \frac{f(b)}{2} \right] h \]

so

\[ S_{2} = \frac{\left[ 2f(a) + 4f(a+h) + 2f(b) \right] - [f(a) + f(b)]}{3} h, = \frac{f(a) + 4f(a+h) + f(b)}{3} h \]

which is the basic Simpson’s Rule. The subscript “2” is because this uses two intervals, with \(h=(b-a)/2\)

4.5.3. Accuracy and Order of Precision of Simpson’s Rule

Rather than derive this the traditional way — by fitting a quadratic to the function values at \(x=a\), \(a+h\) and \(b\) — this can be confirmed “a postiori” by showing that the degree of precision is at least 2, so that it is exact for all quadratics. And actually we get a bonus, thanks to some symmetry.

For \(f(x) = 1\), the exact integral is \(I = b-a, = 2h\), and also

\[ S_2 = \frac{1 + 4 \times 1 + 1}{3} h, = 2h \]

For \(f(x) = x\), the exact integral is \(I = \int_a^b x \, dx = [x^2/2]_a^b = (b^2 - a^2)/2 = (b-a)(b+a)/2 = (a+b)h\)

and

\[ S_2 = \frac{a + 4(a+b)/2 + b}{3} h = \frac{a + 2(a+b) + b}{3} h = (a+b)h \]

However, it is sufficient to traslate the domain to the symmetric interval \([-h, h]\), so redo the \(f(x)=x\) case this easier way:

The exact integral is \(\int_{-h}^h x \, dx = 0\) (because the function is odd)

\[ S_2 = \frac{-h + 4 \times 0 + h}{3} h = 0 \]

For \(f(x) = x^2\), again do it just on the symmetric interval \([-h, h]\): the exact integral is \(\int_{-h}^h x^2 \, dx = [x^3/3]_{-h}^h = 2h^3/3\) and

\[ S_2 = \frac{(-h)^2 + 4\times 0^2 + h^2}{3} h = 2h^3/3 \]

So the degree of precision is at least 2, as expected.

What about cubics? Check with \(f(x)=x^3\), again on interval \([-h, h]\).

Almost no calculation is needed: symmetry does it all for us:

• on one hand, the exact integral is zero due to the function being odd on a symmetric interval: \(\int_{-h}^h x^3 \, dx = [x^4/4]_{-h}^h = 0\)

• on the other hand,

\[ S_2 = \frac{(-h)^3 + 4 \times 0^3 + h^3}{3} h = 0 \]

The degree of precision is at least 3.

Our luck ends here, but looking at \(f(x)=x^4\) is informative:

For \(f(x) = x^4\),

• the exact integral is \(\int_{-h}^h x^4 \, dx = [x^5/5]_{-h}^h = 2h^5/5\)

• on the other hand

\[ S_2 = \frac{(-h)^4 + 4\times 0^4 + h^4}{3} h = 2h^5/3 \]

So there is a discrepancy of \((4/15) h^5, = O(h^5)\).

This Simpson’s Rule has degree of precision 3: it is exact for all cubics, but not for all quartics.

The last result also indicate the order of error:

\[S_2 - I = O(h^5)\]

Just as for the composite Trapezoid and Midpoint Rules, when we combine multiple simple Simpson’s Rule approximations with \(2n\) intervals each of width \(h= (b-a)/(2n)\), the error is roughly multiplied by \(n\), so \(h^5\) goes to \(n h^5, = (b-a)h^4\), leading to

\[S_{2n} - I = O(h^4)\]

4.5.4. Appendix: Deriving Simpson’s Rule by the Method of Undetermined Coefficients

We wish the determine the most accurate approximation of the form

\[ \int_a^b f(x) \, dx \approx \left[ C_1 f(a) + C_2 f(c) + C_3 f(b) \right] h \]

where \(c\) is the midpoint, \(c = (a+b)/2\)

This wilk be done by the first, “hardest” method: inserting Taylor polynomial and error terms, but to make it a bit less hard, we can consider just the symmetric case \(a=-h\), \(b=h\), \(h= (b-a)/2\) by making the change of variables \(x \to x-c\).

As we now know that this will be exact for cubics, use third order Tayloe polynomials:

\[ f(\pm h) = f(0) \pm f'(0) h + \frac{f''(0)}{2} h^2 \pm \frac{f'''(0)}{6} h^3 + \frac{f''''(\xi_\pm)}{24} h^4 \]

(Note that the special values \(\xi_\pm\) are in general different for the “\(+h\)” and “\(-h\)” cases.

As usual, gather terms with the same power of \(h\):

\[\begin{split} \begin{split} S_2 &= h f(0) (C_1 + C_2 + C_3) \\&+ h^2 f^{(1)}(0) (-C_1 + C_3) \\&+ h^3 f^{(2)}(0) (C_1/2 + C_3/2) \\&+ h^4 f^{(3)}(0) (-C_1/6 + C_3/6) \\&+ h^5 (C_1 f^{(4)}(\xi_-) + C_3 f^{(4)}(\xi_+))/24 \end{split} \end{split}\]

The exact integral can also be computed with Taylor’s formula:

\[\begin{split} \begin{split} I = \int_{-h}^h f(x)\, dx &= \int_{-h}^h \left[ f(0) + Df(0) x + \frac{D^2f(0)}{2}x^2 + \frac{D^3f(0)}{6}x^3 + \frac{D^4f(24)}{2}x^4 + \frac{D^5f(\xi_x)}{120}x^5 \right]\, dx \\&= 2 h f(0) + \frac{D^2f(0)}{3} h^3 + \frac{D^3f(0)}{12} h^5 + O(h^6) \end{split} \end{split}\]

(Symmetry causes all the odd power integrals to valish.)

so the error is

\[\begin{split} \begin{split} S_2 - I &= h f(0) (C_1 + C_2 + C_3 - 2) \\&+ h^2 Df(0) (-C_1 + C_3) \\&+ h^3 D^2f(0) (C_1/2 + C_3/2 - 1/3) \\&+ O(h^5) \end{split} \end{split}\]

The best possibility is setting the coeficients of \(h\), \(h^2\) and \(h^3\) to zero:

\[\begin{split} \begin{split} C_1 + C_2 + C_3 &= 2 \\-C_1 + C_3 &= 0 \\C_1/2 + C_3/2 &= 1/3 \end{split} \end{split}\]

Symmetry helps, as the “\(h^2"\) equation \(-C_1 + C_3 = 0\) gives \(C_3 = C_1\), leaving

\[C_1 = 1/3, \quad 2C_1 + C_2 = 2\]

and thus

\[ C_1 = C_3 = 1/3, \quad C_2 = 4/3 \]

as claimed above.

4.5.5. Exercises

Exercise 4.11

Consider the goal of approximating \(\displaystyle I = \int_0^1 f(x)\, dx\) with three evaluations of the function \(f\), by

\[c_0 f(0) + c_1 f(1/2) + c_2 f(1).\]

Determine the coefficients that give the highest possible degree of precision.

Then rescale this formula to give an approximation

\[\int_a^b f(x) dx \approx [ C_0 f(a) + C_1 f((a+b)/2) + C_2 f(b) ] (b-a)\]

What is its degree of precision?

Where have you seen this before?