Kurs:Numerik I/8 Numerische Integration

8.1 Einführung

In Anwendungen der Mathematik müssen häufig Riemann-Integrale für stückweise stetige Funktionen berechnet werden. In vielen Fällen ist eine geschlossene Lösung eines solchen Integrals nicht bekannt, so dass es näherungsweise numerisch gelöst werden muss. Die numerische Lösung eines Integrals bezeichnet man auch als numerische Quadratur. In diesem Abschnitt sollen eine Reihe von Formeln zur numerischen Integration hergeleitet und untersucht werden.

Das Integral über eine stückweise stetige Funktion kann bekanntlich als Summe von Integralen über stetige Funktionen beschrieben werden, so dass wir uns auf die Betrachtung stetiger Funktionen beschränken können. Dazu definieren wir den Operator ${\mathcal {I}}:C[a,b]\to \mathbb {R}$ mit

${\mathcal {I}}(f):=\int \limits _{a}^{b}f(x)\,dx$

für $f\in C[a,b]$. Dieser ist linear, da für alle $f,g\in C[a,b]$ und $\alpha ,\beta \in \mathbb {R}$

${\mathcal {I}}(\alpha f+\beta g)=\int \limits _{a}^{b}[\alpha f(x)+\beta g(x)]\,dx=\alpha \int \limits _{a}^{b}f(x)\,dx+\beta \int \limits _{a}^{b}g(x)\,dx=\alpha {\mathcal {I}}(f)+\beta {\mathcal {I}}(g)$

gilt und er ist positiv, d. h. man hat

$f\in C[a,b],f\geq 0\Rightarrow {\mathcal {I}}(f)\geq 0,$

wobei $f\geq 0$ bedeutet, dass $f(x)\geq 0,x\in [a,b]$ ist. Gesucht sind nun einfach auszuwertende Formeln, die jedem $f\in C[a,b]$ einen Näherungswert ${\hat {\mathcal {I}}}(f)$ für den Wert des Integrals zuordnen und zwar so, dass der Quadraturfehler ${\mathcal {I}}(f)-{\hat {\mathcal {I}}}(f)$ möglichst klein ist.

Definition 8.1

Unter einer Quadraturformel ${\mathcal {I}}_{n}:C[a,b]\to \mathbb {R}$ zur Berechnung des bestimmten Integrals ${\mathcal {I}}(f)$ versteht man eine Summe

(8.1) ${\mathcal {I}}_{n}(f):=(b-a)\sum _{i=0}^{n}\sigma _{i}f(x_{i})$

für $f\in C[a,b]$ mit bekannten Gewichten $\sigma _{i}\in \mathbb {R}$ $(i=0,1,\ldots ,n)$ und Stützstellen bzw. Knoten $x_{i}\in [a,b]$ $(i=0,1,\ldots ,n)$, wobei $x_{i}\neq x_{j}$ $(i\neq j)$ sei.

Wenn wir die Abhängigkeit der Gewichte und Stützstellen von der Wahl von $n$ darstellen wollen, schreiben wir statt $\sigma _{i}$ und $x_{i}$ auch $\sigma _{i}^{(n)}$ und $x_{i}^{(n)}$. Wie ${\mathcal {I}}$ ist auch eine Quadraturformel ${\mathcal {I}}_{n}$ ein linearer Operator, denn man hat

${\mathcal {I}}_{n}(\alpha f+\beta g)=(b-a)\sum _{i=0}^{n}\sigma _{i}[\alpha f(x_{i})+\beta g(x_{i})]=\alpha (b-a)\sum _{i=0}^{n}\sigma _{i}f(x_{i})+\beta (b-a)\sum _{i=0}^{n}\sigma _{i}g(x_{i})$

(8.2) $=\alpha {\mathcal {I}}_{n}(f)+\beta {\mathcal {I}}_{n}(g)$

für alle $f,g\in C[a,b]$ und $\alpha ,\beta \in \mathbb {R}$. Sind die Gewichte nichtnegativ, d. h. hat man $\sigma _{i}\geq 0$ $(i=0,1,\ldots ,n)$, so ist ferner mit ${\mathcal {I}}$ auch ${\mathcal {I}}_{n}$ positiv und gilt also

$f\in C[a,b],f\geq 0\Rightarrow {\mathcal {I}}_{n}(f)\geq 0.$

Wir definieren weiter:

Definition 8.2

Eine Quadraturformel ${\mathcal {I}}_{n}$ hat mindestens den Genauigkeitsgrad $r\in \mathbb {N} _{0}$, wenn

(8.3) ${\mathcal {I}}_{n}(x^{j})={\mathcal {I}}(x^{j}),\quad j=0,1,\ldots ,r$

gilt. Im Fall, dass zusätzlich ${\mathcal {I}}_{n}(x^{r+1})\neq {\mathcal {I}}(x^{r+1})$ richtig ist, sagt man, dass ${\mathcal {I}}_{n}$ den Genauigkeitsgrad $r\in \mathbb {N} _{0}$ hat.

Da ${\mathcal {I}}$ und ${\mathcal {I}}_{n}$ lineare Operatoren sind, folgt aus (8.3)

(8.4) ${\mathcal {I}}_{n}\left(\sum _{j=0}^{r}a_{j}x^{j}\right)=\sum _{j=0}^{r}a_{j}{\mathcal {I}}_{n}(x^{j})=\sum _{j=0}^{r}a_{j}{\mathcal {I}}(x^{j})={\mathcal {I}}\left(\sum _{j=0}^{r}a_{j}x^{j}\right)$

für alle $a_{j}\in \mathbb {R}$ $(j=0,1,\ldots ,r)$, und damit der folgende Satz, wobei $\Pi _{r}$ wieder den Raum aller Polynome vom Grad $\leq r$ bezeichnet:

Satz 8.3

${\mathcal {I}}_{n}$ ist genau dann eine Quadraturformel von mindestens dem Genauigkeitsgrad $r$, wenn gilt:

${\mathcal {I}}_{n}(p)={\mathcal {I}}(p),\quad p\in \Pi _{r}.$

Wir bemerken in diesem Zusammenhang:

Satz 8.4

Zu $n+1$ Stützstellen $x_{i}$ $(i=0,1,\ldots ,n)$ mit $x_{i}\neq x_{k}$ $(i\neq k)$ gibt es genau eine Quadraturformel ${\mathcal {I}}_{n}$, welche mindestens den Genauigkeitsgrad $n\in \mathbb {N} _{0}$ hat, d. h. für die gilt:

(8.5) ${\mathcal {I}}_{n}(p)={\mathcal {I}}(p),\quad p\in \Pi _{n}.$

Diese hat die Gewichte

(8.6) $\sigma _{i}:={\frac {1}{b-a}}{\mathcal {I}}(L_{i}),\quad i=0,1,\ldots ,n,$

wobei $L_{i}\in \Pi _{n}$ $(i=0,1,\ldots ,n)$ die zu den $x_{i}$ gehörenden Lagrangeschen Basispolynome sind (vgl. Definition 6.2).

Beweis.

Für die durch die Stützstellen $x_{i}$ und Gewichte $\sigma _{i}$ in (8.6) definierte Quadraturformel ${\mathcal {I}}_{n}$ gilt für $k=0,1,\ldots ,n$

(8.7) ${\mathcal {I}}_{n}(L_{k})=(b-a)\sum _{i=0}^{n}\sigma _{i}L_{k}(x_{i})=(b-a)\sum _{i=0}^{n}\sigma _{i}\delta _{ki}=(b-a)\sigma _{k}={\mathcal {I}}(L_{k}).$

Da sich jedes Polynom vom Grad $\leq n$ auf eindeutige Weise als Linearkombination der $L_{k}$ $(k=0,1,\ldots ,n)$ darstellen lässt (vgl. (6.6)) und ${\mathcal {I}}$ sowie ${\mathcal {I}}_{n}$ lineare Operatoren sind, folgt damit analog zu (8.4) die Beziehung (8.5). Die so definierte Quadraturformel ${\mathcal {I}}_{n}$ ist eindeutig. Denn für jede andere Quadraturformel ${\mathcal {J}}_{n}$ mit Gewichten ${\hat {\sigma }}_{k}$ und Genauigkeitsgrad $n\geq 0$ hat man wegen $L_{k}\in \Pi _{n}$ die Identität ${\mathcal {J}}_{n}(L_{k})={\mathcal {I}}(L_{k})$ und bekommt man analog zu (8.7) ${\mathcal {J}}_{n}(L_{k})=(b-a){\hat {\sigma }}_{k}$, so dass $(b-a){\hat {\sigma }}_{k}={\mathcal {I}}(L_{k})$ und demnach ${\hat {\sigma }}_{k}=\sigma _{k}$ folgt.

q.e.d.

Weiter stellen wir fest:

Satz 8.5

Ist ${\mathcal {I}}_{n}$ eine Quadraturformel, die einen Genauigkeitsgrad $r\in \mathbb {N} _{0}$ hat, so folgt für ihre Gewichte $\sigma _{i}$

$\sum _{i=0}^{n}\sigma _{i}=1.$

Beweis.

Da ${\mathcal {I}}_{n}$ einen Genauigkeitsgrad $r\geq 0$ hat, folgt

$(b-a)\sum _{i=0}^{n}\sigma _{i}={\mathcal {I}}_{n}(1)={\mathcal {I}}(1)=\int \limits _{a}^{b}1\,dx=b-a.$

q.e.d.

Bezüglich der Konvergenz der durch eine Quadraturformel

${\mathcal {I}}_{n}(f):=(b-a)\sum _{i=0}^{n}\sigma _{i}^{(n)}f(x_{i}^{(n)}$

erzeugten Näherungswerte ${\mathcal {I}}_{n}(f)$ gegen den exakten Wert des Integrals ${\mathcal {I}}(f)$ für $n\to \infty$ kann man allgemein den folgenden Satz angeben, den wir hier jedoch nicht beweisen können (für einen Beweis siehe H. Heuser: Funktionalanalysis, Teubner, Stuttgart, 1992, S. 268).

Satz 8.6 (Szegö)

Man hat

$\lim _{n\to \infty }{\mathcal {I}}_{n}(f)={\mathcal {I}}(f),\quad f\in C[a,b]$

genau dann, wenn gilt:

(a) $\sum _{i=1}^{n}\left|\sigma _{i}^{(n)}\right|\leq M,\quad n\in \mathbb {N} ,$

(b) $\lim _{n\to \infty }{\mathcal {I}}_{n}(x^{j})={\mathcal {I}}(x^{j}),\quad j\in \mathbb {N} _{0}.$

Mit Hilfe von Satz 8.5 erschließt man ferner:

Korollar 8.7

Es sei ${\mathcal {I}}_{n}$ eine Quadraturformel mit Gewichten $\sigma _{i}^{(n)}\geq 0\ (i=0,1,\ldots ,n)$ für alle $n\in \mathbb {N}$ und einem Genauigkeitsgrad $n\geq 0$. Dann hat man

$\lim _{n\to \infty }{\mathcal {I}}_{n}(f)={\mathcal {I}}(f),\quad f\in C[a,b]$

genau dann, wenn gilt:

$\lim _{n\to \infty }{\mathcal {I}}_{n}(x^{j})={\mathcal {I}}(x^{j}),\quad j\in \mathbb {N} _{0}.$

8.2 Interpolatorische Quadraturformeln

8.2.1 Allgemeines

Es seien nun $f\in C[a,b]$ und $x_{i}\in [a,b],i=0,1,\ldots ,n$ mit $x_{i}\neq x_{k}$ $(i\neq k)$ gegeben und $Q_{n}\in \Pi _{n}$ bezeichne das (eindeutige) Interpolationspolynom zu den Stützpunkten $(x_{i},f(x_{i})),i=0,1,\ldots ,n$. Sind $L_{i}\in \Pi _{n}$ $(i=0,1,\ldots ,n)$ wieder die zu den $n+1$ Stützstellen $x_{i}$ gehörenden Lagrangeschen Basispolynome, so kann das Interpolationspolynom $Q_{n}$ damit gemäß (6.7)in der Form

$Q_{n}(x)=\sum _{i=0}^{n}f(x_{i})L_{i}(x)$

geschrieben werden. Wir definieren nun:

Denition 8.8

Eine Quadraturformel ${\mathcal {I}}_{n}$ mit

${\mathcal {I}}_{n}(f):={\mathcal {I}}(Q_{n})=\int \limits _{a}^{b}Q_{n}(x)\,dx=(b-a)\sum _{i=0}^{n}\left\{{\frac {1}{b-a}}\int \limits _{a}^{b}L_{i}(x)\,dx\right\}f(x_{i}),$

d. h. mit Gewichten

(8.8) $\sigma _{i}:={\frac {1}{b-a}}{\mathcal {I}}(L_{i}),\quad i=0,1,\ldots ,n$

heißt interpolatorische Quadraturformel.

Wegen der Übereinstimmung der Gewichte in (8.8) und (8.6) können wir mit Satz 8.4 schließen:

Korollar 8.9

Eine interpolatorische Quadraturformel ${\mathcal {I}}_{n}$ hat mindestens den Genauigkeitsgrad $n\in \mathbb {N} _{0}$ und ist zu den gegebenen Stützstellen die einzige Quadraturformel mit einem Genauigkeitsgrad $\geq n$.

Ferner können wir zeigen:

Satz 8.10

Eine interpolatorische Quadraturformel In besitzt die Gestalt

${\mathcal {I}}_{n}(f)=(b-a)\sum _{i=0}^{n}\sigma _{i}f(x_{i})$

mit

(8.9) $\sigma _{i}:=\int \limits _{0}^{1}\prod _{k=0 \atop k\neq i}^{n}{\frac {t-t_{k}}{t_{i}-t_{k}}}\,dt$ mit $t_{k}:={\frac {x_{k}-a}{b-a}}.$

Beweis.

Mit $t_{k}$ wie in (8.9) lassen sich die Gewichte $\sigma _{i}$ $(i=0,1,\ldots ,n)$ aus (8.8) mit Hilfe der Substitution $x:=(b-a)t+a$ umschreiben in

$\sigma _{i}={\frac {1}{b-a}}{\mathcal {I}}(L_{i})={\frac {1}{b-a}}\int \limits _{a}^{b}L_{i}(x)\,dx={\frac {1}{b-a}}\int \limits _{a}^{b}\prod _{k=0 \atop k\neq i}^{n}{\frac {x-x_{k}}{x_{i}-x_{k}}}\,dx$

(8.10) $=\int \limits _{0}^{1}\prod _{k=0 \atop k\neq i}^{n}{\frac {(b-a)t-(b-a)t_{k}}{(b-a)t_{i}-(b-a)t_{k}}}\,dt=\int \limits _{0}^{1}\prod _{k=0 \atop k\neq i}^{n}{\frac {t-t_{k}}{t_{i}-t_{k}}}\,dt.$

q.e.d.

Die Transformation von $x$ nach $t$ in (8.9) ist sinnvoll, da damit die Gewichte $\sigma _{i}$ in einer interpolatorischen Quadraturformel von den Intervallgrenzen $a$ und $b$ unabhängig werden und nur von der relativen Verteilung der Stützstellen in $[a,b]$ abhängen.

8.2.2 Newton-Cotes-Formeln

Wir wollen nun auf spezielle interpolatorische Quadraturformeln, die Newton-Cotes-Formeln, eingehen. Diese ergeben sich durch äquidistante Wahl der Stützstellen in $[a,b]$. Insbesondere erhält man die abgeschlossenen Newton-Cotes-Formeln, wenn die Randpunkte des Intervalls $[a,b]$ selbst Stützstellen sind, wenn also für $n\geq 1gilt::(8.11)<math>x_{i}:=a+ih\quad (i=0.1,\ldots ,n),\qquad h:=(b-a)/n.$ Bei den offenen Newton-Cotes-Formeln sind die Randpunkte von $[a,b]$ selbst keine Stützstellen, so dass man

$x_{i}:=a+(i+1)h\quad (i=0,1,\ldots ,n),\qquad h:=(b-a)/(n+2)$

hat. Wir wollen hier nur die abgeschlossenen Newton-Cotes-Formeln genauer untersuchen.

Lemma 8.11

Für die Gewichte $\sigma _{i}$ $(i=0,1,\ldots ,n)$ der abgeschlossenen Newton-Cotes-Formeln gilt:

(8.12) $\sigma _{n-i}=\sigma _{i}={\frac {1}{n}}\int \limits _{0}^{n}\prod _{k=0 \atop k\neq i}^{n}{\frac {s-k}{i-k}}\,ds.$

Beweis.

Die zweite Identität in (8.12) folgt mit

$t_{k}:={\frac {x_{k}-a}{b-a}}={\frac {kh}{b-a}}={\frac {k}{n}}$

aus (8.9) mit der Substitution $t:=s/n$, denn man hat

$\sigma _{i}=\int \limits _{0}^{1}\prod _{k=0 \atop k\neq i}^{n}{\frac {t-{\frac {k}{n}}}{{\frac {i}{n}}-{\frac {k}{n}}}}\,dt={\frac {1}{n}}\int \limits _{0}^{n}\prod _{k=0 \atop k\neq i}^{n}{\frac {s-k}{i-k}}\,ds.$

Somit müssen wir noch die erste Identität in (8.12) zeigen. Dazu sei $i\in \{0,\ldots ,n\}$. Sind $L_{i}\in \Pi _{n}$ die Lagrangeschen Basispolynome, so ist $L_{n-i}\in \Pi _{n}$ und $Q(x):=L_{i}(b+a-x)\in \Pi _{n}$ sowie

$Q(x_{n-j})=L_{i}(b+a-x_{n-j})=L_{i}\left(b+a-\left[a+(n-j){\frac {b-a}{n}}\right]\right)=L_{i}\left(a+j{\frac {b-a}{n}}\right)=L_{i}(x_{j})$

$=\delta _{ij}=\delta _{n-i,n-j}=L_{n-i}(x_{n-j})$

für $j=0,1,\ldots ,n$. Da $L_{n-i}$ und $Q$ demnach offenbar Interpolationspolynome zu den Punkten $(x_{n-j},\delta _{n-i,n-j}),j=0,1,\ldots ,n$ sind, muss wegen der Eindeutigkeit des Interpolationspolynoms $L_{n-i}\equiv Q$ gelten, so dass wir schließlich mit der Substitution $t:=b+a-x$ Folgendes erhalten (vgl. (8.10)):

$\sigma _{n-i}={\frac {1}{b-a}}\int \limits _{a}^{b}L_{n-i}(x)\,dx={\frac {1}{b-a}}\int \limits _{a}^{b}Q(x)\,dx={\frac {1}{b-a}}\int \limits _{a}^{b}L_{i}(b+a-x)\,dx={\frac {1}{b-a}}\int \limits _{a}^{b}L_{i}(t)\,dt=\sigma _{i}.$

q.e.d.

Wir geben nun einige Spezialfälle der abgeschlossenen Newton-Cotes-Formeln an.

Beispiel 8.12

(1) Für $n=1$ hat eine interpolatorische Quadraturformel nach Satz 8.10 die Gestalt

${\mathcal {I}}_{1}(f)=(b-a)[\sigma _{0}f(x_{0})+\sigma _{1}f(x_{1})].$

Dabei ergeben sich für die zugehörige abgeschlossene Newton-Cotes-Formel mit $h=b-a$ die Stützstellen $x_{0}=a$ und $x_{1}=b$ und wegen $\sigma _{1}=\sigma _{0}$ (Lemma 8.11) und $\sigma _{0}+\sigma _{1}=1$ (Satz 8.5) die Gewichte

$\sigma _{1}=\sigma _{0}={\frac {1}{2}}.$

Man erhält so die (Sehnen-) Trapezregel

(8.13) ${\mathcal {I}}_{1}(f):={\frac {b-a}{2}}[f(a)+f(b)].$

(2) Für $n=2$ hat man mit $h=(b-a)/2$ die Stützstellen $x_{0}=a,x_{1}=(a+b)/2$ und $x_{2}=b$ und unter Verwendung von Lemma 8.11 und anschließend Satz 8.5 die Gewichte

$\sigma _{2}=\sigma _{0}={\frac {1}{2}}_{0}^{2}{\frac {(s-1)(s-2)}{(0-1)(0-2)}}\,ds={\frac {1}{4}}_{0}^{2}(s^{2}-3s+2)\,ds={\frac {1}{4}}\left[{\frac {8}{3}}-6+4\right]={\frac {1}{6}},$

$\sigma _{1}=1-\sigma _{0}-\sigma _{2}=1-{\frac {2}{6}}={\frac {2}{3}}.$

Unter Verwendung von Satz 8.10 ergibt sich so die Simpson-Regel bzw. Keplersche Fassregel

(8.14) ${\mathcal {I}}_{2}(f):={\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right]$

(3) Der Fall $n=3$ führt auf die Newtonsche 3/8-Regel

${\mathcal {I}}_{3}(f):={\frac {b-a}{8}}\left[f(a)+3f\left({\frac {2a+b}{3}}\right)+3f\left({\frac {a+2b}{3}}\right)+f(b)\right].$

(4) Für $n=4$ bekommt man die Milne-Regel

${\mathcal {I}}_{4}(f):={\frac {b-a}{90}}\left[7f(a)+32f\left({\frac {3a+b}{4}}\right)+12f\left({\frac {2a+2b}{4}}\right)+32f\left({\frac {a+3b}{4}}\right)+7f(b)\right].$

Als Beispiel berechnen wir ein Integral näherungsweise mit der Simpson-Regel.

Beispiel 8.13

Es seien $f(x):=1/(1+x^{2}),a=0$ und $b=1$, so dass

${\mathcal {I}}(f)=\int \limits _{0}^{1}{\frac {1}{1+x^{2}}}\,dx$

ist. Die Simpson-Regel liefert dafür den Näherungswert

${\mathcal {I}}_{2}(f)={\frac {1}{6}}\left[f(0)+4f\left({\frac {1}{2}}\right)+f(1)\right]={\frac {1}{6}}\left(1+4\cdot {\frac {4}{5}}+{\frac {1}{2}}\right)={\frac {47}{60}}=0.783\,33.$

Der exakte Wert des Integrals lautet hier

${\mathcal {I}}(f)=\arctan(x)|_{x=0}^{x=1}=\arctan(1)=0.785\,40.$

Für $n\leq 7$ sind die Gewichte in den abgeschlossenen Newton-Cotes-Formeln nichtnegativ und sind diese Quadraturformeln demzufolge positiv. Für $n=8$ und $n\geq 10$ treten negative Gewichte auf und ist damit als Folge von Satz 8.5

$\sum _{i=0}^{n}\left|\sigma _{i}^{(n)}\right|>1,$

was zu einer Verstärkung von Rundungsfehlern bei den Funktionswerten $f(x_{i})$ führt. Die Verwendung der abgeschlossenen Newton-Cotes-Formeln für $n\geq 8$ ist daher nicht zu empfehlen. Für die (abgeschlossenen) Newton-Cotes-Formeln lässt sich sogar

$\lim _{n\to \infty }\sum _{i=0}^{n}\left|\sigma _{i}^{(n)}\right|=1$

beweisen (Satz von Kusmin), so dass man aus dem Satz 8.6 von Szegö die Existenz eines $f\in C[a,b]$ schließen kann, für das die Konvergenz $\lim _{n\to \infty }{\mathcal {I}}_{n}(f)={\mathcal {I}}(f)$ nicht gilt. Letzteres lässt ja auch der Satz 6.24 von Faber generell für interpolatorische Quadraturformeln vermuten. Eine Erhöhung von $n$ bei den (abgeschlossenen) Newton-Cotes-Formeln muss also nicht zwangsläufig zu einer genaueren Näherung ${\mathcal {I}}_{n}(f)$ von ${\mathcal {I}}(f)$ führen.

Wir geben hier noch einige weitere interpolatorische Quadraturformeln an.

Beispiel 8.14

(1) Für $n=0$ und $x_{0}:=a$ oder $x_{0}:=b$ muss wegen Satz 8.5 $\sigma _{0}=1$ gelten, so dass man alternativ folgende beiden Rechteckregeln erhält:

(8.15) ${\mathcal {I}}_{0}(f):=(b-a)f(a),\quad {\mathcal {I}}_{0}(f):=(b-a)f(b).$

(2) Für $n=0$ bekommt man im Fall der offenen Newton-Cotes-Formeln $h:=(b-a)/2$ und

$x_{0}:=a+{\frac {b-a}{2}}={\frac {a+b}{2}},\quad \sigma _{0}=\sum _{i=0}^{n}\sigma _{i}=1$

und damit eine weitere Rechteckregel, die Mittelpunktregel

${\mathcal {I}}_{0}(f):=(b-a)f\left({\frac {a+b}{2}}\right).$

(3) Die offene Newton-Cotes-Formel für $n=1$ lautet mit $h:=(b-a)/3$,

$x_{0}:=a+{\frac {1}{3}}(b-a),\quad x_{1}:=a+{\frac {2}{3}}(b-a)$

und den Gewichten $\sigma _{0}=\sigma _{1}=1/2$, die man aus der Formel (8.9) errechnet, wie folgt:

${\mathcal {I}}_{1}(f):={\frac {b-a}{2}}\left[f\left({\frac {2a+b}{3}}\right)+f\left({\frac {a+2b}{3}}\right)\right].$

(4) Die offene Newton-Cotes-Formel für $n=2$ lautet mit $h:=(b-a)/4$,

$x_{0}:=a+{\frac {1}{4}}(b-a),\quad x_{1}:=a+{\frac {1}{2}}(b-a),\quad x_{2}:=a+{\frac {3}{4}}(b-a)$

und den mit Hilfe von (8.9) zu berechnenden Gewichten wie folgt:

${\mathcal {I}}_{2}(f):={\frac {b-a}{3}}\left[2f\left({\frac {3a+b}{4}}\right)-f\left({\frac {a+b}{2}}\right)+2f\left({\frac {3a+b}{4}}\right)\right].$

Man beachte, dass sie ein negatives Gewicht beinhaltet.

8.2.3 Quadraturfehler und Genauigkeitsgrad

Für den durch eine beliebige interpolatorische Quadraturformel in Bezug auf den exakten Wert des Integrals entstehenden Fehler, kann man die im folgenden Satz angegebene Abschätzung beweisen.

Satz 8.15

Es sei ${\mathcal {I}}_{n}$ eine interpolatorische Quadraturformel mit Stützstellen $x_{i}$ $(i=0,1,\ldots ,n)$, welche mindestens den Genauigkeitsgrad $r\geq n$ besitze und es sei $f\in C^{(r+1)}[a,b]$. Dann gilt

(8.16) $|{\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)|\leq \gamma _{r}{\frac {(b-a)^{r+2}}{(r+1)!}}\max _{\xi \in [a,b]}\left|f^{(r+1)}(\xi )\right|$

für

$\gamma _{r}:=\min _{t_{n+1},\ldots ,t_{r}\in [0,1]}\int \limits _{0}^{1}\prod _{k=0}^{r}|t-t_{k}|\,dt$

mit

(8.17) $t_{k}:={\frac {x_{k}-a}{b-a}}\quad (k=0,1,\ldots ,n).$

Hat man insbesondere für die $t_{k}$ $(k=0,\ldots ,n)$ aus (8.17) und frei wählbare $t_{k}$ $(k=n+1,\ldots ,r)$ mit

$s(t):=\prod _{k=0}^{r}(t-t_{k})$

die Beziehung $s(t)\geq 0,t\in [0,1]$ oder $s(t)\leq 0,t\in [0,1]$, dann folgt mit

${\hat {\gamma }}_{r}:=\int \limits _{0}^{1}s(t)\,dt$

und einem $\xi \in [a,b]$ die Fehlerdarstellung

(8.18) ${\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)={\hat {\gamma }}_{r}{\frac {(b-a)^{r+2}}{(r+1)!}}f^{(r+1)}(\xi ).$

Beweis.

Seien $x_{i}\in [a,b],i=n+1,\ldots ,r$ zunächst beliebig gewählt, so dass die $x_{i}$ $(i=0,\ldots ,n,n+1,\ldots ,r)$ paarweise verschieden sind und sei $Q_{r}\in \Pi _{r}$ das Interpolationspolynom zur den Stützpunkten $(x_{i},f(x_{i})),i=0,1,\ldots ,r$. Da ${\mathcal {I}}_{n}$ den Genauigkeitsgrad $r$ hat, gilt dann

${\mathcal {I}}_{n}(f)=(b-a)\sum _{i=0}^{n}\sigma _{i}f(x_{i})=(b-a)\sum _{i=0}^{n}\sigma _{i}Q_{r}(x_{i})={\mathcal {I}}_{n}(Q_{r})=I(Qr)$

und demnach

${\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)={\mathcal {I}}(f)-I(Q_{r})=\int \limits _{a}^{b}[f(x)-Q_{r}(x)]\,dx.$

Mit

$\omega (x):=(x-x_{0})\cdots (x-x_{n}),\quad \varphi (x):=(x-x_{n+1})\cdots (x-x_{r})$

und unter Verwendung von Satz 6.11 hat man für ein $\xi (x)\in [a,b]$

$f(x)-Q_{r}(x)={\frac {1}{(r+1)!}}\omega (x)\varphi (x)f^{(r+1)}(\xi (x)).$

Da die linke Seite der letzten Gleichung stetig in $x$ ist, ist es auch die rechte Seite und darum hat man

(8.19) ${\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)={\frac {1}{(r+1)!}}\int \limits _{a}^{b}\left[\omega (x)\varphi (x)f^{(r+1)}(\xi (x))\right]\,dx.$

Man beachte nun, dass die $x_{i}$ $(i=0,1,\ldots ,n)$ durch die Quadraturregel festgelegt sind. Wir wollen abschließend zeigen, dass für die Stützstellen $x_{i}$ $(i=n+1,\ldots ,r)$ die anfangs gemachte Voraussetzung hinsichtlich der paarweisen Unterschiedlichkeit fallen gelassen werden kann. Es seien daher jetzt letztere Punkte vollkommen beliebig aus [a, b] gewählt. Für jedes $m\in \mathbb {N}$ können wir dann Punkte $x_{i}^{(m)}$ $(i=n+1,\ldots ,r)$ finden, die zusammen mit den $x_{i}$ $(i=0,1,\ldots ,n)$ paarweise unterschiedlich sind und für die

$\lim _{m\to \infty }x_{i}^{(m)}=x_{i}\quad (i=n+1,\ldots ,r)$

gilt. Setzen wir

$\varphi _{m}(x):=(x-x_{n+1}^{(m)})\cdots (x-x_{r}^{(m)}),$

so hat man unter Verwendung des ersten Teils des Beweises

$|{\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)|\leq {\frac {1}{(r+1)!}}\max _{\xi \in [a,b]}\left|f^{(r+1)}(\xi )\right|\int \limits _{a}^{b}|\omega (x)\varphi _{m}(x)|\,dx$

$\leq {\frac {1}{(r+1)!}}\max _{\xi \in [a,b]}\left|f^{(r+1)}(\xi )\right|\{\int \limits _{a}^{b}|\omega (x)\varphi (x)|\,dx+\underbrace {\int \limits _{a}^{b}|\omega (x)||\varphi _{m}(x)-\varphi (x)|\,dx} _{\to 0\quad (m\to \infty )}\}$

(8.20) $\leq {\frac {1}{(r+1)!}}\max _{\xi \in [a,b]}\left|f^{(r+1)}(\xi )\right|\int \limits _{a}^{b}|\omega (x)\varphi (x)|\,dx.$

Wählt man nun $x_{i}$ $(i=n+1,\ldots ,r)$ so, dass der Wert des letzten Integrals minimal wird und wendet man die Substitution $x:=(b-a)t+a$ an, so gelangt man schließlich zu

$\gamma _{r}:=\min _{x_{n+1},\ldots ,x_{r}\in [a,b]}\int \limits _{a}^{b}\prod _{i=0}^{r}|x-x_{i}|\,dx=(b-a)^{r+2}\min _{t_{n+1},\ldots ,t_{r}\in [0,1]}\int \limits _{0}^{1}\prod _{i=0}^{r}|t-t_{i}|\,dt,$

womit die Abschätzung (8.16) gezeigt ist.
Ist nun mit gewissen Punkten $x_{i}$ $(i=n+1,\ldots ,r)$

(8.21) $\omega (x)\varphi (x)\geq 0,\quad x\in [a,b],$

so erhält man aus (8.20)

${\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)\leq {\frac {1}{(r+1)!}}\max _{\eta \in [a,b]}\left|f^{(r+1)}(\eta )\right|\int \limits _{a}^{b}[\omega (x)\varphi (x)]\,dx.$

Weiter gewinnt man mit (8.19)

${\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)\geq {\frac {1}{(r+1)!}}\min _{\eta \in [a,b]}\left|f^{(r+1)}(\eta )\right|\int \limits _{a}^{b}[\omega (x)\varphi (x)]\,dx.$

Der Zwischenwertsatz, angewandt auf die Funktion $f^{(r+1)}$, liefert somit für ein $\xi \in [a,b]$

${\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)={\frac {1}{(r+1)!}}f^{(r+1)}(\xi )\int \limits _{a}^{b}[\omega (x)\varphi (x)]\,dx,$

so dass die Substitution $x:=(b-a)t+a$ in diesem Fall zu der Formel (8.18) führt. Analog schließt man im Fall, dass „$\leq$“ statt „$\geq$“ in (8.21) vorliegt.

q.e.d.

Beispiel 8.16

Wir nutzen im Folgenden aus, dass nach Korollar 8.9 der Genauigkeitsgrad einer interpolatorischen Quadraturformel ${\mathcal {I}}_{n}mindestens<math>r:=nist.(1)Sei<math>f\in C^{1}[a,b]$ und ${\mathcal {I}}_{0}(f):=(b-a)f(a)$ die Rechteckregel aus (8.15). Aus (8.18) gewinnt man für ${\mathcal {I}}_{0}$ mit $r=n=0$, mit $x_{0}:=a$ bzw. $t_{0}:=0$ sowie mit

$\prod _{k=0}^{0}(t-t_{k})=t\geq 0,\quad t\in [0,1]$

und ${\hat {\gamma }}_{0}:=\int \limits _{0}^{1}t\,dt={\frac {1}{2}}$ die Fehlerdarstellung

${\mathcal {I}}(f)-{\mathcal {I}}_{0}(f)={\frac {(b-a)^{2}}{2}}f'(\xi ),$

wobei $\xi$ ein Punkt aus $[a,b]$ ist. Entsprechend erhält man für die Rechteckregel ${\mathcal {I}}_{0}(f):=(b-a)f(b)$ mit $r=n=0,x_{0}:=b$ bzw. $t_{0}:=1$ sowie mit

$\prod _{k=0}^{0}(t-t_{k})=t-1\leq 0,\quad t\in [0,1]$

und ${\hat {\gamma }}_{0}:=\int \limits _{0}^{1}(t-1)\,dt=-{\frac {1}{2}}$ die Fehlerdarstellung

${\mathcal {I}}(f)-{\mathcal {I}}_{0}(f)=-{\frac {(b-a)^{2}}{2}}f'(\xi ).$

(2) Im Fall der Trapezregel

${\mathcal {I}}_{1}(f):={\frac {b-a}{2}}[f(a)+f(b)]$

gilt für $f\in C^{2}[a,b]$ mit einem $\xi \in [a,b]$ die Fehlerdarstellung

${\mathcal {I}}(f)-{\mathcal {I}}_{1}(f)=-{\frac {(b-a)^{3}}{12}}f''(\xi ).$

Denn mit $r=n=1,x_{0}:=a,x_{1}:=b$ bzw. $t_{0}:=0,t_{1}:=1$ hat man

$\prod _{k=0}^{1}(t-t_{k})=t(t-1)\leq 0,\quad t\in [0,1]$

sowie

${\hat {\gamma }}_{1}:=\int \limits _{0}^{1}t(t-1)\,dt={\frac {1}{3}}-{\frac {1}{2}}=-{\frac {1}{6}}.$

Der Genauigkeitsgrad einer interpolatorischen Quadraturformel ${\mathcal {I}}_{n}$ ist mindestens $r:=n$. Für gerade $n$ hat man im Fall der abgeschlossenen Newton-Cotes-Formeln sogar das folgende Resultat (für den Beweis siehe Plato, S. 103):

Satz 8.17

Die abgeschlossene Newton-Cotes-Formel ${\mathcal {I}}_{n}$ besitzt für gerades $n\geq 2$ den (exakten) Genauigkeitsgrad $r:=n+1$.

Letzteres Ergebnis können wir z. B. für die Fehlerdarstellung der Simpson-Regel verwenden.

Beispiel 8.18

Es sei $f\in C^{4}[a,b]$. Dann hat man für $n=2$ und $r=3$ mit $x_{0}:=a,x_{1}:=(a+b)/2,x_{2}:=b$ bzw. $t_{0}:=0,t_{1}:=1/2,t_{2}:=1$ und mit dem gewählten Punkt $t_{3}:=1/2$

$\prod _{k=0}^{3}(t-t_{k})=t\left(t-{\frac {1}{2}}\right)^{2}(t-1)\leq 0,\quad t\in [0,1]$

sowie

${\hat {\gamma }}_{3}:=\int \limits _{0}^{1}\left[t\left(t-{\frac {1}{2}}\right)^{2}(t-1)\right]\,dt=-{\frac {1}{120}}.$

Also ergibt sich für die Simpson-Regel

${\mathcal {I}}_{2}(f):={\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right]$

mit einem $\xi \in [a,b]$ der Quadraturfehler

${\mathcal {I}}(f)-{\mathcal {I}}_{2}(f)=-{\frac {(b-a)^{5}}{4!\cdot 120}}f^{(4)}(\xi )=-{\frac {(b-a)^{5}}{2880}}f^{(4)}(\xi ).$

8.3 Summierte abgeschlossene Newton-Cotes-Formel

Wie bereits in Abschnitt 8.2.2 erläutert wurde, garantiert eine Erhöhung von $n$ keineswegs, dass die Newton-Cotes-Formeln Näherungswerte zunehmender Genauigkeit für ${\mathcal {I}}(f)$ liefern. Um Letzteres zu erreichen, müssen wir daher anders vorgehen. Und zwar teilen wir zunächst das Intervall $[a,b]$ mittels Stützstellen

$x_{k}:=a+kh\quad (k=0,1,\ldots ,N),\qquad h:={\frac {b-a}{N}}$

in $N$ gleiche Stücke auf, so dass sich insbesondere

$h=x_{k+1}-x_{k}$

für alle $k\in \{0,1,\ldots ,N-1\}$ ergibt. Dann nähern wir das Integral

${\mathcal {I}}(f)=\int \limits _{a}^{b}f(x)\,dx=\sum _{k=0}^{N-1}\int \limits _{x_{k}}^{x_{k+1}}f(x)\,dx$

durch

${\mathcal {J}}_{n}(f):=\sum _{k=0}^{N-1}\int \limits _{x_{k}}^{x_{k+1}}Q_{n}(x)\,dx$

an, wobei $Q_{n}\in \Pi _{n}$ das (eindeutige) Interpolationspolynom zu $n+1$ paarweise verschiedenen, in jedem Intervall $[x_{k},x_{k+1}]$ in gleichen Abständen gewählten Stützpunkten ist (vgl. Definition 8.8). Wir wählen also eine interpolatorische Quadraturformel und ersetzen jedes der Integrale $\int \limits _{x_{k}}^{x_{k+1}}f(x)\,dx$ durch den sich damit ergebenden Wert. Eine so gewonnene Quadraturformel bezeichnet man als summierte Quadraturformel. Wir wollen solche Formeln nun genauer betrachten, wobei wir uns hier auf die abgeschlossenen Newton-Cotes-Formeln zu deren Generierung beschränken wollen. Letztere Wahl legt die Stützpunkte in jedem Intervall $[x_{k},x_{k+1}]$ durch (8.11) fest, wobei dort $a:=x_{k}$ und $b:=x_{k+1}$ zu wählen ist.

Wir beginnen mit den beiden Rechteckregeln aus (8.15). Für diese erhält man

$\int \limits _{x_{k}}^{x_{k+1}}Q_{0}(x)\,dx=hf(x_{k})$ bzw. $\int \limits _{x_{k}}^{x_{k+1}}Q_{0}(x)\,dx=hf(x_{k+1})$

so dass Summation über $k$ die folgenden summierten Rechteckregeln liefert:

${\mathcal {J}}_{0}(h):=h\sum _{k=0}^{N-1}f(x_{k}),\quad {\hat {\mathcal {J}}}_{0}(h):=h\sum _{k=0}^{N-1}f(x_{k+1}).$

Für diese gelten die nachstehenden Fehlerabschätzungen.

Satz 8.19

Es sei $f\in C^{1}[a,b]$. Dann gibt es $\xi ,{\hat {\xi }}\in [a,b]$, so dass gilt:

(8.22) ${\mathcal {I}}(f)-{\mathcal {J}}_{0}(h)={\frac {b-a}{2}}hf'(\xi ),\quad {\mathcal {I}}(f)-{\hat {\mathcal {J}}}_{0}(h)={\frac {b-a}{2}}hf'({\hat {\xi }}).$

Beweis.

Aus Beispiel 8.16 (1) ergibt sich für $k=0,1,\ldots ,N-1$ die Existenz eines $\xi _{k}\in [a,b]$ mit

$\int \limits _{x_{k}}^{x_{k+1}}f(x)\,dx-hf(x_{k})={\frac {h^{2}}{2}}f'(\xi _{k}).$

Summation über $k$ führt auf

${\mathcal {I}}(f)-{\mathcal {J}}_{0}(h)={\frac {h^{2}}{2}}\sum _{k=0}^{N-1}f'(\xi _{k})={\frac {b-a}{2}}h{\frac {1}{N}}\sum _{k=0}^{N-1}f'(\xi _{k}).$

Aufgrund von

$N\min _{x\in [a,b]}f'(x)\leq \sum _{k=0}^{N-1}f'(\xi _{k})\leq N\max _{x\in [a,b]}f'(x)$

bzw.

$\min _{x\in [a,b]}f'(x)\leq {\frac {1}{N}}\sum _{k=0}^{N-1}f'(\xi _{k})\leq \max _{x\in [a,b]}f'(x)$

existiert nach dem Zwischenwertsatz ein $\xi \in [a,b]$ mit

$f'(\xi )={\frac {1}{N}}\sum _{k=0}^{N-1}f'(\xi _{k}),$

so dass die erste Fehlerdarstellung in (8.22) folgt. Die zweite zeigt man analog.

q.e.d.

Im Fall der Trapezregel (8.13) hat man

$\int \limits _{x_{k}}^{x_{k+1}}Q_{1}(x)\,dx={\frac {h}{2}}[f(x_{k})+f(x_{k+1})].$

Summation über $k$ führt auf die summierte Trapezregel

${\mathcal {J}}_{1}(h):={\frac {h}{2}}\left(f(a)+2\sum _{k=1}^{N-1}f(x_{k})+f(b)\right)$

mit der im folgenden Satz angegebenen Fehlerdarstellung.

Satz 8.20

Es sei $f\in C^{2}[a,b]$. Dann existiert ein $\xi \in [a,b]$ mit

${\mathcal {I}}(f)-{\mathcal {J}}_{1}(h)=-{\frac {b-a}{12}}h^{2}f''(\xi ).$

Beweis.

Der Beweis verläuft analog zu dem von Satz 8.19. Nach Beispiel 8.16 (2) gibt es für $k=0,1,\ldots ,N-1$ ein $\xi _{k}\in [a,b]$ mit

$\int \limits _{x_{k}}^{x_{k+1}}f(x)\,dx-{\frac {h}{2}}[f(x_{k})+f(x_{k+1})]=-{\frac {h^{3}}{12}}f''(\xi _{k}).$

Summation über $k$ liefert mit einem $\xi \in [a,b]$

${\mathcal {I}}(f)-{\mathcal {J}}_{1}(h)=-{\frac {b-a}{12}}h^{2}{\frac {1}{N}}\sum _{k=0}^{N-1}f''(\xi _{k})=-{\frac {b-a}{12}}h^{2}f''(\xi ),$

wobei die Existenz eines solchen $\xi$ aus der Anwendung des Zwischenwertsatzes auf $f''$ geschlossen werden kann.

q.e.d.

Schließlich betrachten wir noch die summierte Simpson-Regel, wobei wir die Darstellung $x_{k}:=a+kh$ mit $h:=(b-a)/N$ für jedes $k\geq 0$ verwenden, so dass insbesondere

${\frac {x_{k}+x_{k+1}}{2}}=x_{k}+{\frac {1}{2}}(x_{k+1}-x_{k})=a+kh+{\frac {1}{2}}h=x_{k+1/2},\quad k=0,1,\ldots ,N-1$

folgt. Die Simpson-Regel, angewandt auf das Intervall $[x_{k},x_{k+1}]$, lässt sich somit in der Form

$\int \limits _{x_{k}}^{x_{k+1}}Q_{2}(x)={\frac {h}{6}}\left[f(x_{k})+4f(x_{k+1/2})+f(x_{k+1})\right]$

schreiben. Summation über $k$ führt auf die summierte Simpson-Regel

${\mathcal {J}}_{2}(h):={\frac {h}{6}}\left(f(a)+4\sum _{k=0}^{N-1}f(x_{k+1/2})+2\sum _{k=1}^{N-1}f(x_{k})+f(b)\right)$

Für diese hat man die im folgenden Satz angegebene Fehlerdarstellung.

Satz 8.21

Es sei $f\in C^{4}[a,b]$. Dann existiert ein $\xi \in [a,b]$ mit

${\mathcal {I}}(f)-{\mathcal {J}}_{2}(h)=-{\frac {b-a}{2880}}h^{4}f^{(4)}(\xi ).$

Beweis.

Der Beweis verläuft wiederum analog zu dem von Satz 8.19. Nach Beispiel 8.18 gibt es für $k=0,1,\ldots ,N-1$ ein $\xi _{k}\in [a,b]$ mit

$\int \limits _{x_{k}}^{x_{k+1}}f(x)\,dx-{\frac {h}{6}}\left[f(x_{k})+4f(x_{k+1/2})+f(x_{k+1})\right]=-{\frac {h^{5}}{2880}}f^{(4)}(\xi _{k}).$

Summation über $k$ liefert mit einem $\xi \in [a,b]$

${\mathcal {I}}(f)-{\mathcal {J}}_{2}(h)=-{\frac {b-a}{2880}}h^{4}{\frac {1}{N}}\sum _{k=0}^{N-1}f^{(4)}(\xi _{k})=-{\frac {b-a}{2880}}h^{4}f^{(4)}(\xi ),$

wobei die Existenz von $\xi$ aus der Anwendung des Zwischenwertsatzes auf $f^{(4)}$ folgt.

q.e.d.

Zur Auswertung der summierten Rechteckregeln müssen $N$, für die der summierten Trapezregel $N+1$ und für die der summierte Simpson-Regel $2N+1$ Funktionswerte bestimmt werden. Der Rechenaufwand bei Verwendung der summierten Simpson-Regel ist damit etwa doppelt so hoch wie der bei Verwendung einer der drei anderen Regeln. Dennoch ist die summierte Simpson-Regel diesen für hinreichend glatte Funktionen wegen der höheren Fehlerordnung in h vorzuziehen. Denn der Quadraturfehler verhält sich bei ihr wie ${\mathcal {O}}(h^{4})$, während er bei den summierten Rechteckregeln und der summierten Trapezregel proportional zu $h$ bzw. $h^{2}$ abnimmt.

Da man die in der jeweiligen Fehlerformel vorkommende Ableitung durch das Maximum des Betrages dieser Ableitung bezüglich aller $x\in [a,b]$ nach oben abschätzen kann, implizieren die angegebenen Fehlerdarstellungen insbesondere, dass die hier angegebenen summierten Quadraturformeln für $h\to 0$ gegen den exakten Wert des Integrals ${\mathcal {I}}(f)$ konvergieren, wobei mit „$h\to 0$“ hier „$h_{k}=(b-a)/N_{k}$ mit $N_{k}\in \mathbb {N}$ und $N_{k}\to \infty$“ gemeint ist.

Wir greifen abschließend nochmals Beispiel 8.13 auf.

Beispiel 8.22

Es seien wieder $f(x):=1/(1+x^{2}),a=0$ und $b=1$, so dass ein Näherungswert für das Integral

${\mathcal {I}}(f)=\int \limits _{0}^{1}{\frac {1}{1+x^{2}}}\,dx$

gesucht ist. Weiter wählen wir $N=3$ und somit $h=1/3$. Der exakte Wert des Integrals lautet $\arctan(1)=0.785\,398\,1$. Mit der summierten Simpson-Regel ergibt sich der Wert

${\mathcal {J}}_{2}(h)={\frac {1}{18}}\left[f(0)+4\left(f(1/6)+f(3/6)+f(5/6)\right)+2\left(f(1/3)+f(2/3)\right)+f(1)\right]$

$={\frac {1}{18}}\left[1+4\left({\frac {36}{37}}+{\frac {36}{45}}+{\frac {36}{61}}\right)+2\left({\frac {9}{10}}+{\frac {9}{13}}\right)+{\frac {1}{2}}\right]=0.785\,397\,94.$

8.4 Extrapolationsverfahren

8.4.1 Einführung

Für die summierte Trapezregel ${\mathcal {J}}_{1}(h)$ gibt der folgende Satz eine asymptotische Entwicklung nach Potenzen von $h^{2}$ an, welche dazu genutzt werden soll, aus einer endlichen Zahl von Auswertungen der summierten Trapezregel eine im Hinblick auf diese Werte genauere Näherung des Integrals ${\mathcal {I}}(f)$ zu berechnen. (Der Satz ist z. B. bei Plato bewiesen.)

Satz 8.23

Für ein $r\geq 0$ sei $f\in C^{2r+2}[a,b]$. Die summierte Trapezregel

${\mathcal {J}}_{1}(h):={\frac {h}{2}}\left(f(a)+2\sum _{k=1}^{N-1}f(x_{k})+f(b)\right)$

mit $h:=(b-a)/N$ für ein $N\in \mathbb {N}$ besitzt die asymptotische Entwicklung

(8.23) ${\mathcal {J}}_{1}(h)=\alpha _{0}+\alpha _{1}h^{2}+\alpha _{2}h^{2\cdot 2}+\ldots +\alpha _{r}h^{2r}+{\mathcal {O}}(h^{2r+2})$ für $h\to 0$

mit $\alpha _{0}:=I(f)$ und gewissen Koeffizienten $\alpha _{i}\in \mathbb {R}$ $(i=1,\ldots ,r)$.

Für periodische Funktionen mit Periode $b-a$ kann man sogar zeigen, dass $\alpha _{i}=0$ $(i=1,\ldots ,r)$ gilt. In einem solchen Fall kann mit dem in diesem Abschnitt beschriebenen Verfahren keine Verbesserung erzielt werden.

Man beachte, dass man ${\mathcal {J}}_{1}(h)$ nur für $h>0$ mit $h:=(b-a)/N$ für eine natürliche Zahl $N$ auswerten kann. Aufgrund von (8.23) (wie auch wegen Satz 8.20) gilt ferner

$\lim _{h\to 0}{\mathcal {J}}_{1}(h)=I(f),$

wobei wir mit „$h\to 0$“ hier „$h_{k}=(b-a)/N_{k}$ mit $N_{k}\in \mathbb {N}$ und $N_{k}\to \infty$“ meinen. Die Entwicklung (8.23) soll nun numerisch dazu ausgenutzt werden, von einer endlichen Zahl berechneter Werte ${\mathcal {J}}_{1}(h_{k}),k=0,1,\ldots ,n$ mit $0<h_{n}<h_{n-1}<\ldots <h_{0}$ auf einen noch genaueren Wert von $I(f)$ als ${\mathcal {J}}_{1}(h_{n})$ zu schließen.

Wir gehen dabei allgemeiner von einer beliebigen Funktion $T(h)$ mit $h>0$ aus, die mit gewissen Koeffizienten $\alpha _{i}\in \mathbb {R}$ $(i=0,1,\ldots ,r)$ und einer Zahl $\gamma >0$ die asymptotische Entwicklung der Ordnung $r$

(8.24) $T(h)=\alpha _{0}+\alpha _{1}h^{\gamma }+\alpha _{2}h^{2\gamma }+\ldots +\alpha _{r}h^{r\gamma }+{\mathcal {O}}(h^{(r+1)\gamma })$ für $h\to 0$

besitzt und für die der Wert

$\lim _{h\to 0}T(h)=\alpha _{0}$

gesucht ist. Typischerweise steht $T(h)$ für ein numerisches Verfahren, das für einen gewählten Diskretisierungsparameter $h>0$ einen Näherungswert für die gesuchte Größe $\alpha _{0}$ liefert. Es sei also angenommen, dass $T(h)$ zumindest für gewisse $h>0$ berechnet werden kann, wie dies z. B. im Fall der Tangententrapezregel für $h:=(b-a)/N$ mit $N\in \mathbb {N}$ der Fall ist.

Wegen (8.24) hat man zunächst für $h>0$ nur die Genauigkeit

$T(h)-\alpha _{0}={\mathcal {O}}(h^{\gamma }).$

Es soll nun ein Verfahren vorgestellt werden, welches ohne großen Mehraufwand aus endlich vielen, bereits berechneten Werten $T(h_{k}),k=0,1,\ldots ,n$ mit $0<h_{n}<h_{n-1}<\ldots <h_{0}$ einen genaueren Wert für die gesuchte Größe $\alpha _{0}$ erzeugt. Setzt man $T(0):=\alpha _{0}$, so extrapoliert dieses Verfahren also $T$ auf den Wert $h=0$ hin, so dass man auch von einem Extrapolationsverfahren spricht. Da die Koeffizienten $\alpha _{i}$ in (8.24) oft nicht explizit bekannt sind oder nur unter einigem Aufwand zu berechnen sind, geht man dabei folgendermaßen vor:

• man vernachlässigt den Restterm ${\mathcal {O}}(h^{(r+1)\gamma })$ in (8.24) und geht davon aus, dass sich $T(h)$ ungefähr wie ein Polynom in $h$ verhält,
• man ersetzt das resultierende (i. A. nicht explizit bekannte) Polynom durch das Interpolationspolynom $P_{0,\ldots ,n}\in \Pi _{n}$ zu den Stützpunkten $(h_{k}^{\gamma },T(h_{k})),k=0,1,\ldots ,n$ (schreibt man $T(h^{\gamma })$ statt $T(h)$, so sind dies mit $z_{k}:=h_{k}^{\gamma }$ die Punkte $(z_{k},T(z_{k}))$) und
• man verwendet den Wert $P_{0,\ldots ,n}(0)$ als Näherung für den unbekannten Wert $\alpha _{0}$.

Im Zusammenhang mit der summierten Trapezregel wird diese Vorgehensweise als Romberg-Verfahren bezeichnet.

8.4.2 Das Verfahren

Wir gehen nun von der asymptotischen Entwicklung (8.24) von $T(h)$ aus und es sei $P_{0,\ldots ,n}\in \Pi _{n}$ das Interpolationspolynom zu den Stützpunkten

(8.25) $(h_{k}^{\gamma },T(h_{k}))\quad k=0,1,\ldots ,n.$

Da dieses nur an einer Stelle, der Stelle 0, ausgewertet werden soll, bietet sich das Neville-Schema zur Verwendung an, wobei hier $P_{j,\ldots ,j+m}\in \Pi _{m}$ das Interpolationspolynom mit

(8.26) $P_{j,\ldots ,j+m}(h_{k}^{\gamma })=T(h_{k}),\quad k=j,\ldots ,j+m$

bezeichnet. Wir setzen dazu

(8.27) $T_{j,\ldots ,j+m}:=P_{j,\ldots ,j+m}(0).$

Satz 6.5 liefert damit

(8.28) $T_{j}=P_{j}(0)=T(h_{j})$

sowie

$T_{j,\ldots ,j+m}={\frac {-h_{j}^{\gamma }T_{j+1,\ldots ,j+m}+h_{j+m}^{\gamma }T_{j,\ldots ,j+m-1}}{h_{j+m}^{\gamma }-h_{j}^{\gamma }}}=T_{j+1,\ldots ,j+m}-h_{j+m}^{\gamma }{\frac {T_{j+1,\ldots ,j+m}-T_{j,\ldots ,j+m-1}}{h_{j+m}^{\gamma }-h_{j}^{\gamma }}}$

(8.29) $=T_{j+1,\ldots ,j+m}+{\frac {T_{j+1,\ldots ,j+m}-T_{j,\ldots ,j+m-1}}{\left({\frac {h_{j}}{h_{j+m}}}\right)^{\gamma }-1}}\quad (j,m\in \mathbb {N} _{0}).$

Das Schema von Neville geht damit in das folgende Extrapolationstableau über, welches zeilenweise aufgebaut wird:

${\begin{matrix}T_{0}=T(h_{0})&&&&&&&\\&\searrow &&&&&&\\T_{1}=T(h_{1})&\to &T_{01}&&&&&\\&\searrow &&\searrow &&&&\\T_{2}=T(h_{2})&\to &T_{12}&\to &T_{012}&&&\\&\searrow &&\searrow &&\searrow &&\\T_{3}=T(h_{3})&\to &T_{23}&\to &T_{123}&\to &T_{0123}&\\\vdots &&\vdots &&\vdots &&\vdots &\ddots \end{matrix}}$

Beispiel 8.24

Für die summierte Trapezregel ${\mathcal {J}}_{1}(h)$ gilt gemäß Satz 8.23 eine Entwicklung der Form (8.24) mit $\gamma =2$. Für die Schrittweiten

$h_{0}:=b-a,\quad h_{1}:={\frac {h_{0}}{2}}={\frac {b-a}{2}}$

erhält man die Werte

$T_{0}={\mathcal {J}}_{1}(h_{0})={\frac {b-a}{2}}[f(a)+f(b)],$

$T_{1}={\mathcal {J}}_{1}(h_{0})={\frac {b-a}{2}}\left[f(a)+2f\left({\frac {a+b}{2}}\right)+f(b)\right]$

und damit

$T_{01}=T_{1}+{\frac {T_{1}-T_{0}}{\left({\frac {h_{0}}{h_{1}}}\right)^{2}-1}}={\frac {b-a}{4}}\left[f(a)+2f\left({\frac {a+b}{2}}\right)+f(b)\right]+{\frac {b-a}{4}}\cdot {\frac {2f\left({\frac {a+b}{2}}\right)-f(a)-f(b)}{4-1}}$

$={\frac {b-a}{4}}\left[{\frac {2}{3}}f(a)+{\frac {8}{3}}f\left({\frac {a+b}{2}}\right)+{\frac {2}{3}}f(b)\right]={\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right].$

Der aus den beiden Auswertungen $T_{0}$ und $T_{1}$ der summierten Trapezregel ermittelte Wert $T_{01}$ entspricht somit dem der Simpson-Regel für ${\mathcal {I}}(f)$.

Im folgenden Satz wird die Größenordnung des Fehlers $T_{j,\ldots ,j+m}-\alpha _{0}$ angegeben. Diese Fehlerbetrachtung macht deutlich, dass sich die Anwendung des hier untersuchten Extrapolationsverfahrens lohnt. Als Hilfsmittel verwenden wir das nachstehende Lemma.

Lemma 8.25

Es seien $L_{k}\in \Pi _{m}$ $(k=0,1,\ldots ,m)$ die Lagrangeschen Basispolynome zu Stützstellen $x_{k}$ $(k=0,1,\ldots ,m)$ mit $x_{k}\neq x_{i}$ $(k\neq i)$. Dann gilt

(8.30) $\sum _{k=0}^{m}L_{k}(0)x_{k}^{j}={\begin{cases}1&f{\ddot {u}}r\ j=0,\\0&f{\ddot {u}}r\ 1\leq j\leq m,\\(-1)^{m}x_{0}x_{1}\cdots x_{m}&f{\ddot {u}}r\ j=m+1.\end{cases}}$

Beweis.

Für $0\leq j\leq m$ ist offenbar $p(x):=x^{j}$ das Interpolationspolynom zu den Punkten $(x_{k},x_{k}^{j}),k=0,1,\ldots ,m$ und daher gemäß (6.7)

$p(x)=x^{j}=\sum _{k=0}^{m}x_{k}^{j}L_{k}(x).$

Setzen wir $x=0$, so folgt die Behauptung für die ersten beiden Fälle in (8.30). Für den Fall $j=m+1$ betrachten wir das Polynom

$q(x):=x^{m+1}-\sum _{k=0}^{m}x_{k}^{m+1}L_{k}(x),$

welches wegen $L_{k}\in \Pi _{m}$ den Grad $m+1$, den führenden Koeffizienten 1 und die Nullstellen $x_{i}$ $(i=0,1,\ldots ,m)$ hat, so dass insbesondere

$q(x)=(x-x_{0})(x-x_{1})\cdots (x-x_{m})$

gilt. Speziell hat man somit

$\sum _{k=0}^{m}L_{k}(0)x_{k}^{m+1}=-q(0)=(-1)^{m}x_{0}x_{1}\cdots x_{m}.$

Satz 8.26

Es sei $T(h)$ mit $h>0$ eine Funktion mit der asymptotischen Entwicklung (8.24) für ein $\gamma >0$ und $r\in \mathbb {N}$. Weiter sei $(h_{k})$ eine Folge von Schrittweiten, so dass mit einer Startschrittweite $h_{0}>0$ gilt:

(8.31) $h_{k}:=h_{0}/n_{k}\quad (k\in \mathbb {N} _{0})$ mit $1=n_{0}\leq n_{1}\leq n_{2}\leq \ldots .$

Schließlich sei $P_{j,\ldots ,j+m}\in \Pi _{m}$ das Interpolationspolynom mit (8.26) und $T_{j,\ldots ,j+m}$ wie in (8.27). Dann genügt der Fehler $T_{j,\ldots ,j+m}-\alpha _{0}$ für $0\leq m\leq r-1$ der asymptotischen Entwicklung

(8.32) $T_{j,\ldots ,j+m}-\alpha _{0}=(-1)^{m}{\frac {\alpha _{m+1}}{n_{j}^{\gamma }\cdots n_{j+m}^{\gamma }}}h_{0}^{(m+1)\gamma }+{\mathcal {O}}(h_{0}^{(m+2)\gamma }).$

Beweis.

Da sich die Indizes in (8.32) auf eine Numerierung der Stützpunkte beziehen und wir den $j$-ten als 0-ten bezeichnen können, können wir o. B. d. A. $j=0$ annehmen. Gemäß der Lagrangeschen Darstellung des Interpolationspolynoms $P_{0,\ldots ,m}$ gilt dann

$P_{0,\ldots ,m}(h^{\gamma })=\sum _{k=0}^{m}T(h_{k})L_{k}(h^{\gamma })=\sum _{k=0}^{m}T(h_{k})\left[\prod _{j=0 \atop j\neq k}^{m}{\frac {h-h_{j}^{\gamma }}{h_{k}^{\gamma }-h_{j}^{\gamma }}}\right],\quad h\in \mathbb {R}$

und somit

(8.33) $T_{0,\ldots ,m}=P_{0,\ldots ,m}(0)=\sum _{k=0}^{m}c_{m,k}T(h_{k})$

für

$c_{m,k}:=L_{k}(0)=\prod _{j=0 \atop j\neq k}^{m}{\frac {h_{j}^{\gamma }}{h_{j}^{\gamma }-h_{k}^{\gamma }}}.$

Nun folgt wegen $m\leq r-1$ aus (8.24)

(8.34) $T(h_{k})=\sum _{s=0}^{m+1}\alpha _{s}h_{k}^{s\gamma }+{\mathcal {O}}(h_{k}^{(m+2)\gamma }).$

Des Weiteren schließt man mit Lemma 8.25

(8.35) $\sum _{k=0}^{m}c_{m,k}h_{k}^{s\gamma }=\sum _{k=0}^{m}L_{k}(0)h_{k}^{s\gamma }={\begin{cases}1&{\mbox{für }}s=0,\\0&{\mbox{für }}1\leq s\leq m,\\(-1)^{m}h_{0}^{\gamma }h_{1}^{\gamma }\cdots h_{m}^{\gamma }&{\mbox{für }}s=m+1.\end{cases}}$

Setzt man die beiden Beziehungen (8.34) und (8.35) in (8.33) ein, so bekommt man schließlich, da die $c_{m,k}$ von $h$ unabhängig sind,

$T_{0,\ldots ,m}=\sum _{k=0}^{m}c_{m,k}\left[\sum _{s=0}^{m+1}\alpha _{s}h_{k}^{s\gamma }+{\mathcal {O}}(h_{k}^{(m+2)\gamma })\right]=\sum _{s=0}^{m+1}\alpha _{s}\left[\sum _{k=0}^{m}c_{m,k}h_{k}^{s\gamma }\right]+\sum _{k=0}^{m}c_{m,k}{\mathcal {O}}(h_{k}^{(m+2)\gamma })$

$=\alpha _{0}+(-1)^{m}\alpha _{m+1}h_{0}^{\gamma }h_{1}^{\gamma }\cdots h_{m}^{\gamma }+{\mathcal {O}}(h_{0}^{(m+2)\gamma }).$

q.e.d.

Der Satz besagt, dass man beim Übergang von $m$ zu $m+1$, d. h. bei Erhöhung der Spaltenzahl in dem Extrapolationstableau um 1, im Prinzip die Ordnung $\gamma$ gewinnt. Diese Sichtweise ist allerdings zu optimistisch, da die Restterme der asymptotischen Entwicklung, die sich hinter ${\mathcal {O}}(h_{0}^{(m+2)\gamma })$ verbergen, nicht bekannt sind und groß werden können.

Es bietet sich also der folgende Algorithmus an:

Algorithmus 10 (Extrapolationsverfahren)

(0) Wähle $h_{0}>0$, eine Folge $h_{k},k=1,2,\ldots$ wie in (8.31) und ein $\varepsilon >0$. Setze $j:=0$.

(1) Berechne $T_{j}:=T(h_{j})$.

(2) Berechne $T_{k,\ldots ,j}$ für $k=j-1,j-2,\ldots ,0$ nach der Formel

$T_{k,\ldots ,j}=T_{k+1,\ldots ,j}+{\frac {T_{k+1,\ldots ,j}-T_{k,\ldots ,j-1}}{\left({\frac {h_{k}}{h_{j}}}\right)^{\gamma }-1}}$

(3) Falls „der Aufwand zu groß wird“ oder

${\frac {|T_{0,\ldots ,j}-T_{0,\ldots ,j-1}|}{|T_{0,\ldots ,j}|}}$

gilt, breche ab. ($T_{0,\ldots ,j}$ ist Näherungswert für $\alpha _{0}$.)

(4) Setze $j:=j+1$ und gehe nach (1).

Man bricht das Extrapolationsverfahren also ab, wenn der Aufwand zur Erzeugung einer neuen Zeile im Extrapolationsschema, den man meistens, wie z. B. für das summierte Trapezverfahren, genau angeben kann, zu groß wird oder die relative Abweichung zweier aufeinanderfolgender Diagonalelemente klein genug wird. In der Praxis ist es jedoch auch möglich, dass aufgrund von Rundungsfehlern Divergenz eintritt, so dass auf früher berechnete Werte im Schema zurückgegriffen werden muss.

Häufig angewandte Schrittweitenfolgen $(h_{k})$ für (8.31) in diesem Zusammenhang sind die Romberg-Folge

(8.36) $n_{k}:=2^{k},\quad h_{k}={\frac {h_{k-1}}{2}}={\frac {h_{0}}{2^{k}}}\quad (k\in \mathbb {N} ),$

die durch

$n_{1}:=2,\quad n_{2}:=3,\quad n_{3}:=4,\quad n_{j}:=2n_{j-2}\quad (j\geq 4)$

definierte Bulirsch-Folge

$h_{1}={\frac {h_{0}}{2}},\quad h_{2}={\frac {h_{0}}{3}},\quad h_{3}={\frac {h_{0}}{4}},\quad h_{4}={\frac {h_{0}}{6}},\quad h_{5}={\frac {h_{0}}{8}},\quad h_{6}={\frac {h_{0}}{12}},\quad h_{7}={\frac {h_{0}}{16}},\quad \ldots$

und die harmonische Folge

$n_{k-1}:=k\quad (k=1,2,\ldots ),\quad h_{k}={\frac {h_{0}}{k+1}}\quad (k\in \mathbb {N} ).$

Insbesondere erhält man für die Romberg-Folge $(j:=0)$:

Korollar 8.27

Unter den Voraussetzungen von Satz 8.26 gilt für die Romberg-Folge (8.36) mit $h_{0}:=(b-a)/N$ und $0\leq m\leq r-1$

$T_{0,\ldots ,m}-\alpha _{0}=\left[{\frac {(-1)^{m}}{2^{\gamma m(m+1)/2}}}\alpha _{m+1}\right]h_{0}^{(m+1)\gamma }+{\mathcal {O}}(h_{0}^{(m+2)\gamma }).$

Beweis.

Im Fall (8.36) hat man mit $n_{0}=1$

$n_{0}^{\gamma }\cdots n_{m}^{\gamma }=2^{0\gamma }\cdot 2^{1\gamma }\cdot 2^{2\gamma }\cdots 2^{m\gamma }=2^{\gamma \sum _{k=0}^{m}k}=2^{\gamma m(m+1)/2},$

so dass die Behauptung unmittelbar aus Satz 8.26 folgt.

q.e.d.

Wir betrachten nun nochmals die summierte Trapezregel als Beispiel.

Beispiel 8.28

(1) Korollar 8.27 wollen wir auf die summierte Trapezregel mit $m:=2$ (und wegen der Forderung $m\leq r-1$ für $r=3$) anwenden. Nach Satz 8.23 ist dann $\gamma =2$. Weiter sei $f\in C^{8}[a,b]$ vorausgesetzt. Korollar 8.27 liefert mit diesen Setzungen

(8.37) $T_{012}-{\mathcal {I}}(f)=\left[{\frac {(-1)^{2}}{2^{2\cdot 2(2+1)/2}}}\alpha _{2+1}\right]h_{0}^{(2+1)2}+{\mathcal {O}}(h_{0}^{(m+2)2})={\frac {\alpha _{3}}{64}}h_{0}^{6}+{\mathcal {O}}(h_{0}^{8}),$

wobei $T_{012}$ mit dem Neville-Schema

${\begin{matrix}T_{0}={\mathcal {J}}_{1}(h_{0})&&&&\\&\searrow &&&\\T_{1}={\mathcal {J}}_{1}(h_{0}/2)&\rightarrow &T_{01}&&\\&\searrow &&\searrow &\\T_{2}={\mathcal {J}}_{1}(h_{0}/4)&\rightarrow &T_{12}&\rightarrow &T_{012}\end{matrix}}$

berechnet wird. Man beachte dabei, dass man bei der Berechnung von $T_{j+1}$ $(j\in \mathbb {N} _{0})$ den zuvor ermittelten Wert $T_{j}$ verwenden kann und nur zusätzlich Funktionsauswertungen für die Mittelpunkte der sich aus der zu $T_{j}$ gehörenden Zerlegung von $[a,b]$ ergebenden Intervalle benötigt. Somit verlangt die Berechnung von $T_{0},T_{1}$ und $T_{2}$ genauso viele Funktionsauswertungen wie die direkte Berechnung von $T_{2}$. Für letzteren Wert alleine hat man aber im Vergleich zu (8.37) gemäß Satz 8.20 für $f\in C^{2}[a,b]$ mit einem $\xi \in [a,b]$ einen Fehler der Größe ${\mathcal {O}}(h_{0}^{2})$:

$T_{2}-{\mathcal {I}}(f)={\mathcal {J}}_{1}(h_{0}/4)-{\mathcal {I}}(f)={\frac {b-a}{12\cdot 4^{2}}}h_{0}^{2}f''(\xi )={\frac {(b-a)f''(\xi )/3}{64}}h_{0}^{2}.$

(2) (Bader) Es soll das Integral

$\int \limits _{0}^{\pi /2}5{\frac {e^{2x}\cos(x)}{e^{\pi }-2}}\,dx=1$

näherungsweise mit der summierten Trapezregel und dem Extrapolationsverfahren mit der Romberg-Folge (8.36) und $h_{0}:=\pi /2$ berechnet werden. Es ergibt sich bei 12-stelliger Rechnung das folgende Extrapolationstableau:

${\begin{matrix}0.185\,755\,068\,924&&&\\0.724\,727\,335\,089&0.904\,384\,757\,145&&\\0.925\,565\,035\,158&0.992\,510\,935\,182&0.999\,386\,013\,717&\\0.981\,021\,630\,069&0.999\,507\,161\,706&0.999\,973\,576\,808&0.999\,998\,776\,222\\0.995\,232\,017\,388&0.999\,968\,813\,161&0.999\,999\,589\,925&1.000\,000\,002\,83\\0.999\,806\,537\,974&0.999\,998\,044\,836&0.999\,999\,993\,614&1.000\,000\,000\,02\end{matrix}}$

Der (in der Tabelle nicht mehr einfügbare) Wert des Diagonalelementes in der 5. Spalte beträgt $1.000\,000\,084\,6$. Er ist offenbar ungenauer als die beiden untersten Werte in der 4. Spalte, wobei für den untersten allerdings auch die summierte Trapezregel einmal mehr ausgewertet werden musste. (Man kann auch zeigen, was für die erste Spalte z. B. aus Satz 8.20 folgt, dass die Werte in jeder einzelnen Spalte des Extrapolationsschemas, d. h. für konstantes $m$ und $j\to \infty$, gegen den gesuchten Wert $\alpha _{0}$ konvergieren.) Für eine Diskussion über ein geeignetes Abbruchkriterium verweisen wir auf Deuflhard/Hohmann.

8.5 Gaußsche Quadraturformeln

8.5.1 Einleitung

In diesem Abschnitt betrachten wir Quadraturformeln für gewichtete Integrale des Typs

(8.38) ${\mathcal {I}}(f):=\int \limits _{a}^{b}w(x)f(x)\,dx,\quad f\in C[a,b],$

wobei das Intervall $[a,b]$ hier halbunendlich oder unendlich, d. h. $a:=-\infty$ und/oder $b:=\infty$ sein darf und $w:(a,b)\to \mathbb {R}$ eine Gewichtsfunktion mit den folgenden Eigenschaften ist:

• $w(x)>0,\quad x\in (a,b),$
• es existieren die Momente

$\mu _{k}:=\int \limits _{a}^{b}w(x)x^{k}\,dx,\quad k=0,1,2,\ldots .$

Häufig in diesem Zusammenhang auftretende Gewichtsfunktionen sind in der folgenden Tabelle wiedergegeben, wobei auch der zuvor betrachtete Fall $w\equiv 1$ von Interesse ist:

Fehler beim Parsen (Konvertierungsfehler. Der Server („https://wikimedia.org/api/rest_“) hat berichtet: „Cannot get mml. TeX parse error: Bracket argument to \\ must be a dimension“): {\displaystyle {\begin{array}{|c|c|}\hline {\text{Intervall}}&{\text{Gewichtsfunktion }}w(x)\\\hline [-1,1]&1\\[-1,1]&1/{\sqrt {1-x^{2}}}\\[0,\infty )&e^{-x}\\(-\infty ,\infty )&e^{-x^{2}}\\[0,\infty )&e^{-x^{2}}x^{\alpha },\ \alpha >-1\\\hline \end{array}}}

Wir definieren in diesem Zusammenhang auf dem Raum aller Polynome $\Pi$ das durch $w$ induzierte Skalarprodukt

(8.39) $\langle f,g\rangle :=\int \limits _{a}^{b}w(x)f(x)g(x)\,dx,\quad f,g\in \Pi .$

Das Integral in (8.39) existiert offenbar unter den Voraussetzungen an $w$. Für alle $f,g,h\in \Pi$ gilt weiter (man verifiziere dies)

$\langle f,f\rangle \geq 0,\quad \langle f,f\rangle =0\Leftrightarrow f\equiv 0,\quad \langle f,g\rangle =\langle g,f\rangle ,$

$\langle \alpha f+\beta g,h\rangle =\alpha \langle f,h\rangle +\beta \langle g,h\rangle =\langle h,\alpha f+\beta g\rangle ,\quad \alpha ,\beta \in \mathbb {R} .$

Insbesondere ist also die Abbildung $\langle \cdot ,\cdot \rangle :\Pi \times \Pi \to \mathbb {R}$ in beiden Eingängen linear. Wir verwenden ferner die durch das Skalarprodukt $\langle \cdot ,\cdot \rangle$ induzierte Norm auf $\Pi$

(8.40) $\|f\|:=\langle f,f\rangle ^{1/2}=\left\{\int \limits _{a}^{b}w(x)f^{2}(x)dx\right\}^{1/2},\quad f\in \Pi .$

Ziel ist es nun wieder, zur numerischen Berechnung von ${\mathcal {I}}(f)$ eine Quadraturformel

(8.41) ${\mathcal {I}}_{n}(f):=\sum _{i=0}^{n}\sigma _{i}f(x_{i})$

herzuleiten. (Man beachte, dass hier der Faktor $b-a$ vor der Summe fehlt.) Und zwar soll eine interpolatorische Quadraturformel entwickelt werden, für welche bei geeigneter Wahl der Stützstellen $x_{i}$ und der Gewichte $\sigma _{i}$ der Genauigkeitsgrad möglichst hoch ist, welche also Polynome bis zu einem möglichst hohen Grad exakt integriert. Man betrachte dazu die Aussagen in Satz 8.15 über den Quadraturfehler. Die Begriffe Genauigkeitsgrad und interpolatorische Quadraturformel sind hierbei analog zu den Definitionen 8.2 und 8.8 auf Integrale mit Gewichten zu übertragen.

Zunächst einmal stellen wir fest, dass man in (8.41) insgesamt 2n+2 Parameter $x_{i}$ und $\sigma _{i}$ zur Verfügung hat, was der Anzahl der Koeffizienten eines Polynoms vom Grad $2n+1$ entspricht. In der Tat werden wir zeigen, dass eine Quadraturformel mit diesem Genauigkeitsgrad existiert. Quadraturformeln mit einem höheren Genauigkeitsgrad kann es nicht geben. Denn wäre ${\mathcal {I}}_{n}$ eine Quadraturformel mit Stützstellen $x_{i}$ $(i=0,1,\ldots ,n)$ und Gewichten $\sigma _{i}$ $(i=0,1,\ldots ,n)$ und hätte ${\mathcal {I}}_{n}$ den Genauigkeitsgrad $2n+2$, so folgte insbesondere für das Polynom

$p(x):=\prod _{i=0}^{n}(x-x_{i})^{2}\in \Pi _{2n+2}$

$p(x_{i})=0$ $(i=0,1,\ldots ,n)$ und daher ${\mathcal {I}}_{n}(p)=0={\mathcal {I}}(p)$. Wegen

$p(x)>0,\quad x\in (x_{i},x_{i+1}),\quad i=0,1,\ldots ,n-1$

ist jedoch ${\mathcal {I}}(p)>0$. Wir können weiter schließen:

Lemma 8.29

Ist ${\mathcal {I}}_{n}$ eine Quadraturformel mit Stützstellen $x_{i}$ $(i=0,1,\ldots ,n)$ und Gewichten $\sigma _{i}$ $(i=0,1,\ldots ,n)$ und hat ${\mathcal {I}}_{n}$ den Genauigkeitsgrad $2n+1$, so gilt

$\langle p,p_{n+1}\rangle =0,\quad p\in \Pi _{n}$

für

$p_{n+1}(x):=a_{n}(x-x_{0})(x-x_{1})\cdots (x-x_{n})$

mit beliebigem $a_{n}\in \mathbb {R} \setminus \{0\}$.

Beweis.

Für $p\in \Pi _{n}$ folgt $pp_{n+1}\in \Pi _{2n+1}$ und somit

$\langle p,p_{n+1}\rangle =\int \limits _{a}^{b}w(x)p(x)p_{n+1}(x)\,dx={\mathcal {I}}_{n}(pp_{n+1})=\sum _{i=0}^{n}\sigma _{i}p(x_{i})\underbrace {p_{n+1}(x_{i})} _{=0}=0.$

q.e.d.

Zwei Funktionen $f$ und $g$, für die $\langle f,g\rangle =0$ gilt, nennt man orthogonal zueinander. Für eine Quadraturformel mit Genauigkeitsgrad $2n+1$ sollten die Stützstellen $x_{i}$ $(i=0,1,\ldots ,n)$ also gerade als Nullstellen eines Polynoms vom Grad $n+1$ gewählt werden, welches bezüglich des Skalarproduktes $\langle \cdot ,\cdot \rangle$ orthogonal zu dem ganzen Raum $\Pi _{n}$ ist. Offenbar kann man ein solches Polynom gewinnen, indem man durch Anwendung des Gram-Schmidt-Orthogonalisierungsverfahren auf die Monome $1,x,\ldots ,x^{n+1}$ orthogonale Polynome $p_{j}\in \Pi _{j}$ hinsichtlich $\langle \cdot ,\cdot \rangle$ erzeugt. Diese Polynome haben nämlich die Eigenschaft

$\Pi _{k}=\operatorname {span} \{p_{0},p_{1},\ldots ,p_{k}\},\quad \langle p_{i},p_{j}\rangle =\delta _{ij}\|p_{i}\|^{2}\quad (i,j,k=0,1,\ldots ,n+1),$

so dass sich insbesondere jedes $p\in \Pi _{n}$ mit gewissen $a_{j}$ $(j=0,1,\ldots ,n)$ in der Form $p(x)=\sum _{j=0}^{n}a_{j}p_{j}(x)$ schreiben lässt und folglich mit den zugehörigen $a_{j}$ gilt:

(8.42) $\langle p,p_{n+1}\rangle =\left\langle \sum _{j=0}^{n}a_{j}p_{j}(x),p_{n+1}\right\rangle =\sum _{j=0}^{n}a_{j}\langle p_{j},p_{n+1}\rangle =0,\quad p\in \Pi _{n}.$

Darüber hinaus haben diese orthogonalen Polynome $p_{j}$ nur reelle und einfache Nullstellen, welche alle im Intervall $(a,b)$ liegen, wie im nächsten Unterabschnitt gezeigt wird. Die Stützstellen $x_{i}$ $(i=0,1,\ldots ,n)$ sollten demzufolge gerade die Nullstellen des $(n+1)$-ten dieser orthogonalen Polynome sein. Die Gewichte $\sigma _{i}$ $(i=0,1,\ldots ,n)$ einer derartigen Quadraturformel sind dann gemäß Satz 8.4, der genauso für gewichtete Integrale gilt, durch

$\sigma _{i}:=I(L_{i})=\int \limits _{a}^{b}w(x)L_{i}(x)\,dx=\langle L_{i},\mathbf {1} \rangle ,\quad i=0,1,\ldots ,n$

festgelegt, wobei $L_{i}$ wieder die Lagrangeschen Basispolynome zu den Stützstellen $x_{k}$ sind. Nach Definition 8.8 (entsprechend für gewichtete Integrale formuliert) handelt es sich bei der so definierten Formel um eine interpolatorische Quadraturformel.

Bevor wir diese sog. Gaußschen Quadraturformeln noch etwas näher betrachten, wollen wir auf ihren wesentlichen Baustein, orthogonale Polynome, näher eingehen.

8.5.2 Orthogonale Polynome

Wie bereits im vorigen Unterabschnitt gesagt wurde, erhält man eine spezielle Folge ${\tilde {p}}_{n}$ paarweise orthonormaler Polynome ${\tilde {p}}_{n}\in \Pi _{n}$ durch Gram-Schmidt-Orthonormalisierung der Monome $1,x,x^{2},\ldots$:

$p_{0}:=1,\quad {\tilde {p}}_{0}:={\frac {p_{0}}{\|p_{0}\|}},\quad p_{n}:=x^{n}-\sum _{j=0}^{n-1}\langle x^{n},{\tilde {p}}_{j}\rangle {\tilde {p}}_{j},\quad {\tilde {p}}_{n}:={\frac {p_{n}}{\|p_{n}\|}}\quad (n=1,2,\ldots ).$

Statt mit den orthonormalen Polynomen ${\tilde {p}}_{n}$ zu arbeiten, deren Hauptkoeffizienten i. a. von 1 verschieden sind, ist es manchmal bequemer, dies mit den orthogonalen Polynomen $p_{n}\in \Pi _{n}$ zu tun, d. h. mit

(8.43) $p_{0}:=1,\quad p_{n}:=x^{n}-\sum _{j=0}^{n-1}{\frac {\langle x^{n},p_{j}\rangle }{\|p_{j}\|^{2}}}p_{j}\quad (n=1,2,\ldots ).$

Diese unterscheiden sich von den ${\tilde {p}}_{n}$ nur durch den Skalar $1/\|p_{n}\|$ und haben offensichtlich den Hauptkoeffizienten 1. Für sie gilt

(8.44) $\Pi _{k}=\operatorname {span} \{p_{0},p_{1},\ldots ,p_{k}\},\quad \langle p_{i},p_{j}\rangle =\delta _{ij}\|p_{i}\|^{2}\quad (i,j,k\in \mathbb {N} _{0})$

und somit (vgl. (8.42))

(8.45) $\langle p,p_{k}\rangle =0,\quad p\in \Pi _{k-1}\quad (k=1,2,\ldots ).$

Nach Konstruktion ist also $p_{n}$ ein Polynom vom genauen Grad $n$ mit Hauptkoeffizienten 1.

Die Polynome $p_{n}$ können statt über die Formel (8.43) auch nach der im folgenden Satz angegebenen Drei-Term-Rekursionsformel berechnet werden.

Satz 8.30

Die Orthogonalpolynome in (8.43) genügen der Drei-Term-Rekursionsformel

$p_{0}(x)=1,\quad p_{1}(x)=x-\beta _{0},\quad p_{n+1}(x)=(x-\beta _{n})p_{n}(x)-\gamma _{n}^{2}p_{n-1}(x),\quad n=1,2,\ldots$

mit den Koeffizienten

$\beta _{n}:={\frac {\langle xp_{n},p_{n}\rangle }{\|p_{n}\|^{2}}}\quad (n=0,1,\ldots ),\quad \gamma _{n}^{2}:={\frac {\|p_{n}\|^{2}}{\|p_{n-1}\|^{2}}}\quad (n=1,2,\ldots ).$

Beweis.

Offenbar ist die behauptete Darstellung richtig für $p_{0}$ und $p_{1}$. Für $n\geq 1$ setzen wir

$q_{n+1}:=(x-\beta _{n})p_{n}-\gamma _{n}^{2}p_{n-1}$

und zeigen im Folgenden $q_{n+1}=p_{n+1}$. Dazu stellen wir fest, dass $p_{n+1}$ sowie $q_{n+1}$ Polynome vom genauen Grad $n+1$ sind und beide den Hauptkoeffizienten 1 haben. Somit gilt

(8.46) $r:=p_{n+1}-q_{n+1}\in \Pi _{n}.$

Wir zeigen nun, dass $q_{n+1}$ wie $p_{n+1}$ orthogonal zu dem ganzen Raum $\Pi _{n}$ ist und damit

(8.47) $\langle p,q_{n+1}\rangle =0,\quad p\in \Pi _{n}$

gilt. Die Beziehungen (8.46) und (8.47) zusammen ergeben dann

$\|r\|^{2}=\langle r,r\rangle =\langle r,p_{n+1}-q_{n+1}\langle =\langle r,p_{n+1}\rangle -\langle r,q_{n+1}\rangle =0$

und folglich $r=0$ bzw., wie behauptet, $p_{n+1}=q_{n+1}$.

Wir wollen nun (8.47) nachweisen. Aufgrund von $\langle p_{n},p_{n-1}\rangle =0$ erhalten wir mit der Definition von $\beta _{n}$

(8.48) $\langle q_{n+1},p_{n}\rangle =\langle (x-\beta _{n})p_{n}-\gamma _{n}^{2}p_{n-1},p_{n}\rangle =\langle xp_{n},p_{n}\rangle -\beta _{n}\|p_{n}\|^{2}=0$

und mit der Definition von $\gamma _{n}$

$\langle q_{n+1},p_{n-1}\rangle =\langle (x-\beta _{n})p_{n}-\gamma _{n}^{2}p_{n-1},p_{n-1}\rangle =\langle xp_{n},p_{n-1}\rangle -\gamma _{n}^{2}\|p_{n-1}\|^{2}=\langle p_{n},xp_{n-1}\rangle -\langle p_{n},p_{n}\rangle$

(8.49) $=\langle p_{n},xp_{n-1}-p_{n}\rangle =0,$

wobei das letzte Gleichheitszeichen wegen $xp_{n-1}-p_{n}\in \Pi _{n-1}$ gilt. Schließlich folgt:

(8.50) $\langle q_{n+1},p_{j}\rangle =\underbrace {\langle p_{n},xp_{j}\rangle } _{=0}-\beta _{n}\underbrace {\langle p_{n},p_{j}\rangle } _{=0}-\gamma _{n}^{2}\underbrace {\langle p_{n-1},p_{j}\rangle } _{=0}=0,\quad j=0,1,\ldots ,n-2.$

Da sich jedes $p\in \Pi _{n}$ gemäß (8.44) als Linearkombination der $p_{j}$ $(j=0,1,\ldots ,n)$ darstellen lässt, folgt aus (8.48), (8.49) und (8.50) für jedes $p\in \Pi _{n}$ mit gewissen $a_{j}$

$\langle q_{n+1},p\rangle =\left\langle q_{n+1},\sum _{j=0}^{n}a_{j}p_{j}(x)\right\rangle =\sum _{j=0}^{n}a_{j}\langle q_{n+1},p_{j}\rangle =0.$

Damit ist alles gezeigt.

q.e.d.

Für die Nullstellen der $p_{j}$ $(j\in \mathbb {N} )$ in (8.43) hat man folgende Aussage:

Satz 8.31

Die Nullstellen $x_{i}$ $(i=0,1,\ldots ,n-1)$ des $n$-ten Orthogonalpolynoms $p_{n}$ in (8.43) sind reell, einfach und liegen alle in $(a,b)$. Sie besitzen die Darstellung

(8.51) $x_{i}={\frac {\langle xL_{i},L_{i}\rangle }{\|L_{i}\|^{2}}}\quad (i=0,1,\ldots ,n-1),$

wobei $L_{i}\in \Pi _{n-1}$ die zu den $x_{k}$ $(k=0,1,\ldots ,n-1)$ gehörenden Lagrangeschen Basispolynome sind.

Beweis.

Es seien die Nullstellen $x_{i}$ von $p_{n}$ so durchnumeriert, dass $a<x_{0}<\ldots <x_{j-1}<b$ $(0\leq j\leq n)$ diejenigen Nullstellen von $p_{n}$ in $(a,b)$ seien, an denen $p_{n}$ sein Vorzeichen wechselt und die somit eine ungerade Vielfachheit haben. Wäre nun $j-1\leq n-2$ bzw. $j\leq n-1$, so hätte das Polynom

$q(x):=\prod _{k=0}^{j-1}(x-x_{k})$

den Grad $j\leq n-1$, so dass wegen (8.45)

(8.52) $\langle p_{n},q\rangle =0$

folgte. Da die $x_{k}$ $(k=j,j+1,\ldots ,n-1)$ Nullstellen von $p_{n}$ mit gerader Vielfachheit wären, wäre dann aber

$p_{n}(x)q(x)=\left(\prod _{k=0}^{j-1}(x-x_{k})^{2}\right)\left(\prod _{k=j}^{n-1}(x-x_{k})\right)\geq 0,\quad x\in [a,b]$

und demzufolge

$\langle p_{n},q\rangle =\int \limits _{a}^{b}w(x)p_{n}(x)q(x)\,dx>0$

im Widerspruch zu (8.52). Also ist $j=n$.

Um zur Darstellung (8.51) zu gelangen, schreibt man $p_{n}$ für $n\geq 1$ in der Form

$p_{n}(x)=(x-x_{i}){\hat {q}}(x),\quad {\hat {q}}(x):=\prod _{k=0 \atop k\neq i}^{n-1}(x-x_{k}).$

Es folgt ${\hat {q}}\in \Pi _{n-1}$ sowie

$0=\langle p_{n},{\hat {q}}\rangle =\langle x{\hat {q}},{\hat {q}}\langle -x_{i}\langle {\hat {q}},{\hat {q}}\rangle .$

Daraus ergibt sich wegen $\langle {\hat {q}},{\hat {q}}\rangle \neq 0$

$x_{i}={\frac {\langle x{\hat {q}},{\hat {q}}\rangle }{\|{\hat {q}}\|^{2}}}={\frac {\langle xL_{i},L_{i}\rangle }{\|L_{i}\|^{2}}}$

wobei sich die letzte Gleichung aus der Tatsache ergibt, dass das Polynom ${\hat {q}}$ bis auf einen konstanten Faktor mit $L_{i}$ übereinstimmt.

q.e.d.

In folgender Tabelle sind für verschiedene Intervalle und Gewichtsfunktionen die Bezeichnungen der zugehörigen orthogonalen Polynome aufgelistet:

Fehler beim Parsen (Konvertierungsfehler. Der Server („https://wikimedia.org/api/rest_“) hat berichtet: „Cannot get mml. TeX parse error: Bracket argument to \\ must be a dimension“): {\displaystyle {\begin{array}{|c|c|c|}\hline {\text{Intervall}}&{\text{Gewichtsfunktion }}w(x)&{\text{Name}}\\\hline [-1,1]&1&{\text{Legendre-Polynome}}\\[-1,1]&(1-x)^{\alpha }(1+x)^{\beta },\alpha ,\beta >-1&{\text{Jacobi-Polynome}}\\[-1,1]&1/{\sqrt {1-x^{2}}}&{\text{Tschebyscheff-Pol. der 1. Art}}\\[-1,1]&{\sqrt {1-x^{2}}}&{\text{Tschebyscheff-Pol. der 2. Art}}\\[0,\infty )&e^{-x^{2}}x^{\alpha },\alpha >-1&{\text{Laguerre-Polynome}}\\(-\infty ,\infty )&e^{-x^{2}}&{\text{Hermite-Polynome}}\\\hline \end{array}}}

Man kann zeigen (siehe z. B. E. W. Cheney: Introduction to Approximation Theory, 2nd ed., Chelsea Publish. Comp., New York, 1982):

Satz 8.32

Für $p_{n}$ aus (8.43) gilt

$\|p_{n}\|=\min _{a_{0},\ldots ,a_{n-1}\in \mathbb {R} }\left\|x^{n}+a_{n-1}x^{n-1}+\ldots +a_{1}x+a_{0}\right\|.$

Unter allen Polynomen vom Grad $n$ mit Hauptkoeffizientem 1 macht also $p_{n}$ die Norm in (8.40) minimal. Im Fall der Tschebyscheff-Polynome 1. Art minimiert $p_{n}$ unter all diesen Polynomen überdies die Maximum-Norm (Satz 6.19) und im Fall der Tschebyscheff-Polynome 2. Art (s. Cheney) die (ungewichtete) $L_{1}$-Norm

$\|f\|:=\int \limits _{a}^{b}|f(x)|\,dx,\quad f\in \Pi .$

Beispiel 8.33

Mit Satz 8.30 sollen die Legendre-Polynome für $n=0,1,2,3$ berechnet werden. Es ist somit $a:=-1,b:=1,w\equiv 1$ und folglich

$\beta _{n}={\frac {\langle xp_{n},p_{n}\rangle }{\|p_{n}\|^{2}}}={\frac {\int _{-1}^{1}xp_{n}^{2}(x)\,dx}{\int _{-1}^{1}p_{n}^{2}(x)\,dx}},\quad \gamma _{n}^{2}={\frac {\|p_{n}\|^{2}}{\|p_{n-1}\|^{2}}}={\frac {\int _{-1}^{1}p_{n}^{2}(x)\,dx}{\int _{-1}^{1}p_{n-1}^{2}(x)\,dx}}.$

Mit $p_{0}(x)=1$ ist

$\beta _{0}={\frac {\int _{-1}^{1}x\,dx}{\int _{-1}^{1}dx}}=0$

und damit weiter $p_{1}(x)=x$. Es ergeben sich ferner

$\beta _{1}={\frac {\int _{-1}^{1}x\cdot x^{2}\,dx}{\int _{-1}^{1}x^{2}\,dx}}=0,\quad \gamma _{1}^{2}={\frac {\int _{-1}^{1}x^{2}\,dx}{\int _{-1}^{1}dx}}={\frac {1}{3}}$

und demnach $p_{2}(x)=x^{2}-{\frac {1}{3}}$ sowie

$\beta _{2}={\frac {\int _{-1}^{1}x\left(x^{2}-{\frac {1}{3}}\right)\,dx}{\int _{-1}^{1}\left(x^{2}-{\frac {1}{3}}\right)^{2}\,dx}}=0,\quad \gamma _{1}^{2}={\frac {\int _{-1}^{1}\left(x^{2}-{\frac {1}{3}}\right)^{2}\,dx}{\int _{-1}^{1}x^{2}\,dx}}={\frac {4}{15}},$

so dass folgt

$p_{3}(x)=x^{3}-{\frac {1}{3}}x-{\frac {4}{15}}x=x^{3}-{\frac {3}{5}}x.$

8.5.3 Die Quadraturformeln

Satz 8.34

Für $n\in \mathbb {N}$ seien $p_{j}$ $(j=0,1,\ldots ,n,n+1)$ die durch (8.43) definierten bezüglich $\langle \cdot ,\cdot \rangle$ orthogonalen Polynome, $x_{i}$ $(i=0,1,\ldots ,n)$ die Nullstellen von $p_{n+1}$ und $\sigma _{i}$ die durch

$\sigma _{i}:=\langle L_{i},1\rangle \quad (i=0,1,\ldots ,n)$

definierten Gewichte. Dann ist die Quadraturformel

(8.53) ${\mathcal {I}}_{n}(f):=\sum _{i=0}^{n}\sigma _{i}f(x_{i})$

interpolatorisch und hat (exakt) den Genauigkeitsgrad $2n+1$.

Beweis.

Nach Definition 8.8 (entsprechend für gewichtete Integrale formuliert) ist ${\mathcal {I}}_{n}$ aufgrund der Wahl der Gewichte eine interpolatorische Quadraturformel. Nach Korollar 8.9 hat eine solche mindestens den Genauigkeitsgrad $n$. Wir wollen nun zeigen, dass sie mindestens den Genauigkeitsgrad $2n+1$ und damit exakt den Genauigkeitsgrad $2n+1$ besitzt, wie aus den Argumenten in Abschnitt 8.5.1 hervorgeht.

Es sei $p\in \Pi _{2n+1}$ beliebig. Dann lässt sich $p$ mit gewissen $q,r\in \Pi _{n}$ nach einer Polynomdivision mit Rest in der Form

$p=qp_{n+1}+r$

schreiben. Wegen $p_{n+1}(x_{i})=0$ gilt dann

$p(x_{i})=r(x_{i}),\quad i=0,1,\ldots ,n.$

Mit der Lagrangeschen Interpolationsformel (6.7), angewandt auf $r$, erhält man weiter

$r(x)=\sum _{i=0}^{n}r(x_{i})L_{i}(x)=\sum _{i=0}^{n}p(x_{i})L_{i}(x).$

Somit schließt man

(8.54) ${\mathcal {I}}(p)=\langle p,\mathbf {1} \rangle =\underbrace {\langle q,p_{n+1}\rangle } _{=0}+\langle r,\mathbf {1} \rangle =\sum _{i=0}^{n}p(x_{i})\langle L_{i},\mathbf {1} \rangle =\sum _{i=0}^{n}\sigma _{i}p(x_{i}),$

womit der Genauigkeitsgrad von mindestens $2n+1$ für ${\mathcal {I}}_{n}$ nachgewiesen ist.

q.e.d.

Die Quadraturformel in (8.53) mit den in Satz 8.34 genannten Stützstellen $x_{i}$ und Gewichten $\sigma _{i}$ bezeichnet man als Gaußsche Quadraturformel. Ihr Genauigkeitsgrad ist optimal, da es keine Quadraturformeln mit Genauigkeitsgrad $2n+2$ gibt (vgl. Abschnitt 8.5.1). Weiter hat man:

Lemma 8.35

Für die Gewichte $\sigma _{i}$ der Gaußschen Quadraturformel ${\mathcal {I}}_{n}$ von Satz 8.34 gilt

$\sigma _{i}=\langle L_{i},L_{i}\rangle >0\quad (i=0,1,\ldots ,n)$

und

(8.55) $\sum _{i=0}^{n}\sigma _{i}=\int \limits _{a}^{b}w(x)\,dx.$

Beweis.

Wendet man die Beziehungen (8.54) auf $p:=L_{j}^{2}\in \Pi _{2n}$ an, so folgt

$\langle L_{j},L_{j}\rangle =\left\langle L_{j}^{2},\mathbf {1} \right\rangle =\sum _{i=0}^{n}\sigma _{i}L_{j}^{2}(x_{i})=\sigma _{j}.$

Weiter gilt $\langle L_{j},L_{j}\rangle >0,daman<math>L_{j}^{2}(x)\geq 0,x\in [a,b]$ sowie $L_{j}^{2}(x)>0$ z. B. für alle $x\in (x_{j-1},x_{j+1})$ hat. Die Beziehung (8.55) folgt nun mit Satz 8.34 aus

$\sum _{i=0}^{n}\sigma _{i}={\mathcal {I}}_{n}(1)={\mathcal {I}}(1)=\int \limits _{a}^{b}w(x)\,dx.$

q.e.d.

Anders als bei den abgeschlossenen Newton-Cotes-Formeln haben also die Gaußschen Quadraturformeln für alle $n\in \mathbb {N} _{0}$ positive Gewichte. Mit dem folgenden Satz geben wir abschließend eine Darstellung für den bei der Gauß-Quadratur entstehenden Quadraturfehler an.

Satz 8.36

Es sei $f\in C^{2n+2}[a,b]$ und In die Gaußsche Quadraturformel mit Stützstellen $x_{k}$ $(k=0,1,\ldots ,n)$. Dann gilt für

${\hat {\gamma }}_{2n+1}:=\int \limits _{0}^{1}w((b-a)t+a)^{2}\prod _{k=0}^{n}(t-t_{k})^{2}\,dt$

mit

$t_{k}:={\frac {x_{k}-a}{b-a}}\quad (k=0,1,\ldots ,n)$

und für ein $\xi \in [a,b]$:

${\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)={\hat {\gamma }}_{2n+1}{\frac {(b-a)^{2n+3}}{(2n+2)!}}f^{(2n+2)}(\xi ).$

Beweis.

Den Satz 8.15 kann man auf den Fall gewichteter Integrale übertragen und dann aufgrund von Satz 8.34 auf die Gaußsche Quadraturformel ${\mathcal {I}}_{n}$ mit $r:=2n+1$ anwenden. Man wählt dort zu den Stützpunkten $x_{k}$ von ${\mathcal {I}}_{n}$ die weiteren Stützpunkte $t_{n+1}:=t_{0},\ldots ,t_{2n+1}:=t_{n}$, so dass insbesondere

$\prod _{k=0}^{2n+1}(t-t_{k})=\prod _{k=0}^{n}(t-t_{k})^{2}=p_{n+1}^{2}(x)\geq 0,\quad t\in [0,1]$

gilt. Weiter folge man dann dem Beweis von Satz 8.15.

q.e.d.

Natürlich ist es auch möglich, summierte Gaußsche Quadraturformeln zu definieren und zu verwenden. Die Resultate in Abschnitt 8.3 lassen sich ganz kanonisch auf solche Formeln übertragen.

8.5.4 Berechnung der Stützstellen und Gewichte

Beispielsweise für die Tschebyscheff-Polynome 1. Art kann man die Nullstellen explizit angeben (vgl. Satz 6.18). Im allgemeinen steht man bei Verwendung der Gaußschen Quadraturformeln für größere Werte von $n$ aber vor dem Problem, die $n+1$ Nullstellen $x_{i}$ des Polynoms $p_{n+1}$ der bezüglich $\langle \cdot ,\cdot \rangle$ orthogonalen Polynome $p_{j}$ $(j=0,1,2,\ldots )$ und/oder die Gewichte $\sigma _{i}:=\langle L_{i},\mathbf {1} \rangle$ zu bestimmen. Wir wollen hier abschließend einen Weg zu ihrer Berechnung aufzeigen. Dazu gehen wir davon aus, dass die Koeffizienten $\beta _{j}$ und $\gamma _{j}$ in der Rekursionsformel (8.43) für die orthogonalen $p_{j}$

(8.56) $p_{0}:=1,\quad p_{1}:=x-\beta _{0},$

(8.57) $p_{k+1}:=(x-\beta _{k})pk-\gamma _{k}^{2}p_{k-1},\quad k=1,2,\ldots$

bereits bekannt sind und somit die symmetrische Tridiagonalmatrix

(8.58) $J:={\begin{pmatrix}\beta _{0}&-\gamma _{1}&0&\ldots &0\\-\gamma _{1}&\beta _{1}&-\gamma _{2}&\ddots &\vdots \\0&-\gamma _{2}&\ddots &\ddots &0\\\vdots &\ddots &\ddots &\ddots &-\gamma _{n}\\0&\ldots &0&-\gamma _{n}&\beta _{n}\end{pmatrix}}\in \mathbb {R} ^{(n+1)\times (n+1)}$

aufgestellt werden kann. Der folgende Satz besagt nun, dass die Stützstellen $x_{i}$ der Gaußschen Quadraturformeln die Eigenwerte von $J$ sind und sich deren Gewichte $\sigma _{i}$ aus zugehörigen Eigenvektoren bestimmen lassen.

Satz 8.37

Für die Stützstellen $x_{i}$ $(i=0,1,\ldots ,n)$ und die Gewichte $\sigma _{i}$ $(i=0,1,\ldots ,n)$ der Gaußschen Quadraturformel ${\mathcal {I}}_{n}$ gilt mit

$v^{(i)}:=(\underbrace {\tau _{0}p_{0}(x_{i})} _{=1},\ldots ,\tau _{n}p_{n}(x_{i}))^{T}$

für

$\tau _{k}:={\begin{cases}1&falls\ k=0,\\(-1)^{k}/(\gamma _{1}\gamma _{2}\cdots \gamma _{k})&falls\ k\in \{1,\ldots ,n\},\end{cases}}$

wobei $p_{k}$ $(k=0,1,\ldots ,n)$ die bezüglich $\langle \cdot ,\cdot \rangle$ orthogonalen Polynome seien:

(8.59) $Jv^{(i)}=x_{i}v^{(i)}\quad (i=0,1,\ldots ,n)$

und

(8.60) $\sigma _{i}={\frac {\langle \mathbf {1} ,\mathbf {1} \rangle }{(v^{(i)})^{T}v^{(i)}}}={\frac {\langle \mathbf {1} ,\mathbf {1} \rangle }{\sum _{k=0}^{n}\tau _{k}^{2}p_{k}^{2}(x_{i})}}\quad (i=0,1,\ldots ,n).$

Beweis.

Aus den Definitionen von $J,v^{(i)},p_{1}$ und $\tau _{1}$ ergibt sich

$\left(Jv^{(i)}\right)_{1}=\beta _{0}-\gamma _{1}\tau _{1}p_{1}(x_{i})=\beta _{0}+p_{1}(x_{i})=\beta _{0}+x_{i}\beta _{0}=x_{i}=x_{i}v_{1}^{(i)}.$

Definiert man $\gamma _{n+1}=\tau _{n+1}:=0$ und berücksichtigt man $p_{n+1}(x_{i})=0$, so erhält man aus den Rekursionsformeln (8.56) und (8.57) mit $x=x_{i}$ weiter

$\left(Jv^{(i)}\right)_{k+1}=-\gamma _{k}\tau _{k-1}p_{k-1}(x_{i})+\beta _{k}\tau _{k}p_{k}(x_{i})-\gamma _{k+1}\tau _{k+1}p_{k+1}(x_{i})$

$=-\gamma _{k}(-\gamma _{k}\tau _{k})p_{k-1}(x_{i})+\beta _{k}\tau _{k}p_{k}(x_{i})-\gamma _{k+1}{\frac {(-1)}{\gamma _{k+1}}}p_{k+1}(x_{i})$

$=\tau _{k}\left[\gamma _{k}^{2}p_{k-1}(x_{i})+\beta _{k}p_{k}(x_{i})+(x_{i}-\beta _{k})p_{k}(x_{i})-\gamma _{k}^{2}p_{k-1}(x_{i})\right]=\tau _{k}x_{i}p_{k}(x_{i})=x_{i}v_{k+1}^{(i)},\quad k=1,\ldots ,n.$

Damit ist (8.59) bewiesen.

Gleichung (8.59) besagt, dass $v^{(i)}$ Eigenvektor zum Eigenwert $x_{i}$ der Matrix $J$ ist. Gemäß Satz 8.31 sind diese Eigenwerte paarweise verschieden. Da für eine reelle symmetrische Matrix Eigenvektoren zu paarweise verschiedenen Eigenwerten orthogonal zueinander sind, folgt

(8.61) $\left(v^{(i)}\right)^{T}v^{(k)}=0\quad (i\neq k).$

Da ferner die Polynome $p_{j}$ $(j=0,1,2,\ldots )$ paarweise orthogonal sind und die Gaußsche Quadraturformel alle $p_{k}$ $(k=0,1,\ldots ,n)$ exakt integriert, folgt weiter mit Satz 8.34

(8.62) $\delta _{k0}\langle \mathbf {1} ,\mathbf {1} \rangle =\langle p_{k},p_{0}\rangle ={\mathcal {I}}(p_{k})=\sum _{j=0}^{n}\sigma _{j}p_{k}(x_{j})\quad (k=0,1,\ldots ,n),$

wobei $\delta _{ij}$ das Kroneckersymbol ist. Multiplikation von (8.62) mit $\tau _{k}^{2}p_{k}(x_{i})$ und anschließende Summation über $k$ liefert schließlich unter Verwendung von (8.61)

$\langle \mathbf {1} ,\mathbf {1} \rangle \sum _{k=0}^{n}\tau _{k}^{2}p_{k}(x_{i})\delta _{k0}=\langle \mathbf {1} ,\mathbf {1} \rangle \tau _{0}^{2}p_{0}(x_{i})=\langle \mathbf {1} ,\mathbf {1} \rangle =\sum _{k=0}^{n}\sum _{j=0}^{n}\sigma _{j}\tau _{k}^{2}p_{k}(x_{i})p_{k}(x_{j})$

$=\sum _{j=0}^{n}\sigma _{j}\sum _{k=0}^{n}\tau _{k}^{2}p_{k}(x_{i})p_{k}(x_{j})=\sum _{j=0}^{n}\sigma _{j}\left(v^{(i)}\right)^{T}v^{(j)}=\sigma _{i}\left(v^{(i)}\right)^{T}v^{(i)}$

Damit ist auch die Gültigkeit von (8.60) bewiesen.

q.e.d.

Beispiel 8.38

Wir verwenden Satz 8.37 für die Herleitung der Gaußschen Quadraturformel mit $n:=2$ zur näherungsweisen Berechnung des Integrals

$\int \limits _{-1}^{1}f(x)\,dx.$

Beispiel 8.33 liefert

$J:={\begin{pmatrix}0&-1/{\sqrt {3}}&0\\-1/{\sqrt {3}}&0&-2/{\sqrt {15}}\\0&-2/{\sqrt {15}}&0\end{pmatrix}}.$

Die Eigenwerte von $J$ berechnen sich aus

$\det(J-\lambda I)=-\lambda \left(\lambda ^{2}-{\frac {4}{15}}\right)+{\frac {1}{3}}\lambda =\lambda \left({\frac {3}{5}}-\lambda ^{2}\right)$

so dass man die Stützstellen $x_{0}=-{\sqrt {3/5}},x_{1}=0$ und $x_{2}={\sqrt {3/5}}$ erhält. Weiter hat man

$\tau _{0}=1,\quad \tau _{1}=-1/\gamma _{1}=-{\sqrt {3}},\quad \tau _{2}=1/(\gamma _{1}\gamma _{2})=3{\sqrt {5}}/2$

sowie

$\langle \mathbf {1} ,\mathbf {1} \rangle =\int \limits _{-1}^{1}dx=2.$

Mit

$p_{0}(x)=1,\quad p_{1}(x)=x,\quad p_{2}(x)=x^{2}-{\frac {1}{3}}$

hat man

${\begin{array}{|c||c|c|c|}\hline i&p_{0}(x_{i})&p_{1}(x_{i})&p_{2}(x_{i})\\\hline 0&1&-{\sqrt {3/5}}&4/15\\1&1&0&-1/3\\2&1&{\sqrt {3/5}}&4/15\\\hline \end{array}}$

und demnach

$\sigma _{0}={\frac {2}{\tau _{0}^{2}p_{0}^{2}(x_{0})+\tau _{1}^{2}p_{1}^{2}(x_{0})+\tau _{2}^{2}p_{2}^{2}(x_{0})}}={\frac {5}{9}},$

$\sigma _{1}={\frac {2}{\tau _{0}^{2}p_{0}^{2}(x_{1})+\tau _{1}^{2}p_{1}^{2}(x_{1})+\tau _{2}^{2}p_{2}^{2}(x_{1})}}={\frac {8}{9}},$

$\sigma _{2}={\frac {2}{\tau _{0}^{2}p_{0}^{2}(x_{2})+\tau _{1}^{2}p_{1}^{2}(x_{2})+\tau _{2}^{2}p_{2}^{2}(x_{2})}}={\frac {5}{9}}.$

Man erhält also die Gaußsche Quadraturformel

${\mathcal {I}}_{2}(f):={\frac {5}{9}}f\left(-{\sqrt {\frac {3}{5}}}\right)+{\frac {8}{9}}f(0)+{\frac {5}{9}}f\left({\sqrt {\frac {3}{5}}}\right)$

Die gesuchten Eigenwerte von $J$ müssen aber, zumindest für größere $n$, normalerweise mit einer numerischen Methode bestimmt werden.