Simpson-Implementation in scipy fehlerhaft?
Verfasst: Sonntag 1. Juni 2008, 20:34
Hallo,
nach Anfängerschwierigkeiten numpy/scipy für win64 aufzutreiben (http://www.python-forum.de/topic-14842.html), habe ich mich daran gemacht einige Funktionen selbst zu implementieren. Ich habe dabei mit der Simpson-Regel (http://de.wikipedia.org/wiki/Simpsonregel) begonnen und bin gleich auf eine Ungereimtheit in scipy gestoßen. Laut Simpson-Regel gilt folgende Reihe (Quelle: alte Vorlesungsmanuskripte):
Als Funktion habe ich Sinus gewählt. Integriert wird von 0 bis pi/2. Analytische Lösung ist 1.
Folgendes Zeitg mein Skript:
calo
nach Anfängerschwierigkeiten numpy/scipy für win64 aufzutreiben (http://www.python-forum.de/topic-14842.html), habe ich mich daran gemacht einige Funktionen selbst zu implementieren. Ich habe dabei mit der Simpson-Regel (http://de.wikipedia.org/wiki/Simpsonregel) begonnen und bin gleich auf eine Ungereimtheit in scipy gestoßen. Laut Simpson-Regel gilt folgende Reihe (Quelle: alte Vorlesungsmanuskripte):
Folgendes Skript:S[n] = h/3 * ( f(a) + f(b) + 4 * Summe( f(a+i*h), für alle ungerade i, n) + 2 * Summe( f(a+i*h), für alle gerade i, n) )
n: Anzahl der Intervalle. Muss gerade sein
h = (b-a)/n ... Interval-Breite
Von a bis b wird integriert
Code: Alles auswählen
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 = 0.
val += y_seq[0]
val += 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)
Folgendes Zeitg mein Skript:
Ist schon eigenartig: in der Scipy Implementation der Simpson-Regel wird das letzte Inkrement weggelassen. Mache ich irgendwas falsch, oder stimmt die scipy/simpson-integration nicht?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
calo