Spline-Interpolation

Bei d​er Spline-Interpolation versucht man, gegebene Stützstellen, a​uch Knoten genannt, m​it Hilfe stückweiser Polynome niedrigen Grades z​u interpolieren. Während d​as Ergebnis e​iner Polynominterpolation d​urch unvorteilhaft festgelegte Stützstellen o​ft bis z​ur Unkenntlichkeit oszilliert, liefert d​ie Splineinterpolation brauchbare Kurvenverläufe u​nd Approximationseigenschaften (Rungephänomen). Die Spline-Interpolation lässt s​ich mit geringem, linearem Aufwand berechnen, liefert a​ber im Vergleich z​ur Polynominterpolation e​ine geringere Konvergenzordnung.

Beispiel eines Splines mit 8 Knoten

Vorlage für d​ie Splineinterpolation (dritten Grades) i​st das traditionelle, biegsame Lineal d​er Schiffbauer, d​ie Straklatte (englisch Spline). Diese w​ird an beliebig vielen, v​om Konstrukteur vorgegebenen Punkten fixiert u​nd verbindet d​ie Punkte d​ann durch e​ine glatte u​nd harmonische Biegelinie. Die Straklatte erzeugt s​o die Linie d​urch alle Punkte m​it minimaler Biegeenergie u​nd kleinsten Krümmungen. Während b​ei der Straklatte d​ie Wendestellen (Orte maximaler Linearität u​nd minimaler Biegeenergie) i​n der Regel zwischen d​en Stützstellen liegen u​nd die Stützstellen selbst Orte maximaler Krümmung s​ind (Orte maximaler Kraft d​urch Fixierung), liegen d​ie Wendestellen b​ei der Polynomeninterpolation n​ahe an d​en Stützstellen, b​ei der polynomialen Bestapproximation s​ogar in d​en Stützstellen.

Die Begriffe Splineinterpolation bzw. Splinefunktion o​hne weitere Zusätze bezeichnen i​mmer die Splineinterpolation bzw. Splinefunktion dritten Grades. Beide Begriffe werden zumeist synonym verwendet. Der Begriff Spline w​ird jedoch zunehmend a​ls Abkürzung für B-Spline, seltener a​uch für andere splineartige Linien w​ie die Bézierkurven, benutzt.

Smoothing Splines s​ind Splines, d​ie nicht d​urch jeden Datenpunkt verlaufen müssen u​nd können z​ur Signalglättung benutzt werden.

Gegebenheiten

Gegeben: eine natürliche Zahl ${\displaystyle n\in \mathbb {N} }$ und ${\displaystyle n+1}$ Stützstellen ${\displaystyle x_{0}<x_{1}<\dotsb <x_{n-1}<x_{n}\in \mathbb {R} }$ sowie ${\displaystyle n+1}$ Funktionswerte ${\displaystyle y_{0},y_{1}\dotsc ,y_{n-1},y_{n}\in \mathbb {R} }$. Gesucht ist eine stückweise polynomiale Funktion, ein Spline,

${\displaystyle S\colon [x_{0},x_{n}]\to \mathbb {R} }$

mit ${\displaystyle S(x_{i})=y_{i}}$ für ${\displaystyle i=0,\dotsc ,n}$, bei der für ${\displaystyle i=0,\dotsc ,n-1}$ die Einschränkungen

${\displaystyle s_{i}:=S|_{[x_{i},x_{i+1}]}\;\colon \;[x_{i},x_{i+1}]\to \mathbb {R} }$

auf die Teilintervalle ${\displaystyle [x_{i},x_{i+1}]}$ Polynome sind.

Linear (einfacher Streckenzug)

Die einfachste Methode i​st die Verwendung v​on Geraden zwischen jeweils z​wei benachbarten Punkten, d​ie Berechnung e​ines einfachen Splines a​ls Streckenzug erfolgt a​uf dieselbe Weise, m​it der m​an auch d​en Graphen zwischen z​wei Punkten ermittelt:

${\displaystyle {\begin{aligned}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{aligned}}}$

oder a​uch

${\displaystyle s(x)={\frac {x_{2}-x}{x_{2}-x_{1}}}\cdot y_{1}+{\frac {x-x_{1}}{x_{2}-x_{1}}}\cdot y_{2}}$

Diese „einfachen“ Spline-Polynome können – w​ie oben angesprochen – s​ehr ungenau sein. Wesentlich bessere Ergebnisse liefern kubische Spline-Polynome.

Kubisch (Polynome 3. Grades)

Bei d​er Verbindung v​on Punkten m​it Polynomen höheren Grades müssen zusätzlich z​u den Stützstellen Eigenschaften definiert werden, w​ie die Polynome ineinander übergehen. Für kubische Splines s​ind in e​iner Dimension 4 Koeffizienten z​u bestimmen u​nd zwei weitere Bedingungen s​ind zu definieren.

Der kubische C²-Spline

C² fordert, dass die zusammengesetzte Funktion ${\displaystyle S}$ aus allen Einschränkungen (Teilintervallen) ${\displaystyle s_{i}}$ zweimal stetig differenzierbar ist. Dafür wird gefordert, dass die erste und zweite Ableitung der Einschränkungen an den Stützstellen

${\displaystyle s'_{i}(x_{i}):=s'_{i-1}(x_{i})}$   und   ${\displaystyle s''_{i}(x_{i}):=s''_{i-1}(x_{i})}$

für ${\displaystyle i=1,\dotsc ,n-1}$ übereinstimmen.

Prinzipiell gilt, dass sich eine Änderung einer Stützstelle ${\displaystyle (x_{i},y_{i})}$ stets global auf den gesamten Spline ${\displaystyle S}$ auswirkt, jedoch wird der Einfluss der Änderung mit zunehmender Distanz zu ${\displaystyle x_{i}}$ – anders als bei Interpolationspolynomen – stark gedämpft. Kubische Splines neigen daher weniger zum Überschwingen.

Der kubische C²-Spline erfüllt e​ine Minimalitätseigenschaft d​er zweiten Ableitung, w​as ihn gegenüber anderen Interpolationen besonders interessant macht.

Konstruktion

Es ist ersichtlich, dass die zweite Ableitung von ${\displaystyle S}$ ein linearer Spline ist. Diese kann daher wie oben beschrieben durch folgende Form beschrieben werden:

${\displaystyle s_{i}''(x)={\frac {x_{i+1}-x}{h_{i}}}\cdot M_{i}+{\frac {x-x_{i}}{h_{i}}}\cdot M_{i+1}}$   mit   ${\displaystyle h_{i}=x_{i+1}-x_{i}}$

für ${\displaystyle i=0,\dotsc ,n-1}$. ${\displaystyle M_{i}}$ sind die sogenannten Momente,[1] welche den Werten von ${\displaystyle S''(x_{i})}$ an den Stützstellen entsprechen und im Folgenden zu berechnen sind. Durch zweifache Integration und geschickte Umformung entstehen aus diesen Gleichungen Polynome dritten Grades mit zwei weiteren Parametern ${\displaystyle c_{i}}$ und ${\displaystyle d_{i}}$ der Form:

${\displaystyle {\frac {1}{6}}\left({\frac {(x_{i+1}-x)^{3}}{h_{i}}}\cdot M_{i}+{\frac {(x-x_{i})^{3}}{h_{i}}}\cdot M_{i+1}\right)+c_{i}\cdot (x-x_{i})+d_{i}}$ .

Um die Stetigkeitsbedingungen ${\displaystyle s_{i}(x_{i})=y_{i}}$ und ${\displaystyle s_{i}(x_{i+1})=y_{i+1}}$ zu erfüllen, wählen wir[2]

${\displaystyle d_{i}=y_{i}-{\frac {h_{i}^{2}}{6}}\cdot M_{i}}$ und ${\displaystyle c_{i}={\frac {y_{i+1}-y_{i}}{h_{i}}}-{\frac {h_{i}}{6}}\cdot (M_{i+1}-M_{i})}$ .

Mit diesem Ansatz stimmen bereits die nullten und die zweiten Ableitungen der Einschränkungen ${\displaystyle s_{i}}$ an den Stützstellen überein. Die Momente sind so zu wählen, dass auch die ersten Ableitungen an den Stützstellen gleich sind. Mit

${\displaystyle s_{i}'(x)={\frac {1}{2}}\left(-{\frac {(x_{i+1}-x)^{2}}{h_{i}}}\cdot M_{i}+{\frac {(x-x_{i})^{2}}{h_{i}}}\cdot M_{i+1}\right)+c_{i}}$ und ${\displaystyle s_{i}'(x_{i})=s_{i-1}'(x_{i})}$

lassen s​ich folgende Gleichungen aufstellen:

${\displaystyle {\frac {h_{i-1}}{6}}\cdot M_{i-1}+{\frac {h_{i-1}+h_{i}}{3}}\cdot M_{i}+{\frac {h_{i}}{6}}\cdot M_{i+1}={\frac {y_{i+1}-y_{i}}{h_{i}}}-{\frac {y_{i}-y_{i-1}}{h_{i-1}}}}$ für ${\displaystyle i=1,\dots ,n-1}$ .

Für ${\displaystyle i=0}$ und ${\displaystyle i=n}$ fehlen hier zwei Gleichungen, welche sich aus den Randbedingungen ergeben.

Dieses lineare Gleichungssystem k​ann auch d​urch folgende, tridiagonale, streng diagonaldominante Matrix ausgedrückt werden:

${\displaystyle {\begin{bmatrix}\mu _{0}&\lambda _{0}&&&&\\{\frac {h_{0}}{6}}&{\frac {h_{0}+h_{1}}{3}}&{\frac {h_{1}}{6}}&&&\\&\ddots &\ddots &\ddots &&\\&&{\frac {h_{i-1}}{6}}&{\frac {h_{i-1}+h_{i}}{3}}&{\frac {h_{i}}{6}}&\\&&&\ddots &\ddots &\ddots \\&&&&\lambda _{n}&\mu _{n}\end{bmatrix}}\cdot {\begin{bmatrix}M_{0}\\M_{1}\\\vdots \\M_{i}\\\vdots \\M_{n}\end{bmatrix}}={\begin{bmatrix}b_{0}\\{\frac {y_{2}-y_{1}}{h_{1}}}-{\frac {y_{1}-y_{0}}{h_{0}}}\\\vdots \\{\frac {y_{i+1}-y_{i}}{h_{i}}}-{\frac {y_{i}-y_{i-1}}{h_{i-1}}}\\\vdots \\b_{n}\end{bmatrix}}}$

Die Werte für die ${\displaystyle \lambda _{i},\mu _{i},b_{i}}$ hängen von den Randbedingungen ab.

Zur Lösung k​ann hier a​uf den komplizierten Gauss-Algorithmus verzichtet werden u​nd z. B. e​in einfacher Vorwärts-Durchlauf z​ur Elimination d​er Elemente u​nter der Hauptdiagonalen m​it anschließender Rückwärtssubstitution verwendet werden (Thomas-Algorithmus).

Randbedingungen

Prinzipiell g​ibt es e​in Interpolationsintervall weniger a​ls Stützstellen (Zaunpfahlproblem). Das heißt, d​ass zwei Gleichungen z​ur Bestimmung a​ller Koeffizienten fehlen. Diese ergeben s​ich aus d​en Randbedingungen. Typische Randbedingungen sind:

• Natürliche Randbedingungen (auch freier Rand)

Bedingung: ${\displaystyle s_{0}''(x_{0})=0}$, ${\displaystyle s_{n-1}''(x_{n})=0}$

Bedeutung: Das Spline schließt mit Wendepunkten ab.

Berechnung: ${\displaystyle \lambda _{0}=\lambda _{n}=b_{0}=b_{n}=0}$ und ${\displaystyle \mu _{0}=\mu _{n}=1}$

• Hermite Randbedingungen (auch eingespannter Rand)

Bedingung: ${\displaystyle s_{0}'(x_{0})=f'(a)}$, ${\displaystyle s_{n-1}'(x_{n})=f'(b)}$

Bedeutung: ${\displaystyle f'(a)}$ und ${\displaystyle f'(b)}$ sind vorgegeben, normalerweise entweder durch die Ableitung einer zu interpolierenden Funktion ${\displaystyle f}$ oder durch eine Approximation derselben.

Berechnung:

${\displaystyle {\begin{matrix}\lambda _{0}={\frac {h_{0}}{6}},&\mu _{0}={\frac {h_{0}}{3}},&b_{0}={\frac {y_{1}-y_{0}}{h_{0}}}-f'(a)\\\lambda _{n}={\frac {h_{n-1}}{6}},&\mu _{n}={\frac {h_{n-1}}{3}},&b_{n}=-{\frac {y_{n}-y_{n-1}}{h_{n-1}}}+f'(b)\end{matrix}}}$

• periodische Randbedingungen

Bedingung: Intervall ${\displaystyle [x_{0},x_{n+1}]}$, ${\displaystyle y_{0}=:y_{n+1}}$, ${\displaystyle s_{0}'(x_{0})=:s_{n}'(x_{n+1})}$, ${\displaystyle s_{0}''(x_{0})=:s_{n}''(x_{n+1})}$

Bedeutung: Nullte, erste und zweite Ableitung von ${\displaystyle S}$ am Anfang und am Ende des Intervalls sind gleich.

Berechnung: Es wird eine zusätzliche Stützstelle ${\displaystyle x_{n+1}}$ eingeführt, welche das Intervall begrenzt. Die Anzahl der Gleichungen zur Berechnung der Momente und die Größe der Matrix bleibt jedoch gleich, da ${\displaystyle M_{n+1}:=M_{0}}$ bereits gegeben ist, damit die zweiten Ableitungen übereinstimmen. Für die erste- und letzte Zeile der Matrix gilt:

${\displaystyle {\begin{matrix}\mu _{0}={\frac {h_{n}+h_{0}}{3}},&\lambda _{0}={\frac {h_{0}}{6}},&b_{0}={\frac {y_{1}-y_{0}}{h_{0}}}-{\frac {y_{0}-y_{n}}{h_{n}}}\\\mu _{n}={\frac {h_{n-1}+h_{n}}{3}},&\lambda _{n}={\frac {h_{n-1}}{6}},&b_{n}={\frac {y_{0}-y_{n}}{h_{n}}}-{\frac {y_{n}-y_{n-1}}{h_{n-1}}}\end{matrix}}}$

Außerdem sind die Ecken der Matrix abseits der Hauptdiagonalen hier nicht Null:

${\displaystyle m_{0,n}=m_{n,0}={\frac {h_{n}}{6}}}$

Die Lösung dieses Systems ist daher komplizierter. Im Falle äquidistanter Stützstellen lässt sich für diesen Fall eine Transformation anwenden.

• not-a-knot Randbedingungen

Bedingung: ${\displaystyle s_{0}'''(x_{1})=s_{1}'''(x_{1})}$, ${\displaystyle s_{n-2}'''(x_{n-1})=s_{n-1}'''(x_{n-1})}$

Bedeutung: Die äußeren drei Punkte werden je durch ein gemeinsames Polynom interpoliert, was zum Beispiel durch Gleichsetzen der dritten Ableitungen erfolgen kann. Für Splines bis einschließlich vier Stützstellen geht der not-a-knot Spline daher in ein gewöhnliches Interpolationspolynom über. Verwendet wird der not-a-knot-Spline zum Beispiel vom Programm Matlab.

Berechnung: Es gelte ${\displaystyle n>4}$. Für die Vorwärtssubstitution wird zunächst die Randbedingung an ${\displaystyle x_{0}}$ betrachtet. Hierfür können folgende Werte verwendet werden:

${\displaystyle {\begin{matrix}\mu _{0}=1,&\lambda _{0}={\frac {h_{1}+2\cdot h_{0}}{h_{1}-h_{0}}},&b_{0}={\frac {-6\cdot h_{0}}{h_{1}^{2}-h_{0}^{2}}}\cdot \left({\frac {y_{2}-y_{1}}{h_{1}}}-{\frac {y_{1}-y_{0}}{h_{0}}}\right).\end{matrix}}}$

Falls ${\displaystyle h_{1}=h_{0}}$ entsteht hier eine Division durch Null und es muss eine Grenzwertbildung bis in die dritte Zeile der Matrix (Null ist hier der Index der ersten Zeile/Spalte) erfolgen:

${\displaystyle {\begin{matrix}m_{2,1}=0,&m_{2,2}={\frac {h_{1}+h_{2}}{3}},&m_{2,3}={\frac {h_{2}}{6}},&b_{2}={\frac {y_{3}-y_{2}}{h_{2}}}-{\frac {y_{2}-y_{1}}{h_{1}}}-{\frac {y_{2}-2\cdot y_{1}+y_{0}}{6\cdot h_{0}}}.\end{matrix}}}$

Gelöst wird nun lediglich die Untermatrix beginnend bei ${\displaystyle m_{2,2}}$ und das Randpolynom ${\displaystyle p=s_{0}=s_{1}}$ auf dem Intervall ${\displaystyle [x_{0},x_{2})}$ wird durch eine zusätzliche Polynominterpolation bestimmt, sodass ${\displaystyle p(x_{1})=y_{1}}$, ${\displaystyle p(x_{2})=y_{2}}$, ${\displaystyle p'(x_{2})=s'_{2}(x_{2})}$, ${\displaystyle p''(x_{2})=s''_{2}(x_{2})}$. ${\displaystyle p(x_{0})=y_{0}}$ ist dadurch bereits erfüllt.

Nach dem Vorwärtsdurchlauf zur Elimination der Elemente unter der Diagonalen wird nun die Randbedingung an ${\displaystyle x_{n}}$ betrachtet:

${\displaystyle {\begin{matrix}\mu _{n}=1,&\lambda _{n}=0,&b_{n}={\frac {(b_{n-1}\cdot m_{n-2,n-1}-b_{n-2}\cdot m_{n-1,n-1}+b_{n-1}\cdot m_{n-2,n-2})\cdot h_{n-1}+b_{n-1}\cdot m_{n-2,n-2}\cdot h_{n-2}}{(m_{n-2,n-1}+m_{n-2,n-2})\cdot m_{n-1,n}\cdot h_{n-1}+(m_{n-2,n-2}\cdot m_{n-1,n}+m_{n-2,n-2}\cdot m_{n-1,n-1})\cdot h_{n-2}}}.\end{matrix}}}$

Weitere Randbedingungen s​ind gebräuchlich, w​ie etwa integrale Randbedingungen[3] o​der die i​m Folgenden vorgestellte, symmetrische Verlängerung.

Optimierung des Rechenaufwands

Bei äquidistanten Stützstellen mit konstantem Abstand ${\displaystyle h}$ vereinfacht sich das Gleichungssystem zu

${\displaystyle M_{i-1}+4\cdot M_{i}+M_{i+1}=6{\frac {y_{i+1}-2\cdot y_{i}+y_{i-1}}{h^{2}}}}$ für ${\displaystyle i=1,\dots ,n-1}$.

Mit der symmetrischen Verlängerung ${\displaystyle \mu _{0}=\mu _{n}=4}$ und ${\displaystyle \lambda _{0}=\lambda _{n}=1}$ entsteht daraus eine sogenannte Tridiagonal-Toeplitz-Matrix, welche besonders effizient, auch parallel gelöst werden kann. Für die rechte Seite kann hier ${\displaystyle b_{0}={\tfrac {y_{1}}{2\cdot h}}}$ und ${\displaystyle b_{n}=-{\tfrac {y_{n-1}}{2\cdot h}}}$ angesetzt werden. Mit ${\displaystyle A={\tfrac {6}{h}}\cdot M}$ lässt sich schreiben:

${\displaystyle 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 \mathbb {R} ^{(n+1)\times (n+1)}}$

Diese h​at die Inverse

${\displaystyle B:=b_{n+1}^{-1}{\begin{bmatrix}b_{n}b_{0}&-b_{n-1}b_{0}&\cdots &(-1)^{n-1}b_{1}b_{0}&(-1)^{n}b_{0}b_{0}\\-b_{n-1}b_{0}&b_{n-1}b_{1}&\cdots &(-1)^{n}b_{1}b_{1}&(-1)^{n-1}b_{0}b_{1}\\\vdots &\vdots &\ddots &\vdots &\vdots \\(-1)^{n-1}b_{1}b_{0}&(-1)^{n}b_{1}b_{1}&\cdots &b_{1}b_{n-1}&-b_{0}b_{n-1}\\(-1)^{n}b_{0}b_{0}&(-1)^{n-1}b_{0}b_{1}&\cdots &-b_{0}b_{n-1}&b_{0}b_{n}\\\end{bmatrix}}}$

mit Koeffizienten ${\displaystyle b_{i}}$, die den Gleichungen ${\displaystyle b_{0}:=1}$, ${\displaystyle b_{1}:=4}$, der Rekursion ${\displaystyle b_{i}:=4b_{i-1}-b_{i-2}}$ und explizit der Formel

${\displaystyle b_{i}=({\sqrt {3}}(2+{\sqrt {3}})^{i+1}-{\sqrt {3}}(2-{\sqrt {3}})^{i+1})/6}$

genügen.

Minimalitätseigenschaft der zweiten Ableitung

Unter allen zweimal stetig differenzierbaren Funktionen, die alle Stützstellen innerhalb eines Intervalls ${\displaystyle [a,b]}$ miteinander verbinden, hat unter Verwendung natürlicher, periodischer oder Hermite-Randbedingungen der kubische Spline die geringste Krümmung: ${\displaystyle \int _{a}^{b}f''(x)^{2}\mathrm {d} x\geq \int _{a}^{b}S''(x)^{2}\mathrm {d} x}$, wobei ${\displaystyle f}$ hier eine beliebige Funktion auf ${\displaystyle [a,b]}$ aus C² ist, die alle Stützstellen schneidet. Anschaulich folgt daraus, dass ein zu volles Glas Wasser, das während einer Zugfahrt auf dem Tisch steht, „am wenigsten überschwappt“, wenn der Streckenzug der Schienenführung mithilfe kubischer Splines parametrisiert wurde.

Diese sogenannte Identität von Holladay wurde im Jahr 1957 von Holladay[4] bewiesen. Sei mit ${\displaystyle C^{2}([a,b])}$ der Raum der zweimal differenzierbaren Funktionen bezeichnet, für welche die nullte und erste Ableitung absolutstetig sind und die zweite Ableitung in ${\displaystyle L^{2}([a,b])}$ liegt. Sei nun ${\displaystyle S}$ eine interpolierende Splinefunktion zu einer beliebigen Funktion ${\displaystyle f\in C^{2}([a,b])}$ und ${\displaystyle \|.\|}$ die ${\displaystyle L^{2}}$-Norm, so gilt: ${\displaystyle \|f''-S''\|^{2}=\|f''\|^{2}-\|S''\|^{2}-2D}$ mit ${\displaystyle D:=\left.(f'(x)-S'(x))S''(x)\right|_{a}^{b}-\sum _{i=1}^{n}\left.(f(x)-S(x))S'''\right|_{x_{i-1}}^{x_{i}}}$.

Erfüllt die Splinefunktion die natürlichen, periodischen oder vollständigen Randbedingungen, so ist ${\displaystyle D=0}$, also: ${\displaystyle \|f''-S''\|^{2}=\|f''\|^{2}-\|S''\|^{2}\geq 0\Rightarrow \|f''\|^{2}\geq \|S''\|^{2}.}$

Damit gilt nun: ${\displaystyle \|f''-S''\|^{2}=\min _{T\in C^{2}}\|f''-T''\|^{2}.}$

Bikubischer C²-Spline

Der bikubische C²-Spline i​st die Verallgemeinerung d​es einfachen, kubischen C²-Splines a​uf zwei Dimensionen, m​an spricht hierbei a​uch von multivariater Interpolation. Dafür müssen d​ie zu interpolierenden Punkte z_ij i​n einem rechteckigen Gitter angeordnet sein.[5] Jede zwischen v​ier Punkten aufgespannte Fläche A_ij w​ird durch e​in zweidimensionales Polynom v​on 16 Koeffizienten charakterisiert:

${\displaystyle z:=A_{ij}(x,y)=\sum _{k=0}^{3}\sum _{l=0}^{3}a_{k,l}^{ij}\cdot x^{k}\cdot y^{l}={\begin{matrix}a_{0,0}^{ij}&+&a_{1,0}^{ij}\cdot x&+&a_{2,0}^{ij}\cdot x^{2}&+&a_{3,0}^{ij}\cdot x^{3}&+\\a_{0,1}^{ij}\cdot y&+&a_{1,1}^{ij}\cdot x\cdot y&+&a_{2,1}^{ij}\cdot x^{2}\cdot y&+&a_{3,1}^{ij}\cdot x^{3}\cdot y&+\\a_{0,2}^{ij}\cdot y^{2}&+&a_{1,2}^{ij}\cdot x\cdot y^{2}&+&a_{2,2}^{ij}\cdot x^{2}\cdot y^{2}&+&a_{3,2}^{ij}\cdot x^{3}\cdot y^{2}&+\\a_{0,3}^{ij}\cdot y^{3}&+&a_{1,3}^{ij}\cdot x\cdot y^{3}&+&a_{2,3}^{ij}\cdot x^{2}\cdot y^{3}&+&a_{3,3}^{ij}\cdot x^{3}\cdot y^{3}&~\end{matrix}}}$

Für einen zweidimensionalen C²-Spline müssen die Koeffizienten so gewählt werden, dass die aus allen Flächen zusammengesetzte Funktion ${\displaystyle S(x,y):\mathbb {R} ^{2}\rightarrow \mathbb {R} }$ zweimal stetig in x- und y-Richtung differenzierbar ist. Das heißt, neben S selbst sind die folgenden Ableitungen stetig auf ganz ${\displaystyle \mathbb {R} ^{2}}$:

${\displaystyle {\begin{matrix}{\frac {\partial S(x,y)}{\partial x}},&{\frac {\partial S(x,y)}{\partial y}},&{\frac {\partial ^{2}S(x,y)}{\partial x^{2}}},&{\frac {\partial ^{2}S(x,y)}{\partial y^{2}}},&{\frac {\partial ^{2}S(x,y)}{\partial x\partial y}},&{\frac {\partial ^{3}S(x,y)}{\partial x^{2}\partial y}},&{\frac {\partial ^{3}S(x,y)}{\partial x\partial y^{2}}},&{\frac {\partial ^{4}S(x,y)}{\partial x^{2}\partial y^{2}}}\end{matrix}}}$

Konstruktion

Veranschaulichung zur Konstruktion eines bikubischen Splines

Es ist ersichtlich, dass jeder Schnitt durch eine Teilfläche ${\displaystyle A_{i,j}}$ parallel zur ${\displaystyle x}$- oder ${\displaystyle y}$-Achse eine eindimensionale Kurve liefert, welche durch ein kubisches Polynom von vier Koeffizienten beschrieben werden kann. Daraus folgt, dass die vier Ränder jeder Teilfläche solche Polynome sind. Um die geforderte zweifach stetige Differenzierbarkeit zu erhalten, können zunächst alle Punkte entlang der Gitterlinien eindimensional interpoliert werden. Das bikubische Spline erbt dabei die verwendeten Randbedingungen der eindimensionalen Splines.

Nun stehen zu jedem zweidimensionalen Polynom ${\displaystyle A_{i,j}}$ mit 16 Koeffizienten vier eindimensionale Randpolynome ${\displaystyle s_{x~i,j}}$, ${\displaystyle s_{x~i,j+1}}$, ${\displaystyle s_{y~i,j}}$, ${\displaystyle s_{y~i+1,j}}$ mit je 4 Koeffizienten zur Verfügung. Hierbei kennzeichnet der erste Index der s, zu welcher Achse sie parallel verlaufen.

Um diese Randpolynome zu ${\displaystyle A_{i,j}}$ „zusammenzurechnen“ ist ein System von 16 linearen Gleichungen aufzustellen. Vier Gleichungen ergeben sich aus der Forderung, dass S an den Gitterpunkten genau die Werte ${\displaystyle z_{i,j}}$ annimmt:

${\displaystyle {\begin{matrix}A_{i,j}(x_{i},y_{j})=z_{i,j},&A_{i,j}(x_{i+1},y_{j})=z_{i+1,j},&A_{i,j}(x_{i},y_{j+1})=z_{i,j+1},&A_{i,j}(x_{i+1},y_{j+1})=z_{i+1,j+1}\end{matrix}}}$

Weitere vier Gleichungen ergeben sich aus der Ableitung der zur ${\displaystyle x}$-Achse parallelen Randpolynome nach ${\displaystyle x}$:

${\displaystyle {\begin{matrix}{\frac {\partial A_{i,j}(x_{i},y_{j})}{\partial x}}={\frac {\partial s_{x~i,j}(x_{i})}{\partial x}},&{\frac {\partial A_{i,j}(x_{i+1},y_{j})}{\partial x}}={\frac {\partial s_{x~i,j}(x_{i+1})}{\partial x}},&{\frac {\partial A_{i,j}(x_{i},y_{j+1})}{\partial x}}={\frac {\partial s_{x~i,j+1}(x_{i})}{\partial x}},&{\frac {\partial A_{i,j}(x_{i+1},y_{j+1})}{\partial x}}={\frac {\partial s_{x~i,j+1}(x_{i+1})}{\partial x}}\end{matrix}}}$

Weitere vier Gleichungen aus der Ableitung der zur ${\displaystyle y}$-Achse parallelen Randpolynome nach ${\displaystyle y}$:

${\displaystyle {\begin{matrix}{\frac {\partial A_{i,j}(x_{i},y_{j})}{\partial y}}={\frac {\partial s_{y~i,j}(y_{j})}{\partial y}},&{\frac {\partial A_{i,j}(x_{i+1},y_{j})}{\partial y}}={\frac {\partial s_{y~i+1,j}(y_{j})}{\partial y}},&{\frac {\partial A_{i,j}(x_{i},y_{j+1})}{\partial y}}={\frac {\partial s_{y~i,j}(y_{j+1})}{\partial y}},&{\frac {\partial A_{i,j}(x_{i+1},y_{j+1})}{\partial y}}={\frac {\partial s_{y~i+1,j}(y_{j+1})}{\partial y}}\end{matrix}}}$

Nun ergibt sich das Problem, dass entlang jedem der vier Ränder zwei Punkte und zwei Ableitungen gegeben sind. Damit ist ein Polynom von Grad 3 entlang des jeweiligen Randes eigentlich vollständig spezifiziert. Das hinzunehmen z. B. der zweiten Ableitungen an den Eckpunkten oder weiterer Funktionswerte entlang der Randpolynome würde eine lineare Abhängigkeit im Gleichungssystem erzeugen. Mit den gegebenen Gleichungen lassen sich jedoch nur 12 der 16 Koeffizienten bestimmen. Ein nicht linear abhängiges System ergibt sich durch Hinzunahme der gemischten Ableitungen nach x und y. Setzt man diese in den Eckpunkten etwa auf null, so ist ${\displaystyle S}$ nur in den Eckpunkten zweimal stetig differenzierbar. Entlang der Ränder ergeben sich in den zweiten Ableitungen jedoch Sprünge. Aus den Randpolynomen lassen sich die gemischten Ableitungen jedoch nicht direkt berechnen.

Um nun korrekte Werte für diese gemischten Ableitungen zu erhalten, welche auch entlang der Ränder stetige, zweite Ableitungen von ${\displaystyle S}$ liefern, kann wie folgt vorgegangen werden:

• Entlang der Gitterlinien, welche parallel zur ${\displaystyle x}$-Achse verlaufen, werden weitere, eindimensionale Splines ${\displaystyle s_{d}}$ gebildet
• Diese interpolieren statt der ${\displaystyle z}$-Werte die Ableitungen von ${\displaystyle s_{y}}$ an deren Schnittpunkten mit den ${\displaystyle x}$-Gitterlinien
• Die ${\displaystyle s_{d}}$-Splines können nun nach ${\displaystyle x}$ abgeleitet werden. Daraus ergeben sich die gemischten Ableitungen

${\displaystyle {\begin{matrix}{\frac {\partial ^{2}A_{i,j}(x_{i},y_{j})}{\partial x\partial y}}={\frac {\partial s_{d~i,j}(x_{i})}{\partial x}},&{\frac {\partial ^{2}A_{i,j}(x_{i+1},y_{j})}{\partial x\partial y}}={\frac {\partial s_{d~i,j}(x_{i+1})}{\partial x}},&{\frac {\partial ^{2}A_{i,j}(x_{i},y_{j+1})}{\partial x\partial y}}={\frac {\partial s_{d~i,j+1}(x_{i})}{\partial x}},&{\frac {\partial ^{2}A_{i,j}(x_{i+1},y_{j+1})}{\partial x\partial y}}={\frac {\partial s_{d~i,j+1}(x_{i+1})}{\partial x}}\end{matrix}}}$

Wie beim eindimensionalen C²-Spline wirkt sich auch beim bikubischen Spline eine Änderung in einer Stützstelle (Datenpunkt) generell auf alle Teilflächen aus. Dies veranschaulicht, warum eine Konstruktion alleine aus den Randpolynomen ${\displaystyle s_{x}}$ und ${\displaystyle s_{y}}$ fehlschlägt: Diese ändern sich nur, wenn ein Datenwert, welcher parallel (in ${\displaystyle x}$- oder ${\displaystyle y}$-Richtung) zur betrachteten Teilfläche liegt, geändert wird. In der Abbildung rechts hat eine Änderung von ${\displaystyle z_{3,0}}$ keinen Einfluss auf die 1D-Splines, welche den Rand von ${\displaystyle A_{1,1}}$ bilden. Erst die erneute Interpolation durch ${\displaystyle s_{d}}$ erzeugt eine Abhängigkeit zwischen diesem Punkt und der Fläche.

Beispiel

Beispiel für bikubische Interpolation:[6] Ein Datenblock aus 6x6 Werten (links) wird bikubisch interpoliert (mitte). Dabei wurden natürliche Randbedingungen angenommen (die zweite Ableitung auf den Randpunkten ist null). Die Farben zwischen linkem und mittlerem Bild wurden synchronisiert und beschreiben die Funktionswerte ähnlich der Farben auf einer Landkarte zur Illustration der Höhe. Um die Korrektheit der Interpolation zu verdeutlichen, wird rechts die zweite Ableitung ${\displaystyle {\frac {\partial ^{4}S}{\partial x^{2}\partial y^{2}}}}$ gezeigt. Diese besteht aus linearen Funktionen und ist immer noch stetig.

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 ${\displaystyle f\colon [a,b]\to \mathbb {R} }$ erbt:

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

Hier z​eigt sich, d​ass klassische Splines e​twas schlechtere Eigenschaften h​aben als Bézierkurven. Zunächst stellt s​ich die Frage, w​ann ein interpolierender Spline konvex ist.

Für klassische Splines gilt, dass die Menge möglicher Splines auf dem Intervall ${\displaystyle [a,b]}$ zum Gitter ${\displaystyle \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 ${\displaystyle a=t_{0}<t_{1}<\dots <t_{m}=b}$ und zugehörige Ordinaten ${\displaystyle f_{0},\dots ,f_{m}\in \mathbb {R} }$ vorgegeben und gefordert, dass der Spline ${\displaystyle s}$ stetig differenzierbar in ${\displaystyle (a,b)}$ ist und darüber hinaus ${\displaystyle \displaystyle s(t_{k})=f_{k}}$ für ${\displaystyle 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 ${\displaystyle Y}$ aller Ordinatentupel ${\displaystyle (f_{0},\dots ,f_{m})}$, für die ein solcher Spline existiert, abgeschlossen ist.[7]

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

Siehe auch

• Akima-Interpolation – Ein Spezialfall der Spline-Interpolation

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

• Online-Tool zur Kubischen Spline-Interpolation mit Visualisierung und JavaScript-Quellcode
• Skript zur Spline-Interpolation, Uni Göttingen

Einzelnachweise

1. Robert Plato: Numerische Mathematik kompakt, 4. Auflage, Seite 27, Bemerkung 2.11
2. Josef Stoer: Numerische Mathematik 1, 9. Auflage, Abschnitt 2.4.2, Seite 109
3. Lebesgue und Birkhoff: Handbook of Splines, Abschnitt 1.9, Seite 56
4. Takeo Akasaki, William F. Donoghue Jr., George S. McCarty: John C. Holladay, 1928–1986. In Memoriam. University of California, 1986, abgerufen am 22. Januar 2021.
5. Lebesgue und Birkhoff: Handbook of Splines, Abschnitt 2.2, Seite 82
6. Die Beispielinterpolation wurde erstellt mit: Programm Numeric Collection von Mitja Stachowiak
7. Jochen W. Schmidt: Staircase Algorithm and Construction of Convex Spline Interpolants up to the Continuity ${\displaystyle C^{3}}$. In: Computers and Mathematics with Applications, Volume 31, Number 4, February 1995, pp. 67–79.