[SciPy-user] Simpson Bug?

[Anne Archibald]

2008/6/3 Calogero Colletto <calogero.colletto at gmx.de>:
> Hi,
>
> recently I tested the simpson-rule from the scipy-package and I figured out an
> inconsistency in the result. Following rule:
>
> S[n] = h/3 * ( f(a) + f(b) + 4 * Sum( f(a+i*h), for all odd i, n) + 2 * Sum(
> f(a+i*h), for all even i, n) )

This is not quite the right formula. The problem is (if there are an
odd number of points) that it counts the first and last points three
times, rather than once. If you start with -(f(a)+f(b)) it works out
all right. As a cross-check, try using some functions whose integral
is known to be exact when computed with Simpson's rule: constants,
cubic polynomials. You will see that, for example, your code returns
1.1333 rather than 1 as the integral of 1 from 0 to 1.

>    x = arange(0., pi/2., pi/(2.*s))

This use of arange is highly error-prone. Depending on round-off
error, you could get either s or s+1 points:

In [10]: len(arange(0,pi/2,pi/(2*61)))
Out[10]: 62

When you care how many points you get, use linspace:

x = linspace(0,pi/2,s)

It's clearer, and it's guaranteed to give you the right number of
points. arange is appropriate when you are working with integers, or
when what you care about is the size of the step.

Anne