# [SciPy-user] Simpson Bug?

Calogero Colletto calogero.colletto@gmx...
Tue Jun 3 13:49:43 CDT 2008

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

```