Spline-Interpolation - LinkFang.de





Spline-Interpolation


Bei der Spline-Interpolation versucht man, gegebene Stützstellen, auch Knoten genannt, mit Hilfe stückweise stetiger Polynome, genauer Splines, zu interpolieren. Während das Ergebnis einer Polynominterpolation durch unvorteilhaft festgelegte Stützstellen oft bis zur Unkenntlichkeit oszilliert, liefert die Splineinterpolation brauchbare Kurvenverläufe und Approximationseigenschaften (Rungephänomen). Die Spline-Interpolation lässt sich mit geringem linearem Aufwand berechnen, liefert aber im Vergleich zur Polynominterpolation eine geringere Konvergenzordnung.

Vorlage für die Splineinterpolation (dritten Grades) ist das traditionelle, biegsame Lineal der Schiffbauer, die Straklatte (englisch Spline). Diese wird an beliebig vielen, vom Konstrukteur vorgegebenen Punkten fixiert und verbindet die Punkte dann durch eine glatte und harmonische Biegelinie. Die Straklatte erzeugt so die Linie durch alle Punkte mit minimaler Biegeenergie und kleinsten Krümmungen. Während bei der Straklatte die Wendestellen (Orte maximaler Linearität und minimaler Biegeenergie) in der Regel zwischen den Stützstellen liegen und die Stützstellen selbst Orte maximaler Krümmung sind (Orte maximaler Kraft durch Fixierung), liegen die Wendestellen bei der Polynomeninterpolation nahe an den Stützstellen, bei der polynomialen Bestapproximation sogar in den Stützstellen.

Die Begriffe Splineinterpolation und Splinefunktion ohne weitere Zusätze bezeichnen immer die Splineinterpolation/Splinefunktion dritten Grades. Beide Begriffe werden zumeist synonym verwendet. Der Begriff Spline wird jedoch zunehmend als Abkürzung für B-Spline, seltener auch für andere splineartige Linien wie die Bézierkurven, benutzt.

Einfacher Ansatz (Streckenzug)

Die einfachste Methode ist die Verwendung von Geraden zwischen jeweils zwei benachbarten Punkten, die Berechnung eines einfachen Splines als Streckenzug erfolgt auf dieselbe Weise, mit der man auch den Graphen zwischen zwei Punkten ermittelt:

[math]\begin{align} s(x) &= m \cdot x + b \\ &= \underbrace{\frac{y_2 - y_1}{x_2 - x_1}}_{=\,m} \cdot x + \underbrace{y_1 - \frac{y_2 - y_1}{x_2 - x_1} \cdot x_1}_{=\,b \,=\, y_{_1}-m\cdot x_{_1}} \end{align}[/math]

Es ist klar, dass diese „einfachen“ Spline-Polynome – wie oben angesprochen – sehr ungenau sein können. Wesentlich bessere Ergebnisse liefern kubische Spline-Polynome.

Der kubische C2-Spline

Ein kubischer Spline [math]S_\Delta[/math] ist ein Spline, dessen Einschränkungen [math]s_j[/math] auf die Teilintervalle [math][x_{j-1},x_j][/math] (also zwischen zwei Stützstellen) kubische Funktionen sind. Kubische Splines sind zweimal stetig differenzierbar ([math]C^2[/math]) und erfüllen eine Minimalitätsbedingung für die zweite Ableitung, was sie gegenüber anderen Splines besonders interessant macht.

Zur Interpolation der Funktion [math]f[/math] fordert man nun [math]S_\Delta(x_j) = y_j = f(x_j),\ j=0,1, \ldots , n[/math]. Die kubischen Splines eignen sich auf Grund ihrer Glattheit gut zur Approximation von „glatten“ Funktionen. Auf Grund ihrer Konstruktion neigen sie im Gegensatz zu Interpolationspolynomen weniger zu Überschwingern.

Auf jedem Teilintervall wählt man nun das Polynom [math]s_j(x)[/math] in Newtondarstellung, um die Spline-Interpolation anzusetzen.

[math]s_j(x) = a_j + b_j \cdot (x - x_j) + c_j \cdot (x-x_j)^2 + d_j \cdot (x-x_j)^3[/math] für [math]x_{j-1}\leq x\leq x_j[/math] und [math]j=1,\ldots,n[/math]

Um das Gleichungssystem eindeutig zu lösen, werden [math]4n[/math] Bedingungen benötigt. Für jedes der [math]n[/math] Intervalle sind zwei Interpolationsbedingungen zu erfüllen:

[math]\begin{align} s_j(x_{j-1}) = y_{j-1} \qquad j=1,\ldots,n\\ s_j(x_j) = y_j \qquad \qquad j=1,\ldots,n \end{align}[/math]

Dadurch entstehen [math]2n[/math] Bedingungen. Weitere [math]2n-2[/math] Bedingungen erhält man dadurch, dass der Spline an allen [math]n-1[/math] inneren Stützstellen zweimal stetig differenzierbar sein muss:

[math]\begin{align} s'_j(x_j) = s'_{j+1}(x_j) \qquad j=1,\ldots,n-1\\ s''_j(x_j) = s''_{j+1}(x_j) \qquad j=1,\ldots,n-1 \end{align}[/math]

Für die restlichen 2 Bedingungen (Randbedingungen) gibt es verschiedene Möglichkeiten, so z. B.:

  • freier Rand oder natürlicher Spline: [math]s_1''(x_0)=0,s_n''(x_n)=0[/math].
  • eingespannter Rand: [math]s_1'(x_0)=y_0',s_n'(x_n)=y_n'[/math], wobei [math]y_0'[/math] und [math]y_n'[/math] vorgegeben, normalerweise entweder durch die Ableitung der zu interpolierenden Funktion [math]f[/math] oder durch eine Approximation derselben.
  • periodische Randbedingung: [math]s_1'(x_0)=s_n'(x_n)[/math], [math]s_1''(x_0)=s_n''(x_n)[/math]
  • not-a-knot (verwendet zum Beispiel vom Programm Matlab): Die äußeren drei Punkte werden je durch ein Polynom interpoliert.

Die erste Ableitung (Steigung) sieht so aus:

[math]s_j'(x) = b_j + {2 \cdot c_j} \cdot (x-x_j) + {{3 \cdot d_j} \cdot (x-x_j)^2}[/math]

Die zweite Ableitung (Krümmung) sieht so aus:

[math]s_j''(x) = 2 \cdot c_j + {6 \cdot d_j} \cdot (x-x_j)[/math]

Obiger naheliegender Ansatz erfordert einige Rechenarbeit.

Baryzentrische Koordinaten

Zur Vereinfachung werden innerhalb des Intervalls [math][x_{i-1},x_i][/math] für [math]i=1,\dots,n[/math] die Hilfsgrößen [math]h_i:=x_i-x_{i-1}[/math] und die baryzentrischen Variablen [math]u:=\tfrac{x-x_{i-1}}{h_i}[/math] und [math]v:=1-u=\tfrac{x_i-x}{h_i}[/math] eingeführt. Jedes Teilintervall wird hierbei auf das Intervall [math][0,1][/math] linear transformiert. Für das resultierende kubische Polynom wird die Notation [math]\sigma_i(u,v)[/math] verwendet. Die Ableitungen nach [math]x[/math] ergeben sich als [math]u'(x)=h_i^{-1}[/math] und [math]v'(x)=-h_i^{-1}[/math]. Es gilt nach der Kettenregel: [math]s'_i(x) = \sigma'_i h_i^{-1}[/math]. Die Stützstelle [math]x_i[/math] ist rechter Rand im Intervall [math][x_{i-1},x_i][/math] und linker Rand in [math][x_i,x_{i+1}][/math]. Das ergibt auf [math][x_{i-1},x_i][/math] den neuen Ansatz

[math]\sigma_i(u,v)=:y_iu^3+p_iu^2v+q_iuv^2+y_{i-1}v^3[/math] mit [math]0\le u,v\le 1[/math],

durch den die Stetigkeit [math]S \in C^0[a,b][/math] bereits hergestellt ist.

Die ersten Ableitungen sind

[math]s_i'(x)=h_i^{-1}((3y_i-p_i)u^2+2(p_i-q_i)uv+(q_i-3y_{i-1})v^2)[/math].

Nun sei zur Abkürzung [math]s_i'(x_i)=:g_i[/math] für [math]i=0,1,\dots,n[/math] mit den noch unbekannten Parametern [math]g_i[/math]. Die Stetigkeit der ersten Ableitung [math]S \in C^1[a,b][/math] wird über die Anschlussbedingungen hergestellt: [math]s_i'(x_i)=s_{i+1}'(x_i)=g_i[/math] für [math]i=1,2,\dots,n-1[/math]. Somit gilt: [math]h_i^{-1}(3y_i-p_i)=g_i=h_{i+1}^{-1}(q_{i+1}-3y_i)[/math]; die unbekannten Parameter [math]p_i[/math] und [math]q_i[/math] können nun durch [math]g_i[/math] ausgedrückt werden:

[math] \begin{align}p_i&=3y_i-h_ig_i \\q_i&=3y_{i-1}+h_ig_{i-1}\end{align}[/math].

Über die zweiten Ableitungen [math]s_i''(x)=h_i^{-2}(6y_i-4p_i+2q_i)u+h_i^{-2}(6y_{i-1}-4q_i+2p_i)v[/math] werden die Anschlussbedingungen für [math]S \in C^2[a,b][/math] hergestellt:

[math]h_i^{-2}(6y_i-4p_i+2q_i)=h_{i+1}^{-2}(6y_i-4q_{i+1}+2p_{i+1})[/math].

Nach Einsetzen von [math]p_i[/math] und [math]q_i[/math] wird daraus:

[math]h_{i+1}^2(6y_i-12y_i+4h_ig_i+6y_{i-1}+2h_ig_{i-1})=h_i^2(6y_i-12y_i-4h_{i+1}g_i+6y_{i+1}-2h_{i+1}g_{i+1})[/math].

Dies führt zu dem tridiagonalen, streng diagonaldominanten, linearen Gleichungssystem

[math]h_{i+1}g_{i-1}+2(h_i+h_{i+1})g_i+h_ig_{i+1}=3\left(\frac{h_{i+1}}{h_i}(y_i-y_{i-1})+\frac{h_i}{h_{i+1}}(y_{i+1}-y_i)\right),\;i=1,2,\dots,n-1[/math]

mit noch zwei freien Parametern, etwa [math]g_0[/math] und [math]g_n[/math], die zur eindeutigen Lösbarkeit noch festgelegt werden müssen.

Bei äquidistanten Stützstellen mit Abstand [math]h[/math] vereinfacht sich das Gleichungssystem zu

[math]g_{i-1}+4g_i+g_{i+1}=3(y_{i+1}-y_{i-1})/h,\;i=1,2,\dots,n-1[/math].

Hier lässt sich das Gleichungssystem folgendermaßen symmetrisch verlängern, um eine leichte Invertierbarkeit der Matrix sicherzustellen. [Quelle] [math]\begin{align} 4g_0+g_1 & = 3y_1/h & (\text{Zeile}~ i=0) \\ g_{n-1}+4g_n & = -3y_{n-1}/h & (\text{Zeile}~ i=n) \end{align}[/math]

was der Parametrisierung [math]y_{-1}=0=y_{n+1}[/math] entspricht und zu der Tridiagonal-Toeplitz-Matrix

[math]A := \begin{bmatrix} 4 & 1 & 0 & \cdots & 0 & 0 & 0\\ 1 & 4 & 1 & \cdots & 0 & 0 & 0\\ 0 & 1 & 4 & \cdots & 0 & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots\\ 0 & 0 & 0 & \cdots & 4 & 1 & 0\\ 0 & 0 & 0 & \cdots & 1 & 4 & 1\\ 0 & 0 & 0 & \cdots & 0 & 1 & 4\\ \end{bmatrix} \in \R^{n+1\times n+1}[/math]

führt. Diese hat die Inverse

[math]B := b_{n+1}^{-1}\begin{bmatrix} b_nb_0 & -b_{n-1}b_0 & \cdots & (-1)^{n-1}b_1b_0 & (-1)^nb_0b_0\\ -b_{n-1}b_0 & b_{n-1}b_1 & \cdots & (-1)^nb_1b_1 & (-1)^{n-1}b_0b_1\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ (-1)^{n-1}b_1b_0 & (-1)^nb_1b_1 & \cdots & b_1b_{n-1} & -b_0b_{n-1}\\ (-1)^nb_0b_0 & (-1)^{n-1}b_0b_1 & \cdots & -b_0b_{n-1} & b_0b_n\\ \end{bmatrix}[/math]

mit Koeffizienten [math]b_i[/math], die den Gleichungen [math]b_0:=1[/math], [math]b_1:=4[/math], der Rekursion [math]b_i:=4b_{i-1}-b_{i-2}[/math] und explizit der Formel

[math]b_i=(\sqrt{3}(2+\sqrt{3})^{i+1}-\sqrt{3}(2-\sqrt{3})^{i+1})/6[/math]

genügen. Die Übertragung auf Räume mit Dimension größer als eins, zum Beispiel auf Rechteckgitter, ist möglich.

Eigenschaften

Im Jahr 1957 bewies Holladay[1] die folgende nach ihm benannte Identität von Holladay. Mit [math]K^2([a,b])[/math] wird der Raum der zweimal differenzierbaren Funktionen bezeichnet, für welche die nullte und erste Ableitung absolutstetig sind und die zweite Ableitung in [math]L^2([a,b])[/math] liegt. Sei [math]S_\Delta [/math] eine interpolierende Splinefunktion zu [math]f \in K^2([a,b])[/math] und [math]\|.\|[/math] die [math] L^2 [/math]-Norm, so gilt

[math] \|f'' - S''_\Delta\|^2 = \| f'' \|^2 - \|S''_\Delta\|^2 - 2 D[/math]

mit [math]D := \left. (f'(x) - S'_\Delta(x))S''_\Delta(x) \right|_a^b - \sum_{i=1}^n \left. (f(x) - S_\Delta(x)) S'''_\Delta \right|_{x_{i-1}}^{x_i} [/math].

Erfüllt die Splinefunktion die natürlichen, periodischen oder vollständigen Randbedingungen, so ist [math]D = 0 [/math], also

[math] \|f'' - S''_\Delta\|^2 = \|f''\|^2 - \|S''_\Delta\|^2 \geq 0 \Rightarrow \|f''\|^2 \geq \|S''_\Delta \|^2.[/math]

Minimalität der Splineinterpolierenden: [math]S_\Delta[/math] sei [math]\in K^2[/math] und erfülle eine der drei Randbedingungen, dann gilt

[math] \|f'' - S''_\Delta\|^2 = \min_{T \in K^2}\|f'' - T''\|^2.[/math]

Interpolation mit Formerhaltung

Splines sind aufgrund ihrer Eigenschaften im CAD weit verbreitet. Es stellt sich die Frage, unter welchen Bedingungen eine Spline-Interpolante eine der folgenden formerhaltenden Eigenschaften der zu interpolierenden Funktion [math]f:[a,b]\to\mathbb R[/math] erbt:

  • Nichtnegativität: [math]f(x)\ge 0[/math] für alle [math]x[/math]
  • Monotonie: [math]f(x)\le f(y)[/math] für [math]a\le x\le y\le b[/math]
  • Konvexität: [math]f(\lambda x+(1-\lambda) y)\le\lambda f(x)+(1-\lambda) f(y)[/math] für alle [math]x,y\in[a,b][/math] und [math]\lambda\in[0,\,1][/math]

Hier zeigt sich, dass klassische Splines etwas schlechtere Eigenschaften haben als Bézierkurven. Zunächst stellt sich die Frage, wann ein interpolierender Spline konvex ist.

Für klassische Splines gilt, dass die Menge möglicher Splines auf dem Intervall [math][a,b][/math] zum Gitter [math]\Delta:a=x_0\ltx_1\lt\dots\ltx_n=b[/math] ein endlichdimensionaler Vektorraum ist. Für die Interpolation werden (nicht notwendig mit dem Gitter zusammenfallenden) Knoten [math]a=t_0\ltt_1\lt\dots\ltt_m=b[/math] und zugehörige Ordinaten [math]f_0,\dots,f_m\in\mathbb R[/math] vorgegeben und gefordert, dass der Spline [math]s[/math] stetig differenzierbar in [math](a,b)[/math] ist und darüber hinaus [math]\displaystyle s(t_k)=f_k[/math] für [math]k=0,1,\dots,m[/math] gilt. Fordert man zusätzlich die Konvexität des interpolierenden Splines und geringe technische Annahmen, so stellt man fest, dass die Menge [math]Y[/math] aller Ordinatentupel [math](f_0,\dots,f_m)[/math], für die ein solcher Spline existiert, abgeschlossen ist[2].

Das hat weitreichende Konsequenzen. [math]Y[/math] ist eine echte Teilmenge des [math]{\mathbb R}^{m+1}[/math], falls [math]m\ge 2[/math], da die Eingangsdaten nicht in konvexer Lage zu sein brauchen. Bei Vorgabe eines Tupels auf dem Rand von [math]Y[/math] kann infolge Rechenungenauigkeiten oder anderer Störungen die Menge [math]Y[/math] verlassen worden sein, so dass trotz Lösbarkeit des Ausgangsproblems keine Lösung gefunden wird. Die andere Folgerung des Satzes ist noch schlimmer. Dazu seien fünf Punkte in Form des Zeichens „[math]\lor[/math]“ so angeordnet, dass der mittlere Punkt genau auf der Spitze liegt. Die einzige konvexe Interpolierende ist dann die Betragsfunktion, und diese ist nicht stetig differenzierbar. Also gehört das 5-Tupel zum Komplement von [math]Y[/math], und dieses ist offen. Somit gibt es eine Umgebung des 5-Tupels, in der es ebenfalls keine konvexe, stetig differenzierbare Interpolierende gibt. Verschiebt man den mittleren Punkt geringfügig nach oben, ohne die Umgebung zu verlassen, dann erhält man folglich fünf Punkte in streng konvexer Lage, zu denen dennoch die Interpolationsaufgabe keine Lösung besitzt. Da dieser Effekt bei Vorgabe vieler Interpolationspunkte zunimmt, bleibt nur ein Ausweg, die Lösbarkeit für Eingangsdaten in streng konvexer Lage zu gewährleisten, nämlich die Voraussetzungen des Satzes zu verletzen. Die Menge, aus der die Splines entnommen werden dürfen, soll kein endlichdimensionaler Vektorraum sein. Dafür bieten sich u. a. an:

  • (gebrochen-)rationale Splines
  • Splines mit frei wählbaren Zwischenknoten
  • Exponentialsplines
  • lakunäre (lückenhafte) Splines

Literatur

  • Stoer, Bulirsch: Numerische Mathematik 1. 10. Auflage. Springer Verlag, Berlin, Heidelberg, New York 2007, ISBN 978-3-540-45389-5, 2.5 Spline-Interpolation, S. 112–148 (mit Beispielen, Beweisen, Übungsaufgaben und umfangreichen Angaben zu weiterer speziellerer Literatur).

Weblinks

Einzelnachweise

  1. John C. Holladay, 1928–1986. http://texts.cdlib.org/view?docId=hb767nb3z6&doc.view=frames&chunk.id=div00051&toc.depth=1&toc.id=
  2. Jochen W. Schmidt: Staircase Algorithm and Construction of Convex Spline Interpolants up to the Continuity [math]C^3[/math] Computers and Mathematics with Applications, Volume 31, Number 4, February 1995, pp. 67-79.

Kategorien: Numerische Mathematik

Quelle: Wikipedia - http://de.wikipedia.org/wiki/Spline-Interpolation (Vollständige Liste der Autoren des Textes [Versionsgeschichte])    Lizenz: CC-by-sa-3.0

Änderungen: Alle Bilder mit den meisten Bildunterschriften wurden entfernt. Ebenso alle zu nicht-existierenden Artikeln/Kategorien gehenden internen Wikipedia-Links (Bsp. Portal-Links, Redlinks, Bearbeiten-Links). Entfernung von Navigationsframes, Geo & Normdaten, Mediadateien, gesprochene Versionen, z.T. ID&Class-Namen, Style von Div-Containern, Metadaten, Vorlagen, wie lesenwerte Artikel. Ansonsten sind keine Inhaltsänderungen vorgenommen worden. Weiterhin kann es durch die maschinelle Bearbeitung des Inhalts zu Fehlern gerade in der Darstellung kommen. Darum würden wir jeden Besucher unserer Seite darum bitten uns diese Fehler über den Support mittels einer Nachricht mit Link zu melden. Vielen Dank!

Stand der Informationen: August 201& - Wichtiger Hinweis: Da die Inhalte maschinell von Wikipedia übernommen wurden, ist eine manuelle Überprüfung nicht möglich. Somit garantiert LinkFang.de nicht die Richtigkeit und Aktualität der übernommenen Inhalte. Sollten die Informationen mittlerweile fehlerhaft sein, bitten wir Sie darum uns per Support oder E-Mail zu kontaktieren. Wir werden uns dann innerhalb von spätestens 10 Tagen um Ihr Anliegen kümmern. Auch ohne Anliegen erfolgt mindestens alle drei Monate ein Update der gesamten Inhalte.