Numerische Quadraturverfahren, Integrationsverfahren

QuadraturverfahrenStreng analytische Verfahren zur Integration ermöglichen nicht immer die Lösung eines bestimmten Integrals. In der Praxis erweist sich sogar, dass die Mehrzahl aller Integrationsaufgaben nicht analytisch gelöst werden kann. In solchen Fällen nutzt man numerische Näherungsverfahren, von denen in diesem Unterprogramm von 27 Verfahren und Gleichungen jeweils drei parallel verwendet werden können.

Dazu geben Sie zuerst unter Funktion y = die gewünschte Funktion und in den Eingabezeilen Intervall x0 und x1 die Intervallgrenzen an. Für die Intervallgrenzen können Sie auch Terme der Form PI oder SIN(1) … verwenden.

Über das grafische Vorschaubild können Sie sich schnell über den Verlauf der Funktion informieren. Dieses Vorschaubild wird bei der Eingabe der Funktionsgleichung automatisch aktualisiert.
Folgende Verfahren sind nutzbar:

Rechteck-Regel
Rechteck-RegelDie Fläche unter einer Funktion lässt sich näherungsweise berechnen, indem man das Gesamtintervall [a,b] in n gleich große Teilintervalle [a(i), b(i)] teilt und in diese Rechtecke einbeschreibt, deren Seiten zum einen von der Teilintervallbreite

    \[ \Delta x=b(i)-a(i) \]

und zum anderen vom Funktionswertminimum bzw. Funktionswertmaximum des Teilintervalls gebildet werden. Zweckmäßig ist dabei eine konstante Teilintervallbreite Δx, die allerdings nicht unbedingt notwendig ist. Die Summe aller „kleinen“ Rechtecke, also die Untersumme, ist kleiner als der gesuchte Flächeninhalt, die Obersumme ist größer.

Ist die Funktion im Riemannschen Sinne integrierbar, streben Untersumme U(n) und Obersumme O(n) für eine wachsende Zahl von Teilintervallen gegen die gesuchte Fläche unter der Funktion.

Im Programm können Sie sowohl die Unter- als auch die Obersumme berechnen lassen. Für die Untersumme gilt (bei gleicher Intervallbreite und monoton steigender Funktion):

    \[ U(n)=\Delta x(f(a)+f(a+\Delta x)+f(a+2\Delta x)+...+f(b-\Delta x)) \]

Mittelpunktsregel
Während bei der Rechteckregel zur Ermittlung der Höhe des eingeschriebenen Rechtecks, je nach Wahl, der Funktionswert der linken bzw. rechten Intervallgrenze genutzt wird, liefert bei der Mittelpunktsregel das arithmetische Mittel der Intervallgrenzen (der Mittelwert) die Stützstelle. Bei einer Vielzahl von Funktionen wird dadurch ein besserer Näherungswert erreicht. Die Güte der Konvergenz wird allerdings nicht verbessert.

Trapezregel / Tangententrapezregel
TrapezregelOffensichtlich konvergieren die Ober- und Untersumme der Rechteckregel nur langsam gegen den wahren Wert der Fläche. Deshalb ersetzt man das Rechteck in jedem Teilintervall durch ein Trapez, d.h., es werden jeweils zwei Punkte des Graphen durch eine Sekante verbunden.

Leitet man die entsprechende Formel (gleiche Intervallbreite Δx) her, ergibt sich

    \[ T(n)=\Delta x(\frac{f(a)}{2}+f(a+\Delta x)+f(a+2\Delta x)+...+\frac{f(b)}{2}) \]

Der Unterschied zur Rechteckregel besteht folglich nur darin, dass in ihr Δx f(b) / 2 mehr und Δx f(a) / 2 weniger addiert werden. Die numerische Auswirkung auf die Genauigkeit ist jedoch beträchtlich.

Während bei der genannten Trapezregel Sekanten das Trapez bilden, benutzt man bei der Tangententrapezregel Tangenten an den jeweiligen Kurvenstücken. Je nach Integrand erreichen Sie so eine bessere oder schlechtere Konvergenz.
Unter der Poncelet-Formel versteht man das arithmetische Mittel aus Trapezregel und Tangententrapez-Regel.

Simpson-Regel
Setzt man den Gedanken fort, die Teilflächen immer besser der Funktion anzupassen, so erhält man die von Newton und Cotes verfolgte Idee:

Die Geradenstücke durch zwei Punkte werden durch Parabelstücke durch drei Punkte oder durch Polynome 3. Grades durch vier Punkte usw. ersetzt. Insbesondere die Verwendung von Parabelstücken erhöht die Konvergenz erheblich und ergibt die sogenannte Simpson-Regel. Die Entwicklung der Teilintervalle führt zur Gleichung:

    \[ S(n)=\frac{\Delta x}{3}(f(a)+f(b)+4f(a+\Delta x)+2f(a+2\Delta x)+...) \]

Betrachtet man nur zwei Teilintervalle, so erhalten Sie die Keplersche Fassregel:

    \[ A=\frac{b-a}{6}(f(a)+f(b)+4f(\frac{a+b}{2})) \]

wobei a und b die begrenzenden Argumente sind.

Wie damals üblich, sorgte Johannes Kepler (1571-1630) durch das jährliche Einlagern einiger Fässer Wein für sich und seine Familie. Bald schon wunderte er sich aber über die Volumenvermessungstechnik der Fassmacher bzw. Weinlieferanten: Es wurde mit einer Rute durch das an der dicksten Fassstelle gelegene Spundloch zum Rand hin gemessen. Da Kepler klar wurde, dass so extrem unterschiedliche Fässer scheinbar das gleiche Volumen hatten, näherte er die Fassbegrenzung durch eine Parabel an und entwickelte so die nach ihm benannte Fassregel.

Zwar konvergiert die Simpson-Regel wesentlich besser als die Rechteck- oder die Trapezregel, jedoch stellt die unterschiedliche Gewichtung der Stützstellen ein rechentechnisches Ärgernis dar, welches durch interne Rundungen den Konvergenzvorteil durchaus wieder aufheben kann. Im Programm wird mit 14 Stellen Genauigkeit gerechnet, sodass dieser Effekt praktisch nicht auftritt.

Newtonsche 3/8-Regel / Newton-Cotes-Formeln
Interpoliert man die gewählten Teilintervalle nicht durch quadratische Parabeln, sondern durch kubische Funktionen ergibt sich die 3. Newton-Cotes-Formel, die sogenannte 3/8-Regel. Für einen Teilbogen gilt hier:

    \[ F(x)=\frac{3}{8}\Delta x(y_0+3y_1+3y_2+y_3) \]

Summiert man über n Teilintervalle, werden die Intervallgrenzen zyklisch mit den Faktoren 3,3 und 2 gewichtet.

Für noch größere Konvergenzgeschwindigkeit können Rundungsfehler in Abhängigkeit von der Funktion immer deutlicher werden. Die Entwicklung weiterer Newton-Cotes-Formeln kompliziert die Berechnung. Das Programm enthält die Newton-Cotes-Formeln 4. bis 7. Grades.
Zu beachten ist hierbei aber, dass nach einem Satz von Kusmin eine Erhöhung der zur Interpolation genutzten Polynome bei einigen Integranden sogar zur Divergenz führen kann. Die 7. Newton-Cotes-Formel garantiert jedoch in nahezu allen Anwendungsfällen eine hervorragende Konvergenz.

Offene Newton-Cotes-Formeln
Die oben genannten Newton-Cotes-Formeln bilden sogenannte geschlossene Formeln, d.h., die Randwerte des Integrationsintervalls gehen in die Berechnung ein. Für uneigentliche Integrale, bei denen die Funktion an den Intervallgrenzen nicht definiert ist, ist dies ein deutlicher Nachteil, da nun Grenzübergänge zu betrachten wären.

Abhilfe bringen die offenen Newton-Cotes-Formeln, auch Steffensen-Formeln genannt. Die Tatsache, dass nun der Anfangs- und Endwert des Integrals nicht mehr berücksichtigt und somit Fehler an nicht definierten Grenzen verhindert werden, führt aber mitunter zu deutlich schlechterer Konvergenz.

Im Programm sind zwei derartige Gleichungen unter der Bezeichnung 3., 4. und 5. offene Newton-Cotes-Formeln implementiert.

Beispiel: Für die Bestimmung des Integrals

    \[ \int_0^1\frac{1}{\sqrt{x}}\mathrm{d}x \]

wird die 4. offene Newton-Cotes-Formel und die 7. geschlossene Newton-Cotes-Formel genutzt. Da die Funktion f(x) = 1/√x für x = 0 nicht definiert ist, konvergiert die offene Formel mit 1,99 besser gegen den korrekten Wert von 2 als die geschlossene Formel mit 1,97.

Romberg-Verfahren , Steklov-Verfahren
Die Simpson- und die Newtonsche 3/8-Regel sind rechentechnisch relativ aufwendig. Der norwegische Mathematiker Romberg schlug deshalb 1955 ein Verfahren vor, das die Trapezformel nutzt und dennoch Genauigkeiten wie bei der Simpson-Regel erzielt.

Ist T(n) der Wert der Trapezregel für n Teilintervalle und T(2n) der Wert für 2n Teilintervalle, so bildet Romberg den Quotienten:

    \[ S_{2n}=\frac{4}{3}(T_{2n}-T_n) \]

und erhält so sehr günstige Näherungen. Während die Trapezformel mit einem Konvergenzfaktor q = 1/4 konvergiert, erreicht das Romberg-Verfahren q = 1/16. Eine Konvergenz von sogar q = 1/64 ergibt sich für

    \[ R_{2n}=\frac{16}{15}(S_{2n}-S_n) \]

Nach Steklov (1916) kann diese Entwicklung weiter fortgesetzt werden. Im Programm wird als Romberg-Verfahren S(2n) ermittelt. Unter dem Listeneintrag Steklov-Verfahren berechnet das Teilprogramm die angegebenen Terme R(2n).

Monte Carlo-Verfahren
Ist die Fläche unter einer Funktion f(x) in einem Intervall [a;b] gesucht, so kann die Berechnung auch mittels Zufallszahlen geschehen. Wählt man für das Argument gleichmäßig verteilte Zufallszahlen x(i) im Intervall [0;1] und berechnet zu jeder Zahl den Funktionswert f(x(i)), so ist der statistische Mittelwert M[f(x)] mit der Intervallbreite multipliziert, ein guter Schätzwert für den Flächeninhalt des gesuchten Intervalls. Da das arithmetische Mittel wiederum ein effektiver Wert für das statistische Mittel ist, gilt für die gesuchte Fläche:

    \[ A=\frac{1}{n}\sum{f(x_i)} \approx \int_0^1{f(x)\mathrm{d}x} \]

Entwickelt man diese Gleichung zu einer rekursiven Formel der Flächeninhalte A(i), ergibt sich

    \[ A_k=\frac{1}{k}(A_{k-1}(k-1)+f(x_k)) \]

Diese Gleichung wird im Programm genutzt. Vergleiche mit anderen Integrationsverfahren demonstriert die zufällige Konvergenzstärke – sie verdeutlicht aber auch, dass diese Anwendung der Monte-Carlo-Methode durchaus ihren praktischen Nutzen besitzt.

Gauß-Legendre-Formeln / Tschebyschow-Verfahren
Rechentechnisch noch anspruchsvoller als die Newton-Cotes-Formeln sind Formeln nach Gauß-Legendre. Der Grundgedanke dieser Gleichungen besteht darin, nicht nur unterschiedliche Gewichte, d.h. Faktoren der Stützstellen, sondern auch nicht äquidistante (unterschiedlicher Abstand) Stützstellen zu nutzen. Diese Stützstellen ergeben sich als die Nullstellen der zugehörigen Legendre-Polynome Pn(x).

Darüber hinaus muss das Integral auf die Grenzen von -1 bis 1 transformiert werden. So ergibt sich als 3. Gleichung (n = 2) nach Gauß-Legendre zum Beispiel:

    \[ I=\frac{h}{9} (5f(\frac{-h}{5}\sqrt{15}) +8f(0) +5f(\frac{h}{5}\sqrt{15}))+\frac{h^7}{15750}f^{(6)}(\xi) \]

wobei ξ eine Zwischenstelle des Intervalls [-h;h] ist.

Die Tragweite besteht darin, dass mit der n.ten Gauß-Legendre-Formel Polynome 2n-1. Grades bei nur n Stützstellen vollkommen exakt berechnet werden können. Damit sind diese Gleichungen weitaus leistungsfähiger als Newton-Cotes-Formeln gleichen Grades.

Sie können die 1. bis 5. Gauß-Legendre-Formel testen. Insbesondere die 5. Gleichung stellt ein sehr mächtiges Werkzeug für die numerische Integration dar. In vielen Beispielen genügen zwei oder fünf Stützstellen, um eine bessere Genauigkeit als 5 · 10-10 zu erreichen.

Pawnuti L. Tschebyschow (1821-1894) nutzte für seine Quadraturformeln ebenfalls den Grundgedanken von Gauß. Allerdings fordert er gleich große Gewichte, womit die Fehler aller Funktionswerte gleich stark eingehen. Die 1. und 2. Tschebyschow-Gleichung unterscheiden sich nicht von denen von Gauß-Legendre. Die 3. Tschebyschow-Formel lautet:

    \[ I=\frac{2h}{3}(f(\frac{-h}{2}\sqrt{12})+f(0)+f(\frac{h}{2}\sqrt{12}) \]

Weiterhin enthält das Programm mit den Einträgen Tschebyschow n = 4 bzw. n = 5 die 4. und 5. derartige Formel. Interessant ist, dass die Tschebyschow-Formeln nur bis n = 7 brauchbar sind, da für n = 8 bzw. n = 10 nicht reelle, komplexe Stützstellen entstehen.

Bessell-Quadratur
Einer Idee von Charles Hermite folgend ist es möglich, für die Quadratur einer Funktion auch deren Ableitungen zu nutzen. Die Problematik der Bestimmung der Ableitung kann man umgehen, indem man Differenzenquotienten einsetzt. Zusätzlich muss jedoch für die Integration der Funktion f(x) im Intervall [a,b] diese auch für a – Δx und b + Δx definiert sein. Dies führt zur Bessellschen Quadraturformel:

    \[ I=\frac{b-a}{n}(\frac{-1}{24}f(a-\Delta x)+\frac{f(a)}{2}+\frac{25}{24}f(a+\Delta x)+... \]

    \[ ...+\frac{25}{24}f(b-\Delta x)+\frac{f(b)}{2}-\frac{1}{24}f(b+\Delta x)) \]

Vergleich der Näherungsverfahren
In diesem Unterprogramm werden diese 22 Verfahren nebeneinander betrachtet. In den aufklappbaren Boxen wählen Sie jeweils eines der genannten Verfahren aus.

Das Programm berechnet die Fläche unter der Funktion für je 2, 5, 10, 20, 50, 75, 100, 200, 500, 750, 1000, 2000, 3000 und 5000 Stützstellen. Da die n-Newton-Cotes-Formel eine jeweils restlos durch n teilbare Stützstellenzahl benötigt, nutzt das Programm in diesen Fällen veränderte Stützstellenzahlen. Verwenden Sie Gauß-Legendre-Formeln, bricht das Programm die Berechnung bei einer erzielten Genauigkeit von 5 · 10-10 automatisch ab.

Zusätzlich kann die angezeigte Genauigkeit des Funktionswertes von 1 bis 9 Kommastellen (Voreinstellung 6) festgelegt werden. Mit dem Schalter Berechnung wird die numerische Integration begonnen. Die Berechnung wird ebenso abgebrochen, wenn die eingestellte Genauigkeit erreicht wurde.

Bemerkenswert ist, dass ab etwa 1000 Stützstellen die Genauigkeit bei sehr komplexen Integranden nicht mehr unbedingt zunimmt – unter Umständen sogar wieder zurückgeht. Die Ursache hierfür liegt in den internen Rundungen der Zwischenergebnisse.