Numerische Mathematik 1 – Approximation

Interpolation mit Polynomen

Lagrange-Form und 4-Punkt-Formel

Gegeben seien \(n + 1\) paarweise verschiedene Stützstellen \(x_0, \dotsc , x_n\) mit Funktionswerten
\(f_0, \dotsc , f_n\). Dann können diese eindeutig durch ein Polynom \(p\) mit Grad \(\le n\) interpoliert werden (polynomiale Interpolation):

\begin{align*} p(x_k) = f_k \quad \text { fÃijr }\quad k = 0, \dotsc , n. \end{align*} Das Interpolationspolynom \(p\) lässt sich in der Lagrange-Form darstellen:

\begin{align*} p(x) = \sum _{k=0}^n f_k q_k(x), \quad q_k(x) = \prod _{j \not = k} \frac {x - x_j}{x_k - x_j}. \end{align*} Die \(q_k\) heißen dabei Lagrange-Polynome. Sie haben in \(x_k\) den Wert \(1\) und verschwinden in allen anderen Punkten \(x_j\), d. h. \(q_k(x_j) = \delta _{kj}\).

Generation von Zwischenwerten mittels kubischer Interpolation (4-Punkt-Formel):
Sind die \(f_k\) an äquidistanten Stützstellen gegeben, d. h. \(x_k = kh\) mit \(h\) Gitterweite, \(k = 0, \dotsc , n\), so kann man Zwischenwerte an den Stützstellen \(x_{k + 1/2} = (k + 1/2)h\) durch kubische Interpolation approximieren. Zur Interpolation verwendet man dazu die vier benachbarten Stützstellen \(x_{k-1}, x_k, x_{k+1}, x_{k+2}\). Dabei ergibt sich die Formel

\begin{align*} f_{k + 1/2} = (-f_{k-1} + 9f_k + 9f_{k+1} - f_{k+2}) / 16. \end{align*} Den Prozess kann man solange wiederholen, bis genügend Daten erzeugt wurden. Die Gewichte \(-\frac {1}{16}, \frac {9}{16}, \frac {9}{16}, -\frac {1}{16}\) sind die Werte der Lagrange-Polynome an der neuen Stützstelle \(x_{k + 1/2}\). Um \(x_{1/2}\) bzw. \(x_{n - 1/2}\) zu approximieren, verwendet man die Stützstellen \(x_0, x_1, x_2, x_3\) bzw.
\(x_{n-3}, x_{n-2}, x_{n-1}, x_n\) (für \(x_{1/2}\) ergibt sich z. B. \(f_{1/2} = (5f_0 + 15f_1- 5f_2 + f_3) / 16\)).

Schätzformel für zweite Ableitung:
Für die erste Ableitung \(f’\) einer Funktion \(f\) gilt \(f’(x) \approx \frac {f(x) - f(x - h)}{h}\). Daraus folgt für die zweite Ableitung durch Taylorentwicklung (\(f(x + h) = f(x) + f’(x)h + \frac {1}{2} f”(x) h^2 + o(h^2)\), \(h \to 0\))

\begin{align*} f”(x) \approx \frac {f(x + h) - f(x) - f’(x)h}{h^2} \approx \frac {f(x + h) - 2f(x) + f(x - h)}{h^2}. \end{align*}

Schema von Aitken-Neville

Gegeben seien wieder \(n + 1\) Datenpunkte \((x_i, f_i)\), \(i = 0, \dotsc , n\) mit \(x_0 < \dotsb < x_n\).
Dann lässt sich der Wert \(p(x)\) des Interpolationspolynoms an der Stelle \(x \in [x_0, x_n]\) mithilfe eines Dreiecksschemas (Schema von Aitken-Neville) berechnen:

\begin{align*} \begin{array}{ccccccc} f_0 = p_0^0 \\ & \searrow \\ f_1 = p_1^0 & \rightarrow & p_0^1 \\ \vdots & & & \ddots \\ f_{n-1} = p_{n-1}^0 & & \cdots & & p_0^{n-1} \\ & \searrow & & & & \searrow \\ f_n = p_n^0 & \rightarrow & p_{n-1}^1 & \cdots & p_1^{n-1} & \rightarrow & p_0^n = p(x) \end {array}\qquad \begin{array}{l} \text {mit} \\ p_i^j := \dfrac {x_{i+j} - x}{x_{i+j} - x_i} p_i^{j-1} + \dfrac {x - x_i}{x_{i+j} - x_i} p_{i+1}^{j-1}, \\ p_i^0 := f_i, \\[3mm] j = 1, \dotsc , n,\quad i = 0, \dotsc , n - j. \end {array} \end{align*}

Dabei ist \(p_i^j = p_i^j(x)\) ein Polynom vom Grad \(\le j\), das an den Punkten \(x_i, \dotsc , x_{i+j}\) interpoliert. Der Vorteil dieses Dreiecksschemas ist, dass zur Verbesserung der Genauigkeit weitere Datenpunkte sehr einfach als neue Zeile am unteren Rand hinzugefügt werden können, ohne alle Werte neu zu berechnen (anders als z. B. mit Lagrange-Polynomen).

Das Aitken-Neville-Schema kann auch dazu benutzt werden, algorithmisch die Koeffizienten des Interpolationspolynoms zu berechnen. Dazu berechnet man das Schema spaltenweise und speichert in einer Matrix die Koeffizienten der Polynome der vorherigen Spalte. Die neuen Koeffizienten können dann mithilfe der Definition von \(p_i^j\) berechnet werden.

Polynome in Newton-Form, Horner-Schema

Seien ein Polynom \(p\) vom Grad \(\le n\) und \(n\) Punkte \(x_0, \dotsc , x_{n-1}\) gegeben.
Dann ist die Newton-Form von \(p\) bezüglich der Punktfolge \(x_0, \dotsc , x_{n-1}\)

\begin{align*} p(x) = a_0 + a_1 (x - x_0) + \dotsb + a_n (x - x_0) \dotsm (x - x_{n-1}). \end{align*} Dies kann als verallgemeinerte Taylor-Darstellung aufgefasst werden, denn für \(x_0 = \dotsb = x_{n-1}\) erhält man insbesondere die Taylor-Entwicklung von \(p\) im Punkt \(x_0\).

Die Auswertung eines Polynoms in Newton-Form kann mittels Horner-Schema erfolgen:

\begin{align*} p(x) = (\dotsb (a_n y_{n-1} + a_{n-1}) y_{n-2} + \dotsb ) y_0 + a_0 \end{align*} mit \(y_k = x - x_k\), \(k = 0, \dotsc , n - 1\).
Durch die geschachtelte Multiplikation benötigt man nur \(3n\) Operationen:
Anfangs setzt man \(p := a_n\) und dann berechnet man \(p \leftarrow p y_k + a_k\) für \(k = n - 1, \dotsc , 0\).

Umwandlung eines Polynoms von Normalform \(\sum _{k=0}^n c_k x^k\) in Newton-Form:
Der Koeffizient \(a_n\) der Newton-Form ist der Koeffizient von \(x^n\), denn \(p(x) = a_n x^n + \O (x^{n-1})\). So kann man die Newton-Form rekursiv berechnen: \(a_n\) ist der höchste Koeffizient von \(p(x)\), dann subtrahiert man den letzten Summanden \(q(x) = p(x) - a_n (x - x_0) \dotsm (x - x_{n-1})\). \(a_{n-1}\) ist wiederum der höchste Koeffizient dieses Restterms \(q(x)\) usw.

Hermite-Interpolation

Seien eine glatte Funktion \(f\) und \(n + 1\) Punkte \(x_0, \dotsc , x_n\) (nicht notwendig verschieden) gegeben.
Dann gibt es genau ein Polynom \(p\) vom Grad \(\le n\) mit

\begin{align*} p^{(j)}(x_k) = f^{(j)}(x_k), \quad 0 \le j < m_k, \end{align*} wobei \(m_k\) die Vielfachheit des Punktes \(x_k\) ist (\(k = 0, \dotsc , n\)). Tritt also ein Punkt mehrfach auf, so werden nicht nur Funktionswerte, sondern auch Ableitungen interpoliert. Dieses Interpolationsverfahren heißt Hermite-Interpolation.

Im Schaubild werden Vielfachheiten durch eng nebeneinander liegende Markierungen auf der \(x\)-Achse oder zusätzliche Kreise um die Interpolationspunkte angedeutet.

Sind Daten in der Form \((x_k, f_k)\) gegeben, so verwendet man meistens die Konvention, dass \(f_k = p^{(j)}(x_k)\), wobei \(j\) die Anzahl der Punkte \(x_i\) mit \(x_i = x_k\) und \(i < k\) ist.

Im Beispiel \((1, 3), (2, 1), (2, 0), (2, 2), (4, 2), (4, 1)\) interpoliert das Polynom der Hermite-Interpolation \(f(1), f(2), f’(2), f”(2), f(4), f’(4)\).

Dividierte Differenzen

Die Dividierte Differenz ist eine Verallgemeinerung des Differenzenquotienten

\begin{align*} \Delta (a, b)f = \frac {f(a) - f(b)}{a - b}. \end{align*} Sie ist rekursiv definiert:

\begin{align*} \Delta (x_0, \dotsc , x_n)f := \frac {\Delta (x_1, \dotsc , x_n)f - \Delta (x_0, \dotsc , x_{n-1})f}{x_n - x_0} \quad \text { fÃijr } x_0 \not = x_n \quad \text { sowie } \end{align*}

\begin{align*} \Delta (\underbrace {x, \dotsc , x}_{n + 1 \text {-mal}})f := \frac {1}{n!} f^{(n)}(x). \end{align*} Insbesondere ist \(\Delta (x)f = f(x)\).

Die Dividierten Differenzen

\begin{align*} \Delta _i^j := \Delta (x_i, \dotsc , x_{i+j}) = \frac {\Delta _{i+1}^{j-1} - \Delta _i^{j-1}}{x_{i+j} - x_i} \end{align*} können in einem Dreiecksschema berechnet werden:

\begin{align*} \begin{array}{c||cccc} x_0 & \Delta _0^0 & \Delta _0^1 & \Delta _0^2 & \Delta _0^3 \\ x_1 & \Delta _1^0 & \Delta _1^1 & \Delta _1^2 \\ x_2 & \Delta _2^0 & \Delta _2^1 \\ x_3 & \Delta _3^0 \end {array} \end{align*} Dabei hängt der Eintrag \(\Delta _i^j\) vom Eintrag links (\(\Delta _i^{j-1}\)) und links unten (\(\Delta _{i+1}^{j-1}\)) ab.
Die Daten müssen vorgegeben sein, falls es keine zwei verschiedene Stellen in \(\Delta _i^j\) gibt. Startwerte sind dabei Funktionswerte oder bei Vielfachheiten Ableitungswerte von \(f\) an den Punkten \(x_k\). Dabei schreibt man in die \(j\)-te Spalte \(\frac {1}{j!} f^{(j)}(x_k)\) (nach Definition der Dividierten Differenz).
Ansonsten, falls es zwei verschiedene Stellen in \(\Delta _i^j\) gibt, kann man den Eintrag mittels der Definition der Dividierten Differenz berechnen.

Integraldarstellung Dividierter Differenzen

Die Dividierte Differenz \(\Delta (x_0, \dotsc , x_n)f\) lässt sich als Integral über den von den Einheitsvektoren aufgespannten Simplex darstellen (Formel von Hermite-Genocchi):

\begin{align*} \Delta (x_0, \dotsc , x_n)f = \int _{s_0 + \dotsb + s_n = 1} f^{(n)}(s_0 x_0 + \dotsb + s_n x_n)\ds . \end{align*} Daraus folgt insbesondere, dass Dividierte Differenzen für glatte Funktionen \(f\) stetig von den Punkten \(x_k\) abhängen und

\begin{align*} \Delta (x_0, \dotsc , x_n)f = \frac {f^{(n)}(\xi )}{n!} \quad \text {mit}\quad \xi \in [\min x_k, \max x_k]. \end{align*}

Ein \(n\)-dimensionaler Simplex ist dabei die konvexe Hülle \(S\) von \(n + 1\) Punkten \(p_0, \dotsc , p_n\), die nicht alle in einem \(n - 1\)-dimensionalen Unterraum liegen:

\begin{align*} S = \left \{\left .\sum _{j=0}^n \alpha _j p_j \;\right |\; \sum _{j=0}^n \alpha _j = 1,\; \alpha _j \ge 0\right \}. \end{align*} Das Volumen eines Simplex lässt sich durch die Vektoren \(p_i - p_0\), \(i = 1, \dotsc , n\) ausdrücken:

\begin{align*} \vol S = \frac {1}{n!} \cdot |\det (p_1 - p_0, \dotsc , p_n - p_0)|. \end{align*}

Newton-Form und Dividerte Differenzen

Die Newton-Form des Polynoms \(p\) vom Grad \(\le n\), das eine Funktion \(f\) an den Punkten \(x_0, \dotsc , x_n\) interpoliert, lässt sich mithilfe Dividerter Differenzen angeben:

\begin{align*} p(x) = a_0 + a_1 (x - x_0) + \dotsb + a_n (x - x_0) \dotsm (x - x_{n-1}), \quad \; a_k = \Delta (x_0, \dotsc , x_k)f,\;\; k = 0, \dotsc , n. \end{align*} Dabei werden an einem Punkt mit Vielfachheit \(m\) zusätzlich auch alle Ableitungen der Ordnung \(< m\) interpoliert.

Die Newton-Form ist insbesondere geeignet, wenn man weitere Interpolationspunkte hinzufügen will. Die Darstellung braucht dann jeweils um nur einen weiteren Term ergänzt zu werden, die vorherigen Terme bleiben inklusive Koeffizienten gleich. Das Schema zur Berechnung des neuen höchsten Koeffizienten als Dividerte Differenz wird um eine neue Diagonale ergänzt.
Außerdem können Ableitungen mit der Newton-Form einfach interpoliert werden.

Fehler bei der Interpolation glatter Funktionen

Der Fehler des Polynoms \(p\) vom Grad \(\le n\), das eine glatte Funktion an den Punkten \(x_0, \dotsc , x_n\) interpoliert, lässt sich darstellen in der Form

\begin{align*} f(x) - p(x) = \frac {f^{(n+1)}(\xi )}{(n + 1)!} (x - x_0) \dotsm (x - x_n) \quad \text {mit}\quad \xi \in \left [\min \{x, x_k\}, \max \{x, x_k\}\right ]. \end{align*} Insbesondere gilt für äquidistante Punkte \(x_k = x_0 + kh\), \(k = 0, \dotsc , n\)

\begin{align*} |f(x) - p(x)| = \O (h^{n+1}),\quad h \to 0, \quad x_0 \le x \le x_n. \end{align*} Für \(x_0 = \dotsb = x_n\) erhält man die Formel für das Taylor-Restglied.

Wie man in der Fehlerformel sieht, hängt der Fehler stark davon ab, welche Stützpunkte man für eine gegebene zu approximierende Funktion wählt. Setzt man die Punkte äquidistant, so weicht die Approximation zu den Rändern hin stark von der Funktion ab und „pendelt“ hin und her. Dies kann man vermeiden, indem man am Rand mehr Punkte (also dichter) wählt, da dort Informationen im Vergleich zur Mitte „fehlen“.

Mit den Tschebyscheff-Polynomen

\begin{align*} T_n(x) = \cos (n \arccos x),\quad x \in [-1, 1] \end{align*} lässt sich besser approxmieren, indem man die Nullstellen von \(T_n\)

\begin{align*} \xi _i^{(n)} = \cos \left (\frac {2i + 1}{2n} \pi \right ),\quad i = 0, \dotsc , n - 1 \end{align*} als Stützpunkte verwendet. Diese heißen Tschebyscheff-Knoten und basieren auf einer äquidistanten Winkelunterteilung im Halbkreis. Sie sind daher am Rand dichter verteilt wie in der Mitte. Verwendet man die Tschebyscheff-Knoten als Stützstellen bei der Interpolation, so verringert sich der Fehler an den Rändern deutlich.

Polynominterpolation mit MATLAB

Die Koeffizienten eines Polynoms \(p(x) = a_1 x^n + \dotsb + a_n x + a_{n+1}\) vom Grad \(\le n\), das die Daten \((x_k, y_k)\) interpoliert, können in M ATLAB mit a = polyfit(x, y, n); ermittelt werden. Wenn \(n\) kleiner ist als die Anzahl der Datenpunkte minus \(1\), so wird das Polynom, das die Fehlerquadratsumme minimiert, bestimmt. Mit p = polyval(a, x); kann das Polynom in den Punkten \(x_k\) ausgewertet werden (d. h. \(p_k = p(x_k)\)).

Orthogonale Polynome

Allgemeines

Zu jeder auf einem Intervall \((a, b)\) positiven Gewichtsfunktion \(w\) existiert eine bzgl. des Skalarprodukts \(\innerproduct {f, g} := \int _a^b f g w\) orthogonale Folge von Polynomen

\begin{align*} p_n(x) = \alpha _n x^n + \O (x^{n-1}),\quad \alpha _n \not = 0. \end{align*} Bis auf die Normierungsfaktoren \(\alpha _n\) sind die orthogonalen Polynome durch die Orthogonalitätsbedingungen \(\innerproduct {p_m, p_n} = 0\) für \(m \not = n\) eindeutig bestimmt und können mit dem Orthongalisierungsverfahren von Gram-Schmidt bestimmt werden.

Orthogonalisierungsverfahren von Gram-Schmidt: Sei \(b_1, \dotsc , b_n\) Basis eines Vektorraums \(V\). Dann kann man eine orthogonale Basis \(u_1, \dotsc , u_n\) durch

\begin{align*} u_j := b_j - \sum _{k=1}^{j-1} \frac {\innerproduct {b_j, u_k}}{\innerproduct {u_k, u_k}} u_k, \quad j = 1, \dotsc , n \end{align*} konstruieren.

Beim analogen Orthonormalisierungsverfahren von Gram-Schmidt vereinfacht sich die Rekursion, da man die Basisvektoren nach jedem Schritt normiert:

\begin{align*} u_j := b_j - \sum _{k=1}^{j-1} \innerproduct {b_j, u_k} u_k, \quad u_j \leftarrow \frac {u_j}{\norm {u_j}}, \quad j = 1, \dotsc , n. \end{align*}

Im Falle der orthogonalen Polynome geht man von der Monom-Basis \(\alpha _0 x^0, \dotsc , \alpha _n x^n\) aus.

Zur Bestimmung der orthogonalen Polynome für \((a, b) = (0, 1)\), \(w(x) \equiv 1\) kann man auch anders vorgehen: Die eine Methode berechnet iterativ das Polynom \(p_n\) aus den vorhergehenden, schon bekannten Polynomen \(p_k\), \(k = 0, \dotsc , n - 1\).
Dabei wird \(p_n(x) = \alpha _n x^n + c_{n-1} x^{n-1} + \dotsb + c_1 x + c_0\) allgemein angesetzt (\(\alpha _n\) bekannter Normierungsfaktor). Aus den Orthogonalitätsbedingungen
\(0 = \innerproduct {p_n, p_k} = \int _0^1 \left (\alpha _n x^n + \sum _{i=0}^{n-1} c_i x^i\right ) p_k(x) \dx \) für \(k = 0, \dotsc , n - 1\) folgen dann \(n\) Gleichungen für die \(n\) Unbekannten \(c_0, \dotsc , c_{n-1}\). Das Lösen des LGS liefert die gesuchten Koeffizienten.

Mit der anderen Art kann das Polynom \(p_n\) auch direkt bestimmt werden, ohne die vorherigen \(p_k\), \(k = 0, \dotsc , n - 1\) zu kennen. Da man weiß, dass \(p_n(x) = \alpha _n x^n + \sum _{k=0}^{n-1} c_k x^k\) zu allen Polynomen vom Grad \(< n\) und damit auch zu den Monomen \(x^j\), \(j = 0, \dotsc , n - 1\) orthogonal ist, setzt man

\begin{align*} \int _0^1 \left (\alpha _n x^n + \sum _{k=0}^{n-1} c_k x^k\right ) x^j \dx = 0 \quad \Leftrightarrow \quad \sum _{k=0}^{n-1} \frac {1}{j + k + 1} c_k = - \frac {\alpha _n}{n + j + 1} \end{align*} und erhält die Koeffizienten durch Lösen des LGS mit der sog. Hilbert-Matrix

\begin{align*} \left (\frac {1}{i + j - 1}\right )_{i,j=1}^n = \begin{pmatrix} \frac {1}{1} & \frac {1}{2} & \dots & \frac {1}{n} \\ \frac {1}{2} & \frac {1}{3} & \dots & \frac {1}{n + 1} \\ \vdots & \vdots & & \vdots \\ \frac {1}{n} & \frac {1}{n + 1} & \dots & \frac {1}{2n - 1} \end {pmatrix}. \end{align*}

Dreigliedrige Rekursion für orthogonale Polynome

Die orthogonalen Polynome \(q_n = x^n + \O (x^{n-1})\) mit Nomierungsfaktor \(1\) zu einer Gewichtsfunktion \(w\) auf einem Intervall \((a, b)\) können rekursiv berechnet werden:

\begin{align*} q_{n+1} = (\xi - \beta _n) q_n - \gamma _n q_{n-1},\quad n \ge 2 \end{align*} mit \(\xi (x) := x\) (dreigliedrige Rekursion). Die Koeffizienten \(\beta _n\) und \(\gamma _n\) lassen sich mithilfe des Skalarprodukts \(\innerproduct {f, g} = \int _a^b f g w\) ausdrücken:

\begin{align*} \beta _n := \frac {\innerproduct {\xi q_n, q_n}}{\varrho _n},\quad \gamma _n := \frac {\varrho _n}{\varrho _{n-1}} \quad \text {mit}\quad \varrho _n := \innerproduct {q_n, q_n}. \end{align*}

Für orthogonale Polynome \(p_n(x) = \alpha _n x^n + \O (x^{n-1})\) mit allgemeinem Normierungsfaktor \(\alpha _n\) gilt die Rekursion

\begin{align*} p_{n+1} = \frac {\alpha _{n+1}}{\alpha _n} (\xi - \beta ’_n) p_n - \frac {\alpha _{n-1} \alpha _{n+1}}{\alpha _n^2} \gamma ’_n p_{n-1}, \end{align*} wobei die Formeln für \(\beta ’_n, \gamma ’_n\) aus denen von \(\beta _n, \gamma _n\) entstehen, wenn man \(q_n\) durch \(p_n\) ersetzt.

Nullstellen orthogonaler Polynome

Das orthogonale Polynom \(p_n\) vom Grad \(n\) zu einer Gewichtsfunktion \(w\) auf \((a, b)\) hat \(n\) einfache Nullstellen in \((a, b)\). (Diese liegen zwischen denen von \(p_{n+1}\).)

Legendre-Polynome

Die Legendre-Polynome

\begin{align*} p_n(x) := \frac {1}{2^n n!} \frac {d^n}{dx^n} (x^2 - 1)^n = \frac {(2n)!}{2^n (n!)^2} x^n + \O (x^{n-1}) \end{align*} sind bzgl. des Skalarprodukts \(\innerproduct {f, g} = \int _{-1}^1 f g\) orthogonal.
Sie sind Lösungen der Differentialgleichung

\begin{align*} ((1 - \xi ^2) p’_n)’ = -n (n + 1) p_n \end{align*} mit \(\xi (x) = x\) und erfüllen die dreigliedrige Rekursion

\begin{align*} (n + 1) p_{n+1} = (2n + 1) \xi p_n - n p_{n-1}. \end{align*}

Tschebyscheff-Polynome

Die Tschbyscheff-Polynome entstehen durch Transformation der Kosinus-Funktionen:

\begin{align*} p_n(x) := \cos (nt), \quad x = \cos (t). \end{align*} Einem Argument \(x \in [-1, 1]\) entspricht der Winkel \(t = \arccos (x) \in [0,\pi ]\), der den Wert des Polynoms als \(\cos (nt)\) bestimmt. Das Polynom \(p_n\) hat in \([0,1]\) \(n\) Nullstellen \(\xi _k\) (\(k = 1, \dotsc , n\)) und \(n + 1\) Extrema \(\eta _k\) (\(k = 0, \dotsc , n\)), nämlich

\begin{align*} \xi _k := \cos \left (\frac {(2k - 1)\pi }{2n}\right ), \quad \eta _k := \cos \left (\frac {k\pi }{n}\right ), \text { wobei } p_n(\eta _k) = (-1)^k. \end{align*} Die Tschebyscheff-Polynome erfüllen die Orthogonalitätsrelation

\begin{align*} \int _{-1}^1 p_m(x) p_n(x) \frac {\dx }{\sqrt {1 - x^2}} = \begin{cases}\pi & m = n = 0 \\ \pi /2 & m = n > 0 \\ 0 & \text {sonst}, \end {cases} \end{align*} d. h. \([a,b] = [-1, 1]\) und \(w(x) = \frac {1}{\sqrt {1 - x^2}}\). Dies impliziert die dreigliedrige Rekursion

\begin{align*} p_{n+1} = 2 \xi p_n - p_{n-1}, \quad n \ge 1 \quad \text { mit } \quad \xi (x) = x. \end{align*}

Die Tschebyscheff-Entwicklung einer Funktion \(f\)

\begin{align*} f(x) \sim \sum _{n=0}^\infty \frac {\innerproduct {f, p_n}}{\varrho _n} p_n(x), \quad \varrho _n = \innerproduct {p_n, p_n} \end{align*} entspricht der Fourier-Reihe der transformierten Funktion \(g(t) = f(x)\), \(x = \cos (t)\), d. h.

\begin{align*} g(t) \sim \sum _{n=0}^\infty \frac {1}{\varrho _n} \left (\int _0^\pi f(\cos t) \cos (nt) \dt \right ) \cos (nt). \end{align*} Damit kann die schnelle Fourier-Transformation zur näherungsweisen Berechnung der
Entwicklungs-Koeffizienten herangezogen werden.

Minimalität der Tschebyscheff-Polynome

Das Tschebyscheff-Polynom \(p_n(x) = \cos (n \arccos (x))\) minimiert

\begin{align*} \max _{x \in [-1,1]} |q(x)| \end{align*} (eindeutig) unter allen Polynomen \(q\) vom Grad \(n\) mit gleichem höchsten Koeffizienten. Anders formuliert ist die Maxiumum-Norm des Produkts

\begin{align*} \prod _{k=1}^n (x - \xi _k) \end{align*} auf dem Intervall \([-1,1]\) für die Nullstellen \(\xi _k\) von \(p_n\) minimal.

Diskrete Fourier-Transformation

Einschub: Fourier-Reihen

Sei \(f\) eine \(2\pi \)-periodische Funktion, d. h. \(\forall _{x \in \real }\; f(x + 2\pi ) = f(x)\). Dann ist die komplexe Fourier-Reihe von \(f\) die Entwicklung nach dem Orthonormalsystem \(e_k(x) = e^{\i kx}\), \(k \in \integer \):

\begin{align*} f(x) \sim \sum _{k \in \integer } c_k e_k(x), \quad c_k = \innerproduct {f, e_k}_{2\pi } := \frac {1}{2\pi } \int _0^{2\pi } f(t) \overline {e_k(t)} \dt . \end{align*} Die Art der Konvergenz der Reihe hängt dabei von der Glattheit von \(f\) bzw. dem Abfallverhalten Fourier-Koeffizienten \(c_k\) ab. Hinreichend für gleichmäßige Konvergenz ist \(\sum _{k \in \integer } |c_k| < \infty \).

Ist \(f\) eine reellwertige \(2\pi \)-periodische Funktion, so ist die reelle Fourier-Reihe von \(f\) die Entwicklung nach dem Orthogonalsystem der Kosinus- und Sinusfunktionen:

\begin{align*} f(x) \sim \frac {a_0}{2} + \sum _{k=1}^\infty (a_k \cos (kx) + b_k \sin (kx)), \end{align*}

\begin{align*} a_k := \frac {1}{\pi } \int _{-\pi }^\pi f(t) \cos (kt) \dt , \quad b_k := \frac {1}{\pi } \int _{-\pi }^\pi f(t) \sin (kt) \dt . \end{align*} Wiederum hängt die Art der Konvergenz der Reihe von der Glattheit von \(f\) ab. Hinreichend für absolute Konvergenz ist, dass die Fourier-Koeffizienten \(a_k\) und \(b_k\) absolut konvergente Reihen bilden.

Auch eine konvergente Fourier-Reihe muss i. A. nicht an allen Stellen den Funktionswert als Grenzwert annehmen. An Unstetigkeitsstellen konvergiert die Reihe meist gegen den Mittelwert aus links- und rechtsseitigem Grenzwert.
Daher schreibt man oft \(f(x) \sim \sum \dotsb \) statt \(f(x) = \sum \dotsb \).

Komplexe Einheitswurzeln

Die Gleichung \(z^n = 1\), \(z \in \complex \) hat genau \(n\) Lösungen

\begin{align*} z_k = w_n^k, \quad w_n := \exp (2 \pi \i / n), \quad k = 0, \dotsc , n - 1, \end{align*} die als Einheitswurzeln bezeichnet werden, wobei \(w_n^n = 1\) und \(w_n^{k + nm} = w_n^k\) für \(m \in \integer \). Die Einheitswurzeln \(w_n^k\) bilden ein dem Einheitskreis einbeschriebenes regelmäßiges \(n\)-Eck.

Fourier-Matrix

Durch Bilden von Potenzen der Einheitswurzel \(w_n = \exp (2 \pi \i / n)\) erhält man die Fourier-Matrix

\begin{align*} W_n := \begin{pmatrix} w_n^{0 \cdot 0} & \dots & w_n^{0 \cdot (n - 1)} \\ \vdots & & \vdots \\ w_n^{(n - 1) \cdot 0} & \dots & w_n^{(n - 1) \cdot (n - 1)} \end {pmatrix} = (w_n^{k\ell })_{k,\ell =0}^{n-1}. \end{align*} Dabei ist \(W_n/\sqrt {n}\) unitär, d. h. \(W_n^\ast W_n/n\) ist die Einheitsmatrix.

Diskrete Fourier-Transformation

Die Multiplikation eines Vektors \(c = (c_0, \dotsc , c_{n-1})^t\) mit der Fourier-Matrix \(W_n\) wird als
diskrete Fourier-Transformation bezeichnet:

\begin{align*} f = W_n c \quad \Leftrightarrow \quad c = \frac {1}{n} W_n^\ast f. \end{align*} Komponentenweise gilt also

\begin{align*} f_j = \sum _{k=0}^{n-1} c_k w_n^{jk} \quad \Leftrightarrow \quad c_k = \frac {1}{n} \sum _{j=0}^{n-1} f_j w_n^{-kj} \end{align*} mit \(k, j = 0, \dotsc , n - 1\), \(w_n = \exp (2 \pi \i / n)\) und \(f = (f_0, \dotsc , f_{n-1})^t\).

Die diskrete Fourier-Transformation entspricht der
Auswertung des trigonometrischen Polynoms

\begin{align*} p(x) = \sum _{k=0}^{n-1} c_k e^{\i k x} \end{align*} an den Punkten \(x_j = 2 \pi j / n\), d. h. \(f_j = p(x_j)\) für \(j = 0, \dotsc , n - 1\).

Die inverse Transformation kann als Riemann-Summe für die Fourier-Koeffizienten
interpretiert werden:

\begin{align*} \innerproduct {f, e_k}_{2\pi } = \frac {1}{2\pi } \int _0^{2\pi } f(x) e^{-\i k x} \dx \approx \frac {1}{n} \sum _{j=0}^{n-1} f(x_j) e^{-\i k x_j}, \quad x_j = 2 \pi j / n. \end{align*} Diese Approximation ist für glatte Funktionen und \(n \gg |k|\) sehr genau.

Schnelle Fourier-Transformation

Die diskrete Fourier-Transformation \(f_j = \sum _{k=0}^{n-1} c_k w_n^{jk}\) eines Vektors \(c = (c_0, \dotsc , c_{n-1})\) mit \(w_n = e^{2\pi \i /n}\) kann für \(n = 2^\ell \) mit der schnellen Fourier-Transformation (FFT) in
\(2n\ell = 2 n \log _2 n\) Operationen berechnet werden.

In der rekursiven Version hat der Algorithmus die folgende Form:

function f = FFT(c)
   n = length(c);
   if n = 1;
      f = c;
   else;
      g = FFT(c_0, c_2, ..., c_{n-2});
      h = FFT(c_1, c_3, ..., c_{n-1});
      p = (1, w_n, w_n^2, ..., w_n^{n/2-1});
      f = (g + p .* h, g - p .* h);
   end;

Die inverse Fourier-Transformation \(c_k = \frac {1}{n} \sum _{j=0}^{n-1} f_j w_n^{-jk}\) eines Vektors \(f = (f_0, \dotsc , f_{n-1})\) kann vollkommen analog berechnet werden. Man bezeichnet den entsprechenden Algorithmus als inverse schnelle Fourier-Transformation (IFFT).

Das Produkt \(r\) zweier Polynome

\begin{align*} p(x) = \sum _{k=0}^{m_p} p_k x^k, \quad q(x) = \sum _{k=0}^{m_q} q_k x^k \end{align*} kann mithilfe der FFT berechnet werden werden. Man wertet die Polynome an den komplexen Einheitswurzeln aus, multipliziert diese Werte und erhält die Koeffizienten \(r_k\) von \(r\) durch Rücktransformation.
Genauer wählt man zunächst \(n = 2^\ell > m_p + m_q\) und ergänzt die Koeffizienten der Polynome mit Nullen zu Vektoren \(\widetilde {p}\) und \(\widetilde {q}\) der Länge \(n\). Dann wird die diskrete Fourier-Transformation der Koeffizientenvektoren gebildet, d. h. \(\widehat {p} = \FFT (\widetilde {p})\), \(\widehat {q} = \FFT (\widetilde {q})\). Schließlich wird der Vektor der Produkte \(\widehat {r}_j = \widehat {p}_j \widehat {q}_j\), \(j = 0, \dotsc , n - 1\) berechnet und zurücktransformiert, d. h. \(\widetilde {r} = \IFFT (\widehat {r})\) ist der Koeffizientenvektor des Produktpolynoms \(r\).
Insgesamt werden \(\O (n \log n)\) Operationen benötigt, während die direkte Berechnung der Koeffizienten \(\O (n^2)\) Operationen erfordert.

Trigonometrische Interpolation

Für \(n = 2^\ell \) können die Koeffizienten des trigonometrischen Polynoms

\begin{align*} p(x) = c_m \cos (mx) + \sum _{|k| < m} c_k e^{\i k x}, \quad m = n/2, \end{align*} das die Daten

\begin{align*} f_j = f(x_j), \quad x_j = \frac {2\pi j}{n}, \quad j = 0, \dotsc , n - 1 \end{align*} interpoliert, mit der inversen schnellen Fourier-Transformation berechnet werden:

\begin{align*} (c_0, \dotsc , c_m, c_{-m+1}, \dotsc , c_{-1}) = \IFFT (f). \end{align*}

Die trigonometrische Interpolation in Verbindung mit der diskreten Fourier-Transformation kann zum Ausblenden hochfrequenter Störungen in Signalen verwendet werden. Man bildet zu den Daten \(f_j \approx f(x_j)\), \(x_j = \frac {2\pi j}{n}\), \(j = 0, \dotsc , n - 1\), \(n = 2^\ell \) zunächst mithilfe der IFFT das trigonometrische Interpolationspolynom \(p(x) = c_m \cos (mx) + \sum _{|j| < m} c_j e^{\i j x}\), \(m = n/2\).
Dann wählt man eine Bandbreite \(k\) und setzt alle Koeffizienten \(c_j\) mit \(|j| > k\) null.
Mit diesem Tiefpass werden für hinreichend kleines \(k\) im Allgemeinen Störungen unterdrückt.
Eine zu kleine Bandbreite führt dabei zu einem unerwünschten Genauigkeitsverlust.

Fourier-Transformation zyklischer Gleichungssysteme

Eine zyklische Matrix

\begin{align*} A = \begin{pmatrix} a_0 & a_{n-1} & \dots & a_1 \\ a_1 & a_0 & \dots & a_2 \\ \vdots & \vdots & & \vdots \\ a_{n-1} & a_{n-2} & \dots & a_0 \end {pmatrix} \end{align*} besitzt die Eigenwerte

\begin{align*} \lambda _j = \sum _{k=0}^{n-1} a_k w_n^{-kj}, \quad w_n = \exp (2 \pi i / n) \end{align*} und kann durch die Fourier-Matrix diagonalisiert werden:

\begin{align*} \frac {1}{n} W_n^\ast A W_n = \diag (\lambda ), \quad \lambda = W_n^\ast a. \end{align*} Folglich lässt sich die Lösung eines zyklischen Gleichungssystems \(Ax = b\) berechnen in der Form

\begin{align*} x = W_n \diag (\lambda )^{-1} (W_n^\ast b / n). \end{align*} Für \(n = 2^\ell \) ist die FFT anwendbar, und man erhält den folgenden Lösungsalgorithmus:

c = IFFT(b);
lambda = n * IFFT(a);
y_j = c_j / lambda_j, j = 0, ..., n - 1;
x = FFT(y);

Das Produkt \(C = AB\) zweier zyklischer Matrizen der Dimension \(n = 2^\ell \) lässt sich mithilfe der FFT berechnen. Zunächst bestimmt man dazu mit der mit \(n\) multiplizierten IFFT die Eigenwerte von \(A\) und \(B\):

\begin{align*} \lambda _j^A = \sum _{k=0}^{n-1} a_k w_n^{-jk}, \quad \lambda _j^B = \sum _{k=0}^{n-1} b_k w_n^{-jk}, \end{align*} wobei \(a\) bzw. \(b\) die erste Spalte von \(A\) bzw. \(B\) ist. Dann sind

\begin{align*} \lambda _j^C = \lambda _j^A \lambda _j^B \end{align*} die Eigenwerte von \(C\), und man erhält durch die mit \(1/n\) multiplizierte FFT von \(\lambda ^C\)

\begin{align*} c_k = \frac {1}{n} \sum _{j=0}^{n-1} \lambda _j^C w_n^{jk} \end{align*} die Elemente der ersten Spalte von \(C\). Damit ist \(C\) berechnet, denn \(C\) ist als Produkt zyklischer Matrizen wieder zyklisch.

Splines

Kubische Hermite-Interpolation

Funktionswerte und Ableitungen an zwei Punkten können durch ein kubisches Polynom interpoliert werden. Der Interpolant (Hermite-Spline) besitzt die Darstellung

\begin{align*} p = f(a) u_a + f(b) u_b + (b - a)(f’(a) v_a + f’(b) v_b) \end{align*} mit den Lagrange-Funktionen

\begin{align*} u_a(x) = (1 + 2s)(1 - s)^2, \quad u_b(x) = (3 - 2s)s^2, \quad v_a(x) = s(1 - s)^2, \quad v_b(x) = -s^2 (1 - s) \end{align*} und \(s = (x - a) / (b - a)\).

Sind Funktionswerte und Ableitungen an mehreren Punkten \(x_0 < \dotsb < x_n\) gegeben, so bilden die kubischen Hermite-Interpolaten einen stetig differenzierbaren kubischen Hermite-Spline \(q\). Nach Konstruktion ist \(q\) dabei eindeutig durch die Daten \(f(x_j), f’(x_j)\), \(j = 0, \dotsc , n\) bestimmt.

Kubische Splines

Ein kubischer Spline \(p\) zu einer Partition \(a = x_0 < \dotsb < x_n = b\) eines Intervalls \([a, b]\) kann (alternativ zur sog. B-Spline-Darstellung) durch seine Werte \(f_{k-1}\), \(f_k\) und Ableitungen \(d_{k-1}^+\), \(d_k^-\) an den Endpunkten der Teilintervalle \([x_{k-1}, x_k]\) festgelegt werden. Aus diesen Daten können die kubischen Polynome \(p_k\) auf den Intervallen \([x_{k-1}, x_k]\) mit kubischer Hermite-Interpolation berechnet werden.

Soll \(p\) an den Stützstellen glatt, d. h. differenzierbar sein, so werden Bedingungen an \(f_k\) und \(d_k^\pm \) gestellt. Stetige Differenzierbarkeit bei \(x_k\) ist äquivalent zu

\begin{align*} d_k^{-} = d_k = d_k^{+}. \end{align*} Soll auch die zweite Ableitung bei \(x_k\) stetig sein, so ist die Bedingung \(p_k”(x_k^-) = p_{k+1}”(x_k^+)\) äquivalent zu einer linearen Gleichung zwischen \(f_{k-1}, f_k, f_{k+1}\) und \(d_{k-1}^+, d_k, d_{k+1}^-\):

\begin{align*} \frac {1}{\Delta _k} d_{k-1}^+ + \left (\frac {2}{\Delta _k} + \frac {2}{\Delta _{k+1}}\right ) d_k + \frac {1}{\Delta _{k+1}} d_{k+1}^- = \frac {3}{\Delta _k^2} (f_k - f_{k-1}) + \frac {3}{\Delta _{k+1}^2} (f_{k+1} - f_k) \end{align*} mit \(\Delta _k := x_k - x_{k-1}\).

Natürliche Spline-Interpolation

Der natürliche Spline-Interpolant der Daten \((x_i, f_i)\), \(a = x_0 < \dotsb < x_n = b\) ist ein kubischer Spline \(p\), der an den Stützstellen \(x_i\) zweifach stetig differenzierbar ist und die Randbedingungen \(p”(x_0) = p”(x_n) = 0\) erfüllt.

Er minimiert unter allen glatten Interpolanten \(f\) das Integral

\begin{align*} \int _a^b |f”(x)|^2 \dx , \end{align*} das als Maß für die Stärke der Oszillation angesehen werden kann.

Die Ableitungen \(d_i = p’(x_i)\), die den Spline zusammen mit den Daten \(f_i\) festlegen, berechnen sich aus den Glattheitsbedingungen für \(i = 1, \dotsc , n - 1\), nämlich \(p_i”(x_i) = p_{i+1}”(x_i)\) bzw.

\begin{align*} \frac {1}{\Delta _i} d_{i-1} + \left (\frac {2}{\Delta _i} + \frac {2}{\Delta _{i+1}}\right ) d_i + \frac {1}{\Delta _{i+1}} d_{i+1} = \frac {3}{\Delta _i^2} (f_i - f_{i-1}) + \frac {3}{\Delta _{i+1}^2} (f_{i+1} - f_i), \end{align*} und den Randbedingungen \(p”(x_0) = p”(x_n) = 0\) bzw.

\begin{align*} 2d_0 + d_1 = \frac {3}{\Delta _1} (f_1 - f_0), \quad d_{n-1} + 2d_n = \frac {3}{\Delta _n} (f_n - f_{n-1}) \end{align*} mit \(\Delta _i = x_i - x_{i-1}\).

Alternativ kann man auch die Randbedingungen \(p’(a) = \alpha \), \(p’(b) = \beta \) stellen. Der resultierende eingespannte natürliche Spline minimiert dann ebenfalls obiges Integral.

Betrachtet man Splines \(p\), die die Lagrange-Daten \((x_k, \delta _{kj})\), \(k = 0, \dotsc , n\), \(j \in \{0, \dotsc , n\}\) interpolieren, so stellt man fest, dass \(p(x)\) schnell mit zunehmender Entfernung von \(x_j\) abklingt. Dieses numerisch günstige Verhalten ist typisch für Splines.

Außerdem können mit Splines auch gut Daten mit nicht-äquidistanten Stützstellen interpoliert werden. Nicht-äquidistante Stützstellen sind bspw. sinnvoll, falls die Daten Bereiche aufweisen, in denen sie unterschiedlich schnell schwanken.

Splineinterpolation mit MATLAB

Ein kubischer Spline-Interpolant zu den Daten \((x_k, y_k)\) kann in M ATLAB mit dem Befehl
p = spline(x, y); berechnet und mit pt = ppval(p, t); an den Punkten \(t_j\) ausgewertet werden. Der Spline wird als Struktur gespeichert, die unter anderem in dem Feld coefs die Koeffizienten der einzelnen Polynomsegmente enthält: Die Polynome werden dabei zeilenweise mit dem höchsten Koeffizienten zuerst gespeichert, wobei statt \(x\) der Term \(x - x_k\) mit \(x_k\) der unteren \(x\)-Stelle eingesetzt wird. Beispielsweise entspricht die Zeile \((1, -2, 0, 1)\) für das Intervall \([2, 3]\) dem Polynom \(p(x) = (x - 2)^3 - 2(x - 2)^2 + 1\).
Das Feld breaks enthält die Stützstellen.

M ATLAB verwendet dabei nicht die Randbedingungen \(p”(x_0) = p”(x_n) = 0\), sondern fordert stattdessen die Stetigkeit der dritten Ableitung an den Punkten \(x_1\) und \(x_{n-1}\) (Not-A-Knot-Bedingung). Dadurch ergibt sich ein genauerer Interpolant, jedoch geht dabei die oben erwähnte Minimal-Eigenschaft verloren. Der Unterschied zwischen den beiden verschiedenen Methoden ist für größere \(n\) allerdings fast zu vernachlässigen.

Alternativ dazu kann man den Datenvektor \(y\) um zwei Werte erweitern. Hat \(y\) genau zwei Einträge mehr als \(x\), so werden der erste und letzte Wert von \(y\) als Randbedingung für die Steigungen an den Enden der Kurve verwendet.

Die Auswertung kann auch unmittelbar mit dem Befehl spline erfolgen. Darüber hinaus ist die simultane Interpolation vektorwertiger Daten möglich.
In curve = spline(t, y, t_plt); stehen dabei in \(y\) spaltenweise die Funktionswerte in den Stützstellen \(x\), und in \(t_\text {plt}\) wird der Spline ausgewertet.

B-Splines

Knotenfolge

Eine Knotenfolge

\begin{align*} \tau \colon \dotsb \le \tau _{-1} \le \tau _0 \le \tau _1 \le \dotsb \end{align*} ist eine bi-infinite monoton wachsende Folge \(\{\tau _k\}_{k \in \integer }\) reeller Zahlen mit \(\lim _{k \to \pm \infty } \tau _k = \pm \infty \). Endliche Teilfolgen von \(\tau \) heißen Knotenvektoren. Die Vielfachheit \(\#\tau _k\) eines Knotens \(\tau _k\) ist die maximale Anzahl der Wiederholungen von \(\tau _k\) in der Folge bzw. Vektor \(\tau \). Man spricht dann von einfachen oder doppelten Knoten usw.

Rekursion für B-Splines

Zu einer Knotenfolge \(\tau \) definiert man die B-Splines \(b_k^n\) vom Grad \(n\) durch die Rekursion

\begin{align*} b_k^n := \gamma _k^n b_k^{n-1} + (1 - \gamma _{k+1}^n)b_{k+1}^{n-1}, \qquad \gamma _k^n(t) = \frac {t - \tau _k}{\tau _{k+n} - \tau _k}, \end{align*} ausgehend von den charakteristischen Funktionen \(b_0^k := \chi _{\left [\tau _k, \tau _{k+1}\right )}\) der Knotenintervalle
\(\left [\tau _k, \tau _{k+1}\right )\), d. h.

\begin{align*} b_0^k(t) := \begin{cases}1 & \tau _k \le t < \tau _{k+1} \\ 0 & \text {sonst}.\end {cases} \end{align*} Terme, für die der Nenner verschwindet, werden dabei nicht berücksichtigt.

Jeder B-Spline \(b_k^n\) wird durch seinen Knotenvektor \((\tau _k, \dotsc , \tau _{k+n+1})\) eindeutig festgelegt und verschwindet außerhalb von \(\left [\tau _k, \tau _{k+n+1}\right )\). Auf jeden nicht-leeren Knotenintervall \(\left [\tau _i, \tau _{i+1}\right )\),
\(k \le i \le k + n\) ist er ein nicht-negatives Polynom vom Grad \(n\).

Stetige Abhängigkeit vom Knotenvektor

Ist der Knotenvektor \(\tau = (\tau _k, \dotsc , \tau _{k+n+1})\) eines B-Splines \(b_k^n\) der Grenzwert einer Folge von Knotenvektoren \(\tau _\ell \), \(\ell \in \natural \) und bezeichnet \(b_{k,\ell }^n\) die zugehörigen B-Splines, so gilt

\begin{align*} \lim _{t \to \infty } b_{k,\ell }^n(t) = b_k^n(t) \end{align*} für alle \(t\), die nicht gleich einem der Knoten \(\tau _i\) sind. Die Konvergenz ist gleichmäßig auf jedem Intervall \([\alpha , \beta ]\), das keinen der Knoten von \(b_k^n\) enthält.

Die stetige Abhängigkeit von den Knoten ist nützlich für das Beweisen von Identitäten für Linearkombinationen von B-Splines. Gilt eine Gleichung \(\sum _k c_k(\tau ) b_k^n(t) = f(t, \tau )\) für einfache Knoten, so lässt sie sich durch ein Approximationsargument auf beliebige Knoten verallgemeinern. Dabei kann der Summationsbereich unendlich sein, da für jedes beschränkte Intervall nur endlich viele B-Splines nicht null sind.

Ableitung von B-Splines

Die Ableitung eines B-Splines vom Grad \(n\) zu einer Knotenfolge \(\tau \) ist eine gewichtete Differenz von zwei B-Splines vom Grad \(n - 1\). Auf jedem Knotenintervall \([\tau _i, \tau _{i+1})\) gilt

\begin{align*} (b_k^n)’ = \alpha _k^n b_k^{n-1} - \alpha _{k+1}^n b_{k+1}^{n-1}, \qquad \alpha _k^n := \frac {n}{\tau _{k+n} - \tau _k}, \end{align*} wobei Terme, die B-Splines mit leerem Träger enthalten, weggelassen werden.

Aus der Rekursion folgt, dass \(b_k^n\) an einem Knoten \(\tau _i\) \(n - m\)-mal stetig differenzierbar ist, falls \(\tau _i\) in der Folge \(\tau _k, \dotsc , \tau _{k+n+1}\) Vielfachheit \(m \le n\) hat. Insbesondere ist \(b_k^n\) stetig auf \(\real \), wenn keiner seiner Knoten Vielfachheit \(n + 1\) hat.

Uniforme B-Splines

Der uniforme B-Spline \(b^n\) vom Grad \(n > 0\) kann ausgehend von der charakteristischen Funktion \(b^0 := \chi _{[0,1)}\) des Intervalls \([0, 1]\) durch die Rekursion

\begin{align*} b^n(x) := \int _0^1 b^{n-1} (x - y) \dy \end{align*} definiert werden. Diese Identität ist äquivalent zu der Ableitungsformel

\begin{align*} \frac {d}{dx} b^n(x) = b^{n-1}(x) - b^{n-1}(x - 1) \end{align*} mit \(b^n(0) = 0\).

Als Spezialfall des allgemeinen B-Splines mit dem Knotenvektor \(\xi = (0, 1, \dotsc , n + 1)\) ist \(b^n\)

• positiv auf \((0, n + 1)\) und null außerhalb des Intervalls,

• ein Polynom vom Grad \(n\) auf jedem Knotenintervall \([k, k + 1]\) und

• \(n - 1\)-stetig differenzierbar.

Darüber hinaus gilt die Rekursionsformel

\begin{align*} n b^n(x) = x b^{n-1}(x) + (n + 1 - x) b^{n-1}(x - 1). \end{align*} Die B-Splines zu einer allgemeinen Knotenfolge \(\xi \colon \dotsc , -h, 0, h, \dotsc \) sind skalierte Translate von \(b^n\), d. h. \(b_k^n(x) = b^n(x/h - k)\), \(k \in \integer \).

Marsden-Identität

Für eine beliebige Knotenfolge \(\tau \) kann jedes Polynom vom Grad \(\le n\) als Linearkombination von B-Splines dargestellt werden. Insbesondere gilt für alle \(s \in \real \) die Marsden-Identität

\begin{align*} (t - s)^n = \sum _{k \in \integer } \psi _k^n(s) b_k^n(t), \qquad \psi _k^n(s) := (\tau _{k+1} - s) \dotsm (\tau _{k+n} - s). \end{align*}

Durch Ableiten der Identität nach \(s\) und Nullsetzen von \(s\) erhält man explizite Formeln für die Monome \(t^m\). Beispielsweise ist

\begin{align*} 1 = \sum _k b_k^n(t), \qquad t = \sum _k \tau _k^n b_k^n(t) \end{align*} mit \(\tau _k^n := (\tau _{k+1} + \dotsb + \tau _{k+n}) / n\) den sogenannten Knotenmitteln.

Splines

Die Splines \(S_\tau ^n\) vom Grad \(\le n\) zu einer Knotenfolge \(\tau \) sind Linearkombinationen von B-Splines:

\begin{align*} S_\tau ^n \ni p := \sum _{k \in \integer } c_k b_k^n. \end{align*} Anders ausgedrückt besteht \(S_\tau ^n\) aus allen Funktionen \(t \mapsto p(t)\), \(t \in \real \), die auf jedem Intervall \([\tau _k, \tau _{k+1})\) Polynome vom Grad \(\le n\) sind und an einem Knoten mit Vielfachheit \(m \le n\) mindestens \(n - m\)-mal stetig differenzierbar sind.

Splines \(S_\tau ^n(D)\) auf beschränkten Intervallen \(D\) erhält man, indem die Variable \(t\) auf \(D\) eingeschränkt wird. Es sind nur die B-Splines \(b_k^n\), die auf einem Teilintervall von \(D\) nicht null sind, und ihre Knotenvektoren \((\tau _k, \dotsc , \tau _{k+n+1})\) relevant. Die entsprechenden Indizes werden mit \(k \sim D\) bezeichnet:

\begin{align*} p(t) = \sum _{k \sim D} c_k b_k^n(t), \quad t \in D. \end{align*} Insbesondere sind \(k = \ell - n, \dotsc , \ell \) die relevanten Indizes für ein nicht-leeres Knotenintervall \(D = [\tau _\ell , \tau _{\ell +1})\).

Auswertung von Splines (de-Boor-Algorithmus)

Ein Spline

\begin{align*} p = \sum _k c_k b_k^n \in S_\tau ^n \end{align*} kann in \(t \in [\tau _\ell , \tau _{\ell +1})\) durch Bilden von Konvexkombinationen der Koeffizienten der relevanten B-Splines \(b_k^n\), \(k \sim t\), ausgewertet werden (de-Boor-Algorithmus).

Beginnend mit

\begin{align*} p_k^0 := c_k, \quad k = \ell - n, \dotsc , \ell \end{align*} berechnet man sukzessive für \(m = 0, \dotsc , n - 1\)

\begin{align*} p_k^{m+1} := \gamma _k^{n-m} p_k^m + (1 - \gamma _k^{n-m}) p_{k-1}^m, \quad k = \ell - n + m + 1, \dotsc , \ell \end{align*} mit

\begin{align*} \gamma _k^{n-m} := \frac {t - \tau _k}{\tau _{k+n-m} - \tau _k} \end{align*} und erhält \(p(t)\) als den letzten Wert \(p_\ell ^n\).

Die \(p_k^m\) können in einem Dreiecksschema berechnet werden. Für \(t = \tau _\ell \) vereinfacht es sich etwas: Hat \(\tau _\ell \) Vielfachheit \(r\), dann ist \(p(t) = p_{\ell -r}^{n-r}\), d. h. nur \(n - r\) Schritte der Rekurstion werden benötigt.

Ableitung von Splines

Sei \(\tau \) eine Knotenfolge mit Vielfachheiten \(\le n\). Die Ableitung eines Splines in \(S_\tau ^n\) ist ein Spline vom Grad \(\le n - 1\) zur gleichen Knotenfolge:

\begin{align*} \left (\sum _{k \in \integer } c_k b_k^n\right )’ = \sum _{k \in \integer } \alpha _k^n \nabla c_k b_k^{n-1} \in S_\tau ^{n-1} \quad \text {mit}\quad \alpha _k^n := \frac {n}{\tau _{k+n} - \tau _k}, \quad \nabla c_k := c_k - c_{k-1}. \end{align*} Enthält \(\tau \) Knoten mit Vielfachheiten \(> n\), an denen die Splines in \(S_\tau ^n\) Sprünge haben können, so behält die Gleichung auf jedem Intervall, auf dem der Spline stetig ist, ihre Gültigkeit. In diesem Fall werden Ausdrücke mit verschwindenen Nennern, die B-Splines mit leerem Träger entsprechen, weggelassen.

Für Splines auf einem beschränkten Parameterintervall \([\alpha , \beta ]\) beschränkt man die Summation auf die relevanten B-Splines. Genauer sind für eine Knotenfolge mit

\begin{align*} \tau _0 \le \tau _1 < \alpha = \tau _n < \tau _{n+1} \le \dotsb \le \tau _{m-1} < \tau _m = \beta < \tau _{m+n-1} \le \tau _{m+n} \end{align*} die B-Splines

\begin{align*} b_k^n \quad \text {mit}\quad k = 0, \dotsc , m - 1 \qquad \text {und}\qquad b_k^{n-1} \quad \text {mit}\quad k = 1, \dotsc , m - 1 \end{align*} relevant. Dies ist konsistent zur Differenzbildung bei den Koeffizienten, die den Bereich der Indizes um Eins reduziert.

Schoenberg-Schema

Schoenbergs Schema benutzt Funktionswerte an den Knotenmitteln
\(\tau _k^n := (\tau _{k+1} + \dotsb + \tau _{k+n})/n\) als Koeffizienten
einer Spline-Approximation einer glatten Funktion \(f\):

\begin{align*} f \mapsto Qf := \sum _{k \in \integer } f(\tau _k^n) b_k^n \in S_\tau ^n. \end{align*} Es hat die Fehlerordnung zwei, d. h.

\begin{align*} |f(t) - Qf(t)| \le \frac {1}{2} \norm {f”}_{\infty ,\; D_t} h(t)^2 \end{align*} mit \(\tau _\ell \le t < \tau _{\ell +1}\) und

\begin{align*} D_t := [\tau _{\ell -n}^n, \tau _\ell ^n], \qquad h(t) := \max _{k=\ell -n,\dotsc ,\ell } |\tau _k^n - t|. \end{align*}

Der Schoenberg-Operator erhält Positivität, Monotonie und Konvexität, d. h.

\begin{align*} f^{(m)} \ge 0 \quad \Rightarrow \quad (Qf)^{(m)} \ge 0 \end{align*} für \(m \le 2\), falls beide Ableitungen existieren. Für eine äquidistante Knotenfolge bleibt das Vorzeichen aller Ableitungen bis zur Ordnung \(n\) erhalten.

Quasi-Interpolant

Ein lineares Approximationsschema

\begin{align*} f \mapsto Qf := \sum _{k \in \integer } (Q_k f) b_k^n \in S_\tau ^n \end{align*} für stetige Funktionen \(f\) bezeichnet man als Quasi-Interpolant, falls

• \(Q_k\) lokale beschränkte lineare Funktionale sind, d. h.

  \(\seteqnumber{0}{}{0}\)

  \begin{align*} |Q_k f| \le \norm {Q} \cdot \norm {f}_{\infty ,\; D_k} \end{align*} mit \(\norm {f}_{\infty ,\; D_k} := \sup _{\tau \in [\tau _k, \tau _{k+n+1})} |f(t)|\), und

• \(Q\) für Polynome \(p\) vom Grad \(\le n\) exakt ist, d. h. \(Qp = p\).

Äquivalent zur zweiten Bedingung ist, dass \(Q_k p = \psi _k(s)\) für alle \(s \in \real \) mit \(p(t) := (t - s)^n\) und \(\psi _k(s) := (\xi _{k+1} - s) \dotsm (\xi _{k+n} - s)\). Diese Identität für Polynome vom Grad \(\le n\) kann man durch Koeffizientenvergleich oder durch Auswertung an \(n + 1\) Punkten prüfen.

Fehler der Quasi-Interpolation

Für den Fehler eines Quasi-Interpolanten

\begin{align*} f \mapsto Qf = \sum _k (Q_k f) b_k^n \in S_\tau ^n \end{align*} gilt

\begin{align*} |f(t) - (Qf)(t)| \le \frac {\norm {Q}}{(n + 1)!} \norm {f^{(n+1)}}_{\infty ,\; D_t} h(t)^{n+1} \end{align*} mit \(D_t\) der Vereinigung der Träger aller für \(t\) relevanten B-Splines und \(h(t) := \max _{s \in D_t} |s - t|\).

Ist das lokale Gitterverhältnis

\begin{align*} r_\tau := \sup _{\tau _{j-1} < \tau _j = \tau _k < \tau _{k+1}} \max \left ( \frac {\tau _{k+1} - \tau _k}{\tau _j - \tau _{j-1}}, \frac {\tau _j - \tau _{j-1}}{\tau _{k+1} - \tau _k}\right ) \end{align*} beschränkt, so lässt sich ebenfalls der Fehler der Ableitungen abschätzen:

\begin{align*} |f^{(j)}(t) - (Qf)^{(j)}(t)| \le \const (n, r) \norm {Q} \norm {f^{(n+1)}}_{\infty ,\; D_t} h(t)^{n+1-j} \end{align*} für alle \(j \le n\), für die die Ableitungen existieren.

Lösbarkeit von Interpolationsproblemen mit B-Splines

Die Koeffizienten eines Splines \(p = \sum _{k=1}^m c_k b_k^n\), der die Daten \(f_i\) an einer monoton wachsenden Folge von Punkten \(t_i\) interpoliert, werden durch das LGS

\begin{align*} Ac = f, \quad a_{i,k} := b_k^n(t_i) \end{align*} bestimmt. Eine eindeutige Lösung existiert für alle Daten \(f\) genau dann, wenn

\begin{align*} b_k^n(t_k) > 0, \quad k = 1, \dotsc , m. \end{align*}