Aus Wikiversity

Inhaltsverzeichnis

1 8.1 Einführung
2 8.2 Interpolatorische Quadraturformeln
3 8.3 Summierte abgeschlossene Newton-Cotes-Formel
4 8.4 Extrapolationsverfahren
- 4.1 8.4.1 Einführung
  - 4.1.1 Satz 8.23
- 4.2 8.4.2 Das Verfahren
5 8.5 Gaußsche Quadraturformeln

[Bearbeiten] 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^b_a 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^b_a [\alpha f(x) + \beta g(x)]\, dx = \alpha \int\limits^b_a f(x)\, dx + \beta \int\limits^b_a 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 \ge 0 \Rightarrow \mathcal I(f) \ge 0,$

wobei $f \ge 0$ bedeutet, dass $f(x) \ge 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.

[Bearbeiten] 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^n_{i=0} \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 $σ i$ und $x i$ auch $\sigma^{(n)}_i$ und $x^{(n)}_i$ . 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^n_{i=0} \sigma_i [\alpha f(x_i) + \beta g(x_i)] = \alpha (b - a) \sum^n_{i=0} \sigma_i f(x_i) + \beta (b - a) \sum^n_{i=0} \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 \ge 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 \ge 0 \Rightarrow \mathcal I_n(f) \ge 0.$

Wir definieren weiter:

[Bearbeiten] 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^r_{j=0} a_jx^j \right) = \sum^r_{j=0} a_j \mathcal I_n(x^j) = \sum^r_{j=0} a_j \mathcal I(x^j) = \mathcal I \left( \sum^r_{j=0} a_jx^j \right)$

für alle $a_j \in \mathbb R$ $(j = 0, 1, \ldots, r)$ , und damit der folgende Satz, wobei $Π r$ wieder den Raum aller Polynome vom Grad $\le r$ bezeichnet:

[Bearbeiten] 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:

[Bearbeiten] 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).

[Bearbeiten] Beweis.

Für die durch die Stützstellen $x i$ und Gewichte $σ 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^n_{i=0} \sigma_i L_k(x_i) = (b - a) \sum^n_{i=0} \sigma_i \delta_{ki} = (b - a) \sigma_k = \mathcal I(L_k).$

Da sich jedes Polynom vom Grad $\le 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 \ge 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:

[Bearbeiten] Satz 8.5

Ist $\mathcal I_n$ eine Quadraturformel, die einen Genauigkeitsgrad $r \in \mathbb N_0$ hat, so folgt für ihre Gewichte $σ i$

$\sum^n_{i=0} \sigma_i = 1.$

[Bearbeiten] Beweis.

Da $\mathcal I_n$ einen Genauigkeitsgrad $r \ge 0$ hat, folgt

$(b - a) \sum^n_{i=0} \sigma_i = \mathcal I_n(1) = \mathcal I(1) = \int\limits^b_a 1\, dx = b - a.$

q.e.d.

Bezüglich der Konvergenz der durch eine Quadraturformel

$\mathcal I_n(f) := (b - a) \sum^n_{i=0} \sigma^{(n)}_i f(x^{(n)}_i$

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).

[Bearbeiten] 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^n_{i=1} \left| \sigma^{(n)}_i \right| \le 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:

[Bearbeiten] Korollar 8.7

Es sei $\mathcal I_n$ eine Quadraturformel mit Gewichten $\sigma^{(n)}_i \ge 0\ (i = 0, 1, \ldots, n)$ für alle $n \in \mathbb N$ und einem Genauigkeitsgrad $n \ge 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.$

[Bearbeiten] 8.2 Interpolatorische Quadraturformeln

[Bearbeiten] 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^n_{i=0} f(x_i)L_i(x)$

geschrieben werden. Wir definieren nun:

[Bearbeiten] Denition 8.8

Eine Quadraturformel $\mathcal I_n$ mit

$\mathcal I_n(f) := \mathcal I(Q_n) = \int\limits^b_a Q_n(x)\, dx = (b - a) \sum^n_{i=0} \left\{ \frac{1}{b - a} \int\limits^b_a 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:

[Bearbeiten] 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 $\ge n$ .

Ferner können wir zeigen:

[Bearbeiten] Satz 8.10

Eine interpolatorische Quadraturformel In besitzt die Gestalt

$\mathcal I_n(f) = (b - a) \sum^n_{i=0} \sigma_i f(x_i)$

mit

(8.9) $\sigma_i := \int\limits^1_0 \prod^n_{k=0 \atop k \neq i} \frac{t - t_k}{t_i - t_k}\, dt$ mit $t_k := \frac{x_k - a}{b - a}.$

[Bearbeiten] Beweis.

Mit $t k$ wie in (8.9) lassen sich die Gewichte $σ 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^b_a L_i(x)\, dx = \frac{1}{b - a} \int\limits^b_a \prod^n_{k=0 \atop k \neq i} \frac{x - x_k}{x_i - x_k}\, dx$

(8.10) $= \int\limits^1_0 \prod^n_{k=0 \atop k \neq i} \frac{(b - a)t - (b - a)t_k}{(b - a)t_i - (b - a)t_k}\, dt = \int\limits^1_0 \prod^n_{k=0 \atop k \neq i} \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 $σ 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.

[Bearbeiten] 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 \ge 1 gilt ::(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.

[Bearbeiten] Lemma 8.11

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

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

[Bearbeiten] 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^1_0 \prod^n_{k=0 \atop k \neq i} \frac{t - \frac{k}{n}}{\frac{i}{n} - \frac{k}{n}}\, dt = \frac{1}{n} \int\limits^n_0 \prod^n_{k=0 \atop k \neq i} \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)$

= δ i j = δ 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^b_a L_{n-i}(x)\, dx = \frac{1}{b - a} \int\limits^b_a Q(x)\, dx = \frac{1}{b - a} \int\limits^b_a L_i(b + a - x)\, dx = \frac{1}{b - a} \int\limits^b_a L_i(t)\, dt = \sigma_i.$

q.e.d.

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

[Bearbeiten] 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 $σ 1 = σ 0$ (Lemma 8.11) und $σ 0 + σ 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} ^2_0 \frac{(s - 1) (s - 2)}{(0 - 1) (0 - 2)}\, ds = \frac{1}{4} ^2_0 (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.

[Bearbeiten] Beispiel 8.13

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

$\mathcal I(f) = \int\limits^1_0 \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=1}_{x=0} = \arctan(1) = 0.785\, 40.$

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

$\sum^n_{i=0} \left| \sigma^{(n)}_i \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 \ge 8$ ist daher nicht zu empfehlen. Für die (abgeschlossenen) Newton-Cotes-Formeln lässt sich sogar

$\lim_{n \to \infty} \sum^n_{i=0} \left| \sigma^{(n)}_i \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.

[Bearbeiten] Beispiel 8.14

(1) Für $n = 0$ und $x 0 : = a$ oder $x 0 : = b$ muss wegen Satz 8.5 $σ 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^n_{i=0} \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 $σ 0 = σ 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.

[Bearbeiten] 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.

[Bearbeiten] 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 \ge n$ besitze und es sei $f \in C^{(r+1)}[a, b]$ . Dann gilt

(8.16) $|\mathcal I(f) - \mathcal I_n(f)| \le \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^1_0 \prod^r_{k=0} |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^r_{k=0} (t - t_k)$

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

$\hat \gamma_r := \int\limits^1_0 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).$

[Bearbeiten] 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^n_{i=0} \sigma_i f(x_i) = (b - a) \sum^n_{i=0} \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^b_a [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^b_a \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 \N$ können wir dann Punkte $x^{(m)}_i$ $(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^{(m)}_i = x_i \quad (i = n + 1, \ldots, r)$

gilt. Setzen wir

$\varphi_m(x) := (x - x^{(m)}_{n+1}) \cdots (x - x^{(m)}_r),$

so hat man unter Verwendung des ersten Teils des Beweises

$|\mathcal I(f) - \mathcal I_n(f)| \le \frac{1}{(r + 1)!} \max_{\xi \in [a, b]} \left| f^{(r+1)}(\xi) \right| \int\limits^b_a |\omega(x) \varphi_m(x)|\, dx$

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

(8.20) $\le \frac{1}{(r + 1)!} \max_{\xi \in [a, b]} \left| f^{(r+1)}(\xi) \right| \int\limits^b_a |\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^b_a \prod^r_{i=0} |x - x_i|\, dx = (b - a)^{r+2} \min_{t_{n+1}, \ldots, t_r \in [0, 1]} \int\limits^1_0 \prod^r_{i=0} |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) \ge 0, \quad x \in [a, b],$

so erhält man aus (8.20)

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

Weiter gewinnt man mit (8.19)

$\mathcal I(f) - \mathcal I_n(f) \ge \frac{1}{(r + 1)!} \min_{\eta \in [a, b]} \left| f^{(r+1)}(\eta) \right| \int\limits^b_a [\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^b_a [\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 „ $\le$ “ statt „ $\ge$ “ in (8.21) vorliegt.

q.e.d.

[Bearbeiten] Beispiel 8.16

Wir nutzen im Folgenden aus, dass nach Korollar 8.9 der Genauigkeitsgrad einer interpolatorischen Quadraturformel $\mathcal I_n mindestens <math>r := n ist. (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^0_{k=0} (t - t_k) = t \ge 0, \quad t \in [0, 1]$

und $\hat \gamma_0 := \int\limits^1_0 t\, dt = \frac{1}{2}$ die Fehlerdarstellung

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

wobei $ξ$ 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^0_{k=0} (t - t_k) = t - 1 \le 0, \quad t \in [0, 1]$

und $\hat \gamma_0 := \int\limits^1_0 (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^1_{k=0} (t - t_k) = t (t - 1) \le 0, \quad t \in [0, 1]$

sowie

$\hat \gamma_1 := \int\limits^1_0 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):

[Bearbeiten] Satz 8.17

Die abgeschlossene Newton-Cotes-Formel $\mathcal I_n$ besitzt für gerades $n \ge 2$ den (exakten) Genauigkeitsgrad $r : = n + 1$ .

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

[Bearbeiten] 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^3_{k=0} (t - t_k) = t \left( t - \frac{1}{2} \right)^2 (t - 1) \le 0, \quad t \in [0, 1]$

sowie

$\hat \gamma_3 := \int\limits^1_0 \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).$

[Bearbeiten] 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^b_a f(x)\, dx = \sum^{N-1}_{k=0} \int\limits^{x_{k+1}}_{x_k} f(x)\, dx$

durch

$\mathcal J_n(f) := \sum^{N-1}_{k=0} \int\limits^{x_{k+1}}_{x_k} 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+1}}_{x_k} 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+1}}_{x_k} Q_0(x)\, dx = hf(x_k)$ bzw. $\int\limits^{x_{k+1}}_{x_k} Q_0(x)\, dx = hf(x_{k+1})$

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

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

Für diese gelten die nachstehenden Fehlerabschätzungen.

[Bearbeiten] 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).$

[Bearbeiten] 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+1}}_{x_k} 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^{N-1}_{k=0} f'(\xi_k) = \frac{b - a}{2} h \frac{1}{N} \sum^{N-1}_{k=0} f'(\xi_k).$

Aufgrund von

$N \min_{x \in [a, b]} f'(x) \le \sum^{N-1}_{k=0} f'(\xi_k) \le N \max_{x \in [a, b]} f'(x)$

bzw.

$\min_{x \in [a, b]} f'(x) \le \frac{1}{N} \sum^{N-1}_{k=0} f'(\xi_k) \le \max_{x \in [a, b]} f'(x)$

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

$f'(\xi) = \frac{1}{N} \sum^{N-1}_{k=0} 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+1}}_{x_k} 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^{N-1}_{k=1} f(x_k) + f(b) \right)$

mit der im folgenden Satz angegebenen Fehlerdarstellung.

[Bearbeiten] 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).$

[Bearbeiten] 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+1}}_{x_k} 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^{N-1}_{k=0} f''(\xi_k) = - \frac{b - a}{12} h^2 f''(\xi),$

wobei die Existenz eines solchen $ξ$ 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 + k h$ mit $h : = (b - a) / N$ für jedes $k \ge 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+1}}_{x_k} 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^{N-1}_{k=0} f(x_{k+1/2}) + 2 \sum^{N-1}_{k=1} f(x_k) + f(b) \right)$

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

[Bearbeiten] 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).$

[Bearbeiten] 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+1}}_{x_k} 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^{N-1}_{k=0} f^{(4)}(\xi_k) = - \frac{b - a}{2880} h^4 f^{(4)}(\xi),$

wobei die Existenz von $ξ$ 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 $2 N + 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 \N$ und $N_k \to \infty$ “ gemeint ist.

Wir greifen abschließend nochmals Beispiel 8.13 auf.

[Bearbeiten] 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^1_0 \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.$

[Bearbeiten] 8.4 Extrapolationsverfahren

[Bearbeiten] 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.)

[Bearbeiten] Satz 8.23

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

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

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

(8.23) $\mathcal J_1(h) = \alpha_0 + \alpha_1h^2 + \alpha_2h^{2 \cdot 2} + \ldots + \alpha_rh^{2r} + \mathcal O(h^{2r+2})$ für $h \to 0$

mit $α 0 : = I (f)$ und gewissen Koeffizienten $\alpha_i \in \R$ $(i = 1, \ldots, r)$ .

Für periodische Funktionen mit Periode $b - a$ kann man sogar zeigen, dass $α 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 \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 \R$ $(i = 0, 1, \ldots, r)$ und einer Zahl $γ > 0$ die asymptotische Entwicklung der Ordnung $r$

(8.24) $T(h) = \alpha_0 + \alpha_1h^\gamma + \alpha_2h^{2\gamma} + \ldots + \alpha_rh^{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 $α 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 \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 $α 0$ erzeugt. Setzt man $T (0): = α 0$ , so extrapoliert dieses Verfahren also $T$ auf den Wert $h = 0$ hin, so dass man auch von einem Extrapolationsverfahren spricht. Da die Koeffizienten $α 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^\gamma_k, T(h_k)), k = 0, 1, \ldots, n$ (schreibt man $T (h γ)$ statt $T (h)$ , so sind dies mit $z_k := h^\gamma_k$ 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 $α 0$ .

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

[Bearbeiten] 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^\gamma_k, 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^\gamma_k) = 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^\gamma_j T_{j+1, \ldots, j+m} + h^\gamma_{j+m} T_{j, \ldots, j+m-1}}{h^\gamma_{j+m} - h^\gamma_j} = T_{j+1, \ldots, j+m} - h^\gamma_{j+m} \frac{T_{j+1, \ldots, j+m} - T_{j, \ldots, j+m-1}}{h^\gamma_{j+m} - h^\gamma_j}$

(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 \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}$

[Bearbeiten] Beispiel 8.24

Für die summierte Trapezregel $\mathcal J_1(h)$ gilt gemäß Satz 8.23 eine Entwicklung der Form (8.24) mit $γ = 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.

[Bearbeiten] 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^m_{k=0} L_k(0)x^j_k = \begin{cases} 1 & f\ddot ur\ j = 0, \\ 0 & f\ddot ur\ 1 \le j \le m, \\ (-1)^mx_0x_1 \cdots x_m & f\ddot ur\ j = m + 1. \end{cases}$

[Bearbeiten] Beweis.

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

$p(x) = x^j = \sum^m_{k=0} x^j_k 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^m_{k=0} x^{m+1}_k 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^m_{k=0} L_k(0) x^{m+1}_k = - q(0) = (-1)^m x_0 x_1 \cdots x_m.$

[Bearbeiten] Satz 8.26

Es sei $T (h)$ mit $h > 0$ eine Funktion mit der asymptotischen Entwicklung (8.24) für ein $γ > 0$ und $r \in \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 \N_0)$ mit $1 = n_0 \le n_1 \le n_2 \le \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 \le m \le r - 1$ der asymptotischen Entwicklung

(8.32) $T_{j, \ldots, j+m} - \alpha_0 = (-1)^m \frac{\alpha_{m+1}}{n^\gamma_j \cdots n^\gamma_{j+m}} h^{(m+1)\gamma}_0 + \mathcal O(h^{(m+2)\gamma}_0).$

[Bearbeiten] 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^m_{k=0} T(h_k) L_k(h^\gamma) = \sum^m_{k=0} T(h_k) \left[ \prod^m_{j=0 \atop j \neq k} \frac{h - h^\gamma_j}{h^\gamma_k - h^\gamma_j} \right], \quad h \in \R$

und somit

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

für

$c_{m,k} := L_k(0) = \prod^m_{j=0 \atop j \neq k} \frac{h^\gamma_j}{h^\gamma_j - h^\gamma_k}.$

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

(8.34) $T(h_k) = \sum^{m+1}_{s=0} \alpha_sh^{s \gamma}_k + \mathcal O(h^{(m+2)\gamma}_k).$

Desweiteren schließt man mit Lemma 8.25

(8.35) $\sum^m_{k=0} c_{m, k} h^{s \gamma}_k = \sum^m_{k=0} L_k(0) h^{s \gamma}_k = \begin{cases} 1 & \mbox{für } s = 0, \\ 0 & \mbox{für } 1 \le s \le m, \\ (-1)^m h^\gamma_0 h^\gamma_1 \cdots h^\gamma_m & \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^m_{k=0} c_{m, k} \left[ \sum^{m+1}_{s=0} \alpha_s h^{s \gamma}_k + \mathcal O(h^{(m+2)\gamma}_k) \right] = \sum^{m+1}_{s=0} \alpha_s \left[ \sum^m_{k=0} c_{m, k} h^{s \gamma}_k \right] + \sum^m_{k=0} c_{m, k} \mathcal O(h^{(m+2)\gamma}_k)$

$= \alpha_0 + (-1)^m \alpha_{m+1} h^\gamma_0 h^\gamma_1 \cdots h^\gamma_m + \mathcal O(h^{(m+2)\gamma}_0).$

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 $γ$ gewinnt. Diese Sichtweise ist allerdings zu optimistisch, da die Restterme der asymptotischen Entwicklung, die sich hinter $\mathcal O(h^{(m+2)\gamma}_0)$ verbergen, nicht bekannt sind und groß werden können.

Es bietet sich also der folgende Algorithmus an:

[Bearbeiten] 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

α 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 \N),$

die durch

$n_1 := 2, \quad n_2 := 3, \quad n_3 := 4, \quad n_j := 2n_{j-2} \quad (j \ge 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 \N).$

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

[Bearbeiten] 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 \le m \le r - 1$

$T_{0, \ldots, m} - \alpha_0 = \left[ \frac{(-1)^m}{2^{\gamma m(m+1)/2}} \alpha_{m+1} \right] h^{(m+1)\gamma}_0 + \mathcal O(h^{(m+2)\gamma}_0).$

[Bearbeiten] Beweis.

Im Fall (8.36) hat man mit $n 0 = 1$

$n^\gamma_0 \cdots n^\gamma_m = 2^{0\gamma} \cdot 2^{1\gamma} \cdot 2^{2\gamma} \cdots 2^{m\gamma} = 2^{\gamma \sum^m_{k=0} 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.

[Bearbeiten] Beispiel 8.28

(1) Korollar 8.27 wollen wir auf die summierte Trapezregel mit $m : = 2$ (und wegen der Forderung $m \le r - 1$ für $r = 3$ ) anwenden. Nach Satz 8.23 ist dann $γ = 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\cdot2(2+1)/2}} \alpha_{2+1} \right] h^{(2+1)2}_0 + \mathcal O(h^{(m+2)2}_0) = \frac{\alpha_3}{64} h^6_0 + \mathcal O(h^8_0),$

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 \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^2_0)$ :

$T_2 - \mathcal I(f) = \mathcal J_1(h_0/4) - \mathcal I(f) = \frac{b - a}{12 \cdot 4^2} h^2_0 f''(\xi) = \frac{(b - a) f''(\xi)/3}{64} h^2_0.$

(2) (Bader) Es soll das Integral

$\int\limits^{\pi/2}_0 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 : = π / 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 $α 0$ konvergieren.) Für eine Diskussion über ein geeignetes Abbruchkriterium verweisen wir auf Deuflhard/Hohmann.

[Bearbeiten] 8.5 Gaußsche Quadraturformeln

[Bearbeiten] 8.5.1 Einleitung

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

(8.38) $\mathcal I(f) := \int\limits^b_a 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 \R$ eine Gewichtsfunktion mit den folgenden Eigenschaften ist:

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

$\mu_k := \int\limits^b_a 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:

Parser-Fehler (die PNG-Konvertierung schlug fehl): \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 $Π$ das durch $w$ induzierte Skalarprodukt

(8.39) $\langle f, g \rangle := \int\limits^b_a 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 \ge 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 \R.$

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

(8.40) $\|f\| := \langle f, f \rangle^{1/2} = \left\{ \int\limits^b_a 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^n_{i=0} \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 $σ 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 $σ i$ zur Verfügung hat, was der Anzahl der Koeffizienten eines Polynoms vom Grad $2 n + 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 $σ i$ $(i = 0, 1, \ldots, n)$ und hätte $\mathcal I_n$ den Genauigkeitsgrad $2 n + 2$ , so folgte insbesondere für das Polynom

$p(x) := \prod^n_{i=0} (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:

[Bearbeiten] Lemma 8.29

Ist $\mathcal I_n$ eine Quadraturformel mit Stützstellen $x i$ $(i = 0, 1, \ldots, n)$ und Gewichten $σ i$ $(i = 0, 1, \ldots, n)$ und hat $\mathcal I_n$ den Genauigkeitsgrad $2 n + 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 \R \setminus \{0\}$ .

[Bearbeiten] 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^b_a w(x)p(x)p_{n+1}(x)\, dx = \mathcal I_n(pp_{n+1}) = \sum^n_{i=0} \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 $2 n + 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 $Π 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^n_{j=0} a_jp_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^n_{j=0} a_jp_j(x), p_{n+1} \right\rangle = \sum^n_{j=0} 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 $σ 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^b_a 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.

[Bearbeiten] 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^{n-1}_{j=0} \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^{n-1}_{j=0} \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 \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.

[Bearbeiten] 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^2_n 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^2_n := \frac{\|p_n\|^2}{\|p_{n-1}\|^2} \quad (n = 1, 2, \ldots).$

[Bearbeiten] Beweis.

Offenbar ist die behauptete Darstellung richtig für $p 0$ und $p 1$ . Für $n \ge 1$ setzen wir

$q_{n+1} := (x - \beta_n)p_n - \gamma^2_n 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 $Π 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 $β n$

(8.48) $\langle q_{n+1}, p_n \rangle = \langle (x - \beta_n)p_n - \gamma^2_n p_{n-1}, p_n \rangle = \langle xp_n, p_n \rangle - \beta_n \|p_n\|^2 = 0$

und mit der Definition von $γ n$

$\langle q_{n+1}, p_{n-1} \rangle = \langle (x - \beta_n)p_n - \gamma^2_n p_{n-1}, p_{n-1} \rangle = \langle xp_n, p_{n-1} \rangle - \gamma^2_n \|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^2_n \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^n_{j=0} a_j p_j(x) \right\rangle = \sum^n_{j=0} 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 \N)$ in (8.43) hat man folgende Aussage:

[Bearbeiten] 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.

[Bearbeiten] Beweis.

Es seien die Nullstellen $x i$ von $p n$ so durchnumeriert, dass $a < x_0 < \ldots < x_{j-1} < b$ $(0 \le j \le 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 \le n - 2$ bzw. $j \le n - 1$ , so hätte das Polynom

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

den Grad $j \le 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^{j-1}_{k=0} (x - x_k)^2 \right) \left( \prod^{n-1}_{k=j} (x - x_k) \right) \ge 0, \quad x \in [a, b]$

und demzufolge

$\langle p_n, q \rangle = \int\limits^b_a 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 \ge 1$ in der Form

$p_n(x) = (x - x_i)\hat q(x), \quad \hat q(x) := \prod^{n-1}_{k=0 \atop k \neq i} (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 x L_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:

Parser-Fehler (die PNG-Konvertierung schlug fehl): \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):

[Bearbeiten] Satz 8.32

Für $p n$ aus (8.43) gilt

$\|p_n\| = \min_{a_0, \ldots, a_{n-1} \in \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^b_a |f(x)|\, dx, \quad f \in \Pi.$

[Bearbeiten] 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^2_n(x)\, dx}{\int^1_{-1} p^2_n(x)\, dx}, \quad \gamma^2_n = \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^2_1 = \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^2_1 = \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.$

[Bearbeiten] 8.5.3 Die Quadraturformeln

[Bearbeiten] Satz 8.34

Für $n \in \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 $σ 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^n_{i=0} \sigma_i f(x_i)$

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

[Bearbeiten] 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 $2 n + 1$ und damit exakt den Genauigkeitsgrad $2 n + 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 = q p 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^n_{i=0} r(x_i)L_i(x) = \sum^n_{i=0} 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^n_{i=0} p(x_i) \langle L_i, \mathbf 1 \rangle = \sum^n_{i=0} \sigma_i p(x_i),$

womit der Genauigkeitsgrad von mindestens $2 n + 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 $σ i$ bezeichnet man als Gaußsche Quadraturformel. Ihr Genauigkeitsgrad ist optimal, da es keine Quadraturformeln mit Genauigkeitsgrad $2 n + 2$ gibt (vgl. Abschnitt 8.5.1). Weiter hat man:

[Bearbeiten] Lemma 8.35

Für die Gewichte $σ 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^n_{i=0} \sigma_i = \int\limits^b_a w(x)\, dx.$

[Bearbeiten] Beweis.

Wendet man die Beziehungen (8.54) auf $p := L^2_j \in \Pi_{2n}$ an, so folgt

$\langle L_j, L_j \rangle = \left\langle L^2_j, \mathbf 1 \right\rangle = \sum^n_{i=0} \sigma_i L^2_j(x_i) = \sigma_j.$

Weiter gilt $\langle L_j, L_j \rangle > 0, da man <math>L^2_j(x) \ge 0, x \in [a, b]$ sowie $L^2_j(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^n_{i=0} \sigma_i = \mathcal I_n(1) = \mathcal I(1) = \int\limits^b_a 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 \N_0$ positive Gewichte. Mit dem folgenden Satz geben wir abschließend eine Darstellung für den bei der Gauß-Quadratur entstehenden Quadraturfehler an.

[Bearbeiten] 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^1_0 w((b - a)t + a)^2 \prod^n_{k=0} (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).$

[Bearbeiten] 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 : = 2 n + 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^{2n+1}_{k=0} (t - t_k) = \prod^n_{k=0} (t - t_k)^2 = p^2_{n+1}(x) \ge 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.

[Bearbeiten] 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 $β j$ und $γ 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^2_k 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 \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 $σ i$ aus zugehörigen Eigenvektoren bestimmen lassen.

[Bearbeiten] Satz 8.37

Für die Stützstellen $x i$ $(i = 0, 1, \ldots, n)$ und die Gewichte $σ 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_iv^{(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^n_{k=0} \tau^2_k p^2_k(x_i)} \quad (i = 0, 1, \ldots, n).$

[Bearbeiten] Beweis.

Aus den Definitionen von $J, v (i), p 1$ und $τ 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_iv^{(i)}_1.$

Definiert man $γ n + 1 = τ 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^{(i)}_{k+1}, \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^n_{j=0} \sigma_j p_k(x_j) \quad (k = 0, 1, \ldots, n),$

wobei $δ i j$ das Kroneckersymbol ist. Multiplikation von (8.62) mit $\tau^2_k 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^n_{k=0} \tau^2_k p_k(x_i) \delta_{k0} = \langle \mathbf 1, \mathbf 1 \rangle \tau^2_0 p_0(x_i) = \langle \mathbf 1, \mathbf 1 \rangle = \sum^n_{k=0} \sum^n_{j=0} \sigma_j \tau^2_k p_k(x_i) p_k(x_j)$

$= \sum^n_{j=0} \sigma_j \sum^n_{k=0} \tau^2_k p_k(x_i) p_k(x_j) = \sum^n_{j=0} \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.

[Bearbeiten] 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^2_0 p^2_0(x_0) + \tau^2_1 p^2_1(x_0) + \tau^2_2 p^2_2(x_0)} = \frac{5}{9},$

$\sigma_1 = \frac{2}{\tau^2_0 p^2_0(x_1) + \tau^2_1 p^2_1(x_1) + \tau^2_2 p^2_2(x_1)} = \frac{8}{9},$

$\sigma_2 = \frac{2}{\tau^2_0 p^2_0(x_2) + \tau^2_1 p^2_1(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.

Kurs:Numerik I/8 Numerische Integration

Aus Wikiversity

Inhaltsverzeichnis

[Bearbeiten] 8.1 Einführung

[Bearbeiten] Definition 8.1

[Bearbeiten] Definition 8.2

[Bearbeiten] Satz 8.3

[Bearbeiten] Satz 8.4

[Bearbeiten] Beweis.

[Bearbeiten] Satz 8.5

[Bearbeiten] Beweis.

[Bearbeiten] Satz 8.6 (Szegö)

[Bearbeiten] Korollar 8.7

[Bearbeiten] 8.2 Interpolatorische Quadraturformeln

[Bearbeiten] 8.2.1 Allgemeines

[Bearbeiten] Denition 8.8

[Bearbeiten] Korollar 8.9

[Bearbeiten] Satz 8.10

[Bearbeiten] Beweis.

[Bearbeiten] 8.2.2 Newton-Cotes-Formeln

[Bearbeiten] Lemma 8.11

[Bearbeiten] Beweis.

[Bearbeiten] Beispiel 8.12

[Bearbeiten] Beispiel 8.13

[Bearbeiten] Beispiel 8.14

[Bearbeiten] 8.2.3 Quadraturfehler und Genauigkeitsgrad

[Bearbeiten] Satz 8.15

[Bearbeiten] Beweis.

[Bearbeiten] Beispiel 8.16

[Bearbeiten] Satz 8.17

[Bearbeiten] Beispiel 8.18

[Bearbeiten] 8.3 Summierte abgeschlossene Newton-Cotes-Formel

[Bearbeiten] Satz 8.19

[Bearbeiten] Beweis.

[Bearbeiten] Satz 8.20

[Bearbeiten] Beweis.

[Bearbeiten] Satz 8.21

[Bearbeiten] Beweis.

[Bearbeiten] Beispiel 8.22

[Bearbeiten] 8.4 Extrapolationsverfahren

[Bearbeiten] 8.4.1 Einführung

[Bearbeiten] Satz 8.23

[Bearbeiten] 8.4.2 Das Verfahren

[Bearbeiten] Beispiel 8.24

[Bearbeiten] Lemma 8.25

[Bearbeiten] Beweis.

[Bearbeiten] Satz 8.26

[Bearbeiten] Beweis.

[Bearbeiten] Algorithmus 10 (Extrapolationsverfahren)

[Bearbeiten] Korollar 8.27

[Bearbeiten] Beweis.

[Bearbeiten] Beispiel 8.28

[Bearbeiten] 8.5 Gaußsche Quadraturformeln

[Bearbeiten] 8.5.1 Einleitung

[Bearbeiten] Lemma 8.29

[Bearbeiten] Beweis.

[Bearbeiten] 8.5.2 Orthogonale Polynome

[Bearbeiten] Satz 8.30

[Bearbeiten] Beweis.

[Bearbeiten] Satz 8.31

[Bearbeiten] Beweis.

[Bearbeiten] Satz 8.32

[Bearbeiten] Beispiel 8.33

[Bearbeiten] 8.5.3 Die Quadraturformeln

[Bearbeiten] Satz 8.34

[Bearbeiten] Beweis.

[Bearbeiten] Lemma 8.35

[Bearbeiten] Beweis.

[Bearbeiten] Satz 8.36

[Bearbeiten] Beweis.

[Bearbeiten] 8.5.4 Berechnung der Stützstellen und Gewichte

[Bearbeiten] Satz 8.37

[Bearbeiten] Beweis.

[Bearbeiten] Beispiel 8.38

Ansichten

Persönliche Werkzeuge

Suche

Navigation

Aktuell

Werkzeuge