Skip to main content

Subsection 1.11.3 Simpson's Rule

When we use the trapezoidal rule we approximate the area \(\int_{x_{j-1}}^{x_j}f(x)\dee{x}\) by the area between the \(x\)-axis and a straight line that runs from \((x_{j-1},f(x_{j-1}))\) to \((x_j, f(x_j))\) — that is, we approximate the function \(f(x)\) on this interval by a linear function that agrees with the function at each endpoint. An obvious way to extend this — just as we did when extending linear approximations to quadratic approximations in our differential calculus course — is to approximate the function with a quadratic. This is precisely what Simpson's  5 Simpson's rule is named after the 18th century English mathematician Thomas Simpson, despite its use a century earlier by the German mathematician and astronomer Johannes Kepler. In many German texts the rule is often called Kepler's rule. rule does.

Simpson's rule approximates the integral over two neighbouring subintervals by the area between a parabola and the \(x\)-axis. In order to describe this parabola we need 3 distinct points (which is why we approximate two subintegrals at a time). That is, we approximate

\begin{align*} \int_{x_0}^{x_1} f(x)\,\dee{x} + \int_{x_1}^{x_2} f(x)\,\dee{x} &= \int_{x_0}^{x_2} f(x)\,\dee{x} \end{align*}

by the area bounded by the parabola that passes through the three points \(\big(x_0,f(x_0)\big)\text{,}\) \(\big(x_1,f(x_1)\big)\) and \(\big(x_2,f(x_2)\big)\text{,}\) the \(x\)-axis and the vertical lines \(x=x_0\) and \(x=x_2\text{.}\)

We repeat this on the next pair of subintervals and approximate \(\int_{x_2}^{x_4} f(x)\,\dee{x}\) by the area between the \(x\)-axis and the part of a parabola with \(x_2\le x\le x_4\text{.}\) This parabola passes through the three points \(\big(x_2,f(x_2)\big)\text{,}\) \(\big(x_3,f(x_3)\big)\) and \(\big(x_4,f(x_4)\big)\text{.}\) And so on. Because Simpson's rule does the approximation two slices at a time, \(n\) must be even.

To derive Simpson's rule formula, we first find the equation of the parabola that passes through the three points \(\big(x_0,f(x_0)\big)\text{,}\) \(\big(x_1,f(x_1)\big)\) and \(\big(x_2,f(x_2)\big)\text{.}\) Then we find the area between the \(x\)-axis and the part of that parabola with \(x_0\le x\le x_2\text{.}\) To simplify this computation consider a parabola passing through the points \((-h,y_{-1}), (0,y_0)\) and \((h,y_1)\text{.}\)

Write the equation of the parabola as

\begin{align*} y &= Ax^2 + Bx +C \end{align*}

Then the area between it and the \(x\)-axis with \(x\) running from \(-h\) to \(h\) is

\begin{align*} \int_{-h}^h \big[ Ax^2 + Bx +C \big] \dee{x} &= \left[ \frac{A}{3}x^3 + \frac{B}{2}x^2 + Cx \right]_{-h}^h\\ &= \frac{2A}{3}h^3 + 2Ch & \text{it is helpful to write it as}\\ &= \frac{h}{3}\left( 2Ah^2 + 6C \right) \end{align*}

Now, the the three points \((-h,y_{-1}), (0,y_0)\) and \((h,y_1)\) lie on this parabola if and only if

\begin{align*} A h^2 - Bh + C &= y_{-1} & \text{at $(-h,y_{-1})$}\\ C &= y_{0} & \text{at $(0,y_{0})$}\\ A h^2 + Bh + C &= y_{1} & \text{at $(h,y_{1})$} \end{align*}

Adding the first and third equations together gives us

\begin{align*} 2Ah^2 + (B-B)h + 2C &= y_{-1}+y_{1} \end{align*}

To this we add four times the middle equation

\begin{align*} 2Ah^2 + 6C &= y_{-1}+4y_0+y_1. \end{align*}

This means that

\begin{align*} \text{area} &= \int_{-h}^h \big[ Ax^2 + Bx +C \big] \dee{x} = \frac{h}{3}\left( 2Ah^2 + 6C \right)\\ &= \frac{h}{3}\left(y_{-1}+4y_0+y_1\right) \end{align*}

Note that here

  • \(h\) is one half of the length of the \(x\)-interval under consideration
  • \(y_{-1}\) is the height of the parabola at the left hand end of the interval under consideration
  • \(y_0\) is the height of the parabola at the middle point of the interval under consideration
  • \(y_{1}\) is the height of the parabola at the right hand end of the interval under consideration

So Simpson's rule approximates

\begin{equation*} \int_{x_0}^{x_2} f(x)\,\dee{x} \approx \tfrac{1}{3}\De x\big[f(x_0)+4f(x_1)+f(x_2)\big] \end{equation*}

and

\begin{equation*} \int_{x_2}^{x_4} f(x)\,\dee{x} \approx \tfrac{1}{3}\De x\big[f(x_2)+4f(x_3)+f(x_4)\big] \end{equation*}

and so on. Summing these all together gives:

\begin{align*} \int_a^b& f(x)\,\dee{x} =\int_{x_0}^{x_2} f(x)\,\dee{x} +\int_{x_2}^{x_4} f(x)\,\dee{x} +\int_{x_4}^{x_6} f(x)\,\dee{x} +\cdots +\int_{x_{n-2}}^{x_n} f(x)\,\dee{x}\\ &\approx\,\tfrac{\De x}{3}\big[f(x_0)+4f(x_1)+f(x_2)\big] +\,\tfrac{\De x}{3}\big[f(x_2)+4f(x_3)+f(x_4)\big] \cr &\ \ \ +\,\tfrac{\De x}{3}\big[f(x_4)+4f(x_5)+f(x_6)\big] +\,\cdots\ +\,\tfrac{\De x}{3}\big[f(x_{n-2})+4f(x_{n-1})+f(x_n)\big]\\ &=\Big[f(x_0)\!+4f(x_1)\!+2f(x_2)\!+4f(x_3)\!+2f(x_4)\! +\cdots+ 2f(x_{n-2})\!+4f(x_{n-1})\!+ f(x_n)\Big]\tfrac{\De x}{3} \end{align*}

In summary

Notice that Simpson's rule requires essentially no more work than the trapezoidal rule. In both rules we must evaluate \(f(x)\) at \(x=x_0,x_1,\cdots,x_n\text{,}\) but we add those terms multiplied by different constants  6 There is an easy generalisation of Simpson's rule that uses cubics instead of parabolas. It is known as Simpson's second rule and Simpson's \(\frac38\) rule. While one can push this approach further (using quartics, quintics etc), it can sometimes lead to larger errors — the interested reader should look up Runge's phenomenon..

Let's put it to work on our two running examples.

Solution: We proceed almost identically to Example 1.11.7 and again use \(n=8\) steps.

  • We have the same \(\De,a,b,x_0,\cdots, x_n\) as Example 1.11.7.
  • Applying Equation 1.11.9 gives
    \begin{align*} &\int_0^1 \frac{4}{1+x^2}\,\dee{x}\\ &\approx \bigg[\frac{4}{1+0^2}+ 4\frac{4}{1+\tfrac{1}{8^2}}+ 2\frac{4}{1+\tfrac{2^2}{8^2}}+ 4\frac{4}{1+\tfrac{3^2}{8^2}}+ 2\frac{4}{1+\tfrac{4^2}{8^2}}\\ &\hskip0.5in +4\frac{4}{1+\tfrac{5^2}{8^2}}+ 2\frac{4}{1+\tfrac{6^2}{8^2}}+ 4\frac{4}{1+\tfrac{7^2}{8^2}}+ \frac{4}{1+\tfrac{8^2}{8^2}}\bigg]\frac{1}{8\times 3}\\ &=\Big[ 4\!+\! 4\times 3.938461538 \!+\! 2\times 3.764705882 \!+\! 4\times 3.506849315 \!+\! 2\times 3.2\\ &\hskip0.5in +4\times 2.876404494 + 2\times 2.56 + 4\times 2.265486726 + 2\Big]\frac{1}{8\times 3}\\ &=3.14159250 \end{align*}
    to eight decimal places.
  • This agrees with \(\pi\) (the exact value of the integral) to six decimal places. So the error in the approximation generated by eight steps of Simpson's rule is \(|3.14159250-\pi|=1.5\times 10^{-7}\text{,}\) which is \(100\tfrac{|3.14159250-\pi|}{\pi}\% =5\times 10^{-6}\%\) of the exact answer.

It is striking that the absolute error approximating with Simpson's rule is so much smaller than the error from the midpoint and trapezoidal rules.

\begin{align*} &\text{midpoint error}\hskip-0.35in &&= 0.0013\\ &\text{trapezoid error}\hskip-0.35in &&= 0.0026\\ &\text{Simpson error}\hskip-0.35in &&= 0.00000015 \end{align*}

Buoyed by this success, we will also redo Example 1.11.8 using Simpson's rule.

Solution: We proceed almost identically to Example 1.11.8 and again use \(n=8\) steps.

  • We have the same \(\De,a,b,x_0,\cdots, x_n\) as Example 1.11.7.
  • Applying Equation 1.11.9 gives
    \begin{align*} &\int_0^\pi\sin x\,\dee{x}\\ &\hskip0.5in\approx \Big[\sin(x_0)+4\sin(x_1)+2\sin(x_2)+\cdots+4\sin(x_7)+\sin(x_8)\Big]\tfrac{\De x}{3}\\ &\hskip0.5in=\Big[\sin(0)+ 4\sin(\tfrac{\pi}{8})+ 2\sin(\tfrac{2\pi}{8})+ 4\sin(\tfrac{3\pi}{8})+ 2\sin(\tfrac{4\pi}{8})\cr &\hskip0.5in\phantom{=\Big[\sin(0)\,} +4\sin(\tfrac{5\pi}{8})+ 2\sin(\tfrac{6\pi}{8})+ 4\sin(\tfrac{7\pi}{8})+ \sin(\tfrac{8\pi}{8})\Big]\tfrac{\pi}{8\times 3}\cr &=\hskip0.5in\Big[ 0+ 4\times 0.382683+ 2\times 0.707107+ 4\times 0.923880+ 2\times 1.0\\ &\hskip0.5in\phantom{=\Big[ 0\,}+ 4\times 0.923880+ 2\times 0.707107+ 4\times 0.382683+0\Big]\tfrac{\pi}{8\times 3}\\ &\hskip0.5in=15.280932\times 0.130900\\ &\hskip0.5in=2.00027 \end{align*}
  • With only eight steps of Simpson's rule we achieved \(100\tfrac{2.00027-2}{2}=0.014\%\) accuracy.

Again we contrast the error we achieved with the other two rules:

\begin{align*} &\text{midpoint error}\hskip-0.35in &&= 0.013\\ &\text{trapezoid error}\hskip-0.35in &&= 0.026\\ &\text{Simpson error}\hskip-0.35in &&= 0.00027 \end{align*}

This completes our derivation of the midpoint, trapezoidal and Simpson's rules for approximating the values of definite integrals. So far we have not attempted to see how efficient and how accurate the algorithms are in general. That's our next task.

login