Simpson’s Rule¶

Consider two consecutive subintervals, \([x_{i-1}, x_i]\) and \([x_i, x_{i+1}]\). Simpson’s Rule approximates the area under \(f(x)\) over these two subintervals by fitting a quadratic polynomial through the points \((x_{i-1}, f(x_{i-1})), (x_i, f(x_i))\), and \((x_{i+1}, f(x_{i+1}))\), which is a unique polynomial, and then integrating the quadratic exactly. The following shows this integral approximation for an arbitrary function.

First, we construct the quadratic polynomial approximation of the function over the two subintervals. The easiest way to do this is to use Lagrange polynomials, which was discussed in Chapter 17. By applying the formula for constructing Lagrange polynomials we get the polynomial

\[\begin{eqnarray*}P_i(x) &=& f(x_{i-1})\frac{(x - x_i)(x - x_{i+1})}{(x_{i-1} - x_i)(x_{i-1} - x_{i+1})} + f(x_i)\frac{(x - x_{i-1})(x - x_{i+1})}{(x_{i} - x_{i-1})(x_{i} - x_{i+1})}\\ &&+ f(x_{i+1})\frac{(x - x_{i-1})(x - x_{i})}{(x_{i+1} - x_{i-1})(x_{i+1} - x_{i})},\end{eqnarray*}\]

and with substitutions for \(h\) results in

\[P_i(x) = \frac{f(x_{i-1})}{2h^2} (x - x_i)(x - x_{i+1}) - \frac{f(x_i)}{h^2} (x - x_{i-1})(x - x_{i+1}) + \frac{f(x_{i+1})}{2h^2} (x - x_{i-1})(x - x_{i}).\]

You can confirm that the polynomial intersects the desired points. With some algebra and manipulation, the integral of \(P_i(x)\) over the two subintervals is

\[\int_{x_{i-1}}^{x_{i+1}} P_i(x) dx = \frac{h}{3}(f(x_{i-1}) + 4f(x_i) + f(x_{i+1}).\]

To approximate the integral over \((a, b)\), we must sum the integrals of \(P_i(x)\) over every two subintervals since \(P_i(x)\) spans two subintervals. Substituting \(\frac{h}{3}(f(x_{i-1}) + 4f(x_i) + f(x_{i+1}))\) for the integral of \(P_i(x)\) and regrouping the terms for efficiency leads to the formula

\[\int_a^b f(x) dx \approx \frac{h}{3} \left[f(x_0)+4 \left(\sum_{i=1, i\ {\text{odd}}}^{n-1}f(x_i)\right)+2 \left(\sum_{i=2, i\ {\text{even}}}^{n-2}f(x_i)\right)+f(x_n)\right].\]

This regrouping is illustrated in the figure below:

WARNING! Note that to use Simpson’s Rule, you must have an even number of intervals and, therefore, an odd number of grid points.

To compute the accuracy of the Simpson’s Rule, we take the Taylor series approximation of \(f(x)\) around \(x_i\), which is

\[f(x) = f(x_i) + f^{\prime}(x_i)(x - x_i) + \frac{f''(x_i)(x-x_i)^2}{2!} + \frac{f'''(x_i)(x-x_i)^3}{3!} + \frac{f''''(x_i)(x-x_i)^4}{4!} + \cdots\]

Computing the Taylor series at \(x_{i-1}\) and \(x_{i+1}\) and substituting for \(h\) where appropriate gives the expressions

\[f(x_{i-1}) = f(x_i) - hf^{\prime}(x_i) + \frac{h^2f''(x_i)}{2!} - \frac{h^3f'''(x_i)}{3!} + \frac{h^4f''''(x_i)}{4!} - \cdots\]

and

\[f(x_{i+1}) = f(x_i) + hf^{\prime}(x_i) + \frac{h''(x_i)}{2!} + \frac{h^3f'''(x_i)}{3!} + \frac{h^4f''''(x_i)}{4!} + \cdots\]

Now consider the expression \(\frac{f(x_{i-1}) + 4f(x_i) + f(x_{i+1})}{6}\). Substituting the Taylor series for the respective numerator values produces the equation

\[\frac{f(x_{i-1}) + 4f(x_i) + f(x_{i+1})}{6} = f(x_i) + \frac{h^2}{6}f''(x_i) + \frac{h^4}{72}f''''(x_i) + \cdots\]

Note that the odd terms cancel out. This implies

\[f(x_i) =\frac{f(x_{i-1}) + 4f(x_i) + f(x_{i+1})}{6} - \frac{h^2}{6}f''(x_i) + O(h^4).\]

By substitution of the Taylor series for \(f(x)\), the integral of \(f(x)\) over two subintervals is then

\[\begin{eqnarray*}\int_{x_{i-1}}^{x_{i+1}} f(x) dx &=& \int_{x_{i-1}}^{x_{i+1}} \left(f(x_i) + f^{\prime}(x_i)(x - x_i) + \frac{f''(x_i)(x-x_i)^2}{2!}\right.\\ &&\qquad\qquad\left. + \frac{f'''(x_i)(x-x_i)^3}{3!}+ \frac{f''''(x_i)(x-x_i)^4}{4!} + \cdots\right) dx.\end{eqnarray*}\]

Again, we distribute the integral and without showing it, we drop the integrals of terms with odd powers because they are zero.

\[\int_{x_{i-1}}^{x_{i+1}} f(x) dx = \int_{x_{i-1}}^{x_{i+1}} f(x_i) dx + \int_{x_{i-1}}^{x_{i+1}}\frac{f''(x_i)(x-x_i)^2}{2!}dx + \int_{x_{i-1}}^{x_{i+1}}\frac{f''''(x_i)(x-x_i)^4}{4!}dx + \cdots.\]

We then carry out the integrations. As will soon be clear, it benefits us to compute the integral of the second term exactly. The resulting equation is

\[\int_{x_{i-1}}^{x_{i+1}} f(x) dx = 2h f(x_i) + \frac{h^3}{3}f''(x_i) + O(h^5).\]

Substituting the expression for \(f(x_i)\) derived earlier, the right-hand side becomes

\[2h \left(\frac{f(x_{i-1}) + 4f(x_i) + f(x_{i+1})}{6} - \frac{h^2}{6}f''(x_i) + O(h^4)\right) + \frac{h^3}{3}f''(x_i) + O(h^5),\]

which can be rearranged to

\[\left[\frac{h}{3}(f(x_{i-1}) + 4f(x_i) + f(x_{i+1})) - \frac{h^3}{3}f''(x_i) + O(h^5)\right] + \frac{h^3}{3}f''(x_i) + O(h^5).\]

Canceling and combining the appropriate terms results in the integral expression

\[\int_{x_{i-1}}^{x_{i+1}} f(x) dx = \frac{h}{3}(f(x_{i-1}) + 4f(x_i) + f(x_{i+1})) + O(h^5).\]

Recognizing that \(\frac{h}{3}(f(x_{i-1}) + 4f(x_i) + f(x_{i+1}))\) is exactly the Simpson’s Rule approximation for the integral over this subinterval, this equation implies that Simpson’s Rule is \(O(h^5)\) over a subinterval and \(O(h^4)\) over the whole interval. Because the \(h^3\) terms cancel out exactly, Simpson’s Rule gains another two orders of accuracy!

TRY IT! Use Simpson’s Rule to approximate \(\int_{0}^{\pi} \text{sin} (x)dx\) with 11 evenly spaced grid points over the whole interval. Compare this value to the exact value of 2.

import numpy as np

a = 0
b = np.pi
n = 11
h = (b - a) / (n - 1)
x = np.linspace(a, b, n)
f = np.sin(x)

I_simp = (h/3) * (f[0] + 2*sum(f[:n-2:2]) \
            + 4*sum(f[1:n-1:2]) + f[n-1])
err_simp = 2 - I_simp

print(I_simp)
print(err_simp)

2.0001095173150043
-0.00010951731500430384

< 21.3 Trapezoid Rule | Contents | 21.5 Computing Integrals in Python >

Python Numerical Methods

Simpson’s Rule¶