Lagrange Interpolation

Die Lagrange Interpolation kann dazu verwendet werden, ein Polynom möglichst kleinen Grades durch eine Menge von Funktionswerten $y_i$ zu gegeben Stützstellen $x_i$ zu legen. Liegen $n$ solche Paare $(x_i,y_i)$ vor, so reicht hierzu ein Polynom vom Grad $d=n-1$ aus. Ziel der Interpolationsaufgabe ist es hierbei, die Koeffizienten $a_0,a_1,\ldots a_d$ für das gesuchte Polynom \[ p(x):=a_0+a_1\cdot x+a_2\cdot x^2+a_3\cdot x^3+\cdots+a_n\cdot x^n \] zu finden. Es muss allerdings gefordert werden, dass die $x_i$ paarweise verschieden sind, andernfalls könnte eine inkonsistente Forderung an einen $y$-Wert entstehen. Eine Möglichkeit diese Koeffizienten zu berechnen besteht darin, explizit das folgende Gleichungssystem zu lösen:

\[ \left( \begin{array}{ccccc} 1&x_1&x_1^2&\ldots& x_1^d \\ 1&x_2&x_2^2&\ldots& x_2^d \\ \vdots&\vdots&\vdots&& \vdots\\ 1&x_n&x_n^2&\ldots& x_n^d \\ \end{array} \right) \left( \begin{array}{c} a_0\\ a_1\\ \vdots\\ a_n \end{array} \right) = \left( \begin{array}{c} y_0\\ y_1\\ \vdots\\y_n \end{array} \right) \]

Die linke Seite wertet das fragliche Polynom an den Stützstellen $x_i$ aus. Die rechte Seite setzt diese Auswertung mit den geforderten $y$-Werten gleich. Erstaunlicherweise (und erfreulicherweise) ist dieses lineare Gleichungssystem tatsächlich immer dann lösbar, wenn die $x_i$ paarweise verschieden sind. Der Grund dafür ist die so genannte Vandermonde Matrix» und deren Determinante (Genaueres hierzu kommt später in der Vorlesung). Zugegebenermaßen ist es nicht die klügste Variante, die Polynomkoeffizienten auf diese Weise zu bestimmen. Das obige Gleichungssystem ist sehr schlecht konditioniert (d.h. im Allgemeinen numerisch instabil) bedingt durch den Umstand, dass in der Matrix zugleich sehr große und sehr kleine Zahlen auftreten können. Eine bessere Methode, das gesuchte Polynom zu bestimmen, besteht aus den folgenden Schritten: Zunächst bestimmt man Polynome \[ q_i(x):=(x-x_1)\cdot(x-x_2)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_{n-1})(x-x_{n}). \] Diese haben folgende drei wichtige Eigenschaften: Alle $q_i(x)$ haben Grad $d$, Für $i\neq j$ gilt $q_j(x_i)\neq 0$, Es gilt $q_i(x_i)= 0$ für alle $i$. Aus diesen bestimmt man normierte Polynome gemäß $p_i(x):=q_i(x)/q_i(x_i)$. Diese Polynome $p_i$ sind an allen Stützstellen außer $x_i$ gleich $0$. Bei $x_i$ ergibt sich $p_i(x_i)=1$. Das gesuchte Polynom $p(x)$ erhält man dann durch eine Linearkombination der $p_i(x)$: \[ p(x):=y_1\cdot p_1(x)+\cdots+y_n\cdot p_n(x). \] Das folgende Applet verdeutlicht das Ergebnis dieser Berechnung für vier frei positionierbare Punkte. Die grünen Punkte sind frei beweglich.

Mit dem Schiebeschalter können die vier gewichteten Basisfunktionen $y_ip_i(x)$ sichtbar gemacht werden. Der algorithmische Kern des Applets besteht aus den foglenden 8 Zeilen in CindyScript

   pts=allpoints();
   xs=apply(pts,#.x);
   ys=apply(pts,#.y);
   n=length(pts);
   q(x,i):=product(1..n--[i],j,x-xs_j);
   p(x,i):=f(x,i)/f(xs_i,i)*ys_i;
   p(x):=sum(1..n,i,p(x,i));
   plot(p(x));