Splineinterpolation

Bei der Spline-Interpolation versucht man, eine Funktion mit Hilfe von Splines zu interpolieren. Im Gegensatz zur Polynominterpolation oszillieren die resultierenden Funktionen deutlich weniger stark, was auf einfache Weise bessere Approximationseigenschaften erlaubt. Der Preis dafür ist, dass die Interpolante nicht mehr unendlich oft differenzierbar ist.

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:

$s(x) = m \cdot x + b$, d.h.

$s(x) = \frac{y_2 - y_1}{x_2 - x_1} \cdot x + y_1 - \frac{y_2 - y_1}{x_2 - x_1} \cdot x_1$

(mit $m = \frac{y_2 - y_1}{x_2 - x_1}$ und $b = y_1 - m \cdot x_1$)

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 C²-Spline

Kubische Splines sind Splines, die auf jedem Teilintervall [x_i,x_{i + 1}] (also zwischen zwei Stützstellen) mit einem kubischen Polynom übereinstimmen. Sie erfüllen eine Minimalitätsbedingung und sind deshalb gegenüber anderen Splines besonders interessant.

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

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

$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$ für $x_j\leq x\leq x_{j+1}$ und $j=0,\ldots,n-1$

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

$s_j(x_j) = f_j \qquad \qquad j=0,1,\ldots,n-1$

$s_j(x_{j+1}) = f_{j+1} \qquad j=0,1,\ldots,n-1$

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

$s'_{j-1}(x_j) = s'_{j}(x_j) \qquad j=1,\ldots,n-1$

$s''_{j-1}(x_j) = s''_{j}(x_j) \qquad j=1,\ldots,n-1$

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

• freier Rand oder natürlicher Spline: s₀''(x₀) = 0,s_{n − 1}''(x_n) = 0
• eingespannter Rand: s₀'(x₀) = f₀',s_{n − 1}'(x_n) = f_n' wobei f₀' und f_n' vorgegeben, normalerweise entweder durch die Ableitung der zu interpolierenden Funktion f oder durch eine Approximation derselben.
• periodische Randbedingung: s₀'(x₀) = s_{n − 1}'(x_n), s₀''(x₀) = s_{n − 1}''(x_n)
• 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:

$s_j'(x) = b_j + {2 \cdot c_j} \cdot (x-x_j) + {{3 \cdot d_j} \cdot (x-x_j)^2}$

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

$s_j''(x) = 2 \cdot c_j + {6 \cdot d_j} \cdot (x-x_j)$

Obiger naheliegender Ansatz erfordert einige Rechenarbeit. Man spart einiges, indem man baryzentrische Koordinaten verwendet. Innerhalb des Intervalls [x_{i − 1},x_i] werden dazu $\displaystyle h_i:=x_i-x_{i-1}$, $u:=\frac{x-x_{i-1}}{h_i}$ und $v:=1-u=\frac{x_i-x}{h_i}$ festgelegt. Das ergibt auf [x_{i − 1},x_i] den Ansatz s(x) = f_{i − 1}v³ + f_iu³ + a_iu²v + b_iuv² mit $0\le u\le 1$, so dass die Stetigkeit bereits gewährleistet ist.

Die C¹-Bedingung liefert zunächst h_is'(x) = (b_i − 3f_{i − 1})v² + (3f_i − a_i)u² + 2(a_i − b_i)uv. Nun sei s'(x_i) = p_i für $i=0,1,\dots,n$ mit den noch unbekannten Parametern p_i. Aus den Anschlussbedingungen s'(x_i − 0) = s'(x_i + 0) = p_i für $i=1,2,\dots,n-1$ folgt somit 3f_i − a_i = b_{i + 1} − 3f_i, also a_i = 3f_i − h_ip_i und b_i = 3f_{i − 1} + h_ip_{i − 1}.

Die zweiten Ableitungen ergeben $h_i^2s''(x)=(6f_{i-1}-4b_i+2a_i)v+(6f_i-4a_i+2b_i)u$ und schließlich das tridiagonale, streng diagonaldominante, lineare Gleichungssystem

$h_ip_{i-1}+2(h_i+h_{i+1})p_i+h_{i+1}p_{i+1}=6(\frac{f_{i+1}-f_i}{h_{i+1}}-\frac{f_i-f_{i-1}}{h_i}),\;i=1,2,\dots,n-1$

mit noch zwei freien Parametern, etwa p₀ und p_n. Eine Übertragung z. B. auf Rechteckgitter bereitet keine Schwierigkeiten.

Bei äquidistanten Stützstellen mit Abstand h vereinfacht sich das Gleichungssystem zu

$p_{i-1}+4p_i+p_{i+1}=\frac{6}{h^2}(f_{i+1}-2f_i+f_{i-1}),\;i=1,2,\dots,n-1$

Eigenschaften

Identität von Halliday: Sei S_Δ eine interpolierende Splinefunktion zu $f \in K^2$ und $\|.\|$ die L²-Norm , so gilt

$\|f''-S''_\Delta\| = \| f'' \|^2 - \|S''_\Delta\|^2 - 2\left( \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} \right).$

Erfüllt die Splinefunktion die natürlichen, periodischen oder vollständigen Randbedingungen so gilt

$\|f''-S''_\Delta\|^2 = \|f''\|^2 - \|S''_\Delta\|^2 \geq 0 \Rightarrow \|f''\|^2 \geq \|S''_\Delta \|^2.$

Minimalität der Splineinterpolierenden: Es sei $S_\Delta \in K^2$ und erfüllt eine der drei Randbedingungen so gilt

$\|f'' - S''_\Delta\|^2 = \min_{T \in K^2}\|f'' - T''_\Delta\|^2.$

Interpolation mit Formerhaltung

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

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

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 [a,b] zum Gitter $\Delta:a=x_0<x_1<\dots<x_n=b$ ein endlichdimensionaler Vektorraum ist. Für die Interpolation werden (nicht notwendig mit dem Gitter zusammenfallenden) Knoten $a=t_0<t_1<\dots<t_m=b$ und zugehörige Ordinaten $f_0,\dots,f_m\in\mathbb R$ vorgegeben und gefordert, dass der Spline s stetig differenzierbar in (a,b) ist und darüber hinaus $\displaystyle s(t_k)=f_k$ für $k=0,1,\dots,m$ gilt. Fordert man zusätzlich die Konvexität des interpolierenden Splines und geringe technische Annahmen, so stellt man fest, dass die Menge Y aller Ordinatentupel $(f_0,\dots,f_m)$, für die ein solcher Spline existiert, abgeschlossen ist^[1].

Das hat weitreichende Konsequenzen. Y ist eine echte Teilmenge des ${\mathbb R}^{m+1}$, falls $m\ge 2$, da die Eingangsdaten nicht in konvexer Lage zu sein brauchen. Bei Vorgabe eines Tupels auf dem Rand von Y kann infolge Rechenungenauigkeiten oder anderer Störungen die Menge Y 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 „$\lor$“ 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 Y, 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

J. Stoer: Numerische Mathematik 1, 9. Auflage, Springer Verlag, 2005, ISBN 3-540-21395-3

Weblinks

• Skript zur Spline-Interpolation, Uni Göttingen
• Online-Tool zum Berechnen von kubischen Splines
• OpenOfficeCalc-Makro zur Berechnung kubischer Splines und deren Anwendung (Wendepunkt-, Maximum- und Minimumsuche)

Einzelnachweis

1. ↑ Jochen W. Schmidt: Staircase Algorithm and Construction of Convex Spline Interpolants up to the Continuity C³ Computers and Mathematics with Applications, Volume 31, Number 4, February 1995, pp. 67-79.