10.8 Modifikation der Darstellung kubischer Splines

Wir wollen hier die kubischen Splines aus Abschnitt 10.7 mit einer Herleitung versehen, die ohne Hermite–Interpolation auskommt. Wir gehen für Interpolationsdaten f_{0}, …, f_{n} in angeordneten Knoten

a ≤ x_{0} < x_{1} < … < x_{n} ≤ b

von den in (10.7.7) benutzten Bezeichnungen h_{j} := x_{j} − x_{j−1} und M_{j} := s′′(x_{j}) aus und schreiben die zweite Ableitung der Interpolierenden s, die ja stückweise linear sein muß, wie in (10.7.8) im Intervall I_{j} := [x_{j−1}, x_{j}] als

s′′(x) = \frac{1}{h_{j}} (M_{j}(x − x_{j−1}) + M_{j−1}(x_{j} − x)).

Zweimalige Integration liefert dann (10.7.9) mit geeigneten Integrationskonstanten a_{j} und b_{j}:

s(x) = \frac{1}{6h_{j}} (M_{j}(x − x_{j−1})^{3} + M_{j−1}(x_{j} − x)^{3}) + b_{j} (x − \frac{x_{j} + x_{j−1}}{2} ) + a_{j}.

Bis hierhin folgen wir dem Buchtext, aber jetzt bilden wir auch die erste Ableitung und setzen Werte ein. Das ergibt

s′(x) = \frac{1}{2h_{j}} (M_{j}(x − x_{j−1})^{2} − M_{j−1}(x_{j} − x)^{2}) + b_{j}

und die Werte

\begin{aligned}s′(x_{j}) & = & \frac{1}{2}M_{j}h_{j} & + & b_{j} \\ s′(x_{j−1}) & = & −\frac{1}{2}M_{j−1}h_{j} & + & b_{j} \\ s(x_{j}) & = & \frac{1}{6}M_{j}h_{j}^{2} & + & a_{j} + b_{j} \frac{h_{j}}{2} \\ s(x_{j−1}) & = & \frac{1}{6}M_{j−1}h_{j}^{2} & + & a_{j} − b_{j} \frac{h_{j}}{2} . \end{aligned}

Das liefert dann eine Darstellung des ersten Differenzenquotienten von s als

Δ^{1}(x_{j−1}, x_{j})s = \frac{s(x_{j}) − s(x_{j−1})}{h_{j}} = \frac{h_{j}}{6} (M_{j} − M_{j−1}) + b_{j}.

Im Prinzip will man nun die Stetigkeit der ersten Ableitung in x_{j} erzwingen, indem man s′′(x_{j}) einmal in [x_{j−1}, x_{j}] und einmal in [x_{j}, x_{j+1}] darstellt und gleichsetzt. Das ergibt zunächst

\begin{aligned}s′(x_{j}) & = & \frac{1}{2}M_{j}h_{j} & + & b_{j} \\ s′(x_{j}) & = & −\frac{1}{2}M_{j+1}h_{j+1} & + & b_{j+1}, \end{aligned}

wobei aber das Auftreten von b_{j} und b_{j+1} stört. Man benutzt den Differenzenquotienten, um diese Variablen zu eliminieren, hat aber darauf zu achten, sich auf das richtige Grundintervall [x_{j−1}, x_{j}] oder [x_{j}, x_{j+1}] zu beziehen. Man erhält

\begin{aligned}s′(x_{j}) & = & \frac{1}{2}M_{j}h_{j} & + & \frac{s(x_{j}) − s(x_{j−1})}{h_{j}} − \frac{h_{j}}{6} (M_{j} − M_{j−1}) \\ s′(x_{j}) & = & −\frac{1}{2}M_{j+1}h_{j+1} & + & \frac{s(x_{j+1}) − s(x_{j})}{h_{j+1}} − \frac{h_{j+1}}{6} (M_{j+1} − M_{j}), \end{aligned}

und damit folgen nun die Gleichungen (10.7.11) bzw. (10.7.12) durch Umsortieren:

Δ^{1}(x_{j}, x_{j+1})s − Δ^{1}(x_{j−1}, x_{j})s = \frac{h_{j+1}}{6} M_{j+1} + \frac{h_{j}}{3} M_{j} + \frac{h_{j+1}}{3} M_{j} + \frac{h_{j}}{6} M_{j−1}
Δ^{2}(x_{j−1}, x_{j}, x_{j+1})s = \frac{h_{j+1}}{6(h_{j} + h_{j+1})}M_{j+1} + \frac{h_{j}}{3(h_{j} + h_{j+1})}M_{j}

+\frac{h_{j+1}}{3(h_{j} + h_{j+1})}M_{j} + \frac{h_{j}}{6(h_{j} + h_{j+1})}M_{j−1}

3Δ^{2}(x_{j−1}, x_{j}, x_{j+1})s = \frac{h_{j+1}}{2(h_{j} + h_{j+1})}M_{j+1} + M_{j} + \frac{h_{j}}{2(h_{j} + h_{j+1})}M_{j−1}
3Δ^{2}(x_{j−1}, x_{j}, x_{j+1})s = μ_{j}M_{j−1} + M_{j} + λ_{j}M_{j+1}

mit den Konstanten

λ_{j} = \frac{h_{j+1}}{2(h_{j} + h_{j+1})}, μ_{j} = \frac{h_{j}}{2(h_{j} + h_{j+1})}, λ_{j} + μ_{j} = \frac{1}{2}.

Sehen wir uns nun an, wie man (abgesehen von den Randbedingungen, um die wir uns später kümmern) das tridiagonale System lösen kann. Zwei Gleichungen der obigen Form sind, nachdem schon in der j-ten Zeile eliminiert wurde, in leicht abgekürzter Schreibweise darstellbar als

\begin{aligned}0M_{j−1} & + & r_{j}M_{j} & + & λ_{j}M_{j+1} & & & = & d_{j} \\ & & μ_{j+1}M_{j} & + & M_{j+1} & + & λ_{j+1}M_{j+2} & = & δ_{j+1}, \end{aligned}

wobei griechische Buchstaben für die ursprünglichen, die lateinischen Buchstaben für bereits veränderte Werte stehen. Die Gauss–Elimination erfordert dann einfach nur

\begin{aligned}r_{j+1} & = & 1 − \frac{μ_{j+1}}{r_{j}} λ_{j} \\ d_{j+1} & = & δ_{j+1} − \frac{μ_{j+1}}{r_{j}} d_{j} \end{aligned}

für steigende j. Dabei ändern sich die Werte λ_{j} in der oberen Nebendiagonale nicht, und man kann die Rückwärtseinsetzung in der Form

M_{j} = (d_{j} − λ_{j}M_{j+1})/r_{j}

für fallende j ausführen. Beginn und Ende dieser Rekursionen hängen aber noch von den Randbedingungen ab.

Ist das Gleichungssystem gelöst, so kann man aus (10.7.9) die Funktionswerte des kubischen interpolierenden Splines berechnen. Voraussetzung ist aber, daß man die Konstanten a_{j} und b_{j} kennt. Der Differenzenquotient liefert nach Einsetzung der Interpolationsdaten

b_{j} = \frac{f_{j} − f_{j−1}}{h_{j}} − \frac{h_{j}}{6} (M_{j} − M_{j−1}),

während man a_{j} aus zwei Mittelwerten

a_{j} = \frac{f_{j} + f_{j−1}}{2} − \frac{h_{j}^{2}}{6} \frac{M_{j} + M_{j−1}}{2}

erhält. Im Prinzip ist hiermit klar, wie die Werte von s zu berechnen sind, aber man kann natürlich noch eine Umformung der Gleichungen vornehmen, um die Vorteile des Hornerschemas oder der Newton–Darstellung auszunutzen.

Die obige Darstellung zeigt, daß ein kubisches Polynomstück s auf [x_{j−1}, x_{j}] durch die vier Werte f_{j−1} = s(x_{j−1}), f_{j} = s(x_{j}), M_{j−1} = s′′(x_{j−1}), M_{j} = s′′(x_{j}) eindeutig bestimmt ist. Wenn man im Zwischenpunkt x̃_{j} = (x_{j} + x_{j−1})/2 die Werte von s und s′′ durch die einfachen Formeln

\begin{aligned}s′′(x̃_{j}) & = & \frac{M_{j} + M_{j−1}}{2} \\ s(x̃_{j}) & = & \frac{f_{j} + f_{j−1}}{2} − \frac{h_{j}^{2}}{8} s′′(x̃_{j}) \end{aligned}

berechnet, kann man eine sukzessive Intervallhalbierung durchführen und sich bei nicht allzu großen Genauigkeitsansprüchen mit einigen Halbierungsschritten begnügen und das Ergebnis einfach durch einen Polygonzug auf den verkleinerten Intervallen darstellen. Diese Technik heißt Subdivision und wird im Computer–Aided Design oft angewendet (vgl. Kapitel 13).