Polynominterpolation

Grundidee und Bedeutung

Polynominterpolation ist ein Verfahren der numerischen Mathematik. Gesucht wird ein Polynom, das exakt durch vorgegebene Punkte verläuft, zum Beispiel durch Messwerte. Dieses Polynom heißt Interpolationspolynom; man sagt, es interpoliert die gegebenen Punkte.

Für n+1 Wertepaare (x_i, f_i) mit paarweise verschiedenen Stützstellen x_i wird ein Polynom P höchstens n-ten Grades gesucht, das P(x_i)=f_i für i=0,...,n erfüllt. Dieses Polynom existiert immer und ist eindeutig bestimmt. Gesucht wird P im Vektorraum Pi_n der Polynome vom Grad n oder kleiner.

Die Wahl der Basis für Pi_n bestimmt, welches lineare Gleichungssystem entsteht. In der Standardbasis P(x)=sum_{k=0}^n a_k x^k erhält man ein Gleichungssystem mit Vandermonde-Matrix. Diese Matrix ist regulär, wenn die Stützstellen paarweise verschieden sind. Theoretisch zeigt das die Eindeutigkeit und Existenz der Lösung. Praktisch wird dieses Gleichungssystem aber meist nicht verwendet, weil seine Lösung aufwendig und im Allgemeinen schlecht konditioniert ist.

Interpolierende Polynome sind wichtig, weil Polynome leicht integriert und abgeleitet werden können. Darum treten sie zum Beispiel bei numerischer Integration und bei Verfahren zur numerischen Lösung gewöhnlicher Differentialgleichungen auf.

Lagrange- und baryzentrische Form

Eine wichtige Darstellung ist die Lagrange-Basis. Die Lagrange-Polynome l_i(x) sind so konstruiert, dass l_i(x_k)=delta_ik gilt. Das Kronecker-Delta delta_ik ist 1 für i=k und 0 für i ungleich k. Dadurch ergibt sich das Interpolationspolynom direkt als P(x)=sum_{i=0}^n f_i l_i(x).

Der Vorteil dieser Form ist: Die Basisfunktionen l_i hängen nur von den Stützstellen x_i ab, nicht von den Stützwerten f_i. Wenn also dieselben Stützstellen mit verschiedenen Werten verwendet werden, kann man die einmal berechneten Basisfunktionen wiederverwenden. Der Nachteil ist: Kommt eine neue Stützstelle hinzu, müssen alle Basisfunktionen neu berechnet werden. Deshalb ist diese Darstellung für viele praktische Zwecke zu aufwendig. In der digitalen Signalverarbeitung wird Lagrange-Interpolation unter dem Namen Farrow Filter für adaptives Resampling eingesetzt.

Eine praktisch wichtigere Umformung ist die baryzentrische Interpolationsformel: P(x)= (sum_{j=1}^n lambda_j f_j/(x-x_j)) / (sum_{j=1}^n lambda_j/(x-x_j)) für x ungleich x_1,...,x_n. Dabei sind die baryzentrischen Gewichte lambda_j:=produkt_{i=1,i ungleich j}^n 1/(x_j-x_i). Für feste Stützstellen können diese Gewichte vorberechnet werden. Dann kostet die Auswertung von P(x) nur O(n). Beim Hinzufügen einer neuen Stützstelle müssen die Gewichte neu bestimmt werden; das kostet O(n), während das Neubestimmen der Lagrangepolynome O(n^2) kostet.

Newton-Verfahren und Neville-Aitken

Beim Newtonschen Algorithmus wird das Polynom in der Newton-Basis dargestellt. Die Basisfunktionen sind N_0(x)=1 und N_i(x)=produkt_{j=0}^{i-1}(x-x_j) für i=1,...,n. Damit lautet die Newtonsche Interpolationsformel: P(x)=sum_{i=0}^n c_i N_i(x). Im Gegensatz zur Standardbasis entsteht dabei eine untere Dreiecksmatrix, die sich einfacher lösen lässt.

Die Koeffizienten c_i werden effizient mit dividierten Differenzen bestimmt. Es gilt c_i=[x_0,...,x_i]f. Die dividierten Differenzen sind rekursiv definiert durch [x_i]f=f_i und [x_i,...,x_j]f=([x_{i+1},...,x_j]f-[x_i,...,x_{j-1}]f)/(x_j-x_i). Fügt man einen weiteren Punkt hinzu, muss nur eine neue Zeile im Schema ergänzt werden; die bisherigen Koeffizienten müssen nicht neu berechnet werden.

Ist das Polynom in Newton-Form bekannt, kann es mit dem Horner-Schema effizient ausgewertet werden. Dafür wird P(x) geschachtelt geschrieben, und man berechnet rekursiv b_n=c_n, b_i=b_{i+1}(x-x_i)+c_i für i=n-1,...,0 und schließlich P(x)=b_0. Der Aufwand beträgt O(n).

Der Neville-Aitken-Algorithmus berechnet Interpolationswerte ebenfalls rekursiv. Er verwendet Teilpolynome p(f|x_i,...,x_j), also die Interpolationspolynome zu Ausschnitten der Stützstellen. Die Rekursionsformel von Aitken kombiniert zwei kleinere Interpolationspolynome zu einem größeren. Der Newtonsche Algorithmus ist besonders geeignet, wenn alle Koeffizienten bestimmt und das Polynom oft ausgewertet werden sollen; sein Aufwand für die Koeffizienten beträgt O(n^2), die Auswertung O(n). Neville-Aitken eignet sich eher, wenn nur wenige Auswertungen gebraucht werden, und ist weniger anfällig gegen Auslöschung.

Beispiel Tangensfunktion

Ein Beispiel interpoliert die Funktion f(x)=tan(x) an fünf Punkten: x_0=-1,5 mit f(x_0)=-14,101420, x_1=-0,75 mit f(x_1)=-0,931596, x_2=0 mit f(x_2)=0, x_3=0,75 mit f(x_3)=0,931596 und x_4=1,5 mit f(x_4)=14,101420.

Mit der Lagrange-Methode werden zuerst die fünf Lagrange-Basisfunktionen l_0 bis l_4 bestimmt. Daraus ergibt sich nach Einsetzen der Funktionswerte das Interpolationspolynom P_Lagrange(x)=-1,477474x+4,834848x^3.

Mit der Newton-Methode werden zuerst die dividierten Differenzen berechnet. Das Polynom wird dann in Newton-Form geschrieben: P_Newton(x)=-14,1014+17,5597(x+1,5)-10,8784(x+1,5)(x+0,75)+4,83484(x+1,5)(x+0,75)x+0(x+1,5)(x+0,75)x(x-0,75). Ausmultipliziert ergibt sich näherungsweise P_Newton(x)=-0,00005-1,4775x-0,00001x^2+4,83484x^3. Verwendet man genauere Startwerte f(x_i), verschwinden der erste und der dritte Koeffizient.

Fehler und Wahl der Stützstellen

Die Interpolationsgüte beschreibt, wie gut das Interpolationspolynom eine Funktion annähert. Sei f eine Funktion, deren n+1 Funktionswerte f_i an Stellen x_i durch P interpoliert werden. Ist I das kleinste Intervall, das die Stützstellen und eine Stelle x enthält, und ist f auf I (n+1)-mal stetig differenzierbar, dann gibt es ein xi in I mit f(x)-P(x)=f^{(n+1)}(xi)/(n+1)! * produkt_{i=0}^n (x-x_i). Daraus folgt eine Abschätzung mit der Maximumsnorm: |f(x)-P(x)| ist höchstens ||f^{(n+1)}||infty/(n+1)! mal ||w_n||infty, wobei w_n(x)=produkt{i=0}^n (x-x_i). Eine weitere Fehlerformel lautet f(x)-P(x)=f[x_0,...,x_n,x] produkt{i=0}^n (x-x_i).

Der Fehler hängt also sowohl von einer Ableitung von f als auch von den Stützstellen ab. Wenn man die Stützstellen frei wählen kann, ist es günstig, ||w_n||infty klein zu machen. Auf [-1,1] gilt für normierte Stützstellen ||w_n||{[-1,1],infty}>=2^{-n}. Tschebyschow zeigte, dass die Nullstellen der Tschebyschow-Polynome optimale Stützstellen sind. Für T_{n+1}(x)=cos((n+1) arccos(x)) lauten die Nullstellen t_k=cos(((2k+1)/(2n+2)) pi) für k=0,1,...,n. Damit gilt scharf ||w_n||{[-1,1],infty}=2^{-n}. Für ein allgemeines Intervall [a,b] erhält man durch Transformation die Abschätzung ||w_n||{[a,b],infty}=2((b-a)/4)^{n+1}.

Grenzen, Konvergenz und Verallgemeinerung

Mehr Stützpunkte verbessern die Interpolation nicht immer. Bei äquidistanten Stützstellen und hohem Polynomgrad kann das Rungesche Phänomen auftreten: Das Polynom ähnelt der zu interpolierenden Funktion kaum noch und zeigt starke Oszillationen nahe den Intervallgrenzen. Das ist besonders problematisch, wenn die Funktion selbst sich anders verhält als Polynome für x gegen plus oder minus unendlich, etwa periodisch oder asymptotisch konstant. Ein klassisches Beispiel ist die Runge-Funktion f(x)=1/(1+x^2) auf [-5;5]. Tschebyschow-Stützstellen können den Gesamtfehler verkleinern; oft ist aber ein anderes Verfahren wie Spline-Interpolation geeigneter.

Es gibt dennoch Konvergenzaussagen. Ist f analytisch auf I=[a,b] und wird das Stützstellengitter immer feiner, dann konvergieren die zugehörigen Interpolationspolynome gleichmäßig gegen f. Formal: Für Intervallteilungen Delta_m mit Norm ||Delta_m||=max_i |x_{i+1}^{(m)}-x_i^{(m)}| gilt aus ||Delta_m|| gegen 0, dass P_{Delta_m} gegen f gleichmäßig konvergiert. Andererseits besagt der Satz von Faber, dass es zu jeder Folge von Intervallteilungen auch eine stetige Funktion f gibt, für die die Interpolationspolynome nicht gleichmäßig gegen f konvergieren.

Die Nähe zur besten Polynomapproximation wird über die Lebesgue-Konstante Lambda_n beschrieben. Für eine stetige Funktion f, ein Interpolationspolynom P und eine Bestapproximation Q vom Grad höchstens n gilt ||f-P||max <= (1+Lambda_n)||f-Q||max. Die Lebesgue-Konstante ist die Operatornorm des Interpolationsoperators phi_n, oft bezüglich der Maximumsnorm, und kann berechnet werden als Lambda_n=max{x in [a,b]} sum{i=0}^n |l_i(x)|.

Eine Verallgemeinerung ist die Hermiteinterpolation. Dort müssen die Stützstellen nicht paarweise verschieden sein. An mehrfach vorkommenden Stützstellen werden nicht nur Funktionswerte, sondern auch Werte der Ableitungen des Interpolationspolynoms vorgegeben.