Python Riemann Summe ergibt keine Null für gleiche positive und negative Flächen

Question

Ich habe ein Programm geschrieben, um ein Integral mit einer Riemann-Summe zu approximieren und es mit Matplotlib in Python graphisch darzustellen. Für Funktionen mit gleichen Bereichen oberhalb und unterhalb der X-Achse sollte der resultierende Bereich Null sein, aber mein Programm gibt stattdessen eine sehr kleine Zahl aus.

Der folgende Code zeigt die ungerade Funktion f (x) = x^3 von -1 nach 1, so dass der Bereich Null sein sollte. Mein Code entspricht stattdessen 1.68065561477562 e^-15.

Was verursacht das? Ist es ein Rundungsfehler in delta_x, x oder y? Ich weiß, ich könnte den Wert einfach auf Null runden, aber ich frage mich, ob es ein anderes Problem oder eine andere Lösung gibt.

Ich habe versucht, mit der Decimal.decimal-Klasse für delta_x, aber ich habe gerade eine noch kleinere Zahl.

Der Python-Code:

import matplotlib.pyplot as plt 
import numpy as np 

# Approximates and graphs integral using Riemann Sum 


# example function: f(x) = x^3 
def f_x(x): 
    return x**3 

# integration range from a to b with n rectangles 
a, b, n = -1, 1, 1000 

# calculate delta x, list of x-values, list of y-values, and approximate area under curve 
delta_x = (b - a)/n 

x = np.arange(a, b+delta_x, delta_x) 

y = [f_x(i) for i in x] 

area = sum(y) * delta_x 

# graph using matplotlib 
fig = plt.figure() 
ax = fig.add_subplot(111) 
ax.plot(x, y) 
ax.bar(x, y, delta_x, alpha=.5) 
plt.title('a={}, b={}, n={}'.format(a, b, n)) 
plt.xlabel('A = {}'.format(area)) 
plt.show() 

Answer 1

Sie müssen beachten, dass das, was Sie berechnen, kein Riemann-Integral im ursprünglichen Sinne ist. Sie teilen das Intervall in n Bins, aber dann summieren Sie über Bins (hier n = 1000 aber len(x) == 1001). Das Ergebnis kann also nahe an dem liegen, was Sie erwarten, aber es ist sicherlich kein guter Weg dorthin.

Mit der Riemann sum teilen Sie Ihr Intervall in n Bins, und summieren Sie dann über die Werte dieser n Bins. Sie haben die Wahl, ob Sie die linke Riemann-Summe, die rechte Riemann-Summe oder eventuell die Mittelpunkte berechnen wollen.

import numpy as np 

def f_x(x): 
    return x**3 

# integration range from a to b with n rectangles 
a, b, n = -1, 1, 1000 

delta_x = (b - a)/float(n) 

x_left = np.arange(a, b, delta_x) 
x_right = np.arange(a+delta_x, b+delta_x, delta_x) 
x_mid = np.arange(a+delta_x/2., b+delta_x/2., delta_x) 

print len(x_left), len(x_right), len(x_mid)   ### 1000 1000 1000 


area_left = f_x(x_left).sum() * delta_x 
area_right = f_x(x_right).sum() * delta_x 
area_mid = f_x(x_mid).sum() * delta_x 

print area_left  # -0.002 
print area_right # 0.002 
print area_mid  # 1.81898940355e-15 

Während der Mittelpunkt Summe bereits ein gutes Ergebnis ergibt, für symmetrische Funktionen wäre es eine gute Idee sein, n selbst zu wählen, und nehmen Sie den Mittelwert der linken und rechten Summe,

print 0.5*(area_right+area_left) # 1.76204537072e-15 

Diese ist gleich nahe 0.

Jetzt ist es erwähnenswert, dass numpy.arange einige Fehler selbst erzeugt.wäre eine bessere Wahl numpy.linspace

x_left = np.linspace(a, b-delta_x, n) 
x_right = np.linspace(a+delta_x, b, n) 
x_mid = np.linspace(a+delta_x/2., b-delta_x/2., n) 

Nachgeben

print area_left  # -0.002 
print area_right # 0.002 
print area_mid  # 8.52651282912e-17 
print 0.5*(area_right+area_left) # 5.68121938382e-17 

5.68121938382e-17 auf 0. Der Grund ziemlich nah dran ist, warum ist es nicht ganz 0 ist in der Tat floating point inaccuracies verwenden.

Das berühmte Beispiel dafür wäre 0.1 + 0.2 - 0.3 , die statt 0 in 5.5e-17 führt sein Dies, dass dies einfach Betrieb die gleichen Fehler in der Größenordnung von 1e-17 als Riemann Integration zeigen einführt.

Answer 2

Ihr Code scheint mir, um gültig zu sein. Es summiert die y, die von f_x zurückgegeben wird, und berücksichtigt den Approximationsfehler von nur 1000 Teilintervallen, indem delta_x zu dem zweiten Argument von arange hinzugefügt wird. Leider denke ich, dass der Rundungsfehler hier eine Rolle spielt.

Answer 3

Ja, dies liegt an der Gleitkomma-Ungenauigkeit.

Zum Erstellen der Partition kann np.linspace() geeigneter sein, da np.arange() möglicherweise den Endpunkt enthält oder nicht, je nachdem, wie es gerundet wird, wenn nicht ganzzahlige Schrittgrößen verwendet werden.

Von numpy.arange() Dokumentation:

Wenn ein nicht-ganzzahligen Schritt verwendet wird, wie zum Beispiel 0,1, werden die Ergebnisse oft nicht konsistent sein. Es ist besser, in diesen Fällen linspace zu verwenden.