Simpson's rule

Section 14.1 Simpson's rule

Consider the trapezoidal rule on the interval \([-1,1]\) with \(n=2\text{.}\) It has \(h=1\text{,}\) so

\begin{equation} \int_{-1}^{1} f(x)\,dx \approx f(0) + \frac{1}{2}(f(-1)+f(1))\label{eq-trap1}\tag{14.1.1} \end{equation}

Surprisingly, the approximation (14.1.1) can be made much more precise without doing more computations. As in Section 12.4, we use the Richardson extrapolation to improve the accuracy. The same trapezoidal rule with doubled step size \(h=2\) yields

\begin{equation} \int_{-1}^{1} f(x)\,dx \approx \frac{2}{2} (f(-1)+f(1)) \label{eq-trap2}\tag{14.1.2} \end{equation}

Since the error of trapezoidal rule is proportional to \(h^2\) (the rule has second order of accuracy), we expect the second approximation (14.1.2) to have error 4 times greater than the first approximation (14.1.1). To cancel out most of the error, we divide (14.1.2) by 4 and subtract it from (14.1.1). The result is

\begin{equation*} \frac34 \int_{-1}^{1} f(x)\,dx \approx f(0) + \frac{1}{2}(f(-1)+f(1)) - \frac14 (f(-1)+f(1)) \end{equation*}

Hence

\begin{equation} \int_{-1}^{1} f(x)\,dx \approx \frac43 f(0) + \frac{1}{3}(f(-1)+f(1))\label{eq-simps1}\tag{14.1.3} \end{equation}

The approximation (14.1.3) uses the same information about function \(f\) as the trapezoidal rule (14.1.1): the values at -1, 0, 1. But the result is much more accurate. Consider, for example, the integral \(\int_{-1}^1 e^x\,dx = e - 1/e \approx 2.3504\text{.}\) The trapezoidal rule (14.1.1) approximates it by

\begin{equation*} e^0 + \frac{1}{2}(e^{-1} + e^1) \approx 2.5431 \end{equation*}

with the error of 0.1927. The formula (14.1.3) produces

\begin{equation*} \frac43 e^0 + \frac{1}{3}(e^{-1}+ e^1) \approx 2.3621 \end{equation*}

with the error of 0.0117. This is 16 times more accurate. This improvement was achieved by changing the weights (coefficients) attached to the values \(f(-1), f(0), f(1)\text{:}\) for trapezoidal rule they are \(1/2, 1, 1/2\text{,}\) but the formula (14.1.3) uses the weights \(1/3, 4/3, 1/3\text{.}\) The formula (14.1.3) is Simpson's rule (with \(n=2\) subintervals) on the interval \([-1, 1]\text{.}\) On a general interval \([a, b]\) it becomes, via a change of variable,

\begin{equation} \int_{a}^{b} f(x)\,dx \approx \frac{b-a}{6} (f(a) + 4f((a+b)/2) + f(b))\label{eq-simps1a}\tag{14.1.4} \end{equation}

Check that (14.1.4) is consistent with (14.1.3): the midpoint gets greater weight than the endpoints by the factor of 4, and the sum of weights is \(b-a\text{.}\)

Simpson's rule can be used with more subintervals, but their number \(n\) has to be even, since we are combining trapezoidal rule with \(n\) and \(n/2\) subintervals. With \(h=(b-a)/n\) we get

\begin{equation} \int_a^b f(x)\,dx \approx \frac{h}{3}(1,4,2,4,\dots,2,4,1) \cdot (\text{vector of y-values}) \label{eq-simps2}\tag{14.1.5} \end{equation}

A better way to think of this formula is that we first divide the interval into some number of equal parts and then use Simpson's rule on each part. The result is called “composite Simpson's rule” in contrast to the “simple Simpson's rule” that we started with.

To summarize: Simpson's rule is the result of Richardson extrapolation of the Trapezoidal rule. Following the process in Section 13.4 one can find that its order of accuracy is 4, meaning the error is proportional to \(h^4\text{.}\)