Easy Simpson’s rule

Simpson’s rule, as I’m sure all readers of this blog will know, is a simple rule for numerical integration (or quadrature). It can be obtained in two ways, either by determining the interpolating polynomial over three equally spaced points and integrating that polynomial, or by finding the values w_i for which

    \[ \int_a^bf(x),dx = w_0f(a)+w_1f(a+h)+w_2f(b) \]

(where h=(b-a)/2) is accurate for f(x)=1,x,x^2. Over the interval [0,2], Simpson’s rule is

    \[ \int_0^2f(x),dx = \frac{1}{3}\Bigl(f(0)+4f(1)+f(2)\Bigr) \]

and over a general interval

    \[ \int_a^bf(x),dx = \frac{b-a}{6}\Bigl(f(a)+4f(a+h)+f(b)\Bigr). \]

Simpson’s rule has the added advantage that it is also accurate for cubic polynomials.

On its own the rule is not hugely accurate, so one way of increasing the accuracy is to divide the interval into smaller, equally sized regions and apply Simpson’s rule to each of them. This provides the composite Simpson’s rule:

    \begin{align*} \int_a^bf(x),dx &\approx \frac{b-a}{3n}\Bigl(f(a)+4f(a+h)+2f(a+2h)+4f(a+3h)+\cdots\Bigr.\ &\Bigl.+2f(a+(n-2)h)+4f(a+(n-1)h)+f(b)\Bigr) \end{align*}

where n is even and h=(b-a)/n. Sometimes this rule is written

    \begin{align*} \int_a^bf(x),dx \approx \frac{b-a}{6n}\Bigl(f(a)+4f(a+h)+2f(a+2h)+4f(a+3h)+\cdots\Bigr.\ \Bigl.+2f(a+(2n-2)h)+4f(a+(2n-1)h)+f(b)\Bigr) \end{align*}

where there is no restriction on n and h=(b-a)/2n. For a spirited criticism of this rule, both from practical and pedagogical perspectives, see here.

However, in a standard generic mathematics course, Simpson’s rule is as far as most students get with quadrature. The trouble is that the order of the weights (the coefficients of the function values):

1,4,2,4,2,4,\ldots,2,4,1

is tricky to program. Either you have to run two loops, one for the 4’s and the other for the 2’s, or you have to do a little test at each stage to see whether you are at a 2 or 4 point, or you might invoke the modulus function. None of those ways are particularly straightforward, especially if, like many beginning students, you have little programming experience.

So here’s a super easy way. We start with end points a and b and an even integer n. We first create a list of all the values at which the function is to be computed:

xs=[a+kh \mbox{\; for\;} k=0,1,2,\ldots n]

where h=(b-a)/n, and then a list of the weights:

w=[3-(-1)^k \mbox{\; for\;} k=0,1,2,\ldots n].

This last list is [2,4,2,4,\ldots,2,4,2].

We then compute

\sum w_if(xs_i)

However! note that our list of weights starts and finishes with 2’s instead of 1’s. We correct for this by simply subtracting f(a) and f(b) from the sum. Our “easy” Simpson’s rule is then:

    \[ \frac{b-a}{2n}\Bigl(\sum w_if(xs_i)-f(a)-f(b)\Bigr). \]

This is almost trivial to compute. For example, on a TI-nspire CAS calculator, we might have:

a:=0

b:=1

\mbox{Define\;}f(x)=e^{-x^2}

n:=16

h:=(b-a)/n

xs:=\mbox{seq}(a+k\cdot h,k,0,n)

w:=\mbox{seq}(3-(-1)^k,k,0,n)

(\mbox{sum}(w\cdot f(xs))-f(a)-f(b))\cdot (b-a)/n

The TI-nspire has the nice attribute that a function applied to a list will produce a list of all function values, and the product of two lists is a list containing all the products of corresponding elements. The Casio Classpad has the same functionality with lists, and the above commands can be used.

This approach of using 3-(-1)^k to obtaining the weights for Simpson’s rule seems not widely used, but it seems to me to be absurdly simple. I’ve only come across it in one other place.

2 thoughts on “Easy Simpson’s rule

  1. Thanks – in fact the problem is that when I imported all my posts from wordpress.com to this wordpress.org blog, the LaTeX backslashes were all removed – which of course means that most mathematics is not properly typeset. So every now and again I fix up a blog post.

Leave a Reply

Your email address will not be published. Required fields are marked *