Kurs:Numerik I/7 Splines

Aus Wikiversity

Wechseln zu: Navigation, Suche

Inhaltsverzeichnis

[Bearbeiten] 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 xk bezeichnet man im Zusammenhang mit Splines (aus historischen Gründen) meist als Knoten.

[Bearbeiten] Definition 7.1

Eine Spline-Funktion oder kurz ein Spline der Ordnung \ell \in \mathbb N zur Zerlegung Δ von [a,b] mit Knoten xk ist eine Funktion s \in C^{\ell - 1}[a, b], die auf jedem Intervall [xk,xk + 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 ur\ 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 yk (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 pk nur auf dem Teilintervall [xk,xk + 1] der Zerlegung Δ betrachtet wird. In den inneren Knoten xk + 1 (k = 0, \ldots, m - 2) von Δ 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^{(\ell - 1)}_k (x_{k+1}) = p^{(\ell - 1)}_{k+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 yk (k = 0, 1, \ldots, m) kommen gemäß (7.2) dazu noch die m + 1 Interpolationsbedingungen

(7.4) pk(xk) = yk (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 \le \ell sowie die abgebrochenen Potenzen vom Grad \ell

(x - x_k)^\ell_+ := \begin{cases} (x - x_k)^\ell, & \text{falls } x \ge 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 \ge 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.

[Bearbeiten] 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.

[Bearbeiten] Beweis.

Zur Konstruktion eines Splines s \in S_{\Delta, \ell} hat man höchstens m + \ell Freiheitsgrade. Denn auf dem Intervall [x0,x1] kann man jedes Polynom vom Grad \le \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 pk 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}) \le m + \ell und es bleibt zu zeigen, dass die m + \ell Funktionen in (7.7) linear unabhängig sind.

Sei dazu

s(x) := \sum^\ell_{j=0} a_j x_j + \sum^{m - 1}_{j=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 \ge 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 δkj

0 = G_k(s) = \sum^\ell_{j=0} a_j \underbrace{G_k(x^j)}_{= 0} + \sum^{m - 1}_{j=1} b_j \underbrace{G_k \left( (x - x_j)^\ell_+ \right)}_{= \delta_{kj}} = b_k.

Also gilt

s(x) = \sum^\ell_{j=0} a_j x^j = 0, \quad x \in [a, b],

was auch aj = 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.

[Bearbeiten] 7.2 Interpolierende lineare Splines

Wir wollen uns zunächst mit interpolierenden linearen Splines s \in S_{\Delta, 1} für gegebene Knoten xk und Stützwerte yk (k = 0, 1, \ldots, m) beschäftigen. Für jedes k \in \{0, 1, \ldots, m - 1\} besitzt ein solcher Spline auf dem Intervall [xk,xk + 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:

[Bearbeiten] Satz 7.3

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

Sind die Stützwerte yk 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.

[Bearbeiten] Satz 7.4

Zu einer Funktion f \in C^2[a, b], Knoten xk und Stützwerten yk: = f(xk) (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 \le \left( \frac{1}{8} \|f''\|_\infty \right) h^2_\max.

[Bearbeiten] Beweis.

Für jedes k \in \{0, 1, \ldots, m - 1\} stimmt der Spline s auf dem Intervall [xk,xk + 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)| \le \left| \frac{(x - x_k)(x - x_{k+1})}{2} \right| \max_{\xi \in [x_k, x_{k+1}]} |f''(\xi)| \le \frac{h^2_\max}{8} \|f''\|_\infty,

wobei eingeht, dass die Funktion | (xxk)(xxk + 1) | ihr Maximum auf [xk,xk + 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^2_\max)

für den Fehler bei der Interpolation mit linearen Splines. Ist weiter für j \in \mathbb N_0

\Delta_j := \Bigl\{ x^{(j)}_k \Big| a := x^{(j)}_0 < x^{(j)}_1 < \ldots < x^{(j)}_{m_j} := b \Bigr\}

mit einem m_j \in \mathbb N_0 eine Zerlegung von [a,b], ist

h^{(j)}_\max := \max_{1 \le k \le m_j} (x^{(j)}_k - x^{(j)}_{k-1})

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

s_j(x^{(j)}_k) = f(x^{(j)}_k), \quad k = 0, 1, \ldots, m_j,

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

\lim_{j \to \infty} h^{(j)}_\max = 0 \Rightarrow \lim_{j \to \infty} \|f - s_j\|_\infty = 0.

[Bearbeiten] 7.3 Interpolierende kubische Splines

[Bearbeiten] 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^b_a |u(x)|^2\, dx \right)^{1/2}, \quad u \in C[a, b]

die L2-Norm auf dem C[a,b].

[Bearbeiten] Lemma 7.5

Zu einer Funktion f \in C^2[a, b], Knoten xk und Stützwerten yk: = f(xk) (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=b}_{x=a}.

[Bearbeiten] Beweis.

Nach Definition der Norm \|\cdot\|_2 gilt

\|f'' - s''\|^2_2 = \int\limits^b_a |f''(x) - s''(x)|^2\, dx = \|f''\|^2_2 - 2 \int\limits^b_a (f''s'')(x)\, dx + \|s''\|^2_2
(7.12) = \|f''\|^2_2 - 2 \int\limits^b_a ([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+1}}_{x_k} ([f'' - s''] s'')(x)\, dx = ([f' - s'] s'')(x)\big|^{x=x_{k+1}}_{x=x_k} - \int\limits^{x_{k+1}}_{x_k} ([f' - s'] s''')(x)\, dx
= ([f' - s'] s'')(x) \big|^{x=x_{k+1}}_{x=x_k} - \underbrace{([f - s] s''')(x) \big|^{x=x_{k+1}}_{x=x_k}}_{=0} + \underbrace{\int\limits^{x_{k+1}}_{x_k} ([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 (xk,xk + 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^b_a ([f'' - s''] s'')(x)\, dx = \sum^{m-1}_{k=0} \{([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:

[Bearbeiten] Satz 7.6

Gegeben seien eine Funktion f \in C^2[a, b], Knoten xk und die Stützwerte yk: = f(xk) (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).

[Bearbeiten] 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:

[Bearbeiten] Korollar 7.7

Sei s \in S_{\Delta, 3} wie in Satz 7.6 definiert. Dann gilt
(7.14) \|s''\|_2 \le \|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.

[Bearbeiten] Beweis.

Die Beziehung (7.14) ergibt sich für Splines mit der Eigenschaft (a), (b) oder (c) wegen \|g'' - s''\|^2_2 \ge 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^b_a [\kappa(x)]^2\, dx \approx \int\limits^b_a [g''(x)]^2\, dx = \|g''\|^2_2.

Kubische Splines besitzen also nach dem letzten Satz unter allen Kurven in C2[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.)

[Bearbeiten] 7.3.2 Vorüberlegungen

Wir wollen nun auf die Berechnung interpolierender kubischer Splines s \in S_{\Delta, 3} mit Knoten xk und zugehörigen Stützpunkten yk (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) = pk(x): = ak + bk(xxk) + ck(xxk)2 + dk(xxk)3,
s'(x) = p'k(x) = bk + 2ck(xxk) + 3dk(xxk)2,
s''(x) = p''k(x) = 2ck + 6dk(xxk).

Insgesamt ist ein (interpolierender) kubischer Spline also durch die 4m Koeffizienten ak,bk,ck und dk 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.

[Bearbeiten] 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_ks''_{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 xk und Stützwerten yk (k = 0, 1, \ldots, m).

[Bearbeiten] Beweis.

Für pk (k = 0, 1, \ldots, m - 1) wie in (7.15) erhält man mit ak 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_kh_k + c_kh^2_k + d_kh^3_k = y_k + (y_{k+1} - y_k) - \frac{h^2_k}{6} \left( s''_{k+1} + 2s''_k \right) + \frac{s''_k}{2} h^2_k + \frac{s''_{k+1} - s''_k}{6} h^2_k
(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_kh_k + 3d_kh^2_k = \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 C2[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''(x0) = p''0(x0) = 2c0 = 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 xk übereinstimmen.

[Bearbeiten] 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.12 bieten sich für deren Festlegung die folgenden drei Möglichkeiten an:

Natürliche Randbedingungen: s''(x0) = s''(xm) = 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^2_{m-1} = \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_0s''_0 + h_0s''_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'(x0) = y'0 und y'm = s'(xm) 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'(x0) = s'(xm) 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_0s''_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:

[Bearbeiten] Lemma 7.9

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

[Bearbeiten] Beweis.

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

Ax = λ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 | xi | = 1 ist. Dann hat man

\sum^s_{j=1} a_{ij} x_j = a_{ii} x_i + \sum^s_{j=1 \atop j \neq i} a_{ij} x_j = \lambda x_i

und damit

|\lambda - a_{ii}| \underbrace{|x_i|}_{=1} = \left| \sum^s_{j=1 \atop j \neq i} a_{ij} x_j \right| \le \sum^s_{j=1 \atop j \neq i} |a_{ij}| \underbrace{|x_j|}_{\le 1} \le \sum^s_{j=1 \atop j \neq i} |a_{ij}|.

Demzufolge gilt

\pm (\lambda - a_{ii}) \le \sum^s_{j=1 \atop j \neq i} |a_{ij}|

bzw. wegen der strikten Diagonaldominanz von A und aii > 0

0 < a_{ii} - \sum^s_{j=1 \atop j \neq i} |a_{ij}| \le \lambda \le a_{ii} + \sum^s_{j=1 \atop j \neq i} |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:

[Bearbeiten] Korollar 7.10

Zu Knoten xk und Stützwerten yk (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:

[Bearbeiten] Satz 7.11

Zu einer Funktion f \in C^4[a, b], Knoten xk und Stützwerten yk: = f(xk) (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)| \le C \|f^{(4)}\|_\infty h^2_\max

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 \le c \left\| f^{(4)} \right\|_\infty h^4_\max,
\|f' - s'\|_\infty \le c \left\| f^{(4)} \right\|_\infty h^3_\max,
\|f'' - s''\|_\infty \le c \left\| f^{(4)} \right\|_\infty h^2_\max,
\|f'''(x) - s'''(x)\|_\infty \le 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^4_\max).

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

\lim_{j \to \infty} h^{(j)}_\max = 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.

[Bearbeiten] 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 + x2) − 1 herrührenden Stützwerte yk: = f(xk) und gesucht sei dazu der natürliche Spline. Offenbar hat man m = 4 und hk = 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 ak: = yk und ck: = 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)| \le 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 xk: = − 5 + (10k) / n (k = 0, 1, \ldots, n) für n \to \infty gegen \infty streben (s. Schwarz, S. 102).

Persönliche Werkzeuge