Integration

 Allgemeine Integrationsaufgabe

Die Grundaufgabe dieses Kapitels ist zu einer gegebenen Funktion $ f\in C[a, b]$ das Integral
$\displaystyle \int_a^b f(x) dx$
numerisch zu bestimmen. Zur Lösung dieser Aufgabe wurden die verschiedensten Algorithmen entwickelt.

Ein Ansatz zur numerischen Integration einer Funktion $ f$ sind die sogenannten Newton-Cotes-Formeln.

Zunächst werden $ n+1$ paarweise verschiedene Stützstellen $ x_i$ gewählt mit

$\displaystyle x_i = a + i \left(\frac{b-a}{n}\right),$
wobei $ i = 0, \hdots , n$ ist.

Mit diesen Stützstellen wird dann das zur Funktion $ f$ gehörige, eindeutig bestimmte Interpolationspolynoms $ P_n (x)$ vom Grad $ n$ unter Verwendung der Lagrange-Polynome $ L_{i,n} (x), i=0,\hdots,n$ bestimmt. Man erhält also:

$\displaystyle P_n (x) = \sum_ {i=0}^n f( x_i ) L_{i,n} (x)$
Als Approximation an das Integral von $ f$ wird dann das Integral des Interpolationspolynom $ P_n$ genommen. Also:
$\displaystyle \int_a^b f(x) dx \approx \int_a^b P_n (x) dx = \int_a^b \sum ... ...}^n f( x_i ) L_{i,n} (x) dx = \sum_{i=0}^n f( x_i ) \int_a^b L_{i,n} (x) dx$
Der Vorteil dieser Darstellung ist, dass das Integral der Lagrange-Polynome nur abhängig ist von den Parametern $ n$ und $ b-a$. Man erhält:
$\displaystyle \int_a^b f(x) dx \approx \sum _{i=0}^n f( x_i ) \int_a^b L_{i,n} (x) dx = \sum_{i=0}^n f( x_i ) \left(\frac{b-a}{n}\right) w_{i,n}$
Die Werte der Gewichte $ w_{i,n}$ kann man aus der folgenden Tabelle entnehmen. Der Name der entsprechenden Regel ist in der letzten Spalte angegeben.
\begin{displaymath}\begin{array}{cc\vert ccccccc} & & & & i & & & & \ & & 0 ... ... & 7/90 & & \mbox{Milne-Regel} \ & \vdots & \ \end{array}\end{displaymath}
Da bei größeren Werten von $ n$ auch negative Gewichte auftauchen, werden die entsprechenden Formeln numerisch unbrauchbar. Um die Genauigkeit dieser Approximation dennoch zu erhöhen, teilt man zunächst das gesamte Intervall in kleinere Teilintervalle. Auf diesen Teilintervallen wendet man nun die Regeln an und addiert dann die ermittelten Werte. Sind die Teilintervalle gleich groß, so erhält man die sogenannten Summenregeln.

Exemplarisch wird hier die Trapezsumme vorgestellt.

Es seien folgende Stützstellen gegeben:

$\displaystyle x_i = a + i h$
wobei $ i = 0 , \hdots , m$ und $ h:=\frac{b-a}{m}$ ist. Die Trapezsumme ist dann
$\displaystyle \int_a^b f(x) dx \approx T_h (f) := h \left[ \frac{1}{2}f( x_0 ) + \sum_{i=0}^{m-1} f( x_i ) + \frac{1}{2} f( x_m ) \right]$

 Romberg Integration

Die Romberg Integration beruht auf der oben vorgestellten Trapezsumme. Man kann zeigen, dass für die Trapezsumme einer $ (2m+2)$-mal stetig differenzierbaren Funktion folgende Fehlerdarstellung gilt:
$\displaystyle \int_a^b f(x) dx - T_h (f) = \sum _{k=1}^m a_k h^{2k} + O( h^{2m+2} ),$
wobei $ a_k$ von $ h$ unabhängige Konstanten sind. Insbesondere ist der Fehler bei der Trapezsumme vom Typ $ O( h^2 )$. Eine Herleitung des Romberg Verfahrens besteht darin die Faktoren $ a_k$ durch geschickte Addition zu eliminieren und so die Fehlerordnung zu erhöhen.

Dazu werden zunächst die Treppensummen $ T_{h_i} (f)$ für $ i = 0 , \hdots , m$ mit $ h_i =\frac{b-a}{2m}$ bestimmt. Nun wird das folgende Romberg-Schema angewandt:

\begin{displaymath}\begin{array}{cccccccccc} T_{h_0}=T_{0,0} & & & & & & & & & ... ..._{m,2} & \to & T_{m,3} & \cdots & \to & T_{m,m} \ \end{array}\end{displaymath}
Dabei ergibt sich $ T_{i,j}$ durch die folgende Additionsformel:
$\displaystyle T_{i,j} = \frac{4^ j \, T_{i,j-1} - T_{i-1,j-1}}{4^ j - 1}$
Damit lässt sich zeigen, dass in jeder Spalte die Ordnung des Fehlers um 2 erhöht wird, in der $ m$-ten Spalten also ein Fehler vom Typ $ O( h^{2m+2} )$ vorliegt.

Man kann ebenso zeigen, dass die zweite Spalte der Simpsonsumme entspricht, d.h. $ T_{i,1}$ entspricht der Simpsonsumme zur Weite $ h_i$.

Das folgende Programm führt die Rombergintegration aus. Man kann Spalten aus dem obigen Rombergschema auswählen und sich den jeweiligen Fehler zum exakten Integral anzeigen lassen. Außerdem kann man sich den Fehler der Diagonalen, d.h. von $ T_{1,1} , T_{2,2} , \hdots , T_{m,m}$ ausgeben lassen.

 Adaptive Integration

Die Romberg Integration liefert bei glatten Funktionen sehr gute Approximationen an das Integral. Jedoch bei stark variierenden Funktionen ist die Romberg Integration eher ineffektiv. Man würde dann lieber nicht äquidistante Stützstellen je nach dem Verhalten von $ f$ verwenden. Diesen Ansatz verfolgt die adaptive Integration.

Zur Erläuterung des Verfahrens wird das folgende Teilintervall $ [ x_i , x_{i+1} ]$ betrachtet. Auf diesem Intervall wird sowohl die Treppensumme $ T$, als auch die Simpsonregel $ S$ berechnet. Anschließend wird geprüft, ob folgende Bedingung erfüllt ist.

$\displaystyle \frac{\vert T-S\vert}{\vert W\vert} \leq \varepsilon \frac{x_{i+1} - x_i}{b-a}$
Hierbei ist $ W$ ein sehr grober Schätzwert für das exakte Integral, welcher die Größenordnung erfasst, und $ \varepsilon$ die geforderte relative Genauigkeit. Wird die obige Bedingung erfüllt, so wird $ S$ als Approximation auf diesem Teilintervall akzeptiert und das nächste Teilintervall betrachtet.

Wird die obige Bedingung nicht erfüllt, so wird das Intervall $ [ x_i , x_{i+1} ]$ in die zwei Teilintervalle $ [ x_i , \tau ]$ und $ [ \tau , x_{i+1} ]$ aufgeteilt, wobei $ \tau = \frac{x_i + x_{i+1}}{2}$ ist. Anschließend wird diese Methode zunächst auf das Intervall $ [ x_i , \tau ]$ und dann auf das Intervall $ [ \tau , x_{i+1} ]$ angewendet.

Das folgende Programm führt die adaptive Integration aus. Man kann sowohl die relative Genauigkeit $ \varepsilon$, als auch die maximale Anzahl an Stützstellen angeben. Anschließend wird die Bestimmung der Stützstellen als Slideshow angezeigt.