# 9.3: Simpson's Rule

- Page ID
- 34853

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.