Numerische Mathematik, Teil 2: Interpolation
Nun sind wir schon bei der zweiten Folge. Heute befassen wir uns mit dem Interpolationsproblem. Bei diesem Problem geht man von einer Anzahl gegebener Wertepaare aus und bestimmt eine durch die gegebenen Punkte verlaufende, möglichst glatte Kurve.
Formulieren wir das Problem mathematischer: Gegeben sind n+1 Wertepaare (xi,yi) (i=0,1,..., n). Die xi sind paarweise verschieden. Gesucht ist nun eine Funktion p(x) derart, daß p(xi)=yi ist. Anschaulich ist dieses Problem in Bild 1 dargestellt. Da bekanntlich viele Wege nach Rom führen, gibt es auch hier mehrere Möglichkeiten, ein passende Funktion p(x) zu finden. Wir beschränken uns hier jedoch ausschließlich auf ganzrationale Funktionen p(x). p(x) ist also ein Polynom. Es gibt zwar auch andere Möglichkeiten. eine Interpolationsfunktion zu finden. Beispielsweise könnte p(x) auch gebrochen rational sein oder aus trigonometrischen Funktionen bestehen. Auch stückweise Interpolation z.B. mit Splines wollen wir hier nicht berücksichtigen. Ein Beispielprogramm zur Spline-Interpolation findet sich in [6] (s.u.).
p(x) heißt in unserem Fall das Interpolationspolynom, die xi heißen Stützstellen und die yi die zugehörigen Stützwerte. Bei n+1 gegebenen Stützstellen (mit den passenden Stützwerten) hat das Polynom p(x) höchstens den Grad n, d. h. es hat die Form p(x)=a0+a1x+a2x2+...+anxn, mit zu bestimmenden ai, i=0, 1, ..., n. Unter den genannten Voraussetzungen gibt es genau ein Interpolationspolynom. Unsere Aufgabe ist es nun, dieses zu finden. Beginnen wir mit der Interpolationsformel von Joseph Louis Lagrange (geb. 25.1.1738 in Turin, gest. 10.4.1813 in Paris). Die Lagrange’sche Interpolationsformel (Bild 2) besteht aus einer Summe mit den Stützwerten und Gewichten, in die die Stützstellen eingehen. Die Interpolationsformel löst die Interpolationsaufgabe, jedoch ist die Verwendung dieser Formei nicht gerade praktisch, da durch Hinzufügen weiterer Stützstellen die gesamte Formel neu berechnet werden müßte. Dafür ist sie jedoch von theoretischer Bedeutung, da man aus ihr z. B. Formeln zur Integration gewinnen kann. Näheres hierzu jedoch erst in Folge 3 dieser Serie.
Zur Auswertung des Interpolationspolynoms an anderen als den Stützstellen kann man die Rechnung etwas abkürzen. Die Formel, die zur Berechnung des Lagrange’schen Interpolationspolynoms herangezogen wird, heißt baryzentrische Formel der Lagrange-Interpolation. Die li (l=Lambda, hier fehlt mir leider der griechische Buchstabe) werden rekursiv berechnet. Das Listing 1 ist dazu nur geringfügig zu verändern. Die Stützkoeffizienten li hängen nicht von der auszuwertenden Stelle ab, sie sind also nur einmal zu berechnen. Aber von dieser Stelle ab muß für jede neu auszuwertende Stelle x die Rechnung komplett neu durchgeführt werden. Somit hätten wir zwei Formen des Lagrange'schen Interpolationspolynoms besprochen.
Wenden wir uns einer anderen Interpolationsformel zu, der Interpolationsformel von Isaac Newton (geb. 4.1.1643 in Wollsthorpe (Lincolnshire), gest. 31.3.1727 in London). Wie in Bild 3 zu sehen ist, hat diese Formel einen ganz anderen Aufbau als die Formel von Lagrange. In den Koeffizienten der einzelnen Summanden tauchen nicht mehr jeweils alle xi auf, sondern sukzessive immereiner mehr als im Vorgänger. Dies ist in einem Programm besonders einfach auszuwerten. Das einzige Problem besteht in der Bestimmung der ck. Die ck werden über das Schema der sogenannte dividierten Differenzen berechnet. Das Schema erkennt man aus der im Bild 3 dargestellten Formel für die ck sowie dem Schema in Bild 4. Jedes ck baut auf zwei Formeln auf, die auf einer niedrigeren Stufe stehen. Auf der untersten Stufe werden die f[xi]=yi gesetzt. Besondere Vorteile bietet die Formel von Newton, weil man neue Stützstellen hinzufügen kann. Es ist pro hinzugefügter Stützstelle nur ein ck neu zu berechnen, welches man zudem noch teilweise aus alten Werten erhält. Aus den beschriebenen Interpolationsformeln von Lagrange und Newton erhält man leicht das Interpolationspolynom.
Mit dem Algorithmus von Neville (Bild 5) wird die Aufgabe der Bestimmung eines zu interpolierenden Wertes direkt angegangen. Das Neville-Schema (Bild 6) ähnelt etwas dem Schema der dividierten Differenzen. Interpolation mit Polynomen hätten wir nun besprochen. Da mit der Anzahl der Stützstellen auch der Grad des Interpolationspolynoms steigt, kann man sich leicht überlegen, daß diese Art der Interpolation nicht für viele Stützstellen geeignet ist. Und genau so ist es auch. Man greift dann eher zu anderen Interpolationsmethode wie z. B. der Spline-Interpolation.
Spline-Interpolation ist zwar auch Interpolation mit Polynomen, jedoch haben diese Polynome beispielsweise bei der kubischen Spline-Interpolation höchstens den Grad 3, unabhängig von der Anzahl der Stützstellen. Dadurch wird ein vernünftiger Verlauf der Interpolationsfunktion erreicht. Nehmen wir als Beispiel die kubische Spline-Funktion. Sie wird stückweise definiert. Zwischen zwei Stützstellen wird jeweils ein Polynom dritten Grades bestimmt. Der Verlauf der zusammengesetzten Graphen dieser Polynome ergibt dann den Verlauf des Graphen der Interpolationsfunktion. Hier merkt man schon den Unterschied zu den oben beschriebenen Interpolationsmethoden, bei nur ein Polynom gesucht wurde. Hier sind es bei n+1 Stützstellen bereits n Polynome. Eine kubische Spline-Funktion ist charakterisiert durch die gerade beschriebene Eigenschaft, auf jedem Teilintervall zwischen je zwei Stützstellen ein Polynom dritten Grades zu sein. Ferner ist die zweite Ableitung stetig (d. h. der zugehörige Graph läßt sich mit einem Stift zeichnen, ohne den Stift vom Papier abzuheben) an den inneren Stützstellen. Selbstverständlich erfüllt die Spline-Funktion auch die Forderungen des Interpolationsproblems. Diese Art der Interpolation wurde beispielsweise im Schiffsbau eingesetzt, um die Biegung einer Straklatte zu bestimmen. Straklatte sagt Ihnen nichts? Das englische Wort dafür ist Spline! Aber es gibt auch spielerische Anwendungen. Für einige Animationen und Computerspiele wird die Bahn eines fliegenden Objekts benötigt. Da diese Bahn schön rund aussehen soll, wird diese über Splines-interpoliert. Bei der rationalen Interpolation wird eine gebrochenrationale Funktion R(x)=P(x)/Q(x) mit ganzrationalen Funktionen P und Q angesetzt. Diese Interpolationsfunktionen können nur in speziellen Fällen benutzt werden, da sie gegenüber den ganzrationalen Interpolationsfunktionen einen gewaltigen Nachteil haben: Sie erfüllen nicht in jedem Fall die Forderungen des Interpolationsproblems. Einige Punkte sind daher unerreichbar.
Im nächsten Monat werden wir teilweise auf bauend auf diese Folge die Themen Integralrechnung und Differentiation näher betrachten. Bis dahin...
Dipl.-Math. Dietmar Rabich
Ergebnisse der Programme
| Wertetabelle, bestimmt mit Lagrange-Interpolation: | Wertetabelle, bestimmt mit Neville-Interpolation: | Wertetabelle, bestimmt mit Newton-Interpolation: | |||
|---|---|---|---|---|---|
| x | p(x) | x | p(x) | x | p(x) |
| 0.00000000 | 0.10000000 | 0.00000000 | 0.09999999 | 0.00000000 | 0.10000000 |
| 0.20000000 | 0.35348859 | 0.20000000 | 0.35348854 | 0.20000000 | 0.35348862 |
| 0.40000001 | 0.35791311 | 0.40000000 | 0.35791308 | 0.40000000 | 0.35791314 |
| 0.60000002 | 0.23721629 | 0.60000000 | 0.23721624 | 0.60000000 | 0.23721629 |
| 0.80000001 | 0.21728800 | 0.80000000 | 0.21728792 | 0.80000000 | 0.21728798 |
| 1.00000000 | 0.50000000 | 1.00000000 | 0.50000000 | 1.00000000 | 0.50000000 |
| 1.20000005 | 1.13724053 | 1.20000000 | 1.13724041 | 1.20000000 | 1.13724083 |
| 1,40000010 | 1.90494657 | 1.40000000 | 1.90494990 | 1.40000000 | 1.90495031 |
| 1.60000014 | 2.17714524 | 1.60000000 | 2.17715383 | 1.60000000 | 2.17715442 |
| 1.80000019 | 0.79999673 | 1.80000000 | 0.79999995 | 1.80000000 | 0.80000000 |
Literatur:
[1] Einführung in die Numerische Mathematik I.
J. Stoer.
Springer Berlin/ Heidelberg/ New York/ Tokyo.
4. Aufl. 1983. S. 3!ff
[2] Formelsammlung zur Numerischen Mathematik mit BASIC-Programmen.
G. Engeln-Miillges/ F. Reutter,
Bibliographisches Institut Mannheim/ Wien/ Zürich. 1. Aufl. 1983. S. I30ff
[3] Numerische Mathematik.
H. R. Schwarz,
Teubner Stuttgart. 1. Aufl. 1986. S. 88ff
[4] Numerische Methoden.
A. Björck/ G. Dahh/uist.
Oldenbourg München Wien. 2. Aufl. 1979. S. 201 ff
[5] Methode der Numerischen Mathematik.
W. Böhm/ G. Gose/ J. Kahmann.
Vieweg Braunschweig/ Wiesbaden,
- Aufl. 1985. S. 89ff
[6] Rund muß es sein.
D. Rabich,
ST-Computer 1/89, S. 84ff
[7] Erfolgreich programmieren mit C. ,
J. A. lllik.
Sybex Düsseldorf/ San Francisco/ Paris/ London, 4. Aufl. 1987
[8] Programmieren in C,
B. W. Kernighan/ D. M. Ritchie.
Hanser München/ Wien. 1. Aufl. 1983
[9] PASCAL für Anfänger.
H. Schauer,
Oldenbourg Wien/ München, 4. Aufl. 1982
[10] PASCAL für Fortgeschrittene.
H. Schauer,
Oldenbourg Wien/ München, 2. Aufl. 1983
[11] Programmieren in Modula-2.
N. Wirth.
Springer Berlin/ Heidelberg/ New York/ Tokyo.
- Aufl. 1985
(***********************************************)
(* Beispielprogramm zur Lagrange-Interpolation.*)
(* ------------------------------------------- *)
(* Entwickelt mit SPC Modula 2. 15.02.1989 *)
{***********************************************)
(* --------- + ---------------- *)
(* Listing 1 / (c) by D. Rabich *)
(* --------- + ---------------- *)
MODULE LagrangeInterpolation;
(* Importe *)
FROM InOut IMPORT WriteReal,WriteLn,WriteString,Read;
(* Konstanten *)
CONST MaximalStuetz = 10;
FloatFehler = 1.0E-6;
(* Variablen *)
VAR Lambda,
x,y : ARRAY [0..MaximalStuetz] OF REAL;
xw : REAL;
MaxStuetz : CARDINAL;
VoidChar : CHAR;
(* Stützwerte einlesen *)
PROCEDURE ReadStuetzStellen;
BEGIN
MaxStuetz:=5;
x[0]:=-0.1;
x[1]:= 0.0;
x[2]:= 0.5;
x[3]:= 0.7;
x[4]:= 1.0;
x[5]:= 1.8;
y[0]:=-0.1;
y[1]:= 0.1;
y[2]:= 0.3;
y[3]:= 0.2;
y[4]:= 0.5;
y[5]:= 0.8
END ReadStuetzStellen;
(* Lambda-Koeffizienten berechnen *)
PROCEDURE CalcLambda;
VAR i,k : CARDINAL;
Summe : REAL;
BEGIN
Lambda[0]:=1.0;
FOR k:=1 TO MaxStuetz DO
FOR i:=0 TO k-1 DO
Lambda[i]:=Lambda[i]/(x[i]-x[k])
END;
Summe:=0.0;
FOR i:=0 TO k-1 DO
Summe:=Summe+Lambda[i]
END;
Lambda[k]:=-Summe
END
END CalcLambda;
(* Mu-Koeffizienten berechnen *)
PROCEDURE CalcMuSums (BerX : REAL; VAR Sum1,Sum2 : REAL);
VAR i : CARDINAL;
Mu : REAL;
BEGIN
Sum1:=0.0;
Sum2:=0.0;
FOR i:=0 TO MaxStuetz DO
Mu :=Lambda[i]/(BerX-x[i]);
Sum1:=Sum1+Mu*y[i];
Sum2:=Sum2+Mu
END
END CalcMuSums;
(* Interpolationswert berechnen *)
PROCEDURE Interpol (BerX : REAL) : REAL;
VAR i : CARDINAL;
Sum1,Sum2 : REAL;
BEGIN
ReadStuetzStellen;
(* Ist BerX Stützstelle? *)
FOR i:=0 TO MaxStuetz DO
IF ABS(x[i]-BerX)<FloatFehler THEN
RETURN y[i]
END
END;
(* Wert berechnen *)
CalcLambda;
CalcMuSums(BerX,Sum1,Sum2);
RETURN Sum1/Sum2
END Interpol;
(* Hauptprogramm *)
BEGIN
xw:=0.0;
WriteString(' x p(x)');
WriteLn;
REPEAT
WriteReal(xw,12,8);
WriteString(' '');
WriteReal(Interpol(xw),12,8);
WriteLn;
xw:=xw+0.2
UNTIL xw>2.0;
Read(VoidChar)
END LagrangeInterpolation.
(**********************************************)
(* Beispielprogramm zur Newton-Interpolation. *)
(* -------------------------------------------*)
(* Entwickelt mit ST Pascal Plus. 13.02.1989 *)
(**********************************************)
(* --------- + ---------------- *)
(* Listing 2 / (c) by D. Rabich *)
(* --------- + ---------------- *)
program newton_interpolation;
(* Konstante *)
const max_dim = 10;
(* Variablen *)
var max_stuetz : short_integer;
x,y,
c : array [0..max_dim] of real;
xw : real;
(* Werte einlesen *)
procedure read_werte;
begin
max_stuetz:=5;
x[0]:=-0.1;
x[1]:= 0.0;
x[2]:= 0.5;
x[3]:= 0.7;
x[4]:= 1.0;
x[5]:= 1.8;
y[0]:=-0.1;
y[1]:= 0.1;
y[2]:= 0.3;
y[3]:= 0.2;
y[4]:= 0.5;
y[5]:= 0.8
end;
(* Koeffizienten berechnen *)
procedure calc_coeff;
var i,k : short_integer;
begin
for i:=0 to max_stuetz do
c[i]:=y[i];
for k:=1 to max_stuetz do
for i:=max_stuetz downto k do
c[i]:=(c[i]-c[i-1])/(x[i]-x[i-k])
end;
(* Wert interpolieren *)
function interpol (ber_x : real) : real;
var i : short_integer;
p : real;
begin
read_werte;
calc_coeff;
(* Polynom auswerten *)
p:=c[max_stuetz];
for i:=max_stuetz-1 downto 0 do
p:=c[i]+(ber_x-x[i])*p;
interpol:=p
end;
(* Hauptprogramm *)
begin
xw:=0.0;
writeln(' x p(x)');
(* Werte ausgeben *)
repeat
writeln(xw:12:8,' ',interpol(xw):12:8);
xw:=xw+0.2
until xw>2.0;
(* auf Tasten warten... *)
repeat
until keypress
end.
/***********************************************/
/* Beispielprogramm zur Neville-Interpolation. */
/* ------------------------------------------- */
/* Entwickelt mit Turbo C. 13.02.1989 */
/***********************************************/
/* --------- + ---------------- */
/* Listing 3 / (c) by D. Rabich */
/* --------- + ---------------- */
/* Ein-/Ausgabe-Routinen importieren */
# include <stdio.h>
/* Konstante */
# define MAXDIM 10
/* globale Variablen */
int max_stuetz; *
float x[MAXDIM],y[MAXDIM],p[MAXDIM];
/* Werte einlesen */
void read_werte (void)
{
max_stuetz=5;
x[0]=-0.1;
x[1]= 0.0;
x[2]= 0.5;
x[3]= 0.7;
x[4]= 1.0;
x[5]= 1.8;
y[0]=-0.1;
y[1]= 0.1;
y[2]= 0.3;
y[3]= 0.2;
y[4]= 0.5;
y[5]= 0.8;
}
/* Wert interpolieren */
float interpol (float ber_x)
{
int i,k;
read_werte();
for (i=0;i<=max_stuetz;i++)
p[i]=y[i];
for (k=1;k<=max_stuetz;k++)
for(i=max_stuetz;i>=k;i—-)
p[i]+=(ber_x-x[i])*(p[i]-p[i-1])/(x[i]-x(i-k]) ;
return(p[max_stuetz]);
}
/* Hauptprogramm */
void main (void)
{
float xw;
printf(" x p(x)\n");
for(xw=0.0;xw<=2.0;xw+=0.2)
printf("%12.8g %12.8g\n",xw,interpol(xw));
getchar();
}
Dietmar Rabich