In einem anderen Artikel habe ich die Interpolation nach dem Verfahren von Neville-Aitken betrachtet. Es wird verwendet, wenn man das Interpolations-Polynom nur an wenigen Stellen auswerten will. Benötigt man allerdings die Auswertung an vielen Stellen, empfiehlt es sich einen anderen Ansatz zu wählen.
Ein Polynom, das in der Newtonschen Darstellung gegeben ist, kann schnell mittels modifiziertem Hornerschema ausgewertet werden. Die Newtonschen Darstellung sieht folgendermaßen aus:

p(x) = c0 + c1(x - x0) + c2(x - x0)(x - x1) + ... + cn(x - x0)(x - x1)...(x - xn-1) =

n ∑ i=0

(ci ⋅

i-1 ∏ j=0

(x - xj))

Zur Berechnung der ci werden sogenannte Dividiert Differenzen benutzt:

f[xk] := fk

f[xk,xk+j] :=

f[xk+1,xk+j] - f[xk,xk+j-1] xk+j - xk

Die ci ergeben sich dann als c0 = f[x0], ci = f[x0, xi]. Diese Koeffizienten kann man sehr gut in einem Tabellenschema darstellen und damit auch berechnen:

x0	f0 = c0
x1	f1	f[x0, x1] = c1
x2	f2	f[x1, x2]	f[x0, x2] = c2
x3	f3	f[x2, x3]	f[x1, x3]	f[x0, x3] = c3

Die Methode der Dividierten Differenzen zur Berechnung der ci wird auch als Steigungsschema bezeichnet. In Matlab sieht der Algorithmus dann wie folgt aus:

• function c = koeffizienten(f,s)
• % Steigungsschema zur Berechnung der Koeffizienten c_i
• c = f(s);
• n = length(s);
• for(i=2:n)
• for(k=n:-1:i)
• c(k) = (c(k)-c(k-1)) / (s(k)-s(k-i+1));
• end
• end
• end
•  
• function y = polynom(x,s,c)
• % Auswertung des Polynoms mit den Koeffizienten c_i durch modifiziertes Hornerschema
• y = c(n);
• for(k=n-1:-1:1)
• y = c(k)+(x-s(k)).*y;
• end
• end