Autor: Stefan Huber

[ < ] [ globale Übersicht ] [ Kapitelübersicht ] [ Stichwortsuche ] [ > ]


6.7.4 Wertberechnung

Eine der wichtigsten Manipulationen mit Interpolationspolynomen ist deren Auswertung (die Bestimmung von Funktionswerten) an einer oder mehreren Stellen.

Horner-Algorithmus

Für ein Polynom der Potenzform

P_d(x) = c_0 + c_1 x + \cdots + c_d x^d (Formel 9.28)

gibt es eine bekannte Rekursion- das Horner-Schema - zur Berechnung des Wertes von $P_d$ an der stelle x:

$p:=c_d$
do
$i=1,2,\dots,d$
$p:=x p+c_{d-i}$
(Formel 9.29)
end do

Das Horner-Schema benötigt jeweils d Addidionen und Multiplikationen zur Berechnung des gesuchten Funktionswertes $P_d(x):=p$. Dieser Rechenaufwand für die Auswertung eines Polynoms mit gegebnen Koeffizienten bezüglich der Monombasis $B_d= \{ 1, x, \ldots, x^d\}$ ist optimal. Es gibt kein Verfahren, das mit weniger Rechenoperationen auskommt.
Als Vermutung stammt diese Optimalitätsaussage aus einer Arbeit von Ostrowski [318] aus dem Jahre 1954, die man oft als die Geburtsstunde der albebraischen Komplexitätstheorie ansieht. Der Beweis der Vermutung wurde 1966 von Pan [320] geliefert.

Clenshaw-Algorithmus

Erhält man als Resultat eines Interpolations- oder Approximationsverfahrens ein Polynom in der Tschebyscheff-Form $P_d=a_0/2 + a_1 T_1 + \cdots + a_d T_d$ und will man dieses an der Stelle x auswerten, dann könnte man es zuerst in die Potenzform (9.28) bringen und dann das Horner-Schema (9.29) anwenden.
Die dafür erforderliche Umrechnung der Koeffizienten ist aber mit unnötigem Aufwand verbunden und verschlechtert außerdem die Genauigkeit des Resultats. Günstiger ist es, mit den Koeffizienten $a_0,a_1,\ldots,a_d$ der Tschebyscheff-Entwicklung den Clenshaw-Algorthmus anzuwenden, der als Resultat den gesuchten Funktionswert $P_d(x) := b_0$ liefert:

$b_{d}:=a_d$
$b_{d-1}:=2xb_d+a_{d-1}$
do $i=2,3,\dots,d\!-\!1$
$b_{d-i}:=2xb_{d-i+1}-b_{d-i+2}+a_{d-i}$
end do
$b_0:=xb_1-b_2+a_0$

Eine genaue Herleitung des Clenshaw-Algorthmus, der auf der Rekursion (9.12) beruht, und einen Beweis seiner numerischen Stabilität findet man z.B. bei Fox und Parker [199].

Software (Auswertung eines Polynoms in Tschebyscheff-Darstellung
Ein Polynom, das durch Koeffizienten seiner Tschebyscheff-Darstellung gegeben ist, kiann mittels der Unterprogramme NAG/e02aef bzw. NAG/e02akf an einer gegebenen Stelle ausgewertet werden, NAG/e02akf läßt sich dabei flexibler verwenden als NAG/e02aef, besitzt aber auch eine entsprechend komplizierte Parameterliste.

Algorithmen für Ableitungen und Integrale

Auch wenn man an Ableitungswerten oder Integralen

{P_d}'(x), \quad {P_d}''(x), \quad \ldots \quad
         \int \limits_I P_d(t)\, dt

interessiert ist, gibt es Algorithmen und Programme, die auf speziell für Ableitungen oder Integrale modifizierten Koeffizienten der Tschebyscheff-Darstellung des Interpolationspolynoms beruhen.

Software (Differenzieren und Integrieren von Polynomen in Tschebyscheff-Darstellung)
Die Tschebyscheff-Darstellung der ersten Ableitung eines bereits in Tschebyscheff-Darstellung gegebenen Polynoms berechnet das Unterprogramm NAG/e02ahf. Die Koeffizienten der Tschebyscheff-Darstellung des unbestimmten Integrals für ein solches Polynom können mit Hilfe von NAG/e02ajf ermittelt werden.

Algorithmus von de Casteljau

Ein Polynom, das durch die Koeffizienten $\beta_0,\beta_1,\dots,\beta_d$ seiner Bézier-Darstellung (9.11) gegeben ist, läßt sich mit der Rekursion

do $i=0,1,\dots,d$
$b_i:=\beta_i$
end do
do $k=1,2,\dots,d\!-\!k$
$b_i:=(1-x)b_i+xb_{i+1}$
end do

auswerten, die als Resultat $P_d(x):=b_0$ liefert ( Locher [62]).

Werteberechnung ohne Koeffizientenberechnung

Wenn man nur einen Wert (oder eine geringe Anzahl von Werten) des Interpolationspolynoms benötigt, dann wäre es unzweckmäßig, zuerst die Koeffizienten einer der obigen Darstellungen des Interpolationspolynoms zu berechnen und dieses dann an der gewünschten Stelle auszuwerten.
Für diese Problemstellung gibt es spezielle Algorithmen, die bezüglich des Rechenaufwandes günstiger sind. Diese Verfahren beruhen auf dem Prinzip der sukzessiven linearen Interpolation, deren Grundlage folgender Satz ist:

Satz 9.3.2 (Lemma von Aitken)
An der Knotenmenge $\{x_0,x_1,\dots,x_d\}$ seien die Werte $f_0,f_1,\dots,f_d$ vorgegeben. Hat man für $i \ne j$

  1. das Interplationspolynom $P_{d-1}^1$ zu den Knoten $\{x_0,\dots,x_{i-1},x_{i+1},\dots,x_d\}$,
  2. das Interpolationspolynom $P_{d-1}^2$ zu den Knoten $\{x_0,\dots,x_{j-1},x_{j+1},\dots,x_d\}$,

so erhält man das Interpolationspolynom $P_d$ zu der gesamten Knotenmenge durch folgende Linearkombinationen:

P_d(x) := \frac{x_j - x}{x_i-x_j} \cdot P_{d-1}^1(x) -
                   \frac{x_i - x}{x_i-x_j} \cdot P_{d-1}^2(x).

Beweis: Henrici [54].

Ausgehend vom Lemma von Aitken kann man eine Vielfalt von iteraktiven Interpolationsverfahren konstruieren, die sich hauptsächlich in der Reihenfolge der Verwendung der Wertepaare $(x_i,f_i)$ unterscheiden. Die verbreitetsten Verfahren sind jene von Neville und von Aitken. Im folgenden soll die iterierte Interpolation nach E.H. Neville - das Neville-Schema - beschrieben werden. Dabei wird das Lemma von Aitken sukzessive auf die Polynome $P_{k-1}^{j-1}$ zu den Knoten $\{x_{j-k}, x_{j-k+1},\dots,x_{j-1}\}$ sowie $P_{k-1}^j$ zu den Knoten $\{x_{j-k+1}, x_{j-k+2},\dots,x_j\}$ angewendet; das resultierende Polynom $P_k^{j}$ interpoliert dann an den Knoten $\{x_{j-k},x_{j-k+1},\dots,x_j\}$:

do $j=0,1,\dots,d$
$P_0^j:=f_j$
end do
do $k=1,2,\dots,d$
do $j=k,k\!+\!1,\dots,d$
$\displaystyle{P_k^{j}(x)\,:=\,\frac{x_j-x}
        {x_j-x_{j-k}}\cdot P_{k-1}^{j-1}(x)\, - \,
        \frac{x_{j-k}-x}{x_j-x_{j-k}}\cdot P_{k-1}^j(x)}$
(Formel 9.30)
end do
end do

Um diese Rekursion zur Berechnung des Wertes von $P_d$ an einer konkreten Stelle $\xi$ zu verwenden, muß man lediglich im Neville-Schema (9.30) x durch $\xi$ ersetzen; $P_d^{d}(\xi)$ ist dann das gesuchte Resultat.

Software (Direkte Interpolation einzelner Funktionswerte)
Wird für eine Menge gegebener Datenpunkte keine expliziete Darstellung des interpolierenden Polynoms benötigt, sondern nur ein einzelner Wert dieser Funktion, so kann dieser mit dem Unterprogramm NAG/e01aaf nach dem Aitken-Verfahren durch sukzessive lineare Interpolation berechnet werden. Die Datenpunkte können dabei beliebig verteilt sein. Dabei wird nicht nur das Endresultat $P_d^d(x)$, sondern auch alle Resultate $P_k^j(x)$ von Teilinterpolationen zurückgeliefert. Duch Vergleich dieser Werte kann man sich Informationen über die vermutliche Größe des absoluten Interpolationsvehlers $P_d^d(x)-f(x)$ verschaffen.
Soferne äquidistante (gleichabständige) Datenpunkte gegeben sind, kann man das Unterprogramm NAG/e01abf verwenden. Dieses berechnet einen einzelnen Wert des Interpolationspolynoms mit Hilfe der Everett-Formel ( Hildebrand [55]).


Autor: Stefan Huber

[ < ] [ globale Übersicht ] [ Kapitelübersicht ] [ Stichwortsuche ] [ > ]