6.3. Gauss Quadratur

Zuletzt wollen wir uns mit einem weiteren Verfahren für die numerische Integration beschäftigen, bei der die Stützstellen nicht äquidistant gewählt werden, sondern so, dass man einen möglichst hohen Exaktheitsgrad erzielt (siehe Bemerkung 6.2). Hierbei verzichtet man bei der Diskretisierung des Intervalls \([a,b] \subset \R\) also inbesondere auf eine fixe Schrittweite \(h = \frac{b-a}{n}\).

In diesem Abschnitt beschäftigen wir uns, anders als zuvor, mit einem gewichteten Integral der Form

\[I(f) \ \coloneqq \ \int_a^b f(x) \sigma(x)\,\mathrm{d}x,\]

wobei \(\sigma \colon [a,b] \rightarrow \R_{\geq 0}\) eine gegebene nichtnegative Gewichtsfunktion ist. Wir lassen hierbei auch unendliche Integrale wie \([0, \infty]\) oder \([-\infty, \infty]\) zu. Die Gewichtsfunktion \(\sigma\) soll folgende Voraussetzungen erfüllen:

1. \(\sigma(x)\) ist für alle \(x \in [a,b]\) nichtnegativ und messbar,

2. Alle Momente der Form

   \[\mu_k \ \coloneqq \ \int_a^b x^k \sigma(x) \, \mathrm{d}x, \qquad k=0,1,\ldots\]

   existieren und sind endlich,

3. Für jedes Polynom \(p(x)\) mit \(\int_a^b \sigma(x) p(x) \mathrm{d}x = 0\) und \(p(x) \geq 0\) für alle \(x \in [a,b]\) gilt schon \(p(x) \equiv 0\).

Man kann zeigen, dass jede Gewichtsfunktion \(\sigma\), die auf dem Intervall \([a,b]\) positiv und stetig ist, die obigen Voraussetzungen erfüllt. Für den Spezialfall \(\sigma(x) \equiv 1\) für alle \(x\in [a,b]\) erhalten wir das herkömmliche bestimmte Integral.

Wir betrachten zunächst wieder Quadraturformeln aus Definition 6.1 der Form

(6.15)\[G_n(f) \ = \ \sum_{i=1}^n w_i f(x_i) \ \approx \ \int_a^b f(x) \sigma(x)\, \mathrm{d}x \ = \ I(f),\]

mit Gewichten \(w_i \in \R, i=1,\ldots,n\). Im Gegensatz zu den bisher betrachteten Quadraturformeln nehmen wir hier jedoch keine äquidistanten Stützstellen \(x_i \in [a,b], i=1, \ldots, n\) an, so dass diese frei gewählt werden können. Somit hat die Quadraturformel (6.15) insgesamt \(2n\) Freiheitsgrade.

Nach Bemerkung 6.2 wissen wir bereits, dass die interpolatorischen Quadraturformeln nach Newton-Cotes mit \(n \in \N_+\) Stützstellen ein Polynom \(n\)ten Grades exakt integrieren. Hierdurch motiviert wollen wir in diesem Abschnitt spezielle Quadraturformeln \(G_n\), genannt Gauss-Quadraturen, herleiten, die einen Exaktheitsgrad von \(2n-1\) besitzen, d.h., für alle Polynome \(p \in \Pi_{2n-1}\) soll folgende Bedingung gelten:

\[G_n(p) \ = \ I(p).\]

Dies führt zu genau \(2n\) Bedingungen für unsere \(2n\) Freiheitsgrade in (6.15).

Der folgende Satz besagt, dass diese Forderung schon maximal ist.

Satz 6.4 (Maximaler Exaktheitsgrad für Quadraturformeln)

Es gibt keine Quadraturformel \(G_n\) wie in (6.15) mit \(2n\) Freiheitsgraden, die vom Exaktheitsgrad \(2n\) ist.

Proof. Nehmen wir an \(G_n\) sei eine Quadraturformel mit \(2n\) Freiheitsgraden und Exaktheitsgrad \(2n\). Dann gilt offensichtlich \(G_n(p) = I(p)\) für alle Polynome \(p \in \Pi_{2n}\). Dies würde also auch gelten für das spezielle Polynom

\[p(x) \ \coloneqq \ \prod_{i=1}^n (x - x_i)^2,\]

für die Stützstellen \(x_i \in [a,b], i=1,\ldots,n\).

Im Allgemeinen gilt offensichtlich \(I(p) \neq 0\), jedoch ist

\[G_n(p) \ = \ \sum_{i=1}^n w_i p(x_i) \ = \ \sum_{i=1}^n w_i \prod_{j=1}^n (x_i - x_j)^2 \ = \ 0.\]

Damit ist also \(G_n(p) \neq I(p)\) und somit kann es keine Quadraturformel der Form (6.15) mit \(2n\) Freiheitsgraden geben, die Exaktheitsgrad \(2n\) besitzt. ◻

Es stellt sich nun die Frage, wie man eine Quadraturformel \(G_n\) finden kann, die den Exaktheitsgrad \(2n-1\) besitzt. Eine theoretische Grundlage für die Beantwortung dieser Frage liefert uns das folgende Theorem.

Satz 6.5 (Exaktheitsgrad einer Quadraturformel)

Seien \(n, m \in \N_+\). Eine Quadraturformel \(G_n\) der Form (6.15) für das gewichtete Integral hat den Exaktheitsgrad \((n+m - 1) \in \N_+\), genau dann wenn sie die folgenden Bedingungen erfüllt:

1. \(G_n\) ist eine interpolatorische Quadraturformel,

2. das zu den \(n\) Stützstellen \(x_i \in [a,b], i=1, \ldots, n\) gehörende Polynom \(\omega_n \in \Pi_n\) der Form \(\omega_{n}(x) = \prod_{i=1}^n (x-x_i)\) erfüllt die folgende Orthogonalitätsbedingung

   (6.16)\[\int_a^b \omega_{n}(x) p(x) \sigma(x)\,\mathrm{d} x \ = \ 0\]

   für alle Polynome \(p \in \Pi_{m-1}\) vom Grad kleiner oder gleich \(m-1\).

Proof. Wir zeigen beide Richtungen der Aussage des Theorems getrennt.
\(\Rightarrow\): Sei zunächst \(G_n(f) = \sum_{i=1}^n w_i f(x_i)\) eine Quadraturformel vom Exaktheitsgrad \(n+m-1\) mit den Stützstellen \(x_i \in [a,b], i=1,\ldots,n\). Wir betrachten nun ein Polynom \(P(x) \coloneqq \omega_n(x) p(x)\) für \(\omega_n \in \Pi_n\) mit \(w_n(x) \coloneqq \prod_{i=1}^n (x - x_i)\) und \(p \in \Pi_{m-1}\) sei ein beliebiges Polynom vom Grad kleiner oder gleich \(m-1\). Dann ist klar, dass \(P(x)\) ein Polynom vom Grad kleiner gleich \(n+m-1\) ist. Da wir angenommen hatten, dass die Quadraturformel den Exaktheitsgrad \(n+m-1\) besitzt gilt per Definition \(I(P) = G_n(P)\) und somit

\[G_n(P) \ = \ I(P) \ = \ \int_a^b P(x) \sigma(x) \,\mathrm{d}x.\]

Wir sehen also, dass \(G_n\) eine interpolatorische Quadratur ist. Außerdem gilt:

\[\int_a^b \omega_n(x) p(x) \sigma(x) \, \mathrm{d}x \ = \ I(P) \ = \ G_n(P) \ = \ \sum_{i=1}^n w_i P(x_i) \ = \ \sum_{i=1}^n w_i \cdot \omega_{n}(x_i) p(x_i) \ = \ 0,\]

da die \(x_i \in [a,b], i=1, \ldots, n\) gerade die Nullstellen von \(\omega_{n}\) sind. Also erfüllt das Polynom \(\omega_n \in \Pi_n\) die Orthogonalitätsbedingungen für alle Polynome \(p \in \Pi_{m-1}\).  
\(\Leftarrow\): Wir nehmen nun an, dass \(G_n\) eine interpolatorische Quadraturformel ist und das Polynom \(\omega_n \in \Pi_n\) die Orthogonalitätsbedingung (6.16) für alle Polynome \(p \in \Pi_{m-1}\) erfüllt.

Sei \(p \in \Pi_n\) ein beliebiges Polynom vom Grad \(n+m-1\), so können wir \(p\) mittels Polynomdivision in der Form

\[p \ = \ \omega_{n} q_1 + q_2\]

schreiben, mit einem Polynom \(q_1 \in \Pi_{m-1}\) vom Grad kleiner oder gleich \(m-1\) und einem Restpolynom \(q_2 \in \Pi_{n-1}\) vom Grad kleiner oder gleich \(n-1\). Da die interpolatorische Quadraturformel \(G_n\) Exaktheitsgrad \(n\) besitzt und somit Polynome vom Grad kleiner oder gleich \(n\) exakt integriert, folgt

\[G_n(q_2) \ = \ \sum_{i=1}^n w_i q_2(x_i) \ = \ \int_a^b q_2(x) \sigma(x)\,\mathrm{d} x \ = \ \int_a^b (p(x) - \omega_{n} q_1(x)) \sigma(x) \,\mathrm{d}x.\]

Da wir die Orthogonalitätsbedingung (6.16) angenommen haben folgt

\[\begin{split} G_n(p) \ &= \ \sum_{i=1}^n w_i p(x_i) \ = \ \sum_{i=1}^n w_i (\omega_n(x_i) q_1(x_i) + q_2(x_i)) \ = \ \underbrace{\sum_{i=1}^n w_i \omega_n(x_i) q_1(x_i)}_{= \, 0} + \sum_{i=1}^n w_i q_2(x_i)\\ &= \ \sum_{i=1}^n w_i q_2(x_i) \ = \ \int_a^b (p(x) - \omega_{n}(x) q_1(x)) \sigma(x) \,\mathrm{d}x \\ &= \ \int_a^b p(x) \sigma(x) \,\mathrm{d}x - \underbrace{\int_a^b \omega_{n}(x) q_1(x) \sigma(x) \,\mathrm{d}x}_{= \, 0} \ = \ \int_a^b p(x) \sigma(x) \,\mathrm{d}x. \end{split}\]

Dies bedeutet, dass die Quadraturformel \(G_n\) den Exaktheitsgrad \(n+m-1\) hat. ◻

Das Satz 6.5 sagt uns also, dass wir für einen Exaktheitsgrad \(n+m-1\) die Nullstellen \(x_i \in [a,b], i=1,\ldots,n\) des Polynoms \(\omega_n\) (die die Stützstellen der Quadraturformel \(G_n\) sind) so geschickt wählen müssen, dass die Orthogonalitätsbedingung für alle Polynome \(p \in \Pi_{m-1}\) vom Grad kleiner oder gleich \(m-1\) erfüllt wird. Wir wollen natürlich den nach Satz 6.4 maximal möglichen Exaktheitsgrad von \(2n-1\) erreichen, d.h. wir setzen \(m=n\) in Satz 6.5. Dann muss die Orthogonalitätsbedingung (6.16) entsprechend für alle Polynome \(p \in \Pi_n\) vom Grad \(n-1\) gelten.

Um die Gausschen Quadraturformeln herzuleiten wollen wir zunächst noch einige grundlegende Eigenschaften von Orthogonalpolynomen wiederholen.

Definition 6.3 (Menge der normierten Polynome)

Sei \(\Pi_n\) der Vektorraum aller reellen Polynome vom Grad kleiner oder gleich \(n \in \N\). Dann definieren wir durch

\[\overline{\Pi}_n \ \coloneqq \ \lbrace p \in \Pi_n \, | \, p(x) = x^n + a_1 x^{n-1} + \ldots + a^{n-1}x + a_n \rbrace\]

die Menge aller normierten, reellen Polynome vom Grad \(n \in \N\).

Weiterhin betrachten wir den Raum \(L^2([a,b])\) aller quadratisch integrierbaren Funktionen versehen mit dem Skalarprodukt

\[\langle f, g \rangle_\sigma \ \coloneqq \ \int_a^b f(x) g(x) \sigma(x) \, \mathrm{d}x,\]

wobei \(\sigma \colon [a,b] \rightarrow \R_+\) eine stetige, positive Gewichtsfunktion ist.

Mit Hilfe dieses Skalarprodukts können wir nun orthogonale Polynome auf \(L^2[a,b]\) definieren.

Definition 6.4 (Orthogonale Funktionen)

Seien \(f,g \in L^2([a,b])\) quadratisch integrierbare Funktionen für ein Intervall \([a,b] \subset \R\). Wir nennen \(f\) und \(g\) orthogonal, falls \(\langle f, g \rangle_\sigma = 0\) gilt.

Das folgende Theorem beweist, dass es eindeutig bestimmte Orthogonalpolynome bezüglich einer gewählten Gewichtsfunktion \(\sigma\) gibt.

Satz 6.6 (Existenz und Eindeutigkeit von Orthogonalpolynomen)

Es gibt für \(n=0,1,\ldots\) eindeutig bestimmte Polynome \(p_n \in \overline{\Pi}_n\) für die die folgende Orthogonalitätseigenschaft gilt:

\[\langle p_n, p_m \rangle_\sigma \ = \ 0, \qquad \text{ für } \quad n \neq m.\]

Diese Orthogonalpolynome genügen der Rekursionformel

(6.17)\[\begin{split} p_0(x) \ &\equiv \ 1, \\ p_{n+1}(x) \ &= \ (x - \delta_{n+1})p_n(x) - \gamma^2_{n+1} p_{n-1}(x) \qquad \text{ für } \quad n \in \N, \end{split}\]

wobei wir \(p_{-1}(x) \equiv 0\) definieren und

\[\begin{split} \delta_{n+1} \ &\coloneqq \ \langle x p_n, p_n \rangle_\sigma \, / \, \langle p_n, p_n \rangle_\sigma \qquad \quad \, \ \text{ für } n \in \N,\\ \gamma^2_{n+1} \ &\coloneqq \ \begin{cases} \ 1, \ &\text{ für } n = 0,\\ \ \langle p_n, p_n \rangle_\sigma \, / \, \langle p_{n-1}, p_{n-1} \rangle_\sigma \ &\text{ für } n \in \N_+. \end{cases} \end{split}\]

Proof. Siehe [FH07]. ◻

Die Orthogonalpolynome \(p_n \in \Pi_n\) erfüllen bereits die gewünschte Orthogonalitätseigenschaft aus (6.16), wie das folgende Korollar aussagt.

Korollar 6.4

Sei \(p_n \in \overline{\Pi}_n\) ein Orthogonalpolynom. Dann gilt \(\langle p_n, p \rangle_\sigma = 0\) für alle Polynome \(p \in \Pi_{n-1}\).

Proof. Da die Orthogonalpolynome \(p_n \in \overline{\Pi}_n\) eine Basis des Vektorraums der Polynome vom Grad kleiner oder gleich \(n \in \N\) bilden, lässt sich jedes Polynom \(p \in \Pi_{n-1}\) als Linearkombination von Orthogonalpolynomen \(p_0, \ldots, p_{n-1} \in \overline{\Pi}_{n-1}\) darstellen. Daher gilt schon:

\[\langle p_n, p \rangle_\sigma \ = \ \langle p_n, \sum_{k=1}^{n-1} \alpha_k p_k \rangle_\sigma \ = \ \sum_{k=1}^{n-1} \alpha_k \underbrace{\langle p_n, p_k \rangle_\sigma}_{= \, 0} \ = \ 0.\]

◻

Folgendes Resultat beschreibt die Eigenschaften der Nullstellen der eindeutig bestimmten Orthogonalpolynome \(p_n \in \overline{\Pi}_n\).

Satz 6.7 (Nullstellen der Orthogonalpolynome)

Die Nullstellen \(x_i, i=1,\ldots, n\) des Orthogonalpolynoms \(p_n \in \overline{\Pi}_n\) sind reell, einfach und liegen alle im offenen Intervall \((a,b)\).

Proof. Siehe [FH07]. ◻

Nach Satz 6.6 kann ein System von Orthogonalpolynomen bezüglich der Gewichtungsfunktion \(\sigma\) schrittweise durch die Rekursionsformel (6.17) mittels des Gram-Schmidt-Verfahrens berechnet werden. Mit Hilfe der Nullstellen des \(n\)-ten Orthogonalpolynoms \(\omega_n \in \overline{\Pi}_n\) lassen sich so die optimalen Stützstellen für die Gauss-Quadratur \(G_n\) finden. Das Vorgehen ist im folgenden Algorithmus nochmal zusammengefasst.

Algorithmus 6.1 (Gausssches Integrationsverfahren)

Sei zunächst ein Intervall \([a,b] \subset \R\) und eine stetige, positive Gewichtungsfunktion \(\sigma \colon [a,b] \rightarrow \R_+\) gegeben. Für die numerische Approximation des bestimmten, gewichteten Integrals \(I(f) = \int_a^b f(x) \sigma(x) \mathrm{d}x\) führen wir folgende Schritte durch:

1. Berechne entsprechend zur Gewichtsfunktion \(\sigma\) das Orthogonalpolynom \(\omega_{n} \in \overline{\Pi}_n\) vom Grad \(n\).

2. Berechne die Nullstellen \(x_i \in (a,b), i=1, \ldots, n\) von \(\omega_{n}\).

3. Verwende die bestimmten Nullstellen von \(\omega_n\) als Stützstellen für die Integrationsformeln nach Newton-Cotes (siehe Integrationsformeln nach Newton-Cotes)

   \[G_n(f) \ = \ \sum_{i=1}^n w_i f(x_i) , \qquad w_k = \int_a^b L_k(x) \sigma(x) \,\mathrm{d} x.\]

Je nach Wahl der Gewichtsfunktion \(\sigma \colon [a,b] \rightarrow \R_{\geq 0}\) und des Intervalls \([a,b] \subset \R\) erhält man ein anderes Orthogonalsystem von normierten Polynomen, wie die folgende Bemerkung zusammenfasst.

Bemerkung 6.9 (Unterschiedliche Systeme von Orthogonalpolynomen)

Die folgende Tabelle fasst die wichtigsten Gewichtsfunktionen \(\sigma\), die zugehörigen Intervalle \([a,b] \subset \R\) und die daraus resultierenden Systeme von Orthogonalpolynomen zusammen.

\([a,b]\)	\(\sigma\)	Bezeichnung der Orthogonalpolynome
\([-1,1]\)	\(1\)	Legendre-Polynome
\([-1,1]\)	\((1-x^2)^{-\frac{1}{2}}\)	Tschebyscheff-Polynome 1. Art
\([-1,1]\)	\((1-x^2)^{\frac{1}{2}}\)	Tschebyscheff-Polynome 2. Art
\([-1,1]\)	\((1-x)^\alpha(1+x)^\beta\)	Jacobi-Polynome
\((-\infty, \infty)\)	\(e^{-\frac{x^2}{2}}\)	Hermite’sche Polynome
\((0, \infty)\)	\(e^{-x}\)	Laguerre’sche Polynome

Das folgende Beispiel bestimmt die optimalen Stützstellen für eine Gauss-Quadratur mit \(n=2\).

Beispiel 6.8 (Stützstellen einer Gauss-Quadratur)

Wir betrachten im folgenden Beispiel das Intervall \([0,1]\) und die Gewichtsfunktion \(\sigma(x) \equiv 1\) für \(x \in [0,1]\). Unser Ziel ist es die Nullstellen des Orthogonalpolynoms \(\omega_n \in \overline{\Pi}_n\) für \(n=2\) zu berechnen.

Wegen der Rekursionsformel (6.17) gilt für das \(0\)-te Orthogonalpolynom

\[\omega_0(x) \ = \ 1.\]

Da wir wissen, dass die Orthogonalpolynome normiert sind, können wir für \(\omega_1\) den Ansatz \(\omega_1(x) = x + a\) machen Aus der Orthogonalitätsbedingung

\[\int_0^1 \underbrace{\omega_0(x)}_{\equiv 1} \omega_1(x) \underbrace{\sigma(x)}_{\equiv 1} \,\mathrm{d}x \ = \ 0\]

erhalten wir mit dem Hauptsatz der Integral- und Differentialrechnung

\[0 \ = \ \int_0^1 x + a \,\mathrm{d}x \ = \ \left[ \frac{1}{2}x^2 + ax + c \right]^1_0 \ = \ \frac{1}{2} + a.\]

Daher muss \(a = -\frac{1}{2}\) gelten und wir erhalten die Stützstelle \(x_1 = \frac{1}{2}\) als Nullstelle des Orthogonalpolynoms \(\omega_1(x) = (x - \frac{1}{2})\). Die entsprechende Gauss-Quadratur entspricht in diesem Fall der Mittelpunktsregel aus Beispiel 6.1.

Um das nächste Orthogonalpolynom \(\omega_2 \in \overline{\Pi}_2\) zu bestimmen machen wir den Ansatz \(\omega_2(x) = x^2 + bx +c\). Die unbekannten Parameter bestimmen wir aus den beiden Orthogonalitätsbedingungen

\[\int_0^1 \omega_0(x) \omega_2(x) \sigma(x) \,\mathrm{d}x \ = \ 0, \qquad \int_0^1 \omega_1(x) \omega_2(x) \sigma(x) \,\mathrm{d}x \ = \ 0.\]

Für die erste Orthogonalitätsbedingung erhalten wir

\[0 \ = \ \int_0^1 x^2 + bx + c \,\mathrm{d}x \ = \ \left[ \frac{1}{3}x^3 + \frac{b}{2}x^2 + cx +d \right]^1_0 \ = \ \frac{1}3 + \frac{b}2 + c.\]

Für die zweite Orthogonalitätsbedingung erhalten wir

\[\begin{split} 0 \ &= \ \int_0^1 (x-\frac{1}{2})(x^2 + bx + c) \,\mathrm{d}x \ = \ \int_0^1 x^3 + (b-\frac{1}{2})x^2 + (c - \frac{b}{2})x - \frac{c}{2} \,\mathrm{d}x\\ &= \ \left[ \frac{1}{4}x^4 + (\frac{b}{3} - \frac{1}{6})x^3 + (\frac{c}{2} - \frac{b}{4})x^2 - \frac{c}{2}x + d\right]^1_0\\ &= \frac{1}{4} + \frac{b}{3} - \frac{1}{6} + \frac{c}{2} - \frac{b}{4} - \frac{c}{2} \ = \ \frac{1}{12} + \frac{b}{12} \ = \ \frac{1}{12}(1+b). \end{split}\]

Hieraus können wir direkt folgern, dass \(b=-1\) sein muss. Setzen wir dies in die erste Orthogonalitätsbedingung wieder ein, so erhalten wir \(c = \frac{1}{2} - \frac{1}{3} = \frac{1}{6}\). Damit ist das zweite Orthogonalpolynom \(\omega_2 \in \overline{\Pi}_2\) eindeutig bestimmt mit

\[\omega_2(x) \ = \ x^2 - x + \frac{1}{6}.\]

Mit Hilfe der \(p\)-\(q\)-Formel lassen sich nun die Nullstellen von \(\omega_2\) bestimmen als:

\[x_1 \ = \ \frac{1}{2} + \frac{1}{2 \sqrt{3}}, \qquad x_2 \ = \ \frac{1}{2} - \frac{1}{2 \sqrt{3}}.\]

Bemerkung 6.10 (Fehlerabschätzung für die Gauss-Quadratur)

Da die Gauss-Quadratur einen Exaktheitsgrad von \(2n-1\) hat lässt sich mit der Fehlerabschätzung für den Approximationsfehler der numerischen Interpolation in Satz 5.2 für eine Funktion \(f \in C^{2n}([a,b])\) zeigen, dass gilt

\[\left\vert G_n(f) - \int_a^b f(x) \sigma(x)\,\mathrm{d}x \right\vert \ \leq \ \frac{\Vert f^{(2n)} \Vert_\infty}{(2n)!} \frac{(b-a)^{2n+1}}{2^{n}}.\]

Natürlich kann die Gauss-Quadratur auch als summierte Quadraturformel angewandt werden (siehe Stückweise interpolatorische Quadratur) indem man das Intervall \([a,b]\) in mehrere Teilintervalle zerlegt und dort die Nullstellen der Orthogonalpolynome verwendet.