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.