Kurs:Numerik I/7 Splines

7.1 Einleitung

Bei der Polynominterpolation stellt sich mit wachsender Stützstellenzahl häufig ein oszillierendes Verhalten der Polynome ein und kann nach dem Satz von Faber im Fall, dass die Stützwerte von einer gegebenen Funktion herrühren, nicht notwendig mit gleichmäßiger Konvergenz der Interpolationspolynome bei feiner werdenden Gittern gerechnet werden. Diese Tatsachen motivieren die Einführung der in diesem Abschnitt betrachteten Splinefunktionen zur Interpolation, welche in vielen mathematischen Bereichen Anwendung finden. Für deren Definition sei

(7.1) $\Delta :=\{x_{i}|a:=x_{0}<x_{1}<\ldots <x_{m}=:b\}$

eine fest gewählte Zerlegung des Intervalls $[a,b]$. Ihre Elemente $x_{k}$ bezeichnet man im Zusammenhang mit Splines (aus historischen Gründen) meist als Knoten.

Definition 7.1

Eine Spline-Funktion oder kurz ein Spline der Ordnung $\ell \in \mathbb {N}$ zur Zerlegung $\Delta$ von $[a,b]$ mit Knoten $x_{k}$ ist eine Funktion $s\in C^{\ell -1}[a,b]$, die auf jedem Intervall $[x_{k},x_{k+1}]$ der Zerlegung mit einem Polynom $\ell$-ten Grades übereinstimmt. Den Raum solcher Splines bezeichnet man mit

$S_{\Delta ,\ell }:=\left\{s\in C^{\ell -1}[a,b]{\Big |}s{\Bigl |}_{[x_{k},x_{k+1}]}=p_{k}{\Bigr |}_{[x_{k},x_{k+1}]}\ f{\ddot {u}}r\ ein\ p_{k}\in \Pi _{\ell }\ (k=0,\ldots ,m-1)\right\}.$

Man spricht von einem interpolierenden Spline $s\in S_{\Delta ,\ell }$ mit (gegebenen) Stützwerten $y_{k}$ $(k=0,1,\ldots ,m)$, wenn gilt:

(7.2) $s(x_{k})=y_{k},\quad k=0,1,\ldots ,m.$

Ein Spline $s\in S_{\Delta ,\ell }$ ist also durch $m$ Polynome

$p_{k}\in \Pi _{\ell },\quad k=0,1,\ldots ,m-1$

gegeben, wobei $p_{k}$ nur auf dem Teilintervall $[x_{k},x_{k+1}]$ der Zerlegung $\Delta$ betrachtet wird. In den inneren Knoten $x_{k+1}$ $(k=0,\ldots ,m-2)$ von $\Delta$ gelten wegen $s\in C^{\ell -1}[a,b]$ die Glattheitsbedingungen

(7.3) ${\begin{matrix}p_{k}(x_{k+1})=p_{k+1}(x_{k+1}),&k=0,\ldots ,m-2,\\p'_{k}(x_{k+1})=p'_{k+1}(x_{k+1}),&k=0,\ldots ,m-2,\\\vdots &\vdots \\p_{k}^{(\ell -1)}(x_{k+1})=p_{k+1}^{(\ell -1)}(x_{k+1}),&k=0,\ldots ,m-2.\end{matrix}}$

Dies sind insgesamt $\ell (m-1)$ Bedingungen. Für einen interpolierenden Spline $s\in S_{\Delta ,\ell }$ mit Stützwerten $y_{k}$ $(k=0,1,\ldots ,m)$ kommen gemäß (7.2) dazu noch die $m+1$ Interpolationsbedingungen

(7.4) $p_{k}(x_{k})=y_{k}$ $(k=0,\ldots ,m-1),\quad p_{m-1}(x_{m})=y_{m}.$

Für die Konstruktion eines interpolierenden Splines sind also insgesamt

(7.5) $\ell (m-1)+(m+1)$

Glattheits- und Interpolationsbedingungen zu erfüllen.

Offenbar ist $S_{\Delta ,\ell }$ mit den üblichen Verknüpfungen ein linearer Vektorraum. Dieser enthält alle Polynome vom Grad $\leq \ell$ sowie die abgebrochenen Potenzen vom Grad $\ell$

$(x-x_{k})_{+}^{\ell }:={\begin{cases}(x-x_{k})^{\ell },&{\text{falls }}x\geq x_{k},\\0,&{\text{falls }}x<x_{k},\end{cases}}$

für $k=1,\ldots ,m-1$. Denn wegen

(7.6) ${\frac {d^{\ell -1}}{dx^{\ell -1}}}(x-x_{k})_{+}^{\ell }={\begin{cases}\ell !(x-x_{k}),&{\text{falls }}x\geq x_{k},\\0,&{\text{falls }}x<x_{k},\end{cases}}$

sind letztere Funktionen insbesondere $(\ell -1)$-mal stetig auf $[a,b]$ differenzierbar. Der folgende Satz besagt, dass die abgebrochenen Potenzen vom Grad $\ell$ zusammen mit den Monomen $1,x,\ldots ,x^{\ell }$ eine Basis des Spline-Raumes $S_{\Delta ,\ell }$ bilden.

Satz 7.2

Es gilt

(7.7) $S_{\Delta ,\ell }=\operatorname {span} \left\{1,x,\ldots ,x^{\ell },(x-x_{1})_{+}^{\ell },\ldots ,(x-x_{m-1})_{+}^{\ell }\right\}$

sowie

$\dim(S_{\Delta ,\ell })=m+\ell .$

Beweis.

Zur Konstruktion eines Splines $s\in S_{\Delta ,\ell }$ hat man höchstens $m+\ell$ Freiheitsgrade. Denn auf dem Intervall $[x_{0},x_{1}]$ kann man jedes Polynom vom Grad $\leq \ell$, also $\ell +1$ Koeffizienten bzw. Parameter frei wählen. Die zu den folgenden $m-1$ Intervallen $[x_{1},x_{2}],\ldots ,[x_{m-1},x_{m}]$ gehörenden Polynome $p_{k}$ haben insgesamt $(m-1)(\ell +1)$ Koeffizienten, von denen aber $\ell (m-1)$ durch die Glattheitsforderungen (7.3) festgelegt sind, so dass man höchstens

$(m-1)(\ell +1)-\ell (m-1)=m-1$

weitere Freiheitsgrade hat. Daher ist $\dim(S_{\Delta ,\ell })\leq m+\ell$ und es bleibt zu zeigen, dass die $m+\ell$ Funktionen in (7.7) linear unabhängig sind.

Sei dazu

$s(x):=\sum _{j=0}^{\ell }a_{j}x_{j}+\sum _{j=1}^{m-1}b_{j}(x-x_{j})_{+}^{\ell }=0,\quad x\in [a,b].$

Für $k=1,\ldots ,m-1$ schließt man aus (7.6)

${\frac {d^{\ell }}{dx^{\ell }}}(x-x_{k})_{+}^{\ell }={\begin{cases}\ell !,&{\text{falls }}x\geq x_{k},\\0,&{\text{falls }}x<x_{k},\end{cases}}$

und daher

$\lim _{x\to x_{k}+}{\frac {d^{\ell }}{dx^{\ell }}}(x-x_{k})_{+}^{\ell }=\ell !,\quad \lim _{x\to x_{k}-}{\frac {d^{\ell }}{dx^{\ell }}}(x-x_{k})_{+}^{\ell }=0.$

Wenden wir für $k=1,\ldots ,m-1$ die linearen Funktionale

$G_{k}(f):={\frac {1}{\ell !}}\left(\lim _{x\to x_{k}+}f^{(\ell )}(x)-\lim _{x\to x_{k}-}f^{(\ell )}(x)\right)$

auf $s$ an, so folgt demzufolge mit dem Kroneckersymbol $\delta _{kj}$

$0=G_{k}(s)=\sum _{j=0}^{\ell }a_{j}\underbrace {G_{k}(x^{j})} _{=0}+\sum _{j=1}^{m-1}b_{j}\underbrace {G_{k}\left((x-x_{j})_{+}^{\ell }\right)} _{=\delta _{kj}}=b_{k}.$

Also gilt

$s(x)=\sum _{j=0}^{\ell }a_{j}x^{j}=0,\quad x\in [a,b],$

was auch $a_{j}=0$ $(j=0,1,\ldots ,\ell )$ impliziert.

q.e.d.

Die in (7.7) angegebene Basis von $S_{\Delta ,\ell }$ ist für praktische Zwecke jedoch nicht geeignet. So erweist es sich als ungünstig, dass die Träger der Funktionen, welche die Basis bilden, d. h. die Bereiche, auf denen diese Funktionen verschieden von Null sind, nicht endlich sind. Ferner sind die Monome für große $j$ sowie die abgebrochenen Potenzfunktionen für dicht beieinander liegende Knoten nahezu linear abhängig, was die Auswertung eines Splines in einem Punkt mittels dieser Basis zu einer numerisch schlecht konditionierten Aufgabe macht. Für die Herleitung einer numerisch günstigeren Basis von $S_{\Delta ,\ell }$, welche aus Funktionen mit kompaktem Träger, d. h. Funktionen, die außerhalb eines abgeschlossenen Intervalls identisch Null sind, gebildet wird, sei auf Deuflhard/Hohmann, S. 245 ff., verwiesen.

Im Folgenden werden wir für interpolierende Splines der Ordnung $\ell =1$ (lineare Splines) und der Ordnung $\ell =3$ (kubische Splines) Algorithmen zu ihrer Berechnung sowie Fehlerabschätzungen angeben. Splines der Ordnung $\ell =2$ (quadratische Splines) spielen in der Praxis eine geringere Rolle und werden daher hier nicht behandelt.

7.2 Interpolierende lineare Splines

Wir wollen uns zunächst mit interpolierenden linearen Splines $s\in S_{\Delta ,1}$ für gegebene Knoten $x_{k}$ und Stützwerte $y_{k}$ $(k=0,1,\ldots ,m)$ beschäftigen. Für jedes $k\in \{0,1,\ldots ,m-1\}$ besitzt ein solcher Spline auf dem Intervall $[x_{k},x_{k+1}]$ mit einem Polynom $p_{k}\in \Pi _{1}$ offenbar die Darstellung

(7.8) $s(x)=p_{k}(x):=a_{k}+b_{k}(x-x_{k}),\quad x\in [x_{k},x_{k+1}],$

wobei die Glattheitsbedingungen (7.3)

$a_{k}+b_{k}(x_{k+1}-x_{k})=a_{k+1}+b_{k+1}(x_{k+1}-x_{k+1}),\quad k=0,1,\ldots ,m-2$

bzw.

$b_{k}={\frac {a_{k+1}-a_{k}}{x_{k+1}-x_{k}}},\quad k=0,1,\ldots ,m-2$

und die Interpolationsbedingungen (7.4)

$a_{k}+b_{k}(x_{k}-x_{k})=y_{k},\quad k=0,1,\ldots ,m-1,\quad a_{m-1}+b_{m-1}(x_{m}-x_{m-1})=y_{m}$

zu erfüllen sind. Die Glattheits- und Interpolationsbedingungen legen also die Koeffizienten in dem allgemeinen Ansatz (7.8) in eindeutiger Weise durch

(7.9) $a_{k}=y_{k},\quad b_{k}={\frac {y_{k+1}-y_{k}}{x_{k+1}-x_{k}}},\quad (k=0,1,\ldots ,m-2)$

fest und liefern den interpolierenden linearen Spline. Wir haben also gezeigt:

Satz 7.3

Zu Knoten $x_{k}$ und Stützwerten $y_{k}$ $(k=0,1,\ldots ,m)$ gibt es genau einen interpolierenden linearen Spline $s\in S_{\Delta ,1}$. Er besitzt auf dem Intervall $[x_{k},x_{k+1}]$ die Darstellung (7.8) mit Koeffizienten (7.9).

Sind die Stützwerte $y_{k}$ Funktionswerte einer gegebenen Funktion $f\in C[a,b]$, fordert man also

$s(x_{k})=f(x_{k}),\quad k=0,1,\ldots ,m,$

so kann man für den Fehler bei der Spline-Interpolation Folgendes schließen, wobei

(7.10) $\|u\|_{\infty }:=\max _{x\in C[a,b]}|u(x)|,\quad u\in C[a,b]$

die Maximumnorm auf $C[a,b]$ bezeichne.

Satz 7.4

Zu einer Funktion $f\in C^{2}[a,b]$, Knoten $x_{k}$ und Stützwerten $y_{k}:=f(x_{k})$ $(k=0,1,\ldots ,m)$ sei $s\in S_{\Delta ,1}$ der interpolierende lineare Spline. Mit

$h_{\max }:=\max _{k=0,1,\ldots ,m-1}(x_{k+1}-x_{k})$

gilt dann

$\|f-s\|_{\infty }\leq \left({\frac {1}{8}}\|f''\|_{\infty }\right)h_{\max }^{2}.$

Beweis.

Für jedes $k\in \{0,1,\ldots ,m-1\}$ stimmt der Spline $s$ auf dem Intervall $[x_{k},x_{k+1}]$ mit dem Polynom $p_{k}\in \Pi _{1}$ überein, welches den Interpolationsbedingungen

$p_{k}(x_{k})=f(x_{k}),\quad p_{k}(x_{k+1})=f(x_{k+1})$

genügt. Satz 6.11 über den Fehler bei der Polynominterpolation liefert daher für alle $x\in [x_{k},x_{k+1}]$ die Abschätzung

$|f(x)-s(x)|=|f(x)-p_{k}(x)|\leq \left|{\frac {(x-x_{k})(x-x_{k+1})}{2}}\right|\max _{\xi \in [x_{k},x_{k+1}]}|f''(\xi )|\leq {\frac {h_{\max }^{2}}{8}}\|f''\|_{\infty },$

wobei eingeht, dass die Funktion $|(x-x_{k})(x-x_{k+1})|$ ihr Maximum auf $[x_{k},x_{k+1}]$ bei

${\frac {x_{k}+x_{k+1}}{2}}=x_{k}+{\frac {x_{k+1}-x_{k}}{2}}$

annimmt. Damit ist die behauptete Fehlerabschätzung bewiesen.

q.e.d.

Nach Satz 7.4 hat man für $f\in C^{2}[a,b]$ die wichtige Aussage

$\|f-s\|_{\infty }={\mathcal {O}}(h_{\max }^{2})$

für den Fehler bei der Interpolation mit linearen Splines. Ist weiter für $j\in \mathbb {N} _{0}$

$\Delta _{j}:={\Bigl \{}x_{k}^{(j)}{\Big |}a:=x_{0}^{(j)}<x_{1}^{(j)}<\ldots <x_{m_{j}}^{(j)}:=b{\Bigr \}}$

mit einem $m_{j}\in \mathbb {N} _{0}$ eine Zerlegung von $[a,b]$, ist

$h_{\max }^{(j)}:=\max _{1\leq k\leq m_{j}}(x_{k}^{(j)}-x_{k-1}^{(j)})$

und $s_{j}\in S_{\Delta ,1}$ der zugehörige interpolierende lineare Spline mit Stützwerten

$s_{j}(x_{k}^{(j)})=f(x_{k}^{(j)}),\quad k=0,1,\ldots ,m_{j},$

so kann man aufgrund von Satz 7.4 für $\lim _{j\to \infty }h_{\max }^{(j)}=0$, anders als im Fall der gewöhnlichen Polynominterpolation (vgl. Satz 6.16), immer die gleichmäßige Konvergenz der $s_{j}$ gegen $f$ schließen, d. h.

$\lim _{j\to \infty }h_{\max }^{(j)}=0\Rightarrow \lim _{j\to \infty }\|f-s_{j}\|_{\infty }=0.$

7.3 Interpolierende kubische Splines

7.3.1 Minimaleigenschaften

Wir wollen als nächstes interpolierende kubische Splines und deren Berechnung ausführlicher betrachten. Wir beginnen damit, eine für die Anwendungen wichtige Minimaleigenschaft solcher Splines vorzustellen. Hierzu bezeichne im Folgenden

$\|u\|_{2}:=\left(\int \limits _{a}^{b}|u(x)|^{2}\,dx\right)^{1/2},\quad u\in C[a,b]$

die $L_{2}$-Norm auf dem $C[a,b]$.

Lemma 7.5

Zu einer Funktion $f\in C^{2}[a,b]$, Knoten $x_{k}$ und Stützwerten $y_{k}:=f(x_{k})$ $(k=0,1,\ldots ,m)$ sei $s\in S_{\Delta ,3}$ ein zugehöriger interpolierender kubischer Spline. Man hat dann

(7.11) $\|f''-s''\|_{2}^{2}=\|f''\|_{2}^{2}-\|s\|_{2}^{2}-2([f'-s']s'')(x){\big |}_{x=a}^{x=b}.$

Beweis.

Nach Definition der Norm $\|\cdot \|_{2}$ gilt

$\|f''-s''\|_{2}^{2}=\int \limits _{a}^{b}|f''(x)-s''(x)|^{2}\,dx=\|f''\|_{2}^{2}-2\int \limits _{a}^{b}(f''s'')(x)\,dx+\|s''\|_{2}^{2}$

(7.12) $=\|f''\|_{2}^{2}-2\int \limits _{a}^{b}([f'-s']s'')(x)\,dx-\|s''\|_{2}^{2},$

so dass wir uns nur noch mit dem mittleren Ausdruck in (7.12) befassen müssen. Für $k=0,1,\ldots ,m-1$ liefert für diesen zweimalige partielle Integration

$\int \limits _{x_{k}}^{x_{k+1}}([f''-s'']s'')(x)\,dx=([f'-s']s'')(x){\big |}_{x=x_{k}}^{x=x_{k+1}}-\int \limits _{x_{k}}^{x_{k+1}}([f'-s']s''')(x)\,dx$

$=([f'-s']s'')(x){\big |}_{x=x_{k}}^{x=x_{k+1}}-\underbrace {([f-s]s''')(x){\big |}_{x=x_{k}}^{x=x_{k+1}}} _{=0}+\underbrace {\int \limits _{x_{k}}^{x_{k+1}}([f-s]s^{(4)})(x)\,dx} _{=0},$

wobei der vorletzte Term aufgrund der Interpolationsforderungen $s(x_{k})=f(x_{k}),k=0,1,\ldots ,m$ identisch Null ist und das letzte Integral verschwindet, da $s^{(4)}\equiv 0$ auf den Teilintervallen $(x_{k},x_{k+1})$ gilt. Summation über $k=0,1,\ldots ,m-1$ liefert aufgrund der Stetigkeit der Funktionen $f',s',s''$ auf dem Intervall $[a,b]$ die folgende Teleskopsumme und damit die Aussage des Lemmas:

$\int \limits _{a}^{b}([f''-s'']s'')(x)\,dx=\sum _{k=0}^{m-1}\{([f'-s']s'')(x_{k+1})-([f'-s']s'')(x_{k})\}=([f'-s']s'')(b)-([f'-s']s'')(a).$

q.e.d.

Unter gewissen zusätzlichen Bedingungen vereinfacht sich die Aussage von Lemma 7.5:

Satz 7.6

Gegeben seien eine Funktion $f\in C^{2}[a,b]$, Knoten $x_{k}$ und die Stützwerte $y_{k}:=f(x_{k})$ $(k=0,1,\ldots ,m)$. Für einen zugehörigen interpolierenden kubischen Spline $s\in S_{\Delta ,3}$ gilt dann

(7.13) $\|f''-s''\|_{2}^{2}=\|f''\|_{2}^{2}-\|s''\|_{2}^{2},$

sofern eine der drei folgenden Randbedingungen erfüllt ist:

(a) $s''(a)=s''(b)=0,$

(b) $s'(a)=f'(a),\quad s'(b)=f'(b),$

(c) $s'(a)=s'(b),\quad s''(a)=s''(b)$, falls $f'(a)=f'(b)$.

Beweis.

In jedem der Fälle (a), (b) und (c) verschwindet in (7.11) der letzte Ausdruck, so dass dann die Identität (7.11) in die Identität (7.13) übergeht.

Aus Satz 7.6 schließt man:

Korollar 7.7

Sei $s\in S_{\Delta ,3}$ wie in Satz 7.6 definiert. Dann gilt

(7.14) $\|s''\|_{2}\leq \|g''\|_{2}$

für alle Funktionen $g\in C^{2}[a,b]$, welche den selben Bedingungen genügen, wie sie für $f$ in Satz 7.6 gefordert sind.

Beweis.

Die Beziehung (7.14) ergibt sich für Splines mit der Eigenschaft (a), (b) oder (c) wegen $\|g''-s''\|_{2}^{2}\geq 0$ unmittelbar aus Satz 7.6.

q.e.d.

Die Krümmung einer Kurve $g\in C^{2}[a,b]$ in der Ebene an der Stelle $x\in [a,b]$ ist durch

$\kappa (x):={\frac {g''(x)}{[1+g'(x)]^{3/2}}}$

definiert. Für kleine Auslenkungen $g'(x)$, wie sie z.B. bei einer dünnen Holzlatte auftreten, gilt näherungsweise $\kappa (x)\approx g''(x)$ und damit für die gesamte Krümmung der Kurve näherungsweise

$\int \limits _{a}^{b}[\kappa (x)]^{2}\,dx\approx \int \limits _{a}^{b}[g''(x)]^{2}\,dx=\|g''\|_{2}^{2}.$

Kubische Splines besitzen also nach dem letzten Satz unter allen Kurven in $C^{2}[a,b]$, welche gewissen Interpolations- und Randbedingungen genügen, in dem genannten genäherten Sinne minimale Krümmung. Diese Minimaleigenschaft stellt den Grund dafür dar, dass in der Praxis, wie beispielsweise bei der Konstruktion von Schiffsrümpfen oder der Festlegung von Schienenwegen, häufig kubische Splinefunktionen für die Interpolation verwendet werden. (Ein „spline“ ist ein englischer Name für eine dünne Holzlatte, die beim Zeichnen benutzt wurde.)

7.3.2 Vorüberlegungen

Wir wollen nun auf die Berechnung interpolierender kubischer Splines $s\in S_{\Delta ,3}$ mit Knoten $x_{k}$ und zugehörigen Stützpunkten $y_{k}$ $(k=0,1,\ldots ,m)$ eingehen. Da ein solcher Spline auf jedem Intervall $[x_{k},x_{k+1}],k=0,1,\ldots ,m-1$ mit einem Polynom 3. Grades $p_{k}\in \Pi _{3}$ identisch ist, gilt dort

(7.15) $s(x)=p_{k}(x):=a_{k}+b_{k}(x-x_{k})+c_{k}(x-x_{k})^{2}+d_{k}(x-x_{k})^{3},$

$s'(x)=p'_{k}(x)=b_{k}+2c_{k}(x-x_{k})+3d_{k}(x-x_{k})^{2},$

$s''(x)=p''_{k}(x)=2c_{k}+6d_{k}(x-x_{k}).$

Insgesamt ist ein (interpolierender) kubischer Spline also durch die $4m$ Koeffizienten $a_{k},b_{k},c_{k}$ und $d_{k}$ für $k=0,1,\ldots ,m-1$ bestimmt. Für deren Berechnung hat man zunächst die $3(m-1)$ Glattheitsbedingungen (7.3) und die $m+1$ Interpolationsbedingungen (7.4), also, wie schon allgemeiner in (7.5) festgestellt wurde, insgesamt $4m-2$ Gleichungen zur Verfügung. (Es deutet sich also schon an, dass man zur eindeutigen Bestimmung eines interpolierenden kubischen Splines 2 weitere Festlegungen benötigt.) Stellt man diese $4m-2$ Gleichungen auf, so kann man diese anschließend durch geschicktes Ineinandereinsetzen auf die in dem folgenden Lemma genannten $m-1$ linearen Gleichungen (7.17) reduzieren. Statt so vorzugehen, gehen wir hier von dem Ergebnis aus und zeigen wir umgekehrt, dass die Lösung der genannten $m-1$ linearen Gleichungen auf einen interpolierenden kubischen Spline führt. Dazu definieren wir

(7.16) $h_{k}:=x_{k+1}-x_{k},\quad k=0,1,\ldots ,m-1.$

Lemma 7.8

Falls $m+1$ reelle Zahlen $s''_{k}$ $(k=0,1,\ldots ,m)$ den $m-1$ gekoppelten Gleichungen

(7.17) $h_{k-1}s''_{k-1}+2(h_{k-1}+h_{k})s''_{k}+h_{k}s''_{k+1}=g_{k},\quad k=1,\ldots ,m-1$

mit rechten Seiten

(7.18) $g_{k}:=6{\frac {y_{k+1}-y_{k}}{h_{k}}}-6{\frac {y_{k}-y_{k-1}}{h_{k-1}}}$

genügen, so liefert der lokale Ansatz (7.15) mit den Setzungen

(7.19) $a_{k}:=y_{k},\quad b_{k}:={\frac {y_{k+1}-y_{k}}{h_{k}}}-{\frac {h_{k}}{6}}(s''_{k+1}+2s''_{k}),$

(7.20) $c_{k}:={\frac {s''_{k}}{2}},\quad d_{k}:={\frac {s''_{k+1}-s''_{k}}{6h_{k}}}$

für $k=0,1,\ldots ,m-1$ einen interpolierenden kubischen Spline $s\in S_{\Delta ,3}$ mit Knoten $x_{k}$ und Stützwerten $y_{k}$ $(k=0,1,\ldots ,m)$.

Beweis.

Für $p_{k}$ $(k=0,1,\ldots ,m-1)$ wie in (7.15) erhält man mit $a_{k}$ aus (7.19)

(7.21) $p_{k}(x_{k})=a_{k}=y_{k},\quad k=0,1,\ldots ,m-1.$

Weiter hat man mit (7.19) und (7.20)

$p_{k}(x_{k+1})=a_{k}+b_{k}h_{k}+c_{k}h_{k}^{2}+d_{k}h_{k}^{3}=y_{k}+(y_{k+1}-y_{k})-{\frac {h_{k}^{2}}{6}}\left(s''_{k+1}+2s''_{k}\right)+{\frac {s''_{k}}{2}}h_{k}^{2}+{\frac {s''_{k+1}-s''_{k}}{6}}h_{k}^{2}$

(7.22) $=y_{k+1},\quad k=0,1,\ldots ,m-1.$

Die Beziehungen (7.21) und (7.22) zusammen liefern die Interpolationsbedingungen (7.4) sowie die Stetigkeitsbedingungen

$p_{k}(x_{k+1})=p_{k+1}(x_{k+1}),\quad k=0,1,\ldots ,m-2.$

Weiter folgt für $k=0,1,\ldots ,m-2$ aus (7.17)

${\frac {y_{k+1}-y_{k}}{h_{k}}}={\frac {y_{k+2}-y_{k+1}}{h_{k+1}}}-{\frac {1}{6}}h_{k}s''_{k}-{\frac {1}{3}}(h_{k}+h_{k+1})s''_{k+1}-{\frac {1}{6}}h_{k+1}s''_{k+2}$

und damit wiederum unter Verwendung von (7.19) und (7.20)

$p'_{k}(x_{k+1})=b_{k}+2c_{k}h_{k}+3d_{k}h_{k}^{2}={\frac {y_{k+1}-y_{k}}{h_{k}}}-{\frac {h_{k}}{6}}(s''_{k+1}+2s''_{k})+h_{k}s''_{k}+{\frac {h_{k}}{2}}(s''_{k+1}-s''_{k})$

$={\frac {y_{k+1}-y_{k}}{h_{k}}}+{\frac {1}{3}}h_{k}s''_{k+1}+{\frac {1}{6}}h_{k}s''_{k}={\frac {y_{k+2}-y_{k+1}}{h_{k+1}}}-{\frac {h_{k+1}}{6}}(s''_{k+2}+2s''_{k+1})$

$=b_{k+1}=p'_{k+1}(x_{k+1}),\quad k=0,1,\ldots ,m-2.$

Schließlich gilt mit (7.20)

(7.23) $p''_{k}(x_{k+1})=2c_{k}+6d_{k}h_{k}=s''_{k}+6d_{k}h_{k}=s''_{k+1},\quad k=0,1,\ldots ,m-1$

und folglich

(7.24) $p''_{k}(x_{k+1})=2c_{k+1}=p''_{k+1}(x_{k+1}),\quad k=0,1,\ldots ,m-2.$

Damit sind auch die Glattheitsbedingungen (7.3) nachgewiesen und ist folglich $s$ mit der lokalen Darstellung (7.15) und Koeffizienten wie in (7.19) und (7.20) Element von $C^{2}[a,b]$.

q.e.d.

In der in Lemma 7.8 beschriebenen Situation bezeichnet man die $m+1$ reellen Zahlen $s''_{k}$ $(k=0,1,\ldots ,m)$ als Momente. Da insbesondere

$s''(x_{0})=p''_{0}(x_{0})=2c_{0}=s''_{0}$

gilt und gemäß (7.24)

$s''(x_{k})=p''_{k-1}(x_{k})=s''_{k},\quad k=1,\ldots ,m$

ist, hat man

(7.25) $s''_{k}=s''(x_{k}),\quad k=0,1,\ldots ,m$

Das heißt, dass die Momente $s''_{k}$ mit den zweiten Ableitungen des Splines $s$ in den Knoten $x_{k}$ übereinstimmen.

7.3.3 Natürliche, vollständige und periodische Splines

Lemma 7.8 zeigt, dass die Koeffizienten in der lokalen Darstellung (7.15) unmittelbar aus den $m+1$ Momenten $s''_{k}$ berechnet werden können. Diese $m+1$ Momente wiederum ergeben sich aus den $m-1$ linearen Gleichungen (7.17), so dass also noch zwei Freiheitsgrade vorliegen. Aufgrund der Bedingungen (a), (b) und (c) in Satz 7.6 bieten sich für deren Festlegung die folgenden drei Möglichkeiten an:

Natürliche Randbedingungen: $s''(x_{0})=s''(x_{m})=0,$

Vollständige Randbedingungen: $s'(x_{0}):=y'_{0},\quad s'(x_{m}):=y'_{m}$ für gegebene $y'_{0},y'_{m}$,

Periodische Randbedingungen: $s'(x_{0})=s'(x_{m}),\quad s''(x_{0})=s''(x_{m}).$

Im Folgenden wollen wir für diese drei Fälle die Gleichungen (7.17) zusammen mit den beiden zusätzlichen Randbedingungen in Matrix-Vektor-Form angeben.

Im Fall der natürlichen Randbedingungen lassen sich die Gleichungen (7.17) in folgender Form schreiben, wobei wegen $s''_{0}=s''_{m}=0$ (vgl. (7.25)) in diesem Fall $s''_{0}$ aus der ersten und $s''_{m}$ aus der letzten dieser Gleichungen gestrichen werden kann:

(7.26) ${\begin{pmatrix}2(h_{0}+h_{1})&h_{1}&0&\ldots &0\\h_{1}&2(h_{1}+h_{2})&h_{2}&\ddots &\vdots \\0&h_{2}&\ddots &\ddots &0\\\vdots &\ddots &\ddots &\ddots &h_{m-2}\\0&\ldots &0&h_{m-2}&2(h_{m-2}+h_{m-1})\end{pmatrix}}{\begin{pmatrix}s''_{1}\\\vdots \\s''_{m-1}\end{pmatrix}}={\begin{pmatrix}g_{1}\\\vdots \\g_{m-1}\end{pmatrix}}.$

Im Fall der vollständigen Randbedingungen verwenden wir die Beziehungen

(7.27) $s'(x_{0})=p'_{0}(x_{0})=b_{0}={\frac {y_{1}-y_{0}}{h_{0}}}-{\frac {h_{0}}{6}}(s''_{1}+2s''_{0})$

und

$s'(x_{m})=p'_{m-1}(x_{m})=b_{m-1}+2c_{m-1}h_{m-1}+3d_{m-1}h_{m-1}^{2}={\frac {y_{m}-y_{m-1}}{h_{m-1}}}-{\frac {h_{m-1}}{6}}(s''_{m}+2s''_{m-1})+s''_{m-1}h_{m-1}+{\frac {h_{m-1}}{2}}(s''_{m}-s''_{m-1})$

(7.28) $={\frac {y_{m}-y_{m-1}}{h_{m-1}}}+{\frac {2}{6}}h_{m-1}s''_{m}+{\frac {1}{6}}h_{m-1}s''_{m-1},$

welche die beiden folgenden zusätzlichen Gleichungen ergeben:

(7.29) $2h_{0}s''_{0}+h_{0}s''_{1}=-6s'(x_{0})+6{\frac {y_{1}-y_{0}}{h_{0}}},$

(7.30) $h_{m-1}s''_{m-1}+2h_{m-1}s''_{m}=6s'(x_{m})-6{\frac {y_{m}-y_{m-1}}{h_{m-1}}}.$

Mit den Randbedingungen $s'(x_{0})=y'_{0}$ und $y'_{m}=s'(x_{m})$ bezeichnen wir die rechten Seiten dieser Gleichungen mit

$g_{0}:=-6y'_{0}+6{\frac {y_{1}-y_{0}}{h_{0}}},\quad g_{m}:=6y'_{m}-6{\frac {y_{m}-y_{m-1}}{h_{m-1}}}.$

Fügt man damit die Gleichungen (7.29) und (7.30) an erster bzw. letzter Stelle den Gleichungen (7.17) hinzu, so gelangt man zu dem Gleichungssystem

${\begin{pmatrix}2h_{0}&h_{0}&0&\ldots &\ldots &0\\h_{0}&2(h_{0}+h_{1})&h_{1}&\ddots &&\vdots \\0&h_{1}&\ddots &\ddots &\ddots &\vdots \\\vdots &\ddots &\ddots &\ddots &\ddots &0\\\vdots &&\ddots &\ddots &2(h_{m-2}+h_{m-1})&h_{m-1}\\0&\ldots &\ldots &0&h_{m-1}&2h_{m-1}\end{pmatrix}}{\begin{pmatrix}s''_{0}\\\vdots \\s''_{m}\end{pmatrix}}={\begin{pmatrix}g_{0}\\\vdots \\g_{m}\end{pmatrix}}.$

Schließlich gewinnt man für die periodischen Randbedingungen wegen $s'(x_{0})=s'(x_{m})$ und $s''_{m}=s''_{0}$ durch Gleichsetzung von (7.27) und (7.28) die zusätzliche Gleichung

$2(h_{m-1}+h_{0})s''_{0}+h_{0}s''_{1}+h_{m-1}s''_{m-1}=6{\frac {y_{1}-y_{0}}{h_{0}}}-6{\frac {y_{m}-y_{m-1}}{h_{m-1}}}=:g_{0}.$

Zusammen mit dieser Gleichung an erster Position und Ersetzung von $s''_{m}$ durch $s''_{0}$ in der letzten der Gleichungen (7.17) gelangt man in diesem Fall zu dem linearen Gleichungssystem

${\begin{pmatrix}2(h_{m-1}+h_{0})&h_{0}&0&\ldots &0&h_{m-1}\\h_{0}&2(h_{0}+h_{1})&h_{1}&\ddots &&0\\0&h_{1}&\ddots &\ddots &\ddots &\vdots \\\vdots &\ddots &\ddots &\ddots &\ddots &0\\\vdots &&\ddots &\ddots &\ddots &h_{m-2}\\h_{m-1}&0&\ldots &0&h_{m-2}&2(h_{m-2}+h_{m-1})\end{pmatrix}}{\begin{pmatrix}s''_{0}\\\vdots \\s''_{m-1}\end{pmatrix}}={\begin{pmatrix}g_{0}\\\vdots \\g_{m-1}\end{pmatrix}}.$

Wir stellen nun weiter fest, dass die Matrizen in den obigen Gleichungssystemen, welche zur Bestimmung eines natürlichen, vollständigen und periodischen kubischen Splines gelöst werden müssen, jeweils symmetrisch und strikt diagonaldominant sind und positive Diagonalelemente besitzen. (Für die Definition einer strikt diagonaldominanten Matrix siehe Definition 3.2.) Für solche Matrizen kann man zeigen:

Lemma 7.9

Es sei $A:=(a_{ij})\in \mathbb {R} ^{s\times s}$ eine symmetrische und strikt diagonaldominante Matrix mit $a_{ii}>0$ $(i=1,\ldots ,s)$. Dann ist A positiv definit.

Beweis.

Es sei $\lambda \in \mathbb {R}$ Eigenwert und $x\in \mathbb {R} ^{s}$ Eigenvektor der symmetrischen Matrix $A$, d. h.

$Ax=\lambda x,$

wobei wir $x$ so normieren, dass $\|x\|_{\infty }:=\max _{j=1,\ldots ,s}|x_{j}|=1$ gilt. Sei nun $i\in \{1,\ldots ,s\}$ derart, dass $|x_{i}|=1$ ist. Dann hat man

$\sum _{j=1}^{s}a_{ij}x_{j}=a_{ii}x_{i}+\sum _{j=1 \atop j\neq i}^{s}a_{ij}x_{j}=\lambda x_{i}$

und damit

$|\lambda -a_{ii}|\underbrace {|x_{i}|} _{=1}=\left|\sum _{j=1 \atop j\neq i}^{s}a_{ij}x_{j}\right|\leq \sum _{j=1 \atop j\neq i}^{s}|a_{ij}|\underbrace {|x_{j}|} _{\leq 1}\leq \sum _{j=1 \atop j\neq i}^{s}|a_{ij}|.$

Demzufolge gilt

$\pm (\lambda -a_{ii})\leq \sum _{j=1 \atop j\neq i}^{s}|a_{ij}|$

bzw. wegen der strikten Diagonaldominanz von $A$ und $a_{ii}>0$

$0<a_{ii}-\sum _{j=1 \atop j\neq i}^{s}|a_{ij}|\leq \lambda \leq a_{ii}+\sum _{j=1 \atop j\neq i}^{s}|a_{ij}|.$

q.e.d.

Nach Lemma 3.20 ist eine positiv definite Matrix insbesondere regulär, so dass jedes der obigen drei hergeleiteten linearen Gleichungssysteme eine eindeutige Lösung besitzt. Somit können wir zusammen mit Lemma 7.8 schließen:

Korollar 7.10

Zu Knoten $x_{k}$ und Stützwerten $y_{k}$ $(k=0,1,\ldots ,m)$ gibt es jeweils genau einen interpolierenden kubischen Spline $s\in S_{\Delta ,3}$ mit natürlichen, (für vorgegebene Zahlen $y'_{0},y'_{m}$) vollständigen bzw. periodischen Randbedingungen.

Die Matrizen in den bei der Bestimmung des natürlichen und vollständigen Splines auftretenden Gleichungssystemen sind offenbar strikt diagonal dominante Tridiagonal-Matrizen, für die man eine LR-Zerlegung mit nur ${\mathcal {O}}(m-1)$ bzw. ${\mathcal {O}}(m+1)$ arithmetischen Rechenoperationen berechnen kann. Außerdem sind die Konditionen dieser Matrizen unproblematisch, so dass man die entsprechenden Systeme numerisch stabil lösen kann. Ferner gibt es auch für eine zyklische, positiv definite Tridiagonal-Matrix, wie sie bei einem periodischen Spline vorliegt, eine effiziente Cholesky-Zerlegung (siehe Schwarz und Werner für Details).

Schließlich sind generell für interpolierende kubische Splines bei Verwendung der Maximumnorm (7.10) die folgenden Fehlerabschätzungen gültig:

Satz 7.11

Zu einer Funktion $f\in C^{4}[a,b]$, Knoten $x_{k}$ und Stützwerten $y_{k}:=f(x_{k})$ $(k=0,1,\ldots ,m)$ sei $s\in S_{\Delta ,3}$ ein interpolierender kubischer Spline. Weiter sei

$h_{\max }:=\max _{k=0,1,\ldots ,m-1}(x_{k+1}-x_{k}),\quad h_{\min }:=\min _{k=0,1,\ldots ,m-1}(x_{k+1}-x_{k}).$

Falls mit einer Konstanten $C>0$

(7.31) $\max _{k=0,1,\ldots ,m}|f''(x_{k})-s''(x_{k})|\leq C\|f^{(4)}\|_{\infty }h_{\max }^{2}$

gilt, so folgen mit der Konstanten

(7.32) $c:={\frac {h_{\max }}{h_{\min }}}\left(C+{\frac {1}{4}}\right)$

die nachstehenden Abschätzungen:

$\|f-s\|_{\infty }\leq c\left\|f^{(4)}\right\|_{\infty }h_{\max }^{4},$

$\|f'-s'\|_{\infty }\leq c\left\|f^{(4)}\right\|_{\infty }h_{\max }^{3},$

$\|f''-s''\|_{\infty }\leq c\left\|f^{(4)}\right\|_{\infty }h_{\max }^{2},$

$\|f'''(x)-s'''(x)\|_{\infty }\leq c\left\|f^{(4)}\right\|_{\infty }h_{\max },\quad x\in [a,b]\setminus \{x_{0},x_{1},\ldots ,x_{m}\}.$

Der Satz ist z.B. bei Plato bewiesen. Eine Bedingung der Form (7.31) lässt sich gerade für den natürlichen, vollständigen und periodischen kubischen Spline verifizieren (siehe ebenfalls Plato), wobei beispielsweise für den natürlichen Spline $C:=3/4$ gezeigt werden kann. Nach Satz 7.11 hat man für $f\in C^{4}[a,b]$ somit für die drei untersuchten Splinetypen die Aussage

$\|f-s\|_{\infty }={\mathcal {O}}(h_{\max }^{4}).$

Ferner hat man damit, ähnlich wie für lineare Splines, mit den in Abschnitt 7.2 eingeführten Definitionen

$\lim _{j\to \infty }h_{\max }^{(j)}=0\Rightarrow \lim _{j\to \infty }\|f-s_{j}\|_{\infty }=0,$

wenn $s_{j}\in S_{\Delta ,3}$ hier den entsprechenden natürlichen, vollständigen oder periodischen kubischen Spline bezeichnet.

Beispiel 7.12

Gegeben seien die in der Tabelle

${\begin{array}{|c|c c c c c|}\hline k&0&1&2&3&4\\\hline x_{k}&-1&-0.5&0&0.5&1\\y_{k}&0.5&0.8&1&0.8&0.5\\\hline \end{array}}$

von der Funktion $f(x):=(1+x^{2})^{-1}$ herrührenden Stützwerte $y_{k}:=f(x_{k})$ und gesucht sei dazu der natürliche Spline. Offenbar hat man $m=4$ und $h_{k}=0.5$ $(k=0,1,2,3)$. Aus (7.18) errechnet man

$g_{1}=6{\frac {1-0.8}{0.5}}-6{\frac {0.8-0.5}{0.5}}=-1.2,$

$g_{2}=6{\frac {0.8-1}{0.5}}-6{\frac {1-0.8}{0.5}}=-4.8,$

$g_{3}=6{\frac {0.5-0.8}{0.5}}-6{\frac {0.8-1}{0.5}}=-1.2.$

Damit lautet das Gleichungssystem (7.26)

${\begin{pmatrix}2&0.5&0\\0.5&2&0.5\\0&0.5&2\end{pmatrix}}{\begin{pmatrix}s''_{1}\\s''_{2}\\s''_{3}\end{pmatrix}}={\begin{pmatrix}-1.2\\-4.8\\-1.2\end{pmatrix}}.$

Seine Lösung ist $s''_{1}=0,s''_{2}=-2.4,s''_{3}=0$, so dass wir aus (7.19) und (7.20) und mit den Randvorgaben $s''_{0}=s''_{4}=0$ des natürlichen Splines folgende Größen erhalten:

$b_{0}={\frac {0.3}{0.5}}=0.6,\quad b_{1}={\frac {0.2}{0.5}}+{\frac {0.5}{6}}2.4=0.6,$

$b_{2}={\frac {0.2}{0.5}}+{\frac {0.5}{6}}4.8=0,\quad b_{3}={\frac {-0.3}{0.5}}=-0.6$

und

$d_{0}=0,\quad d_{1}={\frac {-2.4}{3}}=-0.8,\quad d_{2}={\frac {2.4}{3}}=0.8,\quad d_{3}=0.$

Mit $a_{k}:=y_{k}$ und $c_{k}:=s''_{k}/2$ für $k=0,1,2,3$ gelangen wir somit zu der folgenden Darstellung des gesuchten Splines $s\in S_{\Delta ,3}$:

$s(x):={\begin{cases}0.5+0.6(x+1),&x\in [-1,-0.5],\\0.8+0.6(x+0.5)-0.8(x+0.5)^{3},&x\in [-0.5,0],\\1-1.2x^{2}+0.8x^{3},&x\in [0,0.5],\\0.8-0.6(x-0.5),&x\in [0.5,1]\end{cases}}.$

Maximale Approximationsfehler werden in den Intervallen $[-1,-0.5]$ und $[0.5,1]$ und zwar genau bei $\pm 0.787\,553$ angenommen und der Betrag beider Fehler ist $0.012\,756$. Es gilt hier

$\max _{x\in [-1,1]}\left|f^{(4)}(x)\right|=\max _{x\in [-1,1]}\left|{\frac {24}{(1+x^{2})^{3}}}-288{\frac {x^{2}}{(1+x^{2})^{4}}}+384{\frac {x^{4}}{(1+x^{2})^{5}}}\right|=f^{(4)}(0)=24.$

Weiter erhält man mit $C=3/4$ für $c$ in (7.32) den Wert $c=1$ und damit

$\max _{x\in [-1,1]}|f(x)-s(x)|\leq 24\cdot 0.5^{4}=1.5.$

Für praktische Zwecke sind also die Abschätzungen in Satz 7.11 häufig nicht brauchbar. Für die vorliegende Funktion $f$ hatte Runge gezeigt, dass die Maximumnormen der Interpolationsfehler im Fall der üblichen Polynominterpolation auf dem Intervall $[-5,5]$ für die äquidistanten Stützstellen $x_{k}:=-5+(10k)/n$ $(k=0,1,\ldots ,n)$ für $n\to \infty$ gegen $\infty$ streben (s. Schwarz, S. 102).