Spline-Interpolation

Grundidee und mathematische Aufgabe

Bei der Spline-Interpolation werden gegebene Stützstellen oder Knoten durch stückweise Polynome niedrigen Grades verbunden. Im Gegensatz zur Polynominterpolation entstehen auch bei ungünstig verteilten Stützstellen brauchbare, glatte Kurvenverläufe; starke Runge-Oszillationen werden vermieden. Die Berechnung ist mit geringem, linearem Aufwand möglich, allerdings ist die Konvergenzordnung im Vergleich zur Polynominterpolation geringer. Ohne weitere Zusätze bezeichnen Splineinterpolation und Splinefunktion gewöhnlich die kubische Splineinterpolation, also Splines dritten Grades. Smoothing Splines müssen nicht durch jeden Datenpunkt verlaufen und können zur Signalglättung verwendet werden.

Gegeben sind eine natürliche Zahl n ∈ ℕ, geordnete Stützstellen x₀ < x₁ < … < xₙ₋₁ < xₙ ∈ ℝ und Funktionswerte y₀, y₁, …, yₙ ∈ ℝ. Gesucht ist eine Funktion s: [x₀, xₙ] → ℝ mit s(xᵢ) = yᵢ für i = 0, …, n. Ihre Einschränkungen sᵢ := s|[xᵢ,xᵢ₊₁] auf den Teilintervallen [xᵢ, xᵢ₊₁] müssen Polynome sein.

Lineare Interpolation und ihre Genauigkeit

Die einfachste Splinefunktion ist ein Streckenzug: Zwischen jeweils zwei benachbarten Punkten liegt eine Gerade. Für die Punkte (x₁,y₁) und (x₂,y₂) gilt

s(x) = m·x + b = ((y₂ − y₁)/(x₂ − x₁))·x + y₁ − ((y₂ − y₁)/(x₂ − x₁))·x₁.

Äquivalent kann man schreiben:

s(x) = ((x₂ − x)/(x₂ − x₁))·y₁ + ((x − x₁)/(x₂ − x₁))·y₂.

Für n+1 Stützstellen a = x₀ < x₁ < … < xₙ = b seien hᵢ = xᵢ₊₁ − xᵢ und h := maxᵢ hᵢ. Interpoliert der lineare Spline eine Funktion f ∈ C²([a,b]), dann gilt mit der Maximumsnorm

∥f − s∥L∞([a,b]) := maxₓ∈[a,b] |f(x) − s(x)| ≤ (h²/8)∥f″∥L∞([a,b]).

Lineare Splines konvergieren somit quadratisch mit der Gitterweite h gegen f. Sie sind jedoch an den Stützstellen im Allgemeinen nicht differenzierbar. Kubische Splines beseitigen dieses Problem, weil sie zweimal stetig differenzierbar konstruiert werden können.

Kubische C²-Splines und ihre Konstruktion

Ein kubisches Teilpolynom besitzt vier Koeffizienten. Zusätzlich zu den Interpolationsbedingungen müssen daher Übergangsbedingungen festgelegt werden. Ein C²-Spline S ist so zusammengesetzt, dass S selbst sowie die erste und zweite Ableitung an jedem inneren Knoten stetig sind:

sᵢ′(xᵢ) = sᵢ₋₁′(xᵢ) und sᵢ″(xᵢ) = sᵢ₋₁″(xᵢ) für i = 1, …, n−1.

Eine Änderung eines Knotens wirkt sich grundsätzlich global auf den gesamten Spline aus. Der Einfluss nimmt mit wachsender Entfernung jedoch stark ab; kubische Splines neigen deshalb weniger zum Überschwingen als Interpolationspolynome.

Die zweite Ableitung S″ ist selbst ein linearer Spline. Mit hᵢ = xᵢ₊₁ − xᵢ und den Momenten Mᵢ = S″(xᵢ) gilt auf [xᵢ,xᵢ₊₁]:

sᵢ″(x) = ((xᵢ₊₁ − x)/hᵢ)·Mᵢ + ((x − xᵢ)/hᵢ)·Mᵢ₊₁.

Zweifache Integration liefert

sᵢ(x) = 1/6·[((xᵢ₊₁ − x)³/hᵢ)·Mᵢ + ((x − xᵢ)³/hᵢ)·Mᵢ₊₁] + cᵢ·(x − xᵢ) + dᵢ,

wobei die Interpolation an den Endpunkten durch

dᵢ = yᵢ − (hᵢ²/6)·Mᵢ und cᵢ = (yᵢ₊₁ − yᵢ)/hᵢ − (hᵢ/6)·(Mᵢ₊₁ − Mᵢ)

gesichert wird. Die Momente werden so bestimmt, dass auch die ersten Ableitungen übereinstimmen. Es entsteht für i = 1, …, n−1 das tridiagonale Gleichungssystem

(hᵢ₋₁/6)·Mᵢ₋₁ + ((hᵢ₋₁+hᵢ)/3)·Mᵢ + (hᵢ/6)·Mᵢ₊₁ = (yᵢ₊₁−yᵢ)/hᵢ − (yᵢ−yᵢ₋₁)/hᵢ₋₁.

Für i = 0 und i = n kommen zwei Gleichungen aus den Randbedingungen hinzu. Die Matrix ist tridiagonal und streng diagonaldominant. Zur Lösung genügt beispielsweise der Thomas-Algorithmus: ein Vorwärtsdurchlauf zur Elimination unterhalb der Hauptdiagonalen mit anschließender Rückwärtssubstitution.

Randbedingungen, Minimalität und Konvergenz

Da es ein Interpolationsintervall weniger als Stützstellen gibt, fehlen zwei Gleichungen. Typische Randbedingungen sind:

• Natürlich oder freier Rand: s₀″(x₀) = 0 und sₙ₋₁″(xₙ) = 0. Der Spline schließt mit Wendepunkten ab. In der Matrix setzt man λ₀ = λₙ = b₀ = bₙ = 0 sowie μ₀ = μₙ = 1.
• Hermite- oder eingespannter Rand: s₀′(x₀) = f′(a) und sₙ₋₁′(xₙ) = f′(b). Die Randableitungen stammen aus f oder einer Approximation. Es gilt λ₀ = h₀/6, μ₀ = h₀/3, b₀ = (y₁−y₀)/h₀ − f′(a) sowie λₙ = hₙ₋₁/6, μₙ = hₙ₋₁/3, bₙ = −(yₙ−yₙ₋₁)/hₙ₋₁ + f′(b).
• Periodisch: Auf [x₀,xₙ₊₁] gelten y₀ =: yₙ₊₁, s₀′(x₀) =: sₙ′(xₙ₊₁) und s₀″(x₀) =: sₙ″(xₙ₊₁). Anfang und Ende stimmen somit in Funktion, erster und zweiter Ableitung überein. Mₙ₊₁ := M₀ ist bereits gegeben; die Matrix erhält zusätzlich die Einträge m₀,ₙ = mₙ,₀ = hₙ/6.
• Not-a-knot: s₀‴(x₁) = s₁‴(x₁) und sₙ₋₂‴(xₙ₋₁) = sₙ₋₁‴(xₙ₋₁). Die äußeren drei Punkte werden jeweils durch ein gemeinsames Polynom interpoliert. Für höchstens vier Stützstellen wird daraus ein gewöhnliches Interpolationspolynom. Bei h₁ = h₀ tritt in einer angegebenen Matrixform eine Division durch null auf; dann ist eine Grenzwertbildung erforderlich. Diese Randbedingung wird beispielsweise von Matlab verwendet.

Unter natürlichen, periodischen oder Hermite-Randbedingungen besitzt der kubische Spline unter allen zweimal stetig differenzierbaren interpolierenden Funktionen die kleinste Krümmungsenergie:

∫ₐᵇ f″(x)² dx ≥ ∫ₐᵇ S″(x)² dx.

Die von Holladay 1957 bewiesene Identität lautet mit der L²-Norm ∥·∥:

∥f″−S″∥² = ∥f″∥² − ∥S″∥² − 2D,

wobei D := [(f′(x)−S′(x))S″(x)]ₐᵇ − Σᵢ₌₁ⁿ [(f(x)−S(x))S‴]ₓᵢ₋₁ˣᵢ. Für diese drei Randbedingungstypen ist D = 0. Daraus folgen ∥f″∥² ≥ ∥S″∥², die Minimalität und die Eindeutigkeit der Splineinterpolation.

Für h = maxᵢ hᵢ → 0 konvergieren kubische Splines und mindestens ihre ersten beiden Ableitungen. Die Ordnung hängt von Norm und Differenzierbarkeit von f ab. Unter anderem gelten:

∥f−s∥L∞ ≤ √2·h³ᐟ²∥f″∥L², ∥f−s∥L² ≤ 2h²∥f″∥L²,

∥f−s∥L∞ ≤ √8·h⁷ᐟ²∥f⁗∥L², ∥f−s∥L² ≤ 4h⁴∥f⁗∥L²,

∥f−s∥L∞ ≤ h⁴∥f⁗∥L∞,

∥f′−s′∥L∞ ≤ √2·h¹ᐟ²∥f″∥L², ∥f′−s′∥L² ≤ 2h²∥f″∥L²,

∥f′−s′∥L∞ ≤ √8·h⁵ᐟ²∥f⁗∥L², ∥f′−s′∥L² ≤ 4h³∥f⁗∥L²,

∥f″−s″∥L∞ ≤ (h²/2)∥f⁗∥L∞.

Dabei sind ∥f∥L∞([a,b]) = maxₓ∈[a,b]|f(x)| und ∥f∥L²([a,b])² = ∫ₐᵇ f(x)² dx. Für f ∈ C³([a,b]) gibt es analoge Formeln.

Bikubische Interpolation

Der bikubische C²-Spline verallgemeinert den eindimensionalen kubischen C²-Spline auf zwei Dimensionen. Die Datenpunkte zᵢⱼ müssen in einem rechteckigen Gitter liegen. Jede Teilfläche Aᵢⱼ zwischen vier Gitterpunkten wird durch ein Polynom mit 16 Koeffizienten beschrieben:

Aᵢⱼ(x,y) = Σₖ₌₀³Σₗ₌₀³ aₖ,ₗⁱʲ·xᵏ·yˡ.

Die zusammengesetzte Funktion S(x,y) soll zweimal stetig in x- und y-Richtung differenzierbar sein. Stetig sein müssen neben S die Ableitungen ∂S/∂x, ∂S/∂y, ∂²S/∂x², ∂²S/∂y², ∂²S/∂x∂y, ∂³S/∂x²∂y, ∂³S/∂x∂y² und ∂⁴S/∂x²∂y².

Jeder Schnitt einer Teilfläche parallel zu einer Koordinatenachse ist eine eindimensionale kubische Kurve. Deshalb werden zunächst die Gitterlinien eindimensional interpoliert; dabei übernimmt der bikubische Spline die gewählten Randbedingungen. Für jede Teilfläche stehen vier Randpolynome sₓ,ᵢⱼ, sₓ,ᵢⱼ₊₁, sᵧ,ᵢⱼ und sᵧ,ᵢ₊₁,ⱼ mit je vier Koeffizienten zur Verfügung.

Aus den vier Eckwerten, den vier x-Ableitungen und den vier y-Ableitungen entstehen 12 lineare Gleichungen. Die 16 Koeffizienten sind damit noch nicht vollständig bestimmt. Zusätzliche Bedingungen für gemischte Ableitungen ∂²S/∂x∂y sind nötig. Ein bloßes Setzen dieser Werte auf null gewährleistet die zweimalige Differenzierbarkeit nur in den Eckpunkten; entlang der Ränder können Sprünge der zweiten Ableitungen auftreten. Daher werden entlang der x-parallelen Gitterlinien weitere eindimensionale Splines s_d gebildet, die statt der z-Werte die Ableitungen der y-Randpolynome an den Schnittpunkten interpolieren. Durch Ableiten der s_d nach x erhält man die korrekten gemischten Ableitungen an den Ecken. So entsteht die erforderliche Abhängigkeit aller Teilflächen von den Datenpunkten.

Als Beispiel wird ein Datenblock mit 6×6 Werten bikubisch interpoliert. Dabei werden natürliche Randbedingungen angenommen, also eine zweite Ableitung von null an den Randpunkten. Zur Kontrolle wird ∂⁴S/∂x²∂y² dargestellt; diese Ableitung besteht aus linearen Funktionen und ist weiterhin stetig.

Höhere Ordnung und Formerhaltung

Splines lassen sich mit stückweisen Polynomen vom Grad p > 3 erweitern. Um symmetrische Randbedingungen zu ermöglichen, verwendet man ungerade Grade, insbesondere lineare Splines p = 1, kubische Splines p = 3 und quintische Splines p = 5. Aussagen zu Randbedingungen, Eindeutigkeit und Konvergenz sind im Wesentlichen analog. Für p = 5 entsteht weiterhin ein dünn besetztes lineares Gleichungssystem, dessen Matrix jedoch nicht mehr tridiagonal ist. Bei zu großen Graden können wie bei der Polynominterpolation Runge-Oszillationen auftreten; außerdem sind geeignete praktische Randbedingungen schwieriger zu definieren. Splines mit p > 5 sind daher selten.

Wegen ihrer glatten Kurvenverläufe werden Splines häufig im CAD eingesetzt. Wichtige formerhaltende Eigenschaften einer Funktion f: [a,b] → ℝ sind:

• Nichtnegativität: f(x) ≥ 0 für alle x.
• Monotonie: f(x) ≤ f(y) für a ≤ x ≤ y ≤ b.
• Konvexität: f(λx + (1−λ)y) ≤ λf(x) + (1−λ)f(y) für x,y ∈ [a,b] und λ ∈ [0,1].

Klassische Splines besitzen hierbei teilweise schlechtere Eigenschaften als Bézierkurven. Für ein Gitter Δ: a = x₀ < x₁ < … < xₙ = b bildet die Menge der klassischen Splines einen endlichdimensionalen Vektorraum. Werden Knoten a = t₀ < t₁ < … < tₘ = b und Ordinaten f₀, …, fₘ vorgegeben, so sei Y die Menge aller Datentupel, für die ein stetig differenzierbarer, zusätzlich konvexer interpolierender Spline existiert. Unter den genannten technischen Voraussetzungen ist Y abgeschlossen, aber für m ≥ 2 eine echte Teilmenge von ℝᵐ⁺¹. Rechenungenauigkeiten können deshalb Daten auf den Rand von Y aus dem lösbaren Bereich hinausführen.

Ein Beispiel zeigt eine grundsätzliche Einschränkung: Liegen fünf Punkte in Form des Zeichens „∨“, wobei der mittlere Punkt genau auf der Spitze liegt, ist die einzige konvexe Interpolierende die Betragsfunktion; diese ist nicht stetig differenzierbar. Das 5-Tupel liegt daher außerhalb von Y, und auch geringfügig nach oben verschobene Punkte können trotz streng konvexer Lage unlösbar bleiben. Als mögliche Auswege nennt der Artikel gebrochen-rationale Splines, Splines mit frei wählbaren Zwischenknoten, Exponentialsplines und lakunäre beziehungsweise lückenhafte Splines.