Simpson’s Rule

6.3 Simpson’s Rule

If one uses also the function value at the center of the integration interval, one can approximate \(f\) by a quadratic polynomial defined by the following three points: \((a;f(a))\), \(\left(\frac{a+b}{2};f(\frac{a+b}{2})\right)\), \((b;f(b))\). In general, such a polynomial can most conveniently be written as a Lagrange polynomial (cf. Chapter 2):

\begin{equation*} P_2(x)=f(a)L_0(x)+f\left({\textstyle\frac{a+b}{2}}\right)L_1(x)+f(b)L_2(x). \end{equation*}

(6.10)

The approximative value for the integral can be calculated analytically:²

\begin{equation*} \label{eq:S1} \begin{split} S(f)&=\int_a^b P_2(x){\rm d}x \\ &= f(a)\int_a^bL_0(x){\rm d}x+f\left({\textstyle\frac{a+b}{2}}\right)\int_a^bL_1(x){\rm d}x+f(b)\int_a^bL_2(x){\rm d}x \\ &= \frac{b-a}{6}\left[f(a)+4f\left({\textstyle\frac{a+b}{2}}\right)+f(b)\right]. \end{split} \end{equation*}

(6.11)

Example:
Consider again the function \(f(x)=\frac{1}{1+x}\) and its integral over \([0,1]\), now approximated using Simpson’s rule:

\begin{equation*} S(f)=\frac{1-0}{6}\left(\frac{1}{1+0}+4\frac{1}{1+\frac{1}{2}}+\frac{1}{1+1}\right) =\frac{1}{6}\left(1+\frac{8}{3}+\frac{1}{2}\right)=0.69\bar{4}. \end{equation*}

(6.12)

This result has an absolute error of about \(\left|I(f)-S(f)\right|\approx0.0013\) (relative error \(\approx 0.2\) %). Using the trapezoidal rule at the same sampling points it was \(0.015\) (cf. Eq. 6.9), i.e. it was more than 10 times larger. As above, one can improve the approximation by splitting the full interval into \(n\) sub intervals of the same size \(h_n=(b-a)/n\). Since always three sampling points are needed at a time for using Simpson’s rule, the total number of sampling points must be odd, i.e. the number \(n\) of sub intervals must be even. Repeatedly applying Eq. (6.11) to two successive sub intervals so that always the odd-numbered sampling points are at the “center position”, one obtains Simpson’s rule for a total of \(n+1\) sampling points:

\begin{equation*} \label{eq:Sn} \begin{split} S_n(f)&=\frac{2h_n}{6}[f(x_0)+4f(x_1)+f(x_2)]+\frac{2h_n}{6}[f(x_2)+4f(x_3)+f(x_4)]\\ & \quad +\ldots+\frac{2h_n}{6}[f(x_{n-2})+4f(x_{n-1})+f(x_{n})]\\ &=\frac{h_n}{3}[f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+2f(x_4)+\\ & \quad +\ldots+2f(x_{n-2})+4f(x_{n-1})+f(x_{n})]. \end{split} \end{equation*}

(6.13)

Theorem 2: error of Simpson’s rule (without proof)
Let \(f\) be four times continuously differentiable on \([a,b]\). For the approximation of the integral of \(f\) over \([a,b]\), \(I(f)=\int_a^bf(x){\rm d}x\), by Simpson’s rule \(S_n(f)\), \(n \geq 2\) even, it holds that

\begin{equation*} \label{eq:errorSn} |I(f)-S_n(f)|=\frac{n}{2}\,\frac{h_n^5}{90}|f^{(4)}(\epsilon_n)| =\frac{1}{180}\,\frac{(b-a)^5}{n^4}|f^{(4)}(\epsilon_n)|;\ \ \ \epsilon_n\in[a,b]. \end{equation*}

(6.14)

It follows from Theorem 2 that (i) by doubling the number of subintervals in Simpson’s rule the error bound is reduced by a factor of \(\frac{1}{16}\) and (ii) Simpson’s rule, although being based on a quadratic interpolation, provides exact results for polynomials up to degree three, since the error depends on the fourth derivative of the integrand.³
Example:
To estimate the error of the usage of Simpson’s rule in the previous example we need to calculate the fourth derivative of \(f(x)=\frac{1}{1+x}\):

\begin{equation*} f'(x)=-\frac{1}{(1+x)^2}; \ f''(x)=\frac{2}{(1+x)^3}; \ f'''(x)=-\frac{6}{(1+x)^4}; \ f^{(4)}(x)=\frac{24}{(1+x)^5}. \end{equation*}

(6.15)

Applying Eq. (6.14) for the case of \(n = 2\) we obtain

\begin{equation*} |I(f)-S_2(f)|=\frac{1}{180}\,\frac{1}{2^4}\,|f^{(4)}(\epsilon_n)| =\frac{1}{2880}\,\frac{24}{(1+\epsilon_n)^5}. \end{equation*}

(6.16)

To obtain the maximum error, we assume \(\epsilon_n = 0\):

\begin{equation*} |I(f)-S_2(f)| \leq \frac{24}{2880} = \frac{1}{120} = 0.008\bar{3}. \end{equation*}

(6.17)

The true error was considerably smaller, but of the same order of magnitude.
Remark: summary and higher-order rules
All numerical integration formulae discussed so far (and many more) rely on the mean-value theorem of integral calculus, \(\int^b_a f(x){\rm d}x = (b-a)m\) with \(\min_{[a,b]}f(x) \leq m \leq \max_{[a,b]}f(x)\), by providing an estimate for \(m\). This is obtained by a certain average of the function values taken at the sampling points. Therefore, to specify a certain integration rule, one only needs to know the weighting coefficients used in averaging the function values. These weighting coefficients are given in the following table, together with some additional data, for the numerical integration formulae discussed so far and also for some higher-order rules. The normalization constant is just the sum of the weighting coefficients. In the last column, the error scaling behavior for an application on \(n\) sub intervals is given. The full error terms can be found in the “Bronstein”, Table 7.10, and on the web page given in footnote 2.

name	polynomial	normalization	weighting coefficients					error
	degree	constant						scaling

midpoint	0	1	1					\(n^{-2}\)
trapezoid	1	2	1	1				\(n^{-2}\)
Simpson	2	6	1	4	1			\(n^{-4}\)
“\(\frac{3}{8}\) rule”	3	8	1	3	3	1		\(n^{-4}\)
Boole / Bode	4	90	7	32	12	32	7	\(n^{-6}\)

² See, e.g., http://mathworld.wolfram.com/Newton-CotesFormulas.html

³The latter is due to an error cancellation effect that becomes apparent when constructing Simpson’s rule as a certain weighted average of the trapezoidal rule and the midpoint rule; for details see, e.g., http://en.wikipedia.org/wiki/Simpson%27s_rule, section 1.2, “Averaging the midpoint and the trapezoidal rules”.