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

References:

• Sections 5.2.2 and 5.2.3 of Chapter 5 Numerical Differentiation and Integration in [Sauer, 2022].

• Sections 4.3 and 4.4 of Chapter 5 Numerical Differentiation and Integration in [Burden et al., 2016].

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

From the section on The Composite Trapezoid and Midpoint Rules, we have

\[ T_n = \int_a^b f(x) \,dx + \frac{Df(b) - Df(a)}{12} h^2 + O(h^4), = I + c_2 h^2 + O(h^4) \]

where \(I\) is the integral to be approximated (the “Q” in the section on Richardson Extrapolation, and \(c_2 = (Df(b) - Df(a))/12\).

Thus the “n form” of Richardson Extrapolation with \(p=2\) gives a new approximation that I will call \(S_{2n}\):

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

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.