Simpson-Implementation in scipy fehlerhaft?

[calo]

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

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

Folgendes Skript:

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)

Als Funktion habe ich Sinus gewählt. Integriert wird von 0 bis pi/2. Analytische Lösung ist 1.

Folgendes Zeitg mein Skript:

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

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?

calo

[CM]

Hm, Deine Implementierung deckt zwar nicht alle Erweiterungen der Regel ab, scheint aber kurz und elegant! Man kann noch zusammenfassen:

Code: Alles auswählen

val = y_seq[0] + y_seq[-1]

Ich finde es auch seltsam, das scipy soweit daneben liegt. Für "echte" Funktionen magst Du vielleicht integrate.quad nehmen:

Code: Alles auswählen

integrate.quad(sin, 0, pi/2)

integrate.simps kommt näher an das richtige Ergebnis ran, wenn s größer wird (ok, das ist trivial). Offenbar gibt es hier ein Samplingproblem.

Magst Du Deine Frage mal auf der Mailingliste stellen? Ich wäre echt neugierig, was diejenigen sagen, die die Routine implementiert haben.

Gruß,
Christian

[calo]

Hi,

Ich habe die Frage mittlerweile in der Scipy.Mailingliste eingestllt (http://article.gmane.org/gmane.comp.pyt ... user/16572) und hab auch gleich eine Antwort erhalten(http://article.gmane.org/gmane.comp.pyt ... user/16574). Es war wohl doch mein Fehler. Zum einen habe ich die Simpson-Regel nicht ganz korrekt eingesetzt (f(a) und f(b) werden doppelt gezählt) und zum anderen das numpy.arange() falsch verwendet. Da haben sich zwei Fehler fast gegenseitig aufgehoben.

Statt des arange() wurde mir außerdem -aufgrund des genaueren Werte und der sichereren Anzahl an Stützpunkten- das linspace() empfohlen:

Code: Alles auswählen

>>> x1 = arange(0., pi/2, (pi/2)/s)
>>> x1
array([ 0.        ,  0.14279967,  0.28559933,  0.428399  ,  0.57119866,
        0.71399833,  0.856798  ,  0.99959766,  1.14239733,  1.28519699,
        1.42799666])
>>> x2 = linspace(0., pi/2, s)
>>> x2
array([ 0.        ,  0.15707963,  0.31415927,  0.4712389 ,  0.62831853,
        0.78539816,  0.9424778 ,  1.09955743,  1.25663706,  1.41371669,
        1.57079633])
>>> x3 =  arange(0., pi/2+.001, (pi/2)/s)
>>> x3
array([ 0.        ,  0.14279967,  0.28559933,  0.428399  ,  0.57119866,
        0.71399833,  0.856798  ,  0.99959766,  1.14239733,  1.28519699,
        1.42799666,  1.57079633])

Mit Hilfe des linspace() erhalte ich dann beim numpy.integrate.simps() auch ein viel genaueres Ergebnis als bei meiner eigenen Implementierung:

loop 0 0.0523598775598
loop 1 0.0851234353024
loop 2 0.117483619281
loop 3 0.2125671672
loop 4 0.274119894883
loop 5 0.422215992822
loop 6 0.50693605436
loop 7 0.693548024406
loop 8 0.793142429898
loop 9 1.00000339222
loop 10 1.10472314734
My own implementation: 1.10472314734
Scipy implementation: 1.00000339222

Gruß, Calo

[CM]

Arrgh, warum hab' ich das nicht gesehen?

Mit Hilfe des linspace() erhalte ich dann beim numpy.integrate.simps() auch ein viel genaueres Ergebnis als bei meiner eigenen Implementierung

Wie man's nimmt: Schau Dir mal die Ausgabe bei Durchlauf Nr. 9 an und vergleiche mit Annes Antwort, erster Teil.

Gruß,
Christian

[calo]

Hi, ja richtig. 1. und letztes Element müssen in der For-Schleife ausgeklammert werden, da sie ja bereits mit

Code: Alles auswählen

val = y_seq[0] + y_seq[-1]

berücksichtigt wurden.

loop 1 0.0851234353024
loop 2 0.117483619281
loop 3 0.2125671672
loop 4 0.274119894883
loop 5 0.422215992822
loop 6 0.50693605436
loop 7 0.693548024406
loop 8 0.793142429898
loop 9 1.00000339222
My own implementation: 1.00000339222
Scipy implementation: 1.00000339222
Diff: -2.22044604925e-016

Einen größeren Rundungsfehler durch das arange kann ich allerdings nicht erkennen. In beiden Fällen beträgt die Differenz zw meiner integration und der numpy-Integration -2.22044604925e-016

Code: Alles auswählen

x = arange(0., pi/2.+0.001, (pi/2)/(s-1))
#x = linspace(0., pi/2., s)

Zugegebenermaßen sieht aber die linspace-Zeile eleganter aus.

Gruß, Calo