11.2    Verfahren der Flächenberechnung

Für viele Funktionen ƒ existieren keine Stammfunktionen F, oder sie lassen sich nur mit einem hohen Aufwand ermitteln. Bei vielen Fragestellungen aus den Ingenieur- und Naturwissenschaften ist nur der Flächeninhalt von Interesse – die Stammfunktion ist nur Mittel zum Zweck. Numerische Verfahren ermöglichen es, den Flächeninhalt näherungsweise und ohne die aufwendige Bestimmung der Stammfunktion zu ermitteln.

Die Grundidee der numerischen Integration besteht darin, die Gesamtfläche in sehr viele Teilflächen aus Rechtecken, Trapezen oder Parabelbögen aufzuteilen. Die einzelnen Teilflächen ƒ(x) ⋅ Δx werden berechnet und dann aufsummiert.

In diesem Abschnitt behandele ich nur die leicht verständlichen Verfahren der Rechteck- und Trapez- sowie die Simpson-Regel. Die Grundformeln sind in Tabelle 11.1 zusammengestellt.

Bezeichnung	Formel
Mittelpunktregel
Trapezregel
Simpson-Regel

Tabelle 11.1     Formeln für die numerische Integration

Mit diesen Formeln können Sie nur für kleine Intervalle [a,b] die Fläche unter Kurven näherungsweise berechnen. Für größere Intervalle und bessere Näherungen müssen Sie die summierten Regeln anwenden.

11.2.1    Die summierte Mittelpunktregel

Abbildung 11.4 zeigt eine Normalparabel, deren Fläche unterhalb des Kurvenverlaufs in den Grenzen von 0 bis 10 durch Rechteckflächen approximiert wird. In dem Intervall [a,b] von a = 0 bis b = 10 wird die x-Achse in zehn gleiche Abschnitte Δx = h aufgeteilt. Diese Aufteilung wird auch als Zerlegung bezeichnet [Walter: 197]. Die Schrittweite, d. h. der Abschnitt auf der x-Achse, beträgt also h = 1.

Abbildung 11.4     Rechtecksummen für eine Normalparabel

Zerlegung

Die Zerlegung heißt äquidistant, wenn alle Teilintervalle auf der x-Achse dieselbe Länge h = (b – a)/n haben.

Die Schrittweite h wird aus den Integrationsgrenzen und der Anzahl der Teilintervalle n ∈ [a,b] für b > a berechnet:

Je kleiner die Schrittweite h gewählt wird, desto genauer wird die Fläche berechnet. Eine Verkleinerung der Schrittweite verlangt aber auch eine entsprechende Erhöhung der Rechenschritte. Die Schrittweite h repräsentiert die Breite einer Rechteckfläche. Die Höhen der Rechtecke werden durch die Funktionswerte der Mitte zwischen den Abständen x_i und x_i+1, also durch

bestimmt.

Die zu großen Flächenanteile, die oberhalb der Parabel liegen, werden durch die zu kleinen Flächenanteile, die unterhalb der Parabel liegen, näherungsweise wieder ausgeglichen.

Für a = 0 bis b = 10 können die zehn Teilsummen der Rechteckflächen gebildet werden:

Für den Flächeninhalt eines Teilintervalls [x_i,x_i+1] gilt allgemein:

Durch Aufsummieren der Teilrechteckflächen erhält man den gesamten Flächeninhalt:

Diese Berechnungsvorschrift wird als summierte Mittelpunktregel bezeichnet, weil die Stützstellen in der Mitte zwischen x_i und x_i+1 liegen.

Setzt man x_i = a + i ⋅ h als neues Argument ein, erhält man eine Formel, die sich leicht in einen Algorithmus übersetzen lässt:

In Pseudocode übersetzt erhalten wir dann:

h=(b-a)/n
wiederhole von i=0 bis n-1
      summe = summe + f(a + i*h + h/2)
A = summe*h

Listing 11.3 berechnet den Flächeninhalt für stetige Funktionen aus n Rechtecken.

01  #03_rechtecksummen.py
02  from math import *
03  #Funktionsdefinition
04  def f(x):
05      return x**2 #x**3,x**4,sin(x),exp(x)
06  #Stammfunktion
07  def F(x):
08      return x**3/3 #x**4/4,x**5/5,-cos(x),exp(x)
09  #Rechtecksummen
10  def rechteck(f,a,b,n=10):
11      h=(b-a)/n
12      summe=0
13      for i in range(0,n):
14          summe=summe+f(a+i*h+h/2)
15      return h*summe
16  #Grenzen
17  a,b=0,10
18  n=10
19  Aa=rechteck(f,a,b,n) #approximiert
20  Ag=F(b)-F(a) #genau
21  #Ausgabe
22  print("Grenzen: a =",a,"\t b =",b)
23  print("Anzahl der Zerlegungen n =",n)
24  print("Schrittweite h =",(b-a)/n)
25  print("Aa =",round(Aa,6), " approximiert")
26  print("Ag =",round(Ag,6), " genau")
27  print("Fehler E =",round(fabs(Ag-Aa),6))

Listing 11.3     Integration mit der Mittelpunktregel

Ausgabe

Grenzen: a = 0    b = 10
Anzahl der Zerlegungen n = 10
Schrittweite h = 1.0
Aa = 332.5  approximiert
Ag = 333.333333  genau
Fehler E = 0.833333

Analyse

In Zeile 05 können Sie die auskommentierten mathematischen Funktionen testen. Dann müssen Sie in Zeile 08 die zugehörige Stammfunktion auswählen.

In Zeile 14 wird das Funktionsargument um die halbe Schrittweite auf der x-Achse nach rechts verschoben.

In Zeile 17 können Sie die Integrationsgrenzen ändern. Durch Erhöhung von n in Zeile 18 können Sie die Genauigkeit der Flächenberechnung beeinflussen. Mit steigendem n vergrößern Sie die Anzahl der Stützstellen und verkleinern damit die Schrittweite. Versuchen Sie herauszufinden, ab welchem Wert von n eine Genauigkeit von vier Nachkommastellen erreicht wird.

11.2.2    Die summierte Trapezregel

Die Flächenberechnung kann auch mit Sekanten-Trapezsummen durchgeführt werden. Abbildung 11.5 zeigt, wie der Flächeninhalt unter einer Normalparabel mit zehn Sekanten-Trapezflächen approximiert werden kann.

Abbildung 11.5     Trapezsummen

Innerhalb der Grenzen von a = 0 bis b = 10 können Sie die Fläche unterhalb der Normalparabel durch Trapeze nun mit h = (b – a)/n annähern:

Der Faktor 1/2 kann ausgeklammert werden:

Die mittleren Funktionswerte von f(x₁) bis f(x_n-1) kommen doppelt vor:

Die Verallgemeinerung ergibt die summierte Trapezformel:

Mit x_i = a + i ⋅ h erhält man:

Diese Berechnungsvorschrift lässt sich leicht in einen Algorithmus übersetzen. Die Summe ∑f(a + i ⋅ h) wird innerhalb einer Schleife berechnet. Danach wird zu dieser Summe der Mittelwert aus der unteren und oberen Integrationsgrenze addiert. Zum Schluss wird dieses Zwischenergebnis mit der Schrittweite h multipliziert. Der Algorithmus sieht so aus:

h=(b-a)/n
wiederhole von i=1 bis n-1
        summe = summe + f(a+i*h)
A = (summe + (f(a)+f(b))/2)*h

Mit Listing 11.4 können Sie die Trapezsummen berechnen.

01  #04_trapez.py
02  from math import *
03  #Funktionsdefinition
04  def f(x):
05      return x**2 #x**3,x**4,sin(x),exp(x)
06  #Stammfunktion
07  def F(x):
08      return x**3/3 #x**4/4,x**5/5,-cos(x),exp(x)
09  #Trapezsummen
10  def trapez(f,a,b,n=10):
11      h=(b-a)/n
12      summe=0
13      for i in range(1,n):
14          summe=summe+f(a+i*h)
15      return (summe + (f(a)+f(b))/2)*h
16  #Grenzen
17  a,b=0,10
18  n=10
19  Aa=trapez(f,a,b,n)
20  Ag=F(b)-F(a)
21  #Ausgabe
22  print("Grenzen: a =",a,"\t b =",b)
23  print("Anzahl der Zerlegungen n =",n)
24  print("Schrittweite h =",(b-a)/n)
25  print("Aa =",round(Aa,6), " approximiert")
26  print("Ag =",round(Ag,6), " genau")
27  print("Fehler E =",round(fabs(Ag-Aa),6))

Listing 11.4     Integration mit der Trapezregel

Ausgabe

Grenzen: a = 0    b = 10
Anzahl der Zerlegungen n = 10
Schrittweite h = 1.0
Aa = 335.0  approximiert
Ag = 333.333333  genau
Fehler E = 1.666667

Analyse

Das Programm ist im Prinzip so aufgebaut wie bei der Berechnung der Rechtecksummen. In Zeile 14 berechnet der Summenalgorithmus die Funktionswerte (Stützstellen) für die inneren Trapezflächen. In Zeile 15 wird der Mittelwert (f(a)+f(b))/2) des linken und des rechten Randes zur Summe summe addiert und mit der Schrittweite h multipliziert.

Versuchen Sie wieder herauszufinden, ab welcher Anzahl n der Schleifendurchläufe eine Genauigkeit von vier Nachkommastellen erreicht wird (Zeile 18).

Vergleich des Mittelpunktverfahrens mit dem Trapezverfahren

Um den Einfluss der Schrittweite h bzw. der Anzahl der Approximationsschritte n auf die Genauigkeit der numerischen Integration näher quantifizieren zu können, sollen das Mittelpunkt- und das Trapezverfahren miteinander verglichen werden. Die Messwerte (Tabelle 11.2) gelten für die numerische Integration der Parabel y = x² in den Grenzen von 0 bis 10.

n	h	Rechtecke	Fehler	Trapeze	Fehler
100	0,1	333,325000	0,008333	333,350000	0,016667
200	0,05	333,331250	0,002083	333,337500	0,004167
400	0,025	333,332813	0,000521	333,334375	0,001042
800	0,0125	333,333203	0,000130	333,333594	0,000260
1600	0,00625	333,333301	0,000003	333,333398	0,000007

Tabelle 11.2     Vergleich zwischen Mittelpunkt- und Trapezverfahren

Wenn man die Anzahl der Rechenschritte verdoppelt, scheint sich die Genauigkeit um eine signifikante Nachkommastelle zu verbessern. Für n = 400 gilt diese Regel aber nicht.

Der Fehler wird mit jedem neuen Rechenschritt geviertelt. Anders ausgedrückt: Wird die Schrittweite h halbiert, verringert sich der Fehler um den Faktor 1/4.

Dem Augenschein nach würde man durch Vergleich von Abbildung 11.4 und Abbildung 11.5 vermuten, dass die Trapezregel besser approximiert als die Mittelpunktregel. Festzuhalten bleibt, dass beide Verfahren für die Praxis nicht geeignet sind. Sie sind aber nützlich, um das Prinzip der numerischen Integration zu veranschaulichen und besser zu verstehen.

11.2.3    Das Simpson-Verfahren

Beim Simpson-Verfahren wird der zu berechnende Flächeninhalt durch Teilflächen unter Parabelbögen angenähert. Die Parabelbögen werden durch drei Punkte gelegt.

Abbildung 11.6     Approximation durch Parabelbögen (Quelle: Wikipedia)

Das Intervall [a,b] setzt sich aus vier Rechtecken und zwei Trapezen zusammen. Es enthält also insgesamt sechs Elementarflächen. Man kann also die gesamte Fläche durch

annähern.

Mit der Mittelpunktregel

und der Trapezregel

erhält man die Simpson-Regel:

Für die Normalparabel y = x² berechnet diese Formel den exakten Wert des Flächeninhalts. Die Beispielrechnung gilt für a = 0 und b = 10:

Die Simpson-Formel berechnet den Flächeninhalt für Polynome bis zum 3. Grad genau. Sollen die bestimmten Integrale für Polynome höheren Grades oder andere Funktionen berechnet werden, muss aus der Simpson-Formel die summierte Simpson-Formel hergeleitet werden. Dazu wird das Intervall [a,b] in n-Teilintervalle zerlegt. Die einzelnen Teilintervalle können mit

berechnet werden.

Durch Aufsummieren erhält man mit h = (b – a)/n, x_i = a + i ⋅ h und x_i+1 = a + (i + 1) ⋅ h die summierte Simpson-Formel für das Intervall [a,b]:

Aus der Summenformel können Sie direkt den Algorithmus ablesen:

h=(b-a)/n
wiederhole von i=0 bis n-1
        summe = summe + f(a+i*h) + 4*f(a+h/2+i*h) + f(a+(i+1)*h)
 A = h*summe/6

Listing 11.5 berechnet ein bestimmtes Integral mit der Simpson-Regel.

01  #05_simpson.py
02  from math import *
03  #Funktionsdefinition
04  def f(x):
05      return x**2 #x**3,x**4,sin(x),exp(x)
06  #Stammfunktion
07  def F(x):
08      return x**3/3 #x**4/4,x**5/5,-cos(x),exp(x)
09  #Simpson-Regel
10  def simpson(f,a,b,n=10):
11      h=(b-a)/n
12      summe=0
13      for i in range(0,n):
14          summe = summe + f(a+i*h) + 4*f(a+h/2+i*h) + f(a+(i+1)*h)
15      return h*summe/6
16  #Grenzen
17  a,b=0,10
18  n=10
19  Aa=simpson(f,a,b,n) #approximiert
20  Ag=F(b)-F(a) #genau
21  #Ausgabe
22  print("Grenzen: a =",a,"\t b =",b)
23  print("Anzahl der Zerlegungen n =",n)
24  print("Schrittweite h =",(b-a)/n)
25  print("Aa =",Aa," approximiert")
26  print("Ag =",Ag, " genau")
27  print("Fehler E =",fabs(Ag-Aa))

Listing 11.5     Simpson-Regel

Ausgabe

Grenzen: a = 0    b = 10
Anzahl der Zerlegungen n = 10
Schrittweite h = 1.0
Aa = 333.3333333333333  approximiert
Ag = 333.3333333333333  genau
Fehler E = 0.0

Analyse

Der Quelltext (Zeilen 10 bis 15) stimmt formal mit dem Pseudocode überein. Innerhalb der Schleife werden drei Funktionswerte so oft aufaddiert, bis die obere Integrationsgrenze erreicht wird.

Wenn Sie in Zeile 18 der Variablen n den Wert 1 zuweisen, wird schon nach einem Schleifendurchlauf der exakte Wert des Flächeninhalts berechnet.

Testen Sie das Programm für Parabeln vom Grad 3 bis 5, und beobachten Sie den Fehler. Wie groß ist der Fehler bei einer Parabel 4. Grades (n = 10, a = 0, b = 10). Wenn Sie andere Funktionen testen, müssen Sie in Zeile 08 die Stammfunktionen anpassen.

11.2.4    Vergleich der Verfahren

Um die bisher vorgestellten Verfahren zu vergleichen, sollen für die folgenden bestimmten Integrale

die Flächeninhalte mit Rechteck- und Trapezsummen sowie dem Simpson-Verfahren berechnet werden. Die in Tabelle 11.3 dargestellten Ergebnisse wurden mit Listing 11.6 berechnet.

	1	2	3	4
Rechtecksummen	0,333125	0,199583	1,000257	0,693069
Trapezsummen	0,333750	0,200833	0,999485	0,693303
Simpson	0,333333	0,200000	1,000000	0,693147
quad	0,333333	0,2	0,999999	0,693147
genau	1/3	1/5	1	ln(2)

Tabelle 11.3     Vergleich der Verfahren für n = 20

Der Vergleich zeigt, dass das Simpson-Verfahren für die hier betrachteten Beispiele die genaueren Ergebnisse liefert. Die Ergebnisse der Rechteck- und Trapezsummen sind noch akzeptabel. Wenn aber die Integrationsgrenzen weiter auseinanderliegen, dann fallen die Flächenberechnungen, je nach Funktionsmerkmal (konkav, konvex, Steigung), mitunter wesentlich ungenauer aus.

Aus dieser kleinen vergleichenden Untersuchung können wir folgendes Fazit ziehen:

»Grundsätzlich gibt es keine in allen Fällen günstige Allround-Formel. Optimalität kann immer nur in Bezug auf spezielle Merkmale erzielt werden.« [Walter: 298]

Für die folgenden Programme werde ich daher den besonders einfach zu handhabenden Algorithmus der Rechtecksummen verwenden.

Wenn Sie das nächste Programm testen wollen, dann müssen Sie die Funktionen rechteck() aus Listing 11.3, trapez() aus Listing 11.4 und simpson() aus Listing 11.5 in eine Datei integrieren.py speichern und dieses Modul mit der Anweisung from integrieren import * importieren (Zeile 03).

Mit Listing 11.6 werden die bestimmten Integrale Nr. 1, 2, 3 und 4 mit den Algorithmen aus Tabelle 11.3 berechnet.

01  #06_vergleich.py
02  from math import *
03  from integrieren import *
04  from scipy.integrate import quad
05  #Funktionsdefinition
06  def f(x):
07      #y=x**2
08      #y=x**4
09      #y=sin(x)
10      y=1/x
11      return y
12  #Hauptprogramm
13  a,b=1,2 #Grenzen
14  n=20
15  A1=rechteck(f,a,b,n)
16  A2=trapez(f,a,b,n)
17  A3=simpson(f,a,b,n)
18  A4=quad(f,a,b) [0]
19  #Ausgaben
20  print("Grenzen: a =",a,"\t b =",b)
21  print("Anzahl der Zerlegungen n =",n)
22  print(A1,"Rechtecksummen")
23  print(A2,"Trapezsummen")
24  print(A3,"Simpson")
25  print(A4,"quad")

Listing 11.6     Vergleich der Verfahren

Ausgabe

Grenzen: a = 1    b = 2
Anzahl der Zerlegungen n = 20
0.693069098225587 Rechtecksummen
0.6933033817926941 Trapezsummen
0.6931471927479561 Simpson
0.6931471805599454 quad

Analyse

Der natürliche Logarithmus (Logarithmus naturalis) ln(2) hat einen näherungsweise genauen Wert von 0,6931471805599453… Die SciPy-Funktion quad() erreicht eine Genauigkeit von 15 Nachkommastellen. Das Simpson-Verfahren berechnet den ln(2) bis auf fünf Nachkommastellen genau. Die anderen beiden Verfahren (Rechteck- und Trapezsummen) erreichen nur eine Genauigkeit von drei Nachkommastellen.

Die Ergebnisse der übrigen Testfunktionen aus den auskommentierten Zeilen 07 bis 09 werden in Tabelle 11.3 zusammengestellt.

Fehlerabschätzung für summierte Regeln

Mit den Formeln

können die Fehler für die Mittelpunktregel (Rechtecksummen), die Trapezregel und die Simpson-Regel vorausgesagt werden [Knorrenschild: 127]. Der zu erwartende Fehler steht auf der rechten Seite der Ungleichung.

Am Beispiel der Normalparabel y = x² mit der 2. Ableitung y'' = 2 sollen für die Integrationsgrenzen von 0 bis 10 und mit einer Schrittweite von h = 1 die oberen Schranken des zu erwartenden Fehlers berechnet werden.

Für die summierte Mittelpunktregel erhält man den maximalen Fehler:

Die summierte Trapezregel erzeugt den maximalen Fehler:

Die errechneten Werte bestätigen die Ausgaben aus Listing 11.3 und Listing 11.4.

Der Fehler für eine Normalparabel kann mit der Fehlerformel für die summierte Simpson-Regel nicht berechnet werden, weil die 4. Ableitung für diese Funktion nicht existiert. Die Simpson-Regel berechnet den Flächeninhalt unter Parabeln bis zur 3. Ordnung fehlerfrei. Für die Parabel 4. Ordnung y = x⁴ erhält man mit y⁽⁴⁾ = 24 den maximalen Fehler (h = 1):

Aus den Formeln für die Fehlerabschätzung kann für eine vorgegebene Genauigkeit die entsprechende Schrittweite h näherungsweise berechnet werden. Am Beispiel der Mittelpunktregel soll für einen maximalen Fehler von E ⩽ 10^-6 für a = 0, b = 10 und y = x² mit y'' = 2 die Berechnung durchgeführt werden:

Nach h aufgelöst erhält man:

Gewählt wird h = 1 ⋅ 10^-3. Für diese Schrittweite sind

Schleifendurchläufe notwendig. Dieser Wert ist für die Praxis nicht akzeptabel. Wenn Sie diesen Wert in Listing 11.3 eingeben, gibt das Programm den Fehler E=1e-06 aus.

Die Simpson-Formel berechnet schon nach einem Schleifendurchlauf den genauen Wert für die Normalparabel.

Das Problem der Fehlerabschätzung

Das Problem bei der Fehlerabschätzung besteht darin, dass die Ableitungen der Funktion f (Integrand) nicht immer bekannt sind. Wenn man nun vorsichtshalber eine zu kleine Schrittweite wählt, riskiert man zu große Rundungsfehler. Bei jedem Schleifendurchlauf wird in allen Algorithmen der numerischen Integration jeweils eine Addition und eine Multiplikation x=a+i*h durchgeführt. Wird die Schrittweite h zu klein gewählt, dann summieren sich die Rundungsfehler wegen des anwachsenden n, wird sie zu groß gewählt, dann wird die Berechnung zu ungenau.

Eine Lösung dieses Problems besteht darin, die Schrittweite an den Funktionsverlauf anzupassen. Bei kleinen Steigungen wählt der Algorithmus eine größere Schrittweite, und bei großen Steigungen verkleinert er die Schrittweite. Solche Verfahren nennt man adaptive Quadratur.

Mit der SciPy-Funktion quad(f,a,b)

from scipy.integrate import quad
print(quad(f,0,10))
(333.33333333333326, 3.700743415417188e-12) #Ausgabe

erhalten Sie wesentlich genauere Ergebnisse. Das Ergebnis wird als Tupel ausgegeben. Das erste Element des Tupels enthält den Wert des Flächeninhalts. Das zweite Element macht eine Aussage über die Fehlerabschätzung. Das Ergebnis hat eine Genauigkeit von 12 Nachkommastellen (Exponent: –12).