Interpolation - Numerische Mathematik (i282)

Für die Umwandlung von Punktfolgen in Funktionen haben wir bisher nur die Ausgleichsrechnung kennengelernt. Liegen nur wenige und kaum fehlerbehaftete Punkte vor, so ist neben der Approximation der Punkte durch eine Funktion auch die Interpolation möglich, also das Finden einer Funktion, die exakt durch alle gegebenen Punkte verläuft.

Um Funktionen mit dem Computer beschreiben zu können, müssen diese durch endlich viele Werte vollständig festlebar sein. Bei der Ausgleichsrechnung hatten wir deshalb auf gewisse Funktionenklassen beschränkt, in denen jede Funktion durch endlich viele Parameter beschrieben war. Bei der Interpolation beschränkt man sich praktisch immer auf Polynome

a_0+a_1\,x+a_2\,x^2+\cdots+a_{n-1}\,x^{n-1}+a_n\,x^n

(1)

oder sogenannte trigonometrische Polynome (werden hier nicht behandelt). Gründe:

Auswertung eines Polynoms an einer Stelle $x$ erfordert nur wenige Grundoperationen (Multiplikationen, Additionen).
Polynome sind bei theoretischen Betrachtungen sehr gut handhabbar.
Polynome sind beliebig oft differenzierbar, also in vielfältigen Szenarien einsetzbar (Optimierung!).
Polynome sind leicht integrierbar (Flächenberechnungen!).

Beschäftigen uns zunächst mit der klassischen Polynominterpolation (ein Polynom auf dem gesamte Definitionsbereich) und geben dann einen Überblick über Splines (stückweise polynomiale Funktionen).

8.1Formale Problemstellung¶

Gegen ist eine Folge von Punkten $(x_1,y_1),\ldots,(x_n,y_n)\in\bbR^2$ . Finde ein Polynom $p:[a,b]\to\bbR$ (oder einen Spline), sodass

p(x_k)=y_k\quad\text{für }k=1,\ldots,n

(2)

gilt, wobei die Intervallgrenzen $a$ und $b$ so sein sollen, dass $x_k\in[a,b]$ für $k=1,\ldots,n$ .

Die Zahlen $x_1,\ldots,x_n$ werden Stützstellen genannt. Im Folgenden nehmen wir stets an, dass alle Stützstellen verschieden sind, also $x_k\neq x_l$ für $k\neq l$ gilt.

8.2Polynominterpolation¶

8.2.1Direkter Ansatz¶

Sind $n$ zu interpolierende Punkte gegeben, so erscheint ein Polynom $(n-1)$ -ten Grades als Interpolante zweckmäßig, da dieses genau $n$ frei wählbare Parameter enthält. Setzt man die einzelnen Punkte in die Interpolationsbedingung (1) ein, so erhält man das lineare Gleichungssystem

{\begin{bmatrix} 1&x_1&x_1^2&\cdots&x_1^{n-1}\\ \vdots&\vdots&\vdots&&\vdots\\ 1&x_n&x_n^2&\cdots&x_n^{n-1} \end{bmatrix} \begin{bmatrix}a_0\\\vdots\\a_{n-1}\end{bmatrix} =\begin{bmatrix}y_1\\\vdots\\y_n\end{bmatrix}}

(3)

Die Systemmatrix heißt Vandermonde-Matrix und ist stets invertierbar (wenn alle $x_k$ verschieden sind). Das Lösen des Gleichungssystems benötigt etwa $n^3$ Rechenoperationen und ist damit deutlich aufwendiger als weiter unten behandelte Verfahren zur Polynominterpolation.

8.2.2Lagrange-Interpolation¶

Statt das obige Gleichungssystemzu lösen, können wir für das gesuchte Interpolationspolynom auch direkt eine Formel angeben. Benötigen dazu die sogenannten Langrange’schen Basispolynome vom Grad $n-1$ zu den Stützstellen $x_1,\ldots,x_n$ :

L_{n-1,k}(x):=\prod_{l\neq k}\frac{x-x_l}{x_k-x_l},\quad k=1,\ldots,n.

(4)

Offensichtlich gilt

L_{n-1,k}(x_l)=\begin{cases}1,&\text{wenn }l=k,\\0,&\text{wenn }l\neq k.\end{cases}

(5)

import numpy as np
import matplotlib.pyplot as plt

colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red']

nodes = [0.2, 0.5, 0.7, 0.9]
n = len(nodes)
x = np.linspace(0.1, 1, 100)

fig, ax = plt.subplots()

for k in range(n):
    y = np.ones_like(x)
    for l in range(n):
        if l != k:
            y *= (x - nodes[l]) / (nodes[k] - nodes[l])
    ax.plot(x, y, '-', label=f'$L_{{{n-1},{k+1}}}(x)$', color=colors[k])
    ax.plot(nodes[k], 1, 'o', color=colors[k])
ax.plot(nodes, np.zeros(n), 'ok')

ax.set_xticks(nodes)
ax.set_xticklabels([f'$x_{{{k}}}$' for k in range(n)])
ax.set_xlabel('x')
ax.set_xlim(x.min(), x.max())
ax.grid()
ax.legend()
plt.show()

Somit ist

y_1\,L_{n-1,1}(x)+\cdots+y_n\,L_{n-1,n}(x)

(6)

das gesuchte Interpolationspolynom. In dieser Darstellung wird das Interpolationspolynom als Langrange’sches Interpolationspolynom bezeichnet.

Damit ist die Existenz eines interpolierenden Polynoms vom Grad $n-1$ gezeigt. Die Eindeutigkeit folgt leicht aus der Tatsache, dass ein Polynom vom Grad $n-1$ höchstens $n$ Nullstellen hat (IDVID 810).

Beachte: Durch die zu interpolierenden Punkte ist das Interpolationspolynom eindeutig bestimmt. Ja nach Darstellung dieses Polynoms als Formel gibt man ihm aber verschiedenen Namen, z.B. Langrange’sches Interpolationspolynom und weiter unten dann Newton’sches Interpolationspolynom. Das durch die jeweilige Formel beschriebene Polynom ist stets das selbe.

Das Lagrange’sche Interpolationspolynom ist für den praktischen Einsatz kaum von Bedeutung, leistet aber bei theoretischen Überlegungen (Existenz, Eindeutigkeit, Kondition) sehr gute Dienste.

8.2.3Anwendungsfälle¶

Polynominterpolation ist Baustein in zahlreichen numerischen Verfahren. Grundsätzlich sind zwei Anwendungsfälle zu unterscheiden:

Das Interpolationspolynom soll nur an wenigen Stellen ausgewertet werden (z.B. Interpolation von 6-stündlichen Temperaturmessungen um 0 Uhr, 6 Uhr, 12 Uhr, 18 Uhr zum Erhalt einer Temperatur um 8 Uhr).
Das Interpolationspolynom selbst ist von Interesse oder soll an sehr vielen Stellen ausgewertet werden (z.B. Ersatz einer komplizierten Funktion wie $\sin$ durch ein Polynom).

Im ersteren Fall lohnt sich das Aufstellen des Polynoms nicht. Man setzt hier üblicherweise das Aitken-Neville-Verfahren (siehe unten) zur effizienten Auswertung des (nicht explizit aufgestellten) Interpolationspolynoms an einer Stelle ein.

Im zweiten Fall ist man bestrebt, das Interpolationspolynom in einer Form aufzustellen, die einerseits effizient zu bekommen ist und andererseits die effiziente Auswertung des Polynoms ermöglicht (siehe unten).

8.2.4Kondition¶

Üblicherweise sind die Stützstellen $x_1,\ldots,x_n$ nicht oder nur geringfügig mit Fehlern behaftet. Das die zugehörigen $y_1,\ldots,y_n$ jedoch meist aus Messungen entstehen, ist hier mit relevanten Fehlern zu rechnen. Für die praktische Arbeit stehen uns als statt der Punkte $(x_1,y_1),\ldots,(x_n,y_n)$ nur die fehlerbehafteten Punkte $(x_1,\tilde{y_1}),\ldots,(x_n,\tilde{y_n})$ zur Verfügung. Daraus ergeben sich das exakte Interpolationspolynom $p$ und das auf fehlerhaften Eingaben beruhende Interpolationspolynom $\tilde{p}$ .

8.2.4.1Fehlermaße¶

Als Fehlermaß für die Eingaben wählen wir den relativen Fehler bzgl. der $\infty$ -Norm der Vektoren $y$ bzw $\tilde{y}$ , also

\frac{\|y-\tilde{y}\|_\infty}{\|y\|_\infty},

(7)

da wir für den Abstand zweier Polynome mangels Kenntnis anderer Abstandsmaße zwischen Funktionen ebenfalls den maximalen punktweisen Abstand wählen müssen. Für eine stetige Funktion $f:[a,b]\to\bbR$ setzen wir

\|f\|_\infty:=\max_{x\in[a,b]}|f(x)|

(8)

und drücken den Ausgabefehler mittels

\frac{\|p-\tilde{p}\|_\infty}{\|p\|_\infty}.

(9)

aus.

8.2.4.2Fehlerabschätzung¶

Stellen wir $p$ und $\tilde{p}$ als Lagrange’sches Interpolationspolynom dar, also

p(x)=y_1\,L_{n-1,1}(x)+\cdots+y_n\,L_{n-1,n}(x)

(10)

und

\tilde{p}(x)=\tilde{y}_1\,L_{n-1,1}(x)+\cdots+\tilde{y}_n\,L_{n-1,n}(x),

(11)

so erhalten wir

\frac{\|p-\tilde{p}\|_\infty}{\|p\|_\infty}\leq\kappa\,\frac{\|y-\tilde{y}\|_\infty}{\|y\|_\infty}

(12)

mit

\kappa:=\max_{x\in[a,b]}\sum_{k=1}^n|L_{n-1,k}(x)|

(13)

(IDVID 620).

Aus der Definition der Langrange’schen Basispolynome $L_{n-1,k}$ sieht man leicht (schauen wir uns hier aber nicht an), dass $\kappa$ bei festem $n$ nicht von der Lage des Intervalls $[a,b]$ abhängt, sondern nur von der Verteilung der Stützstellen innerhalb des Intervalls. Auch die Werte $y_1,\ldots,y_n$ gehen nicht in $\kappa$ ein. Die Frage nach der Kondition ist also ausschließlich eine Frage der Stützstellenverteilung.

Merke!

Für äquidistante Stützstellen, also $x_k=a+(b-a)\,\frac{k-1}{n-1}$ , kann man zeigen, dass

\kappa\approx\frac{2^n}{(n-1)\,\log(n-1)}

(14)

für große $n$ gilt, also die Kondition schlecht ist (für große $n$ ).

Wählt man hingegen beispielsweise die sogenannten Tschebyscheff-Stützstellen

x_k=a+(b-a)\,\left(\frac{1}{2}-\frac{1}{2}\,\cos\left(\frac{2\,k-1}{2\,n}\,\pi\right)\right),

(15)

so kann man

\kappa\approx\frac{2}{\pi}\,\log(n-1)

(16)

für große $n$ zeigen. Man erhält also eine gute Kondition. Die Tschebyscheff-Stützstellen nicht äquidistant, sondern verdichten sich mit zunehmendem Abstand von der Intervallmitte.

8.2.5Algorithmen¶

8.2.5.1Direkter Ansatz und Horner-Schema¶

Das Aufstellen des Interpolationspolynoms mittels Vandermonde-Matrix benötigt ca. $n^3$ Rechenopertationen (Gauss-Algorithmus). Auch hat die Vandermonde-Matrix oft schon bei niedrigem zweistelligem $n$ eine zu hohe Konditionszahl, sodass dieser Ansatz in der Praxis keine Rolle spielt.

Das anschließende Auswerten des Polynoms an einer Stelle erfordert \frac{n,(n+1)}{2} Multiplikationen und $n$ Additionen, liegt also in der Größenordnung von $n^2$ Grundoperationen.

Liegt das Polynom in der Form

p(x)=a_0+a_1\,x+a_2\,x^2+\cdots+a_{n-2}\,x^{n-2}+a_{n-1}\,x^{n-1}

(17)

vor, so kann man die Auswertung an einer Stelle noch beschleunigen, denn offensichtlich gilt

p(x)=a_0+x\,(a_1+x\,(a_2+\cdots+x\,(a_{n-2}+x\,a_{n-1})\cdots))

(18)

Der daraus abgeleitete Algorithmus zur Polynomauswertung heißt Horner-Schema und benötigt nur $n-1$ Multiplikationen und $n-1$ Additionen, liegt also in der Größenordnung $n$ :

Setze $y:=a_{n-1}$ .
Wiederhole für $k=n-2,\ldots,0$ :
1. Setze $y:=a_k+x\,y$ .

Bei gegebenen $a_0,\ldots,a_{n-1}$ und $x$ liefert dieser Algorithmus $y=p(x)$ .

Das Horner-Schema ist für nicht zu große $n$ stabil wie man mittels Rückwärtsanalyse zeigen kann (tun wir hier nicht). Insbesondere sind sie Zwischenergebnisse kleiner als beim direkten Auswerten (Potenzen von $x$ !), sodass die Gefahr von Exponentenüberläufen geringer ist.

8.2.5.2Aitken-Neville-Schema¶

Soll das Interpolationspolynom nur an wenigen Stellen ausgewertet werden, so muss das Polynom gar nicht explizit aufgestellt werden. Um auch ohne Kenntnis des Polynoms $p$ den Wert $p(x)$ für ein Stelle $x$ zu bekommen, verwenden wir den folgenden Zusammenhang zwischen den Interpolationspolynom vom Grad $n-1$ und zwei Interpolationspolynomen vom Grad $n-2$ zu unterschiedlichen Stützpunkten:

Bezeichnen wir mit $P_{l,k}$ den Wert des Interpolationspolynoms vom Grad $l-1$ zu den Punkten $x_k,\ldots,x_{k+l-1}$ an einer festen Stelle $x$ , so erhalten wir den als Aitken-Neville-Schema bekannten Zusammenhang

P_{l,k}=\frac{x_{k+l-1}-x}{x_{k+l-1}-x_k}\,P_{l-1,k}+\frac{x-x_k}{x_{k+l-1}-x_k}\,P_{l-1,k+1}.

(20)

Ausgehend von

P_{1,k}=y_k,\quad k=1,\ldots,n

(21)

erhalten wir daraus Schritt für Schritt den gesuchten Wert

P_{n,1}=p(x)

(22)

(IDVID 635). Das Aitken-Neville-Schema benötigt etwa $n^2$ Rechenoperationen und lässt sich mit $n$ Speicherplätzen umsetzen, wenn man nich mehr benötigte Zwischenergebnisse wieder überschreibt:

Setze $p_1:=y_1$ .
Wiederhole für $i=2,\ldots,n$ :
1. Setze $p_i:=y_i$ .
2. Wiederhole für $j=i-1,\ldots,1$ :
  1. Setze
    $p_j:=\frac{(x_i-x)\,p_j+(x-x_j)\,p_{j+1}}{x_i-x_j}.$
    (23)
Liefere $p_1$ als Ergebnis.

Kann man die Eingaben $y_1,\ldots,y_n$ überschreiben, so wird kein zusätzlicher Speicher benötigt.

Die Rekursionsvorschrift in Schritt 2.2.1 des Algorithmus ist als stabil anzusehen. Auslöschungseffekte sind im Allgemeinen nicht zu erwarten.

8.2.5.3Newton-Interpolation¶

Ist das Interpolationspolynom explizit gesucht, z.B. um es an vielen verschiedenen Stellen auswerten zu können, kann man ähnlich zum Aitken-Neville-Schema ebenfalls eine rekursive Vorschrift finden, die Schritt für Schritt Interpolationspolynome der Ordnungen $0,1,\ldots,n-1$ zu stets um einen Punkt erweiterten Punktmengen liefert.

Merke!

Mit den Newton’schen Basispolynomen

N_0(x):=1

(24)

und

N_l(x):=(x-x_1)\cdots(x-x_l)

(25)

für $l=1,\ldots,n-1$ kann das Interpolationspolynom $p$ vom Grad $n-1$ zu den Punkten $(x_1,y_1),\ldots,(x_n,y_n)$ als Newton’sches Interpolationspolynom

p(x)=\sum_{l=1}^n b_{l,1}\,N_{l-1}(x)

(26)

dargestellt werden. Dabei sind die Koeffizienten $b_{l,1}$ durch die rekursive Berechnungsvorschrift

b_{1,k}:=y_k

(27)

für $k=1,\ldots,n$ und

b_{l,k}:=\frac{b_{l-1,k+1}-b_{l-1,k}}{x_{k+l-1}-x_k}

(28)

für $l=2,\ldots,n$ und $k=1,\ldots,n-l+1$ gegeben.

Alternativ können die Newton’schen Interpolationspolynome $p_0,p_1,\ldots,p_{n-1}=p$ zu den von ${(x_1,y_1)}$ ausgehenden und stets um einen Punkt wachsenden Punktmengen über die rekursive Vorschrift

p_0(x)=y_0,\quad p_{l-1}(x)=p_{l-2}(x)+b_{l,1}\,N_{l-1}(x)

(29)

erhalten werden.

Herleitung der obigen Formeln: IDVID 650.

import numpy as np
import matplotlib.pyplot as plt

nodes = [0.2, 0.5, 0.7, 0.9]
n = len(nodes)
x = np.linspace(0.1, 1, 100)

fig, ax = plt.subplots()

y = np.ones_like(x)
for k in range(n):
    ax.plot(x, y, '-', label=f'$N_{{{k}}}(x)$')
    y *= (x - nodes[k])
ax.plot(nodes, np.zeros(n), 'ok')

ax.set_xticks(nodes)
ax.set_xticklabels([f'$x_{{{k}}}$' for k in range(n)])
ax.set_xlabel('x')
ax.set_xlim(x.min(), x.max())
ax.grid()
ax.legend()
plt.show()

Die Berechnung der Koeffizienten $b_{1,1},\ldots,b_{n,1}$ benötigt ca. $n^2$ Operationen, ist also deutlich schneller möglich als über das Lösen eines Gleichungssystems. Der Algorithmus dafür ist völlig analog zum Aitken-Neville-Schema.

Das Auswerten des des Newton’schen Interpolationspolynoms an einer Stelle $x$ kann ähnlich zum Horner-Schema in ca. $n$ Operationen erfolgen, denn es gilt

p(x)=b_{1,1}+(x-x_1)\,(b_{2,1}+\cdots+(x-x_{n-2})\,(b_{n-1,1}+(x-x_{n-1})\,b_{n,1})\cdots).

(30)

Algorithmus zum Auswerten:

Setze $P:=b_{n,1}.$
Wiederhole für $l=n-1,\ldots,1$ :
1. Setze $P:=b_{l,1}+(x-x_l)\,P$ .

Am Ende des Algorithmus gilt $P=p(x)$ .

8.2.6Approximationsfehler¶

Nutzt man durch Interpolation erhaltene Polynome als einfach handhabbaren Ersatz für kompliziertere Funktionen, so stellt sich die Frage nach dem Abstand $\|f-p\|_\infty$ des Interpolationspolynoms $p$ zur ursprünglichen Funktion $f$ zwischen den Stützstellen. Diese Frage ist eine andere als die nach der Kondition, da hier Fehler in den Eingaben gar keine Rolle spielen.

Man kann folgendes Resultat zeigen (tun wir hier aber nicht):

Einige Bemerkungen zur Interpretation:

Der Faktor $\frac{1}{n!}$ wird umso kleiner je mehr Stützstellen genutzt werden.
Der Faktor $\|f^{(n)}\|_\infty$ ist unabhängig von der Lage der Stützstellen und wird meistens kleiner je größer $n$ ist. Nur Polynome erfüllen $\|f^{(n)}\|_\infty=0$ für alle hinreichend großen $n$ . Grob formuliert: Je schneller $\|f^{(n)}\|_\infty$ gegen Null geht für große $n$ , desto polynomähnlicher ist die Funktion $f$ .
Der Faktor $\|N_n\|_\infty$ hängt nicht von $f$ ab und auch nicht vom Intervall $[a,b]$ (Verschiebung des Intervall ändert den Wert nicht), sondern nur von der Verteilung der Stützstellen innerhalb des Intervalls. Für äquidistante Stützstellen wird dieser Faktor groß sein, für die Tschebyscheff-Stützstellen hingegen klein (vgl. Betrachtungen zur Kondition).
Bei festem $n$ wird $\|N_n\|_\infty$ umso kleiner je kleiner das Intervall $[a,b]$ ist.

Daraus lassen sich folgende Grundsätze für die Polynominterpolation ableiten:

Nur auf kleinen Intervallen interpolieren. Falls nötig, mehrer Interpolationspolynome “zusammenstückeln” (siehe Splines unten).
Falls kleine Intervalle nicht möglich, dann Stützstellen günstig wählen, z.B. Tschebyscheff-Stützstellen.

Das folgende Beispiel zeigt, dass eine Erhöhung der Stützstellenzahl nicht automatisch den Approximationsfehler verbessert.

Beispiel

Wird die Funktion

f(x)=\frac{1}{1+x^2}

(33)

auf dem Intervall $[-5,5]$ an $n$ äquidistanten Stützstellen durch ein Polynom $p_{n-1}$ interpoliert, so kann man

\lim_{n\to\infty}\|f-p_{n-1}\|_\infty=\infty

(34)

zeigen (tun wir hier aber nicht).

import numpy as np
import matplotlib.pyplot as plt

f = lambda x: 1 / (1 + x ** 2)
x = np.linspace(-5, 5, 300)

fig, ax = plt.subplots()

for n in [5, 9, 17]:
    nodes = np.linspace(-5, 5, n)
    b = f(nodes)
    N = np.ones_like(x)
    y = b[0] * N
    for i in range(1, n):
        for j in range(i - 1, -1, -1):
            b[j] = (b[j+1] - b[j]) / (nodes[i] - nodes[j])
        N *= (x - nodes[i - 1])
        y += b[0] * N
    ax.plot(x, y, '-', label=f'$p_{{{n-1}}}(x)$')

ax.plot(nodes, f(nodes), 'o', color='tab:green')
ax.plot(x, f(x), '-k', label='f(x)')

ax.set_xlabel('x')
ax.set_xlim(x.min(), x.max())
ax.set_ylim(-1.1, 1.5)
ax.grid()
ax.legend()
plt.show()

8.3Splines¶

8.3.1Idee¶

Hatten festgestellt, dass der Approximationsfehler bei Interpolation auf kleinen Intervallen geringer ist als auf großen Intervallen. Auch führt ein hoher Polynomgrad (also viele Stützstellen) nicht automatisch zu besserer Ergebnissen. Aus diesen Beobachtungen sind die so genannten Splines (Englisch für Straklatte) als Ersatz für Interpolationspolynome entstanden:

Zu gegebenen Stützstellen $x_1,\ldots,x_n$ finde für jedes Intervall $[x_k,x_{k+1}]$ ein separates Interpolationspolynom vom Grad $m$ (meist klein, z.B. $m=3$ ).
An den Intervallenden sollen die Ableitungen der benachbarten Polynome bis zur Ordnung $m-1$ übereinstimmen, damit die Verbindungsstellen “glatt” sind.
Bei $x_1$ und $x_n$ können in Summe $m-1$ Ableitungen frei gewählt werden.

Aus diesem Ansatz ergeben sich

$2\,(n-1)$ Bedingungen für die intervallweise Interpolation,
$(m-1)\,(n-2)$ Bedingungen für die Ableitungen an den inneren Intervallgrenzen,
$m-1$ Bedingungen für die Ableitungen an den äußeren Intervallgrenzen,

also insgesamt $(m+1)\,(n-1)$ Bedingungen. Die Anzahl der zu bestimmenden Koeffizienten bei $n-1$ Polynomen vom Grad $m$ ist ebenfalls $(m+1)\,(n-1)$ . Alle Bedingungen liefern lineare Gleichungen, sodass die Spline-Koeffizienten als Lösung eines linearen Gleichungssystems berechnet werden können.

Für Splines existiert analog zur Polynominterpolation eine vollständige, ausgereifte Theorie und auch eine umfangreich Software-Unterstützung.

8.3.2Beispiel¶

Wollen die Funktion

f(x)=\frac{1}{1+x^2}

(35)

auf dem Intervall $[-5,5]$ durch einen kubischen Spline approximieren. Wählen die Stützstellen $\{-5,-2,0,2,5\}$ , also 4 Teilintervalle. An den äußeren Intervallgrenzen soll die erste Ableitung des Splines der ersten Ableitung der Funktion gleichen.

Die vier Polynome $s_1$ , $s_2$ , $s_3$ , $s_4$ , die zusammen den Spline ergeben, sind in der Abbildung dargestellt (IDVID 680):

import numpy as np
import matplotlib.pyplot as plt

a1, b1, c1, d1, a2, b2, c2, d2 = np.linalg.solve(
    np.array([
        [-125, 25, -5, 1, 0, 0, 0, 0],
        [-8, 4, -2, 1, 0, 0, 0, 0],
        [0, 0, 0, 0, -8, 4, -2, 1],
        [0, 0, 0, 0, 0, 0, 0, 1],
        [12, -4, 1, 0, -12, 4, -1, 0],
        [-12, 2, 0, 0, 12, -2, 0, 0],
        [75, -10, 1, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 1, 0]
    ]),
    np.array([1/26, 1/5, 1/5, 1, 0, 0, 5/338, 0])
)
s1 = lambda x: a1 * x ** 3 + b1 * x ** 2 + c1 * x + d1
s2 = lambda x: a2 * x ** 3 + b2 * x ** 2 + c2 * x + d2

f = lambda x: 1 / (1 + x ** 2)
nodes = np.array([-5, -2, 0, 2, 5])

x = np.linspace(-5, 5, 300)
x1 = x[np.logical_and(-5 <= x, x <= -2)]
x2 = x[np.logical_and(-2 <= x, x <= 0)]

fig, ax = plt.subplots()

ax.plot(x, f(x), '-k', label='$f(x)$')
ax.plot(x1, s1(x1), label='$s_1(x)$')
ax.plot(x2, s2(x2), label='$s_2(x)$')
ax.plot(-x1, s1(x1), label='$s_3(x)$')
ax.plot(-x2, s2(x2), label='$s_4(x)$')
ax.plot(nodes, f(nodes), 'ok')

ax.set_xlabel('x')
ax.set_xlim(x.min(), x.max())
ax.set_ylim(-0.1, 1.1)
ax.grid()
ax.legend()
plt.show()

8.3.3Kondition und Algorithmen¶

Die Kondition für das Berechnen der Spline-Koeffizienten ist sehr gut (ohne Beweis).

Für die Berechnung der Koeffizienten muss analog zu den Interpolationspolynomen kein Gleichungssystem gelöst werden. Stattdessen kann man spezielle Basis-Splines einführen, die dann das effiziente und stabile Aufstellen des gesuchten Splines ermöglichen (vgl. auch Langrange’sche und Newton’sche Basispolynome). Mit dieser speziellen Basis dargestellte Splines heißen auch B-Splines (kurz für Basis-Splines).

8.3.4Kurven und Flächen¶

Kurven in der Ebene oder im Raum sind Funktionen, die ein Intervall $[a,b]$ nach $\bbR^2$ oder $\bbR^3$ abbilden. Ist jede Koordinatenfunktion ein Spline, so erhält man leicht zu handhabende Kurven für die grafische Modellierung. Die Koeffizienten bzgl. der B-Spline-Basis können in Grafiksoftware üblicherweise über sogenannte Kontrollpunkte verändert werden.

Ein alternative Basis zur B-Spline-Basis sind die Bernstein-Polynome, die dann sogenannte Bézier-Kurven liefern. Diese können ebenfalls über Kontrollpunkte leicht und intuitiv angepasst werden.

NURBS (Non-Uniform Rational B-Splines) sind eine Verallgemeinerung der Splines, welche auch Quotienten von Polynomen zulassen. Neben Kurven können NURBS auch nahezu beliebige Flächen beschreiben, sodass sie heute das Standardwerkzeug für die Freiformmodellierung sind.