5.3. Interpolationsformel nach Newton

Die in Interpolationsformel nach Lagrange beschriebenen Verfahren, d.h. die Lagrangesche Interpolationsformel und die Vandermonde-Matrix, haben den Nachteil, dass man alle durchgeführten Berechnungen zur Bestimmung des Interpolationspolynoms durch \(n+1\) Stützpunkte nicht erneut verwenden kann, sobald ein weiterer Stützpunkt hinzugenommen wird. Das heißt bei Hinzunahme weiterer Daten ist der bisher betriebene numerische Rechenaufwand umsonst, da man zur Bestimmung der neuen Polynomkoeffizienten wieder von vorne beginnen muss.

Glücklicherweise gibt es eine alternative Herangehensweise, die uns einen signifikanten Teil der Neuberechnung von Koeffizienten erspart. Anstatt die herkömmliche Polynomdarstellung wie in (5.4) zu verwenden verfolgt man den folgenden Ansatz nach Newton.

Definition 5.6 (Newtonsche Interpolationsformel)

Seien \((x_i,f_i) \in \mathbb{R}^2, i=0,\ldots,n\) eine Menge von \(n+1\) gegebenen Datenpunkten. Das unbekannte Interpolationspolynom \(P \in \Pi_n\) vom Grad \(n\) lässt sich mit der Newtonschen Interpolationsformel wie folgt darstellen:

(5.10)\[\begin{split} P(x) \: &= \: a_0 \, + \, a_1(x - x_0) \, + \, a_2(x - x_1)(x - x_0) \, + \, \ldots \, + \, a_n(x - x_{n-1})\ldots (x - x_0) \\ \: &= \: \sum_{k=0}^n a_k N_k(x), \end{split}\]

wobei wir die Newtonschen Basispolynome \(N_k \in \Pi_k, k=0,\ldots,n\) definieren als

\[N_k(x) \, \coloneqq \, \prod_{i=0}^{k-1} (x-x_i) \quad \text{ und } \quad N_0(x) \equiv 1.\]

Bemerkung 5.4

Wir können folgende Beobachtungen zur Newtonschen Interpolationsformel machen:

1. Durch die Wahl der Basispolynome \(g_k(x) \coloneqq N_k(x), k=0,\ldots,n\) lässt sich die Bestimmung der unbekannten Koeffizienten \(a_0,\ldots, a_n \in \R\) eines Polynoms \(P \in \Pi_n\) mit \(P(x_i) = f_i, i=0,\ldots,n\) als ein lineares Gleichungssystem \(Lx = b\) schreiben mit einer linken, unteren Dreiecksmatrix \(L\in \mathbb{R}^{(n+1)\times (n+1)}\) von der Form:

   \[\begin{pmatrix} N_0(x_0) & 0 & & 0 \\ N_0(x_1) & N_1(x_1) & 0 \\ \vdots & \vdots & \ddots & 0 \\ N_0(x_n) & N_1(x_n) & \cdots & N_n(x_n) \end{pmatrix} \begin{pmatrix} a_0 \\ \vdots \\ a_n \end{pmatrix} \ = \ \begin{pmatrix} f_0 \\ \vdots \\ f_n \end{pmatrix}.\]

   Man sieht ein, dass man durch Hinzufügen von weiteren Daten nur eine zusätzliche Zeile und Spalte in \(L\) erhält, die an der Dreiecksgestalt nichts ändern.

2. Für bekannte Koeffizienten \(a_0,\ldots,a_n\) benötigt die Auswertung des Lagrange-Polynoms an einer Stelle \(\xi\) insgesamt \(2n\) Additionen und \(2n\) Multiplikationen. Das effizientere Horner-Schema der Newtonschen Interpolationsformel lässt sich folgendermaßen angeben:

   (5.11)\[P(x) \ = \ ( \ldots ( a_n(x - x_{n-1}) + a_{n-1} ) (x - x_{n-2}) + \ldots + a_1)(x - x_0) + a_0.\]

   In dieser Form benötigt die Berechnung des Polynoms an einer Stelle \(\xi\) insgesamt nur noch \(2n\) Additionen und \(n\) Multiplikationen. Der numerische Rechenaufwand ist also höher als in der herkömmlichen Polynomdarstellung (5.4). Dennoch bietet die Newtonsche Interpolationsformel den Vorteil, dass man weitere Datenpunkte hinzunehmen kann ohne alle Berechnungen erneut durchzuführen, wie wir im weiteren Verlauf noch genauer diskutieren werden.

Offensichtlich handelt es sich in (5.10) um ein Polynom von Grad \(n\) mit \(n+1\) unbekannten Koeffizienten \(a_0,\ldots,a_n\), die nach Satz Satz 5.1 eindeutig bestimmt sind, da wir die Stützstellen \(x_i, i=0,\ldots,n\) als paarweise verschieden angenommen haben.

Zur numerischen Berechnung dieser unbekannten Koeffizienten kann man folgenden sukzessiven Algorithmus verwenden.

Algorithmus 5.1 (Interpolationsschema nach Newton)

Es sei eine Menge von \(n+1\) Datenpunkten \((x_i,f_i) \in \mathbb{R}^2, i=0,\ldots,n\) gegeben, deren Stützstellen paarweise verschieden sind, d.h. es gilt \(x_i \neq x_j\) für \(i \neq j\). Bei der Betrachtung der Newtonschen Interpolationsformel (5.10) fällt auf, dass alle Summanden bis auf den ersten für die Stützstelle \(x = x_0\) wegfallen und wir die folgende Gleichung erhalten:

\[P(x_0) \ = \ a_0.\]

Da \(P \in \Pi_n\) jedoch gleichzeitig Interpolationspolynom bezüglich der gegebenen Daten sein muss, wissen wir, dass auch \(P(x) = f_0\) gelten muss. Also setzen wir im ersten Schritt

\[a_0 \ = \ f_0.\]

Wir können nun wieder die Newtonsche Interpolationsformel betrachten für den Fall \(x = x_1\) und erkennen, dass alle Summanden bis auf die ersten beiden wegfallen. Wir erhalten somit die Gleichung

\[P(x_1) \ = \ a_0 + (x_1 - x_0)a_1 \ = \ f_0 + (x_1 - x_0)a_1.\]

Da \(P(x_1) = f_1\) gelten muss können wir den unbekannten Koeffizienten \(a_1\) also im zweiten Schritt bestimmen als

\[a_1 \ = \ \frac{f_1 - f_0}{x_1 - x_0}.\]

Dieses Vorgehen lässt sich nun sukzessiv anwenden und wir können die Berechnung des Koeffizienten \(a_j, 1\leq j \leq n\) im j-ten Schritt allgemein angeben als:

\[\begin{split} a_j \ &= \ \frac{f_j - a_0}{\prod_{i=0}^{j-1}(x_j - x_i)} \, - \, \frac{a_1}{\prod_{i=1}^{j-1}(x_j - x_i)} \, - \, \ldots \, - \, \frac{a_{j-1}}{(x_j - x_{j-1})} \\ \ &= \ \frac{f_j}{\prod_{i=0}^{j-1}(x_j - x_i)} \, - \, \sum_{k=0}^{j-1}\frac{a_k}{\prod_{i=k}^{j-1}(x_j - x_i)}. \end{split}\]

Bemerkung 5.5 (Rechenaufwand des Newtonschen Interpolationsschemas)

Die Berechnung der unbekannten Koeffizienten \(a_i, i=0,\ldots,n\) benötigt bei der Nutzung des Algorithmus 5.1 einen numerischen Rechenaufwand von \(n\) Divisionen, \(n(n-1)\) Multiplikationen und \(n(n+1)\) Additionen.

Wir stellen bei der Berechnung der unbekannten Koeffizienten \(a_i, i=0,\ldots,n\) in der Newtonschen Interpolationsformel (5.10) fest, dass der \(j\)-te Koeffizient \(a_j\) nur von Koeffizienten und Datenpunkten mit niedrigerem Index abhängt. Dies ist eine sehr nützliche Eigenschaft, denn dies erlaubt es uns zu einem bereits berechneten Interpolationspolynom in \(n+1\) Datenpunkten weitere Daten hinzuzufügen ohne die bisher berechneten Koeffizienten neu bestimmen zu müssen. Man kann sogar zeigen, dass die Anordnung der Stützstellen \(x_i, i=0,\ldots,n\) bei der Berechnung der Koeffizienten keine Rolle spielen. Das bedeutet, dass wir auch neue Daten zwischen zwei bekannten Stützstellen hinzunehmen können, wie im folgenden Beispiel.

Beispiel 5.2 (Newton-Interpolation)

Es seien zunächst die folgenden Datenpunkte \((x_i,f_i) \in \mathbb{R}^2, i=0,1\) gegeben:

\(x_i\)	0	2
\(f_i\)	0	4

Wir wollen das eindeutige Interpolationspolynom \(P \in \Pi_1\) mit dem Newton Verfahren für Polynominterpolation bestimmen. Hierfür bestimmen wir die unbekannten Koeffizienten \(a_0\) und \(a_1\) wie in Algorithmus Algorithmus 5.1 hergeleitet als:

\[a_0 \, = \, f_0 \, = \, 0, \qquad a_1 \, = \, \frac{f_1 - f_0}{x_1 - x_0} \, = \, \frac{4 - 0}{2 - 0} \, = \, 2.\]

Durch Einsetzen in die Newtonsche Interpolationsformel (5.10) erhalten wir das eindeutige, lineare Interpolationspolynom \(P \in \Pi_1\) als:

\[P(x) \ = \ a_0 + (x-x_0)a_1 \ = \ 0 + (x - 0)2 \ = \ 2x.\]

Nehmen wir nun an, dass wir einen weiteren Datenpunkt \((x_2,f_2) = (1,1)\) zusätzlich erhalten, dessen Stützstelle \(x_2=1\) zwischen den vorigen Stützstellen \(x_0=0\) und \(x_1=2\) liegt. Wir wollen nun das eindeutige, quadratische Interpolationspolynom \(P \in \Pi_2\) finden, das alle Daten interpoliert. Hierfür müssen wir nur den zusätzlichen Koeffizienten \(a_2\) mit Algorithmus 5.1 berechnen:

\[a_2 \ = \ \frac{f_2 - a_0}{(x_2 - x_0)(x_2 - x_1)} - \frac{a_1}{(x_2 - x_1)} \ = \ \frac{1}{-1} - \frac{2}{-1} \ = \ 1.\]

Durch Erweitern des bereits bekannten linearen Interpolationspolynoms aus den ursprünglichen Daten erhalten wir:

\[P(x) \ = \ 2x + a_2(x - x_1)(x - x_0) \ = \ 2x + 1(x - 2)(x - 0) \ = \ 2x + x^2 - 2x \ = \ x^2.\]

Wie aus dem obigen Beispiel 5.2 klar werden sollte, handelt es sich bei den Koeffizienten \(a_i\) nicht um die Vorfaktoren der üblichen Polynombasis aus Monomen \(x^i\). Daher muss man die geläufige Darstellung eines Polynoms wie in (5.4) erst aus der Newtonschen Interpolationsformel (5.10) durch Einsetzen berechnen.

5.3.1. Dividierte Differenzen

Anstatt die unbekannten Koeffizienten \(a_i, i=0,\ldots,n\) des Interpolationspolynoms \(P \in \Pi_n\) mit dem oben beschriebenen Algorithmus Algorithmus 5.1 zu berechnen verwendet man eine effizientere Methode, die auf die Idee der dividierten Differenzen basiert. Wir betrachten hierzu noch einmal die Newtonsche Interpolationsformel (5.10) für das eindeutige Interpolationspolynom Polynom \(P_{0,\ldots,n} \in \Pi_n\) vom Grad \(n\) durch \(n+1\) Datenpunkte \((x_i,f_i) \in \mathbb{R}^2, i=0,\ldots,n\). Wir tauschen hierbei die Bezeichnung der unbekannten Koeffizienten \(a_i\) wie folgt:

(5.12)\[P_{0,\ldots,n}(x) \ = \ f_0 \, + \, f_{0,1}N_1(x) \, + \, \ldots \, + \, f_{0,\ldots,n}N_n(x).\]

Die Zahlen \(f_{0,\ldots,n}\) nennt man dividierte Differenzen und es lässt sich zeigen, dass folgende Rekursionsformel für sie gilt:

(5.13)\[f_{i,\ldots,i+k} \ = \ \frac{f_{i+1,\ldots,i+k} \, - \, f_{i,\ldots,i+k-1}}{x_{i+k} \, - \, x_i}.\]

Die rekursive Abhängigkeit der dividierten Differenzen, welche ein Spezialfall der sogenannten Rekursionsformel nach Neville-Aitken ist, lässt sich mit Hilfe des folgenden Differenzenschemas veranschaulichen, welches sich spaltenweise beginnend mit der Spalte \(k=0\) berechnen lässt.

Tab. 5.1 Ein Ausschnitt des dividierte Differenzen Schemas.

	\(k=0\)	\(k=1\)	\(k=2\)	…
\(x_0\)	\(\mathbf{f_0}\)
		\(\mathbf{f_{0,1}}\)
\(x_1\)	\(f_1\)		\(\mathbf{f_{0,1,2}}\)
		\(f_{1,2}\)
\(x_2\)	\(f_2\)
⋮

Mit steigendem Grad \(n\) des zu bestimmenden Polynoms \(P_{0,\ldots,n} \in \Pi_n\) wächst die rekursive Abhängigkeit der neu zu bestimmenden Koeffizienten. Hierbei lässt sich eine zu berechnende dividierte Differenzen einfach aus den beiden linken Nachbarn aus der vorigen Spalte bestimmen. Zum Beispiel gilt folgender Zusammenhang im obigen Schema in Tab. 5.1 für die dividierten Differenzen:

\[f_{0,1} \ = \ \frac{f_1 - f_0}{x_1 - x_0}, \quad f_{1,2} \ = \ \frac{f_2 - f_1}{x_2 - x_1}, \quad f_{0,1,2} \ = \ \frac{f_{1,2} - f_{0,1}}{x_2 - x_0}.\]

Die fett markierten dividierten Differenzen in der obersten Schrägzeile des Schemas ergeben die eindeutig bestimmten Koeffizienten des Interpolationspolynoms \(P \in \Pi_n\) in (5.12).

Bemerkung 5.6 (Rechenaufwand für dividierte Differenzen)

Der numerische Rechenaufwand zur Bestimmung der dividierten Differenzen für ein Interpolationspolynom \(P \in \Pi_n\) vom Grad \(n\) beträgt \(n(n+1)/2\) Divisionen und \(n(n+1)\) Additionen und ist damit deutlich effizienter als das Newtonsche Interpolationsschema (vgl. Bemerkung 5.5).

Folgendes Beispiel erklärt die Berechnung der dividierten Differenzen anschaulich.

Beispiel 5.3 (Dividierte Differenzen)

Es seien die folgenden Datenpunkte \((x_i,f_i) \in \mathbb{R}^2, i=0,1,2\) gegeben:

\(x_i\)	0	1	3
\(f_i\)	1	3	2

Wir wollen das eindeutige Interpolationspolynom \(P_{0,1,2} \in \Pi_2\) mittels dividierter Differenzen bestimmen. Hierfür verwenden wir das Schema Tab. 5.1.

	\(k=0\)	\(k=1\)	\(k=2\)	…
\(x_0=0\)	\(\mathbf{f_0 = 1}\)
		\(\mathbf{f_{0,1} = 2}\)
\(x_1=1\)	\(f_1 = 3\)		\(\mathbf{f_{0,1,2} = -\frac{5}{6}}\)
		\(f_{1,2} = -\frac{1}{2}\)
\(x_2=3\)	\(f_2 = 2\)
⋮

Durch Ablesen der Koeffizienten \(f_0, f_{0,1}, f_{0,1,2}\) des Newtonschen Interpolationspolynoms \(P_{0,1,2}\) entlang der obersten Schrägzeile erhalten wir das Polynom

\[P_{0,1,2}(x) \ = \ f_0 + f_{0,1}(x-x_0) + f_{0,1,2}(x-x_0)(x-x_1) \ = \ 1 + 2(x - 0) - \frac{5}{6}(x-0)(x-1).\]

Wollen wir das Polynom an einer Stelle \(\xi = 2\) auswerten so berechnen wir in dieser Darstellung mittels Horner-Schema in (5.11)

\[P_{0,1,2}(2) \ = \ (-\frac{5}{6}(2-1) + 2)(2 - 0) + 1 \ = \ \frac{10}{3}.\]

5.3.2. Dämpfung von Rundungsfehlern

Es lässt sich zeigen, dass man die Indizierung der Datenpunkte \((x_i, f_i), i=0,\ldots,n\) beliebig durch Permutation ändern kann in der Newtonschen Interpolationsformel (5.10) und man trotzdem das gleiche Interpolationspolynom erhält, wie wir bereits angemerkt haben. Das heißt, wir würden bei der Methode der dividierten Differenzen für eine beliebige Permutation \({i_0,\ldots,i_n} = \Pi({0,\ldots,n})\) das gleiche Interpolationspolynom \(P_{i_0,\ldots,i_n} \equiv P_{0,\ldots,n}\) erhalten. Leider spielen bei der Implementierung des Verfahrens Rundungsfehler eine Rolle, so dass dieses theoretische Resultat numerisch nicht immer zutrifft. Jedoch hilft uns diese Aussage, um bei der Auswertung eines gegebenen Polynoms weniger Rundungsfehler zu machen, wie man im Folgenden sehen wird.

Wir nehmen an, dass die gegebenen Stützstellen \(x_i, i=0,\ldots,n\) monoton steigend angeordnet sind, d.h. es gilt \(x_0 < x_1 < \ldots < x_n\). Das Ziel wird es sein das gegebene Interpolationspolynom \(P_{0,\ldots,n} \in \Pi_n\) vom Grad \(n\) an einer Stelle \(\xi\) so auszuwerten, dass möglichst wenig Rundungsfehler entstehen. Sei nun \(x_{i_0}\) die zu \(\xi\) nächstgelegene Stützstelle, d.h. es gilt

\[| \xi - x_{i_0} | \ = \ \min \lbrace |\xi - x_i| ~ | ~ i=0,\ldots,n \rbrace.\]

Weiter können wir die Stützstellen \(x_i, i=0,\ldots,n\) nach ihrer Entfernung zu \(\xi\) sortieren und dabei die folgende Indizierung wählen:

\[| \xi - x_{i_k} | \ = \ \min \lbrace |\xi - x_i| ~ | ~ i \in I_k \rbrace,\]

mit der Indexmenge \(I_k = \lbrace 0,\ldots,n \rbrace \backslash \lbrace i_0,\ldots, i_{k-1} \rbrace\). Durch diese Umsortierung wird eine Permutation \(\Pi({0,\ldots,n}) = {i_0,\ldots,i_n}\) induziert. Dadurch erhalten wir eine alternative Darstellung des Newtonschen Interpolationspolynoms mit

\[P_{i_0,\ldots,i_n}(x) \ = \ f_{i_0} + f_{i_0,i_1}(x - x_{i_0}) + \ldots + f_{i_0,\ldots,i_n}(x - x_{i_0})\ldots(x - x_{i_{n-1}}).\]

Schreiben wir das Polynom \(P_{i_0,\ldots,i_n}(x)\) bei der Auswertung an der Stelle \(\xi\) in das zugehörige Horner-Schema (5.11) um, so erkennt man den Vorteil dieser Permutation:

\[P_{i_0,\ldots,i_n}(\xi) \ = \ (\ldots ( f_{i_0,\ldots,i_n}(\xi - x_{i_{n-1}}) + f_{i_0,\ldots,i_{n-1}}) (\xi - x_{i_{n-2}}) + \ldots + f_{i_0,i_1})(\xi - x_{i_0}) + f_{i_0}.\]

Es ist davon auszugehen, dass eine Dämpfung der Rundungsfehler bei der Auswertung von \(P_{i_0,\ldots,i_n}(\xi)\) auftritt, da jeder Term mit einem immer kleiner werdenden Faktor multipliziert wird.

Da wir angenommen haben, dass die Stützstellen \(x_0 < x_1 < \ldots < x_n\) monoton angeordnet sind folgt sofort, dass die von uns gewählte Permutation eine zusammenhängende Indexmenge liefert, d.h. es gilt

\[\lbrace i_l,\ldots,i_{l+k}\rbrace \ = \ \lbrace 0 \leq i \leq n ~ | ~ \min_{I_{l,k}} i_j \leq i \leq \max_{I_{l,k}} i_j \rbrace,\]

mit der Indexmenge \(I_{l,k} = \lbrace i_l,\ldots,i_{l+k}\rbrace\). Dadurch lassen sich die dividierten Differenzen aus dem Schema Tab. 5.1 verwenden. Jedoch wird bei der Auswertung nicht die oberste Schrägzeile genutzt, sondern ein Zickzack-Pfad im Schema verwendet, der durch die Permutation induziert wird.

Beispiel 5.4 (Dämpfung von Rundungsfehlern)

Es seien die gleichen Datenpunkte wie in Beispiel Beispiel 5.3 gegeben. Wir wollen die bereits berechneten dividierten Differenzen in der Art umsortieren, so dass bei Auswertungen mittels Horner-Schema ein möglichst geringer Rundungsfehler entsteht. Sei die Auswertestelle nun \(\xi = 2\), so erhalten wir nach dem oben beschriebenen Vorgehen eine Permutation der Indizes mit:

\[i_0 \ = \ 1, \quad i_1 = 2, \quad i_2 \ = \ 0.\]

Der zugehörige Pfad im dividierte Differenzenschema ist fett markiert:

	\(k=0\)	\(k=1\)	\(k=2\)	…
\(x_0=0\)	\(f_0 = 1\)
		\(f_{0,1} = 2\)
\(x_1=1\)	\(\mathbf{f_1 = 3}\)		\(\mathbf{f_{0,1,2} = -\frac{5}{6}}\)
		\(\mathbf{f_{1,2} = -\frac{1}{2}}\)
\(x_2=3\)	\(f_2 = 2\)
⋮

Durch Ablesen der fett markierten Koeffizienten \(f_1, f_{1,2}, f_{0,1,2}\) des Newtonschen Interpolationspolynoms \(P_{0,1,2}\) erhalten wir das Polynom

\[P_{0,1,2}(x) \ = \ f_1 + f_{1,2}(x-x_1) + f_{0,1,2}(x-x_1)(x-x_2) \ = \ 3 - \frac{1}{2}(x - 1) - \frac{5}{6}(x-1)(x-3).\]

Wollen wir das Polynom nun an der Stelle \(\xi = 2\) auswerten, so berechnen wir in dieser Darstellung mittels Horner-Schema (5.11)

\[P_{0,1,2}(2) \ = \ (-\frac{5}{6}(2-3) - \frac{1}{2})(2 - 1) + 3 \ = \ \frac{10}{3}.\]