Interpolation und Quadratur: Lagrange, Newton, Simpson und Gauß

Messungen und Rechnungen liefern Funktionswerte an einzelnen Stellen: eine Tabelle von Temperaturen, Sinuswerte aus einem Handbuch, die Ergebnisse einer teuren Simulation. Was ist der Wert dazwischen? Und wie berechnet man ein Integral von einer Funktion, die man nur punktweise kennt oder für die es keine Stammfunktion gibt? Beide Fragen führen auf dieselbe Idee: Ersetze die Funktion durch ein einfaches Polynom, das durch die bekannten Punkte geht, und rechne mit dem Polynom.

Die Interpolation bestimmt das Polynom, die Quadratur integriert es. Auf dieser Seite lernst du beides mit Fehlerabschätzungen, und du siehst, warum man in der Praxis nicht einfach ein Polynom hohen Grades durch alle Punkte legt. Die Grundlagen zu Rundungsfehlern stehen auf der Seite Numerik-Grundlagen, die Integralrechnung auf Integralrechnung und Das Riemann-Integral.

Das Interpolationspolynom

Beweis. Eindeutigkeit: Sind $p$ und $q$ zwei solche Polynome, hat $p-q$ den Grad höchstens $n$ und mindestens $n+1$ Nullstellen $x_0,\dots ,x_n$. Ein Polynom vom Grad $\le n$ mit $n+1$ Nullstellen ist das Nullpolynom. Existenz: Setze $L_i(x)=\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}$. Das Polynom $L_i$ hat den Grad $n$, und es gilt $L_i(x_k)=1$ für $k=i$ und $L_i(x_k)=0$ sonst. Dann leistet $p(x)=\sum _{i=0}^n{y}_i{L}_i(x)$ das Verlangte.

Der Ansatz $p(x)=a_0+a_{1}x+\dots +a_n{x}^n$ führt auf ein lineares System mit der Vandermonde-Matrix. Es ist lösbar, aber für große $n$ schlecht konditioniert (siehe Kondition einer Matrix), deshalb benutzt man die folgenden Formen.

Lagrange-Form

Die Formel $p(x)=\sum _i{y}_i{L}_i(x)$ aus dem Beweis heißt Lagrange-Form. Sie ist ideal für Theorie und Beweise, weil die Basispolynome nur von den Knoten abhängen, nicht von den Werten.

Newton-Form und dividierte Differenzen

Bei der Lagrange-Form muss man alles neu berechnen, wenn ein Punkt dazukommt. Die Newton-Form vermeidet das. Man schreibt

Die Koeffizienten sind die dividierten Differenzen $c_k=f[x_0,\dots ,x_k]$, rekursiv erklärt durch $f[x_i]=y_i$ und

Kommt ein Punkt hinzu, wird nur ein weiterer Summand angehängt. Die Auswertung geschieht mit dem Horner-Schema.

Der Interpolationsfehler

Beweis. Für einen Knoten $x=x_i$ ist beides gleich 0. Sei $x$ kein Knoten und $\omega (t)=\prod (t-x_i)$. Wähle die Konstante $K$ so, dass $g(t)=f(t)-p(t)-K\omega (t)$ bei $t=x$ verschwindet. Dann hat $g$ die $n+2$ Nullstellen $x_0,\dots ,x_n,x$. Nach dem Satz von Rolle hat $g^{\prime }$ mindestens $n+1$ Nullstellen, $g^{\prime }$ mindestens $n$, und schließlich $g^{(n+1)}$ mindestens eine, $\xi$. Wegen $p^{(n+1)}=0$ und ${\omega }^{(n+1)}=(n+1)!$ folgt $0=f^{(n+1)}(\xi )-K(n+1)!$, also $K=\frac{f^{(n+1)}(\xi )}{(n+1)!}$.

Der Fehler hängt also von zwei Dingen ab: von der Glätte der Funktion (der Ableitung $f^{(n+1)}$) und von der Lage der Knoten (dem Knotenpolynom $\omega$).

Das Runge-Phänomen

Man könnte hoffen, dass mehr Knoten immer bessere Näherungen liefern. Das stimmt nicht. Carl Runge zeigte 1901, dass die Interpolation mit gleichmäßig verteilten Knoten und hohem Grad stark zu schwingen beginnt, am Rand des Intervalls.

Bei gleichmäßigen Knoten wird der Fehler mit wachsendem Grad größer, bei Tschebyschow-Knoten $x_i=\cos \frac{(2i+1)\pi }{2(n+1)}$, die sich zum Rand hin häufen, sinkt er. Der Grund liegt im Knotenpolynom $\omega$: Bei gleichmäßigen Knoten ist es am Rand viel größer als in der Mitte, die Tschebyschow-Knoten minimieren sein Maximum. Deshalb legt man in der Praxis kein Polynom hohen Grades durch alle Punkte, sondern Splines: stückweise Polynome niedrigen Grades (meist kubisch), die an den Knoten glatt aneinanderpassen.

Numerische Integration

Ist $f$ nicht elementar integrierbar (wie $e^{-x^2}$), kennt man nur Werte, oder soll schnell gerechnet werden, ersetzt man $f$ durch das Interpolationspolynom $p$ und integriert dieses. Die entstehenden Formeln haben die Form ${\int }_a^{b}f(x)dx\approx \sum _{i=0}^n{w}_{i}f(x_i)$ mit Knoten $x_i$ und Gewichten $w_i$. Wählt man die Knoten gleichmäßig verteilt, heißen sie Newton-Cotes-Formeln.

Trapezregel (zwei Knoten, Gerade). Die Gerade durch $(a,f(a))$ und $(b,f(b))$ integriert zu $T=\frac{b-a}{2}(f(a)+f(b))$. Nach der Fehlerformel für $n=1$ ist $f(x)-p(x)=\frac{f^{\prime }(\xi )}{2}(x-a)(x-b)$, und mit ${\int }_a^b(x-a)(x-b)dx=-\frac{{(b-a)}^3}{6}$ folgt für den Integralfehler ${\int }_a^{b}f-T=-\frac{{(b-a)}^3}{12}{f}^{\prime }(\eta )$. Die zusammengesetzte Regel mit $n$ Teilintervallen der Länge $h=\frac{b-a}{n}$ hat den Fehler $-\frac{(b-a)h^2}{12}{f}^{\prime }(\eta )$, also $O(h^2)$.

Simpson-Regel (drei Knoten, Parabel). Mit dem Mittelpunkt $m=\frac{a+b}{2}$ ergibt die Parabel durch $a$, $m$, $b$ die Formel

Die Gewichte $\frac{1}{6},\frac{4}{6},\frac{1}{6}$ erhält man, indem man die drei Lagrange-Basispolynome über $[a,b]$ integriert. Überraschend ist die Genauigkeit: Die Regel ist nicht nur für Parabeln, sondern auch für alle kubischen Polynome exakt, weil sich der Fehler des kubischen Anteils wegen der Symmetrie gegenseitig aufhebt. Der Fehler der zusammengesetzten Regel ist $-\frac{(b-a)h^4}{180}{f}^{(4)}(\eta )$, also $O(h^4)$. Dieselbe Formel kannte schon Kepler 1615 für die Volumen von Weinfässern, daher heißt sie im Deutschen auch Keplersche Fassregel.

Aus den beiden Regeln lässt sich eine bessere gewinnen. Die Kombination $\frac{4T_{2n}-T_n}{3}$ ergibt genau die Simpson-Regel $S_{2n}$ (im Beispiel $\frac{4·\mathrm{1,75393}-\mathrm{1,85914}}{3}=\mathrm{1,71886}$), und die Romberg-Extrapolation setzt dieses Prinzip fort, um aus Werten der Trapezregel immer höhere Ordnungen zu gewinnen.

Gauß-Quadratur

Bei den Newton-Cotes-Formeln sind die Knoten fest. Gibt man sie frei, gewinnt man doppelt so viele Freiheitsgrade: Mit $n$ Knoten und $n$ Gewichten kann man $2n$ Bedingungen erfüllen.

Für $n=1$ ist es die Mittelpunktsregel, für $n=2$ gilt ${\int }_{-1}^{1}f(x)dx\approx f(-\frac{1}{\sqrt{3}})+f\left(\frac{1}{\sqrt{3}}\right)$. Diese Regel ist exakt für $1,x,x^2,x^3$: Es ist ${\left(\frac{1}{\sqrt{3}}\right)}^2+{\left(\frac{1}{\sqrt{3}}\right)}^2=\frac{2}{3}={\int }_{-1}^1{x}^{2}dx$, und die ungeraden Potenzen heben sich auf. Bei $x^4$ ist sie nicht mehr exakt ($\frac{2}{9}$ statt $\frac{2}{5}$). Für ein beliebiges Intervall transformiert man mit $x=\frac{a+b}{2}+\frac{b-a}{2}t$.

Typische Fehler

Übungen

Geschichte: von Weinfässern zu Legendre-Polynomen

Die Interpolationsformel trägt den Namen von Joseph-Louis Lagrange, der sie 1795 veröffentlichte. Entdeckt wurde sie schon 1779 von Edward Waring, und sie folgt leicht aus einer Formel, die Leonhard Euler 1783 publizierte. Die Newton-Form mit dividierten Differenzen ist nach Isaac Newton benannt. Die Simpson-Regel ist nach Thomas Simpson (1710 bis 1761) benannt, im Deutschen nach Johannes Kepler, der sie 1615 für Weinfässer herleitete. Carl Friedrich Gauß veröffentlichte 1814 die Quadratur mit optimalen Knoten, und Carl Gustav Jacob Jacobi brachte sie 1826 in die heutige Form mit orthogonalen Polynomen. Das Runge-Phänomen entdeckte Carl Runge 1901.

Häufige Fragen

Quellen und weiterführende Literatur

• Wikipedia: Polynominterpolation
• Wikipedia: Simpsonregel
• Wikipedia: Gauß-Quadratur
• Wikipedia (englisch): Lagrange polynomial
• Peter Deuflhard, Andreas Hohmann: Numerische Mathematik 1, De Gruyter, Berlin.
• Josef Stoer, Roland Bulirsch: Numerische Mathematik 1, Springer, Berlin.