Robert Sedgewick: Algorithmen

[ Inhaltsverzeichnis ] [ vorhergehende Seite ] [ nächste Seite ] [ Stichwort ]


38. Kurvenanpassung



Spline-Interpolation

Dennoch sind Polynome von niedrigem Grad einfache Kurven, mit denen man analytisch leicht arbeiten kann, und sie werden in breitem Umfang für die Kurvenanpassung angewandt. Der Trick besteht darin, daß man sich von der Idee löst, ein Polynom durch alle Punkte legen zu wollen, und stattdessen verschiedene Polynome verwendet, um benachbarte Punkte zu verbinden, wobei man sie auf glatte Weise zusammensetzt. Ein eleganter Spezialfall, welcher auch nur eine relativ einfache Rechnung erfordert, wird Spline-Interpolation genannt.

Das Wort »spline« bezeichnet eine mechanische Vorrichtung, die von Zeichnern benutzt wird, um ästhetisch ansprechende Kurven zu zeichnen: Der Zeichner fixiert eine Menge von Punkten (Knoten) auf seiner Zeichnung, biegt dann ein biegsames Band aus Plastik oder Holz (spline genannt) um sie herum und zeichnet es nach, um die Kurve zu erzeugen. Die Spline-Interpolation ist das mathematische Gegenstück zu diesem Prozeß und führt zu der gleichen Kurve. Abbildung 38.1 zeigt einen durch zehn Knoten verlaufenden Spline.

Abbildung 38.1
Abbildung 38.1 Ein Spline, der durch zehn Knoten verläuft.

Mit Hilfe von Gesetzen der elementaren Mechanik kann gezeigt werden, daß die von dem Band zwischen zwei benachbarten Knoten angenommene Form einem Polynom dritten Grades (das heißt, einem kubischen Polynom) entspricht. In Übertragung auf unser Problem der Datenanpassung bedeutet das, daß wir davon ausgehen sollten, daß sich die Kurve aus N - 1 verschiedenen kubischen Polynomen

si(x) = aix3 + bix2 + cix + di, i = 1, 2, ..., N - 1

zusammensetzt, wobei si(x) als das kubische Polynom definiert ist, das im Intervall zwischen xi und xi+1 zu verwenden ist.

Der Spline kann in offensichtlicher Weise in Form von vier eindimensionalen Feldern (oder von einem zweidimensionalen Feld vom Typ 4 mal (N - 1)) dargestellt werden. Die Erzeugung eines Splines besteht in der Berechnung der erforderlichen Koeffizienten a, b, c, d aus den gegebenen Punkten x und Werten y. Die physikalischen Gesetze, denen das Band gehorchen muß, entsprechen gleichzeitig zu erfüllenden Gleichungen, die gelöst werden können, um die Koeffizienten zu bestimmen.

Zum Beispiel muß offenbar si(xi) = yi und si(xi+1) = yi+1 für i = 1, 2, ..., N - 1 gelten, da das Band die Knoten berühren muß. Doch das Band berührt die Knoten nicht nur, sondern es biegt sich auch glatt um sie herum, ohne scharfe Biegungen oder Knicke. Mathematisch bedeutet das, daß die ersten Ableitungen der Spline-Polynome in den Knoten übereinstimmen müssen (s'i-1(xi) = s'i(xi) für i = 2, 3, ..., N - 1). Tatsächlich erweist es sich, daß die zweiten Ableitungen der Polynome in den Knoten ebenfalls übereinstimmen müssen. Diese Bedingungen ergeben eine Gesamtzahl von 4N - 6 Gleichungen in den 4(N - 1) unbekannten Koeffizienten. Zwei weitere Bedingungen müssen angegeben werden, um die Situation in den Endpunkten des Splines zu beschreiben. Dafür stehen verschiedene Möglichkeiten zur Verfügung; wir wollen den sogenannten »natürlichen« Spline benutzen, der sich aus s''1(x1) = 0 und s''N-1(xN) = 0 ableiten läßt. Diese Bedingungen ergeben ein vollständiges System von 4N - 4 Gleichungen mit 4N - 4 Unbekannten, das unter Benutzung des Gaußschen Eliminationsverfahrens gelöst werden könnte, um alle Koeffizienten zu berechnen, die den Spline beschreiben.

Der gleiche Spline kann jedoch in einer etwas effizienteren Weise berechnet werden, da in Wirklichkeit nur N - 2 »Unbekannte« vorhanden sind: Die meisten der Bedingungen an den Spline sind redundant. Nehmen wir zum Beispiel an, daß pi der Wert der zweiten Ableitung des Splines in xi ist, so daß s''i-1(xi) = s''i(xi) = pi für i = 2, ..., N - 1 gilt, mit p1 = pN = 0. Wenn die Werte von p1, ..., pN bekannt sind, können alle Koeffizienten a, b, c, d für die Abschnitte des Splines berechnet werden, da für jeden Abschnitt des Splines vier Gleichungen mit vier Unbekannten vorliegen: Für i = 1, 2, ..., N - 1 muß gelten

si(xi) = yi
si(xi+1) = yi+1
s''i(xi) = pi
s''i(xi+1) = pi+1.

Die x-Werte und y-Werte sind gegeben; um den Spline vollständig zu bestimmen, müssen wir nur die Werte von p2, ..., pN-1 berechnen. Zu diesem Zweck benutzen wir die Bedingung, daß die ersten Ableitungen übereinstimmen müssen; diese N - 2 Bedingungen liefern genau die N - 2 Gleichungen, die benötigt werden, um sie nach den N - 2 Unbekannten aufzulösen, den Werten der zweiten Ableitung pi.

Wenn wir die Koeffizienten a, b, c und d über die Werte der zweiten Ableitung p ausdrücken und dann diese Ausdrücke in die vier oben aufgeführten Gleichungen für jeden Abschnitt des Splines einsetzen würden, würde dies zu einigen unnötig komplizierten Ausdrücken führen. Stattdessen ist es zweckmäßig, die Gleichungen für die Abschnitte des Splines in einer gewissen kanonischen Form darzustellen, die weniger unbekannte Koeffizienten enthält. Wenn wir Variablen mit Hilfe von t = (x - xi)/(xi+1 - xi) substituieren, kann der Spline wie folgt formuliert werden:

si(t) = tyi+1 + (1 - t)yi + (xi+1 - xi)2((t3 - t)pi+1 - ((1 - t)3 - (1 - t))pi)/6.

Nunmehr ist jeder Spline auf dem Intervall [0,1] definiert. Diese Gleichung ist weniger schwierig, als es den Anschein hat, da für uns hauptsächlich die Endpunkte 0 und 1 von Interesse sind, und in diesen Punkten hat entweder t oder (1 - t) den Wert 0. Durch diese Darstellung kann sehr leicht überprüft werden, daß der Spline durch die gegebenen Punkte verläuft und stetig ist, denn es gilt si-1(1) = si(0) = yi für i = 2,..., N - 1, und es ist nur wenig schwieriger, sich davon zu überzeugen, daß auch die zweite Ableitung stetig ist, da s''i(1) = s''i+1(0) = pi+1 gilt. Dies sind kubische Polynome, die der geforderten Bedingung in den Endpunkten genügen, so daß sie zu den oben beschriebenen Abschnitten des Splines äquivalent sind. Wenn wir für t den zugehörigen Ausdruck einsetzen und die Koeffizienten von x3 usw. bestimmen würden, so würden wir die gleichen Ausdrücke für die a, b, c und d, ausgedrückt über die x, y und p, erhalten, wie wenn wir die im vorangegangenen Absatz beschriebene Methode anwenden würden. Es gibt jedoch keinen Grund, das zu tun, da wir nachgeprüft haben, daß diese Abschnitte des Splines den Randbedingungen genügen, und wir können jeden von ihnen in jedem beliebigen Punkt des zugehörigen Intervalls berechnen, indem wir t berechnen und die obige Formel benutzen (nachdem wir die p kennen).

Um nach den p aufzulösen, müssen wir die ersten Ableitungen der Abschnitte des Splines in den Endpunkten gleichsetzen. Die erste Ableitung (nach x) der obigen Gleichung lautet

s'i(t) = zi + (xi+1 - xi)((3t2 - 1)pi+1 + (3(1 - t)2 - 1)pi)/6,

wobei zi = (yi+1 - yi)/(xi+1 - xi) gilt. Indem wir nun s'i-1(1) = s'i(0) für i = 2, ..., N - 1 setzen, ergibt sich unser System von N - 2 Gleichungen:

(xi - xi-1)pi-1 + 2(xi+1 - xi-1)pi + (xi+1 - xi)pi+1 = 6(z - zi-1).

Dieses Gleichungssystem hat eine einfache tridiagonale Form und kann daher, wie wir in Kapitel 37 gesehen haben, leicht mit einer modifizierten Variante des Gaußschen Eliminationsverfahrens gelöst werden.

Wenn wir die Bezeichnungen ui=xi+1 - xi, di = 2(xi+1 - xi-1) und wi = 6(zi - zi-1) einführen, so erhalten wir zum Beispiel für N = 7 das folgende Gleichungssystem:

Formel

Dabei handelt es sich um ein symmetrisches tridiagonales System, bei dem die Diagonale unter der Hauptdiagonalen mit der Diagonale über der Hauptdiagonalen übereinstimmt. Es erweist sich, daß kein Pivotieren mit Hilfe des größten zur Verfügung stehenden Elements erforderlich ist, um für dieses Gleichungssystem eine exakte Lösung zu erhalten.

Das im obigen Absatz beschriebene Verfahren zur Berechnung eines kubischen Splines läßt sich sehr leicht in Pascal formulieren:

    procedure makespline;
      var i: integer;
      begin
      readln(N);
      for i:=1 to N do readln(x[i],y[i]);
      for i:=2 to N-1 do d[i]:=2*(x[i+1]-x[i-1]);
      for i:=1 to N-1 do u[i]:=x[i+1]-x[i];
      for i:=2 to N-1 do
        w[i]:=6.0*((y[i+1]-y[i])/u[i]-(y[i]-y[i-1])/u[i-1]);
      p[1]:=0.0; p[N]:=0.0;
      for i:=2 to N-2 do
        begin
        w[i+1]:=w[i+1]-w[i]*u[i]/d[i];
        d[i+1]:=d[i+1]-u[i]*u[i]/d[i]
        end;
      for i:=N-1 downto 2 do
        p[i]:=(w[i]-u[i]*p[i+1])/d[i];
      end;

Die Felder d und u dienen zur Darstellung der tridiagonalen Matrix, die unter Verwendung des Programms aus Kapitel 37 gelöst wird. Wir verwenden d[i] dort, wo in jenem Programm a[i,i] benutzt wurde, u[i] dort, wo a[i+1,i] oder a[i,i+1] benutzt wurde, und w[i] dort, wo a[i,N+1] benutzt wurde.

Eigenschaft 38.1 Ein kubischer Spline für N Punkte kann in linearer Zeit berechnet werden.

Diese Tatsache ist aufgrund des Programms offensichtlich, welches einfach eine Folge linearer Durchläufe durch die Daten darstellt. w.z.b.w.

Als Beispiel für die Konstruktion eines kubischen Splines betrachten wir die Anpassung eines Splines an die sechs Datenpunkte

(1,0; 2,0), (2,0; 1,5), (4,0; 1,25), (5,0; 1,2), (8,0; 1,125), (10,0; 1,1)

(Diese Punkte gehören zu der Funktion 1 + 1/x.) Die Parameter des Splines werden durch Lösung des Gleichungssystems

Formel

ermittelt; man erhält p2 = 0,39541, p3 = -0,06123, p4 = 0,02658, p5 = -0,00047.

Um den Wert des Splines für einen beliebigen Wert von x aus dem Intervall [x1,xN] zu berechnen, bestimmen wir einfach das Intervall [xi,xi+1], das x enthält, berechnen dann t und benutzen die obige Formel für si(x) (in der wiederum die für pi und pi+1 berechneten Werte verwendet werden).

    function eval(v: real): real;
      var t: real; i: integer;
      function f(x: real): real;
        begin f:=x*x*x-x end;
      begin
      i:=0; repeat i:=i+1 until v<=x[i+1];
      t:=(v-x[i])/u[i];
      eval:=t*y[i+1]+(1-t)*y[i]+u[i]*u[i]*(f(t)*p[i+1]+f(1-t)*p[i])/6.0
      end;

Dieses Programm überprüft keine Fehlerbedingung, wenn v nicht zwischen x[1] und x[N] liegt. Falls die Anzahl der Abschnitte des Splines groß ist (das heißt, falls N groß ist), könnte eines der effizienteren Suchverfahren aus Kapitel 14 angewandt werden, um das Intervall zu finden, das v enthält.

Es gibt viele Varianten der Idee der Kurvenanpassung durch Zusammensetzen von Polynomen auf eine »glatte« Weise; die Berechnung von Splines ist ein sehr gründlich erforschtes Gebiet. Andere Typen von Splines beruhen auf anderen Arten von Glattheitskriterien sowie auf Änderungen wie der Abschwächung der Forderung, daß der Spline exakt durch jeden Datenpunkt verlaufen muß. Was die Berechnung anbelangt, so erfordern sie genau die gleichen Schritte zur Bestimmung der Koeffizienten für jeden der Abschnitte des Splines durch Lösung des linearen Gleichungssystems, das aus der Formulierung von Bedingungen bezüglich der Art und Weise ihrer Zusammensetzung abgeleitet wird.


[ Inhaltsverzeichnis ] [ vorhergehende Seite ] [ nächste Seite ] [ Stichwort ]