# 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

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}\):

To start, look at the simplest case of this:

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

and

so

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

**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

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)

**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

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,

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

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:

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

## 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

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:

(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\):

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

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

so the error is

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

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

and thus

as claimed above.