[SciPy-user] Simpson Bug?

[Calogero Colletto]

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) )

integration from a to b
n: number of intervals. must be even
h = (b-a)/n ... width of an interval

My Script (simple implementation of the simpson rule):
______________________________________________
from math import pi, sin, cos
from numpy import arange
from scipy import integrate


def simps(y_seq, x_seq):
    h = x_seq[1] - x_seq[0]
    val = y_seq[0] + y_seq[-1]

    for i, y in enumerate(y_seq):
        if (i % 2) == 1:
            val += 4 * y
        elif (i % 2) == 0:
            val += 2 * y
        print "loop ", i, h*val/3.

    return h*val/3.


if __name__ == '__main__':
    s = 11   # Anzahl der Stützstellen: s = n+1
                # Je mehr Stützstellen, umso genauer das Ergebnis
    f = lambda x: sin(x)

    x = arange(0., pi/2., pi/(2.*s))
    y = tuple([f(i) for i in x])

    result = simps(y, x)

    print "My own implementation: ", result
    print "Scipy implementation: ", integrate.simps(y, x) 
_______________________________________________________

I tested the script with the sin-function and integrate this from 0 to pi/2. The
analytical result is 1.

Following results gives me the script:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
loop 0 0.0471153904573
loop 1 0.0742120723007
loop 2 0.101032948993
loop 3 0.180127782511
loop 4 0.231596667976
loop 5 0.356281860151
loop 6 0.428229051385
loop 7 0.588403349479
loop 8 0.675000112936
loop 9 0.85768714791
loop 10 0.951917928825

My own implementation: 0.951917928825
Scipy implementation: 0.85768714791
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Maybe its my fault. But I think that the scipy-implementation is loosing the
last increment.

Calo