9.3: Simpson's Rule

Our analysis of the mid-point rule and the trapezium rule showed that both methods have $O(1/N^2)$ numerical error. In both cases, the error can be traced to the same source: the fact that the integral estimate over each segment differs from the Taylor series result in the second-order term—the term proportional to $f''(x_n)$. This suggests a way to improve on the numerical integration result: we could take a weighted average of the mid-point rule and the trapezium rule, such that the second-order numerical errors from the two schemes cancel each other out! This is the numerical integration method known as Simpson's rule.

To be precise, let's again consider a pair of adjacent segments, which lie between the equally-spaced discretization points $\left\{x_{n-1}, x_n, x_{n+1}\right\}$. As derived above, the integral over these segments can be Taylor expanded as

$$\mathcal{I}_n = \; 2 f(x_n) \Delta x \;+\; \frac{f''(x_n)}{3} \Delta x^3 + O(\Delta x^5) \;+\; \cdots$$

By comparison, the mid-point rule and trapeium rule estimators for the integral are

$$\mathcal{I}_n^{\mathrm{mp}} = 2f(x_n) \Delta x$$

$$\mathcal{I}_n^{\mathrm{trapz}} = 2 f(x_n) \Delta x + \frac{f''(x_n)}{2} \Delta x^3 + O(\Delta x^5).$$

Hence, we could take the following weighted average:

$$\mathcal{I}_n^{\mathrm{simp}} = \frac{1}{3}\mathcal{I}_n^{\mathrm{mp}} + \frac{2}{3}\mathcal{I}_n^{\mathrm{trapz}} = 2 f(x_n) \Delta x + \frac{f''(x_n)}{3} \Delta x^3 + O(\Delta x^5).$$

Such a weighted average would match the Taylor series result up to $O(\Delta x^4)$! (You can check for yourself that the $O(\Delta x^5)$ terms differ.) In summary, Simpson's rule for this set of three points can be written as

$$\begin{align} \mathcal{I}_n^{\mathrm{simp}} &= \frac{1}{3} \Big[2f(x_n) \Delta x\Big] + \frac{2}{3} \, \Delta x \Big[ \frac{f(x_{n-1})}{2} + f(x_n) + \frac{f(x_{n+1})}{2}\Big]\\ &= \frac{\Delta x}{3} \Big[f(x_{n-1}) + 4f(x_n) + f(x_{n+1})\Big]. \end{align}$$

The total numerical error, over a set of $O(N)$ segments, is $O(1/N^4)$. That is an improvement of two powers of $1/N$ over the trapzezium and mid-point rules! What's even better is that it involves exactly the same number of arithmetic operations as the trapezium rule. This is as close to a free lunch as you can get in computational science.

9.3.1 Python Implementation of Simpson's Rule

In Scipy, Simpson's rule is implemented by the scipy.integrate.simps function, which is defined in the scipy.integrate submodule. Similar to the trapz function, this can be called as either simps(y,x) or simps(y,dx=s) to estimate the integral $\int y \;dx$, using the elements of x as the discretization points, with y specifying the set of values for the integrand.

Because Simpson's rule requires dividing the segments into pairs, if you specify an even number of discretization points in x (i.e. an odd number of segments), the function deals with this by doing a trapezium rule estimate on the first and last segments. Usually, the error is negligible, so don't worry about this detail

Here is an example of simps in action:

>>> from scipy import *
>>> from scipy.integrate import simps
>>> x = linspace(0,10,25)
>>> y = exp(-x)
>>> t = simps(y,x)
>>> print(t)
1.00011864276

For the same number of discretization points, the trapezium rule gives $1.01438$; the exact result is $0.9999546\dots$ Clearly, Simpson's rule is more accurate.