"Man muß es wohl einfach mal ausprobieren", dachte ich mir in der VL Mathe für Informatiker III. Die beiden vorgestellten Methoden zur Bestimmung eines Polynoms, das eine Kurve, die durch endlich viele Stützwerte gegeben ist, annähert, sind wirklich einfach in ein Programm "zu gießen".
Im folgenden habe ich für beide Methoden ein C++ Programm geschrieben, in dem zu erkennen ist, wie die angegebenen Formel 1:1 umgesetzt worden sind.
Außerdem habe ich ein Java-Applet "InterpolationApplet" geschrieben, in dem Referenzkurven und Stützwerte ausgewählt werden können und das resultierende Interpolationspolynom graphisch und als Funktionsterm ausgegeben werden. Mathematik interaktiv, sozusagen. (Endlich mal ein nützliches Applet)
Hierbei beziehe ich mich auf die Rekursionsformel, die in der VL als erste unter der Teilüberschrift "Praktische Berechnung" angegeben wurde. Danach ist p0,n das gesuchte, eindeutig bestimmte Interpolationspolynom. Die Rekursionsvorschrift ist in Satz 1.2.2 angegeben. Folgendes Programm berechnet pi,l(x) für einzelne x.
Zunächsteinmal benötigen wir sinnvollerweise eine Datenstruktur für einen Punkt im R2:
typedef long double fzahl;
typedef struct s_point {
fzahl x;
fzahl y;
} point;
Obwohl eine Rekursionsformel angegeben wurde, wird diese Rekursion
beseitigt und die sonst durch Rekursion ermittelten Zwischenwerte werden
in einem 2-dimensionalen Array pa gehalten. Außerdem beötigen
wir ein Array der Stützwerte, pstellen. Die Dimension hängt
von n ab, es gibt n+1 Stützstellen. Als C++-Text:
/** Polynomgrad */
const int n = 3;
fzahl pa[n+1][n+1]; // pa[i][l] entspricht p_{i,l}
point pstellen[n+1]; // die stützstellen
Bevor die Berchnung beginnen kann, müssen natürlich die Stützwerte
bekannt sein, d.h. pstellen muß mit Werten gefüllt werden.
Das Array pa läuft direkt auf ein Dreiecksschema hinaus:
als Initalisierung wird pi,0=yi für i=0,...,n
gesetzt (z.B. in main()):
for(int i = 0;i<=n;i++) pa[i][0] = pstellen[i].y;
und dann verlangt die Berechnung von pi,l die Werte von pi,l-1 und pi+1,l-1 (siehe Formel!). Mit dieser Information können wir die Funktion schreiben:
fzahl y1(fzahl x)
{
for(int l=1;l<=n;l++)
for(int i=0;i<=n-l;i++){
pa[i][l] = ((x-pstellen[i].x)*pa[i+1][l-1]
+ (pstellen[i+l].x - x)*pa[i][l-1])/
(pstellen[i+l].x - pstellen[i].x);
}
return pa[0][n];
}
Das war auch schon alles! Natürlich gibts den kompletten Quelltext.
Die Newton-Darstellung ist lediglich eine andere Formel, um auf das selbe Interpolationspolynom zu kommen, wie bei der Aitken-Form.
Hier zunächst die Formel: p(x) = y[x0] + (x-x0)y[x0,x1] +(x-x0)(x-x1)y[x0,x1,x2] + ... + (x-x0)...(x-xn-1)y[x0,...xn] .
Zunächst stellen wir fest, daß sich hier die Faktoren (x-xi) n-imal wiederholen. Für ein festes x muß dies also nur einmal berechnet werden. y[] bezeichnet man als dividierten Differenzen. Sie ergeben sich durch Rekursion, die sich wieder (ähnlich wie bei der Aitkenform) in einer Dreiecksform darstellen lassen. Hier die Rekursionsformeln:
y[xi] = yi.
y[x0,...,xi+l] = (1/(xi+l-xi)) *(y[xi+1,...,xi+l]-y[xi+1,...,xi+l])
Im Programm brauchen wir dafür wieder ein zweidimensionales Feld der Dimension n+1, daß entsprechend der ersten Zeile initalisiert werden muß:
fzahl divdiff[n+1][n+1];
int main()
{
// ...
for(int i = 0;i<=n;i++)
divdiff[i][0] = pstellen[i].y;
// ...
}
Dabei entspricht divdiff[i][j] dem Term
y[xi,...,xi+j].
In der Funktion wird noch ein Feld der Länge n angelegt, in dem die Werte der Faktoren der Form (x-xi) gespeichert werden. Dann erfolgt die Berechnung. Genau wie bei Aitken gibt es zwei verschachtelte Schleifen, in der Inneren werde die dividierten Differenzen berechnet, in der Äußeren dann der entsprechende Summand des Interpolationpolynoms:
fzahl yNewton(fzahl x)
{
/** xdiff[i] = produkt [j=0 bis i-1] (x-x_j) */
fzahl xdiff[n+1];
// ergebnis
fzahl erg = 0.0;
xdiff[0] = 1.0;
erg = divdiff[0][0]; // = y_0
for(int l=1;l<=n;l++){
for(int i=0;i<=n-l;i++){
divdiff[i][l] = (divdiff[i+1][l-1]-divdiff[i][l-1])/
(pstellen[i+l].x-pstellen[i].x);
}
xdiff[l] = xdiff[l-1] * (x-pstellen[l-1].x);
erg += xdiff[l]*divdiff[0][l];
}
return erg;
}
Wie hierbei ersichtlich ist, sind die dividierten Differenzen gar nicht vom aktuellen Wert von x abhängig, sondern nur von den Stützwerten. Wenn also ein Interpolationspolynom an mehreren Stellen bei den gleichen Stützwerten berechnet wird, reicht die einmalige Berechnung der dividierenden Differenzen. Dies wurde in der Java-Version implementiert.
Die Newton-Darstellung (und auch die Aitken-Form) ist gegenüber der Reihenfolge der Stützwerte invariant, soll heißen, es muß nicht unbedingt x0<...<xn gelten; die Stützwerte können vielmehr in beliebiger Reihenfolge vorliegen.
Dies ermöglicht es, bei einer Hinzunahme weiterer Stützwerte einfach nur den notwendigen Teil der dividierten Differenzen neu zu berechnen, da nicht alle Werte von den neuen abhängen.
Die Ausgabe des Interpolationspolynoms als Zeichenkette habe ich nur in Java implementiert, da hier gegenüber C++ eine wesentlich einfachtes Stringhandling vorliegt. Es unterscheidet sich jedoch (logischerweise) nicht grundsätzlich von der Berechnung des Interpolationspolynoms.
(P)+(C) by
Erik Pischel
Last modified: Mon Nov 23 21:13:22 MET