Der Hobby-Algorithmus

von Philipp Stephani

Einführung

In der zweidimensionalen Computergrafik stehen wir häufig vor dem Problem, durch eine gegebene Anzahl von Punkten eine glatte Kurve zu legen. Meist werden hierfür Splines verwendet. Dabei handelt es sich um miteinander verbundene Kurvenstücke. Besonders beliebt sind kubische Bézierkurven, da sie durch vier Kontrollpunkte mit einfacher geometrischer Bedeutung definiert sind und mit Hilfe des de-Casteljau-Algorithmus effizient und numerisch stabil berechnet werden können. In Abbildung 1 sehen wir zwei verbundene Bézierkurven, die beispielsweise Teil eines Béziersplines sein könnten. Die Bézierkurve auf der rechten Seite schneidet die beiden Knotenpunkte $Latex: z_i$ und $Latex: z_{i+1}$ . Die Tangenten an den Schnittpunkten schneiden die zugehörigen Kontrollpunkte $Latex: p_i$ und $Latex: q_i$ . Der Abstand von $Latex: p_i$ und $Latex: q_i$ von $Latex: z_i$ bzw. $Latex: z_{i+1}$ gibt grob gesprochen an, wie sehr sich die Kurve an die Tangenten anschmiegt.

Abbildung 1: Winkel und Kontrollpunkte von Bézierkurven

Manchmal wollen wir Béziersplines durch eine Anzahl Knotenpunkte zeichnen, ohne die anderen Kontrollpunkte explizit gegeben zu haben. Vielmehr soll das resultierende Spline die Knotenpunkte so natürlich wie möglich, d. h. mit möglichst geringer Gesamtkrümmung, interpolieren. Dies mathematisch zu quantifizieren, ist ziemlich kompliziert. Ich möchte hier lediglich als Ergebnis den Hobby-Algorithmus vorstellen. Dieser wurde von dem amerikanischen Mathematiker John D. Hobby im Jahre 1985 entwickelt und von Donald E. Knuth für sein Schriftdefinitionssystem METAFONT, welches in TeX und LaTeX benutzt wird, eingesetzt. Verwendung findet der Algorithmus ebenfalls in den auf METAFONT basierenden Vektorgrafiksprachen METAPOST und Asymptote.

Definitionen

Sei $Latex: n \in \mathbb N$ mit $Latex: n > 1$ . Gegeben seien $Latex: n+1$ Knotenpunkte $Latex: (x_0, y_0), \ldots, (x_n, y_n)$ , durch die wir ein kubisches Bézierspline mit durch den Hobby-Algorithmus ermittelten Kontrollpunkten legen wollen. Zur notationellen Vereinfachung definieren wir:

$Latex: \vec z_i : \! {{=}} \, (x_i, y_i)^T \text{ fuer } 0 \le i \le n$

$Latex: \delta x_i : \! {{=}} \, x_{i+1} - x_i \text{ fuer } 0 \le i < n$

$Latex: \delta y_i : \! {{=}} \, y_{i+1} - y_i \text{ fuer } 0 \le i < n$

$Latex: \vec\delta_i : \! {{=}} \, \vec z_{i+1} - \vec z_i \text{ fuer } 0 \le i < n$

$Latex: d_i : \! {{=}} \, \| \vec \delta_i \| \text{ fuer } 0 \le i < n$

Dabei bezeichnen $Latex: (\vec v)^T$ die Transposition des Vektors $Latex: \vec v$ (sodass die $Latex: \vec z_i$ Spaltenvektoren im $Latex: \mathbb R^2$ sind) und $Latex: \|\!\cdot\!\|$ die übliche euklidische Norm.
Grundsätzlich lassen sich zwei Typen von Splines unterscheiden: offene und geschlossene. Geschlossene Splines zeichnen sich gegenüber offenen dadurch aus, dass der letzte Punkt $Latex: \vec z_n$ durch ein weiteres Kurvenstück wieder mit dem ersten Punkt $Latex: \vec z_0$ verbunden ist. Geschlossene Splines heißen auch zyklisch. Bezüglich des Hobby-Algorithmus unterscheiden sich offene und geschlossene Splines in einigen wichtigen Details. Das Verhalten an den inneren Knotenpunkten $Latex: \vec z_1, \ldots, \vec z_{n-1}$ ist jedoch genau gleich, sodass wir zunächst diese Punkte betrachten und erst später zwischen den beiden Splinetypen differenzieren.

Die inneren Knotenpunkte

Zunächst führen wir für jeden inneren Knotenpunkt die in Abbildung 1 eingezeichneten Winkel ein. $Latex: \psi_i$ ist der Winkel zwischen der Verbindungsstrecke zwischen $Latex: \vec z_i$ und $Latex: \vec z_{i+1}$ sowie der Verbindungsstrecke zwischen $Latex: \vec z_{i-1}$ und $Latex: \vec z_i$ , d. h. zwischen den Vektoren $Latex: \vec\delta_i$ und $Latex: \vec\delta_{i-1}$ . Die Winkel $Latex: \psi_i$ werden dabei so gewählt, dass $Latex: -\pi < \psi_i \le \pi$ gilt. $Latex: \theta_i$ ist der Winkel zwischen $Latex: \vec\delta_i$ und der Verbindungsstrecke von $Latex: \vec z_i$ und $Latex: \vec p_i$ . Schließlich ist $Latex: \phi_i$ der Winkel zwischen der Verbindungsstrecke von $Latex: \vec z_{i-1}$ und $Latex: \vec z_i$ und der Verbindungsstrecke von $Latex: \vec z_i$ und $Latex: \vec q_{i-1}$ . Zu beachten ist, dass auf diese Weise auch $Latex: \theta_0$ und $Latex: \phi_n$ , nicht aber $Latex: \theta_n$ , $Latex: \phi_0$ , $Latex: \psi_0$ und $Latex: \psi_n$ definiert werden können.
Damit der Spline bei $Latex: \vec z_i$ stetig differenzierbar ist, müssen $Latex: \vec q_{i-1}$ , $Latex: \vec z_i$ und $Latex: \vec p_i$ auf einer Geraden liegen, und zwar derart, dass $Latex: \vec z_i$ zwischen $Latex: \vec q_{i-1}$ und $Latex: \vec p_i$ liegt. Daraus folgt die Bedingung, dass $Latex: \phi_i$ , $Latex: \theta_i$ und $Latex: \psi_i$ sich zum Vollwinkel bzw. Nullwinkel ergänzen. Diese Bedingung werden wir später benötigen.
Durch die Winkel $Latex: \phi_i$ und $Latex: \theta_i$ ist zwar die Richtung der Kontrollpunkte $Latex: \vec p_i$ und $Latex: \vec q_{i-1}$ bezüglich $Latex: \vec z_i$ gegeben, nicht aber ihr Abstand von $Latex: \vec z_i$ . Wir benötigen also pro innerem Knotenpunkt noch zwei weitere Parameter. Diese sollen mit der lokalen Spannungsenergie des Splines an den Kontrollpunkten in Verbindung stehen und werden deshalb Spannungsparameter genannt. $Latex: \tau_i$ soll hierbei den Spannungsparameter des bei $Latex: \vec z_i$ auslaufenden Kurvensegments bezeichnen, $Latex: \bar\tau_i$ entsprechend den des einlaufenden Segments. Im Hobby-Algorithmus werden die Spannungsparameter so definiert, dass sich minimale Gesamtspannungsenergie (und damit natürlichstes Aussehen) für $Latex: \tau_i = \bar\tau_i = 1$ ergibt.
Zunächst betrachten wir lediglich eine einzelne Bézierkurve zwischen den Punkten $Latex: \vec z_i$ und $Latex: \vec z_{i+1}$ . Wir können das Problem vereinfachen, indem wir ausnutzen, dass bei einer gegebenen affinen Transformation $Latex: \mathcal T$ diejenige Bézierkurve, die durch die mit $Latex: \mathcal T$ transformierten Kontrollpunkte definiert ist, identisch mit der durch $Latex: \mathcal T$ transformierten Bézierkurve ist. Mathematisch formuliert bedeutet dies, dass für die Bézierkurve $Latex: \mathfrak B(\vec a, \vec b, \vec c, \vec d)$ mit den Kontrollpunkten $Latex: \vec a$ , $Latex: \vec b$ , $Latex: \vec c$ und $Latex: \vec d$ gilt:

$Latex: \mathcal T(\mathfrak B(\vec a, \vec b, \vec c, \vec d))= \mathfrak B(\mathcal T(\vec a), \mathcal T(\vec b), \mathcal T(\vec c), \mathcal T(\vec d))$

Daraus folgern wir, dass wir eine Bézierkurve mit den Knotenpunkten $Latex: \vec{\tilde z}_i = (0, 0)^T$ und $Latex: \vec{\tilde z}_{i+1} = (1, 0)^T$ betrachten und die ermittelten Kontrollpunkte mit der affinen winkeltreuen Transformation $Latex: \mathcal T_i$ , die $Latex: \vec{\tilde z}_i$ auf $Latex: \vec z_i$ und $Latex: \vec{\tilde z}_{i+1}$ auf $Latex: \vec z_{i+1}$ abbildet, transformieren können. $Latex: \mathcal T_i$ ist eindeutig bestimmt durch

$Latex: \mathcal T_i(\vec v) = \vec z_i + \begin{pmatrix} \delta x_i & -\delta y_i \\ \delta y_i & \delta x_i \end{pmatrix} \vec v.$

Für die zwei untransformierten Kontrollpunkte $Latex: \vec{\tilde p}_i$ und $Latex: \vec{\tilde q}_i$ gilt dann:

$Latex: \vec{\tilde p}_i = \vec{\tilde z}_i + \xi_i \cdot \begin{pmatrix} \cos\theta_i \\ \sin\theta_i \end{pmatrix} = \begin{pmatrix} \xi_i \cos\theta_i \\ \xi_i \sin\theta_i \end{pmatrix},$

$Latex: \vec{\tilde q}_i = \vec{\tilde z}_{i+1} + \zeta_i \cdot \begin{pmatrix} -\cos\phi_{i+1} \\ \sin\phi_{i+1} \end{pmatrix} = \begin{pmatrix} 1 - \zeta_i \cos\phi_{i+1} \\ \zeta_i \sin\phi_{i+1} \end{pmatrix}.$

Dabei sind $Latex: \xi_i$ der Abstand zwischen $Latex: \vec z_i$ und $Latex: \vec p_i$ sowie $Latex: \zeta_i$ der Abstand zwischen $Latex: \vec z_{i+1}$ und $Latex: \vec q_i$ . Hat man $Latex: \vec{\tilde p}_i$ und $Latex: \vec{\tilde q}_i$ erst einmal gefunden, ergeben sich die gesuchten Kontrollpunkte $Latex: \vec p_i$ und $Latex: \vec q_i$ einfach durch Anwendung von $Latex: \mathcal T_i$ :

$Latex: \vec p_i = \mathcal T_i(\vec{\tilde p}_i),$

$Latex: \vec q_i = \mathcal T_i(\vec{\tilde q}_i).$

Die eigentliche Schwierigkeit am Hobby-Algorithmus ist die Ermittlung von $Latex: \xi_i$ , $Latex: \zeta_i$ , $Latex: \phi_i$ und $Latex: \theta_i$ . Hobby entwickelt in seiner Publikation zunächst eine Methode, um $Latex: \xi_i$ und $Latex: \zeta_i$ aus den Winkeln zu berechnen. Es gilt:

$Latex: \xi_i :\!= \frac{\rho(\theta_i, \phi_{i+1})}{3 \tau_i},$

$Latex: \zeta_i :\!= \frac{\sigma(\theta_i, \phi_{i+1})}{3 \bar\tau_{i+1}}.$

Entscheidend für die Qualität des Ergebnisses ist die Auswahl der Funktionen $Latex: \rho$ und $Latex: \sigma$ .
Hobby wählt:

$Latex: \rho(\theta, \phi) :\!= \frac{2 + \alpha(\theta, \phi)}{1 + (1-c) \cos\theta + c \cos\phi},$

$Latex: \sigma(\theta, \phi) :\!= \frac{2 - \alpha(\theta, \phi)}{1 + (1-c) \cos\phi + c \cos\theta}$

mit

$Latex: \alpha(\theta, \phi) :\!= a \cdot (\sin\theta - b \sin\phi) \cdot (\sin\phi - b \cos\theta) \cdot (\cos\theta - \cos\phi).$

Hierbei sind $Latex: a$ , $Latex: b$ und $Latex: c$ experimentell zu ermittelnde Konstanten; METAFONT und Asymptote wählen:

$Latex: a :\!= \sqrt2,$

$Latex: b :\!= \frac{1}{16},$

$Latex: c :\!= \frac{3 - \sqrt5}{2}.$

Damit ist das Problem zumindest für die inneren Knotenpunkte bis auf die Bestimmung der Winkel $Latex: phi_i$ und $Latex: theta_i$ gelöst. Für letzteres führt Hobby eine Krümmungsfunktion k ein:

$Latex: k(\theta, \phi, \tau, \bar\tau) :\!= \frac{2 \tau^2}{\rho(\theta, \phi)^2} \left (\frac{\sigma(\theta, \phi) \sin(\theta + \phi)}{\bar\tau} - 3 \sin\theta \right )$

Diese Funktion ist jedoch nichtlinear in den Variablen $Latex: \theta$ und $Latex: \phi$ , was auf ein schwierig zu lösendes nichtlineares Gleichungssystem für $Latex: \theta_i$ und $Latex: \phi_i$ hinausliefe. Um dies zu umgehen, wird $Latex: k$ nach $Latex: (\theta, \phi)$ entwickelt und nur der lineare Term $Latex: \hat k(\theta, \phi, \tau, \bar\tau)$ benutzt. Man erhält unter Benutzung der linearen Näherungen $Latex: \sin\theta \approx \theta$ und $Latex: \cos\theta \approx 1$ :

$Latex: \alpha(\theta, \phi) \approx a \cdot (\theta - b \phi) \cdot (\phi - b) \cdot (1 - 1) = 0$ ,

$Latex: \rho(\theta, \phi) \approx \frac{2}{1 + 1 - c + c} = 1$ ,

$Latex: \sigma(\theta, \phi) \approx \frac{2}{1 + 1 - c + c} = 1$ ,

$Latex: \hat k(\theta, \phi, \tau, \bar\tau) = 2 \tau^2 \left (\frac{\theta + \phi}{\bar\tau} - 3 \theta \right )$ .

Diese Linearisierung führt dazu, dass die Interpolation bei sehr großen Winkeln $Latex: \theta$ , $Latex: \phi$ nicht mehr optimal ist. Derartig große Winkel kommen jedoch in der Praxis relativ selten vor, sodass die linearisierte Krümmungsfunktion für unsere Zwecke eine mehr als ausreichende Näherung darstellt.
Damit die Krümmung an den Knotenpunkten übereinstimmt, muss

$Latex: \frac{k(\theta_i, \phi_{i+1}, \tau_i, \bar\tau_{i+1})}{d_i} = \frac{k(\phi_i, \theta_{i-1}, \bar\tau_i, \tau_{i-1})}{d_{i-1}}$

gelten. Aus der Stetigkeit der linearisierten Krümmungsfunktion ergibt sich wiederum:

$Latex: \frac{\hat k(\theta_i, \phi_{i+1}, \tau_i, \bar\tau_{i+1})}{d_i} = \frac{\hat k(\phi_i, \theta_{i-1}, \bar\tau_i, \tau_{i-1})}{d_{i-1}}$

$Latex: \frac{\tau_i^2}{d_i} \left (\frac{\theta_i + \phi_{i+1}}{\bar\tau_{i+1}} - 3 \theta_i \right ) = \frac{\bar\tau_i^2}{d_{i-1}} \left (\frac{\phi_i + \theta_{i-1}}{\tau_{i-1}} - 3 \phi_i \right )$

$Latex: \frac{\tau_i^2}{d_i \bar\tau_{i+1}} \phi_{i+1} + \frac{\tau_i^2}{d_i} \left (\frac{1}{\bar\tau_{i+1}} - 3 \right ) \theta_i = \frac{\bar\tau_i^2}{d_{i-1} \tau_{i-1}} \theta_{i-1} + \frac{\bar\tau_i^2}{d_{i-1}} \left (\frac{1}{\tau_{i-1}} - 3 \right ) \phi_i$

$Latex: d_{i-1} \tau_{i-1} \tau_i^2 \phi_{i+1} + d_{i-1} \tau_{i-1} \tau_i^2 \left (1 - 3 \bar\tau_{i+1} \right ) \theta_i = d_i \bar\tau_{i+1} \bar\tau_i^2 \theta_{i-1} + d_i \bar\tau_{i+1} \bar\tau_i^2 \left (1 - 3 \tau_{i-1} \right ) \phi_i$

Mit den Abkürzungen

$Latex: \lambda_i :\!= 3 \tau_i - 1$ ,

$Latex: \mu_i :\!= 3 \bar\tau_i - 1$ ,

$Latex: C_i :\!= d_{i-1} \tau_i^2 \tau_{i-1}$ ,

$Latex: A_i :\!= d_i \bar\tau_i^2 \bar\tau_{i+1}$

vereinfacht sich die letzte Gleichung zu:

$Latex: A_i \theta_{i-1} + C_i \mu_i \theta_i - A_i \lambda_i \phi_i - C_i \phi_{i+1} = 0$

Weiterhin gilt, wie oben erwähnt, aufgrund der stetigen Differenzierbarkeit an den Knotenpunkten:

$Latex: \phi_i = -\psi_i - \theta_i$

Dies alles gilt natürlich nur an den inneren Knotenpunkten. Wir haben aus der Stetigkeit der Krümmungsfunktion $Latex: n-1$ und aus der stetigen Differenzierbarkeit des Splines ebenfalls $Latex: n-1$ , insgesamt also $Latex: 2n-2$ Bedingungen erhalten. Für die $Latex: n$ Kurvensegmente benötigen wir je zwei zusätzliche Kontrollpunkte. Die Abstände $Latex: \xi_i$ und $Latex: \zeta_i$ der Kontrollpunkte von den Knotenpunkten lässt sich wie oben gezeigt aus den Winkeln $Latex: \theta_i$ und $Latex: \phi_i$ berechnen, somit benötigen wir pro Kontrollpunkt nur eine Gleichung (für den jeweils anliegenden Winkel). Dennoch fehlen immer noch zwei Bedingungen, für die wir die äußeren Knotenpunkte $Latex: \vec z_0$ und $Latex: \vec z_n$ betrachten und erstmals zwischen offenen und geschlossenen Splines unterscheiden müssen.

Offene Splines

Hier müssen wir zwei weitere Gleichungen künstlich einführen, da wir keine periodischen Randbedingungen haben. Analog zu den Spannungsparametern definiert Hobby die Wirbelparameter $Latex: chi_0$ und $Latex: \chi_n$ derart, dass sich für $Latex: \chi_0 = \chi_n = 1$ ein natürliches Aussehen ergibt. Die Wirbelparameter beeinflussen das Aussehen des Splines an den Endpunkten. Für $Latex: \vec z_0$ gibt es darüberhinaus einen weiteren Ausgangsspannungsparameter $Latex: \tau_0$ ; analog existiert für $Latex: \vec z_n$ ein Eingangsspannungsparameter $Latex: \bar\tau_n$ . Man wählt die folgenden Randbedingungen:

$Latex: \hat k(\theta_0, \phi_1, \tau_0, \bar\tau_1) = \chi_0 \,\hat k(\phi_1, \theta_0, \bar\tau_1, \tau_0) \iff \theta_0 = F_1 \phi_1$ ,

$Latex: \hat k(\phi_n, \theta_{n-1}, \bar\tau_n, \tau_{n-1}) = \chi_n \,\hat k(\theta_{n-1}, \phi_n, \tau_{n-1}, \bar\tau_n) \iff \phi_n = G_{n-1} \theta_{n-1}$ .

mit

$Latex: F_1 :\!= \frac{\tau_0^3 + \lambda_0 \chi_0 \bar\tau_1^3}{\mu_1 \tau_0^3 + \chi_0 \bar\tau_1^3}$ ,

$Latex: G_{n-1} :\!= \frac{\bar\tau_n^3 + \mu_n \chi_n \tau_{n-1}^3}{\lambda_{n-1} \bar\tau_n^3 + \chi_n \tau_{n-1}^3}$ .

Damit lassen sich $Latex: \theta_0$ und $Latex: \phi_n$ aus den anderen Gleichungen eliminieren:

$Latex: 0 = C_1 \mu_1 \theta_1 + A_1 \left (F_1 - \lambda_1 \right ) \phi_1 - C_1 \phi_2$ $Latex: = \left (A_1 \left (\lambda_1 - F_1 \right ) + C_1 \mu_1 \right ) \theta_1 - C_1 \phi_2 + A_1 \left (\lambda_1 - F_1 \right ) \psi_1$ ,

$Latex: 0 = A_{n-1} \theta_{n-2} + C_{n-1} \left (\mu_{n-1} - G_{n-1} \right ) \theta_{n-1} - A_{n-1} \lambda_{n-1} \phi_{n-1}$ $Latex: = A_{n-1} \theta_{n-2} + \left (C_{n-1} \left (\mu_{n-1} - G_{n-1} \right ) + A_{n-1} \lambda_{n-1} \right ) \theta_{n-1} + A_{n-1} \lambda_{n-1} \psi_{n-1}$ .

Für $Latex: n>2$ lässt sich $Latex: \phi_2$ ebenfalls durch $Latex: \theta_2$ ersetzen:

$Latex: \left (C_1 \mu_1 + A_1 \left (\lambda_1 - F_1 \right ) \right ) \theta_1 + C_1 \theta_2 + A_1 \left (\lambda_1 - F_1 \right ) \psi_1 + C_1 \psi_2 = 0$

Weiterhin gilt für $Latex: 1 < i < n-1$ :

$Latex: A_i \theta_{i-1} + (A_i \lambda_i + C_i \mu_i) \theta_i + C_i \theta_{i+1} + A_i \lambda_i \psi_i + C_i \psi_{i+1} = 0$ .

Definiert man nun

$Latex: F_i :\!= 0 \qquad (1 < i < n)$ ,

$Latex: G_i :\!= 0 \qquad (0 < i < n-1)$ ,

$Latex: \psi_n :\!= 0$ ,

$Latex: Q_i :\!= A_i \left (\lambda_i - F_i \right )$ ,

$Latex: R_i :\!= C_i \left (\mu_i - G_i \right )$ ,

$Latex: B_i :\!= Q_i + R_i$ ,

$Latex: D_i :\!= -Q_i \psi_i - C_i \psi_{i+1}$ ,

so gilt:

$Latex: B_1 \theta_1 + C_1 \theta_2 = D_1$ ,

$Latex: A_i \theta_{i-1} + B_i \theta_i + C_i \theta_{i+1} = D_i \qquad (i < 1 < n-1)$ ,

$Latex: A_n \theta_{n-1} + B_n \theta_n = D_n$ .

Dies ist ein ein lineares tridiagonales Gleichungssystem der Form

$Latex: \begin{pmatrix} B_1 & C_1 & & \cdots & & & \\ A_2 & B_2 & C_2 & \cdots & & & \\ & & & \ddots & & & \\ & & & \cdots & A_{n-2} & B_{n-2} & C_{n-2} \\ & & & \cdots & & A_{n-1} & B_{n-1} \end{pmatrix}$ $Latex: \begin{pmatrix} \theta_1 \\ \theta_2 \\ \vdots \\ \theta_{n-2} \\ \theta_{n-1} \end{pmatrix}$ $Latex: = \begin{pmatrix} D_1 \\ D_2 \\ \vdots \\ D_{n-2} \\ D_{n-1} \end{pmatrix}$

mit Subdiagonale $Latex: (A_i)$ , Diagonale $Latex: (B_i)$ , Superdiagonale $Latex: (C_i)$ und rechtsseitigem Vektor $Latex: (D_i)$ , welches mit dem Gaußalgorithmus in linearer Zeit gelöst werden kann.
Sind alle $Latex: \tau_i$ und $Latex: \bar\tau_i$ größer als $Latex: \frac23$ , so gilt wegen $Latex: C_i > 0$ , $Latex: A_i > 0$ , $Latex: \lambda_i > 1$ und $Latex: \mu_i > 1$ an den inneren Knotenpunkten:

An den Randpunkten folgt analog, falls zusätzlich $Latex: \lambda_1 > F_1$ und $Latex: \mu_{n-1} > G_{n-1}$ gelten:

$Latex: \left | B_1 \right | = A_1 \left (\lambda_1 - F_1 \right ) + C_1 \mu_1 > \left | C_1 \right |$ ,

$Latex: \left | B_{n-1} \right | = A_{n-1} \lambda_{n-1} + C_{n-1} \left (\mu_{n-1} - G_{n-1} \right ) > \left | A_{n-1} \right |$ .

Somit ist die Systemmatrix diagonaldominant. Das führt dazu, dass in diesem (in der Realität meistens angenommenen) Fall keine Pivotisierung nötig ist.

Geschlossene Splines

Hier haben wir periodische Randbedingungen vorliegen, weshalb es sich anbietet, die Indizes zyklisch zu definieren:

$Latex: \vec\delta_n :\!= z_0 - z_n$ ,

$Latex: d_n :\!= \| \vec\delta_n \|$ ,

$Latex: d_{-1} :\!= d_n$ ,

$Latex: d_{n+1} :\!= d_0$ ,

$Latex: \theta_{-1} :\!= \theta_n$ ,

$Latex: \theta_{n+1} :\!= \theta_0$ ,

$Latex: \phi_{-1} :\!= \phi_n$ ,

$Latex: \phi_{n+1} :\!= \phi_0$ ,

$Latex: \tau_{-1} :\!= \tau_n$ ,

$Latex: \bar\tau_{n+1} :\!= \bar\tau_0$ .

Weiterhin definieren wir $Latex: \psi_0$ und $Latex: \psi_{n+1}$ als den Winkel zwischen $Latex: \vec\delta_n$ und $Latex: \vec\delta_0$ sowie $Latex: \psi_n$ als den Winkel zwischen $Latex: \vec\delta_{n-1}$ und $Latex: \vec\delta_n$ . Damit können wir die Gleichungen Kruemmung und Differenzierbarkeit auch für $Latex: i = 0$ und $Latex: i = n$ anwenden. Mit den Definitionen

$Latex: Q_i :\!= A_i \lambda_i$ ,

$Latex: R_i :\!= C_i \mu_i$ ,

$Latex: B_i :\!= Q_i + R_i$ ,

$Latex: D_i :\!= -Q_i \psi_i - C_i \psi_{i+1}$

erhalten wir wieder

$Latex: A_i \theta_{i-1} + B_i \theta_i + C_i \theta_{i+1} = D_i.$

Für $Latex: i = 0$ und $Latex: i = n$ gelten mit obigen Definitionen die Spezialfälle:

$Latex: B_0 \theta_0 + C_0 \theta_1 + A_0 \theta_n = D_0$ ,

$Latex: C_n \theta_0 + A_n \theta_{n-1} + B_n \theta_n = D_n$ .

Aufgrund der Periodizität des Problems ist die Gleichung für $Latex: i = n+1$ identisch zu der Gleichung für $Latex: i = 0$ . Damit haben wir also für die $Latex: n+1$ Winkel $Latex: \theta_0, \ldots, \theta_n$ (ein Winkel pro Knotenpunkt) $Latex: n+1$ Gleichungen gefunden, die nun aber kein einfaches tridiagonales Gleichungssystem mehr, sondern ein zyklisches tridiagonales Gleichungssystem der Form

$Latex: \begin{pmatrix} B_0 & C_0 & & \cdots & & & A_0 \\ A_1 & B_1 & C_1 & \cdots & & & \\ & & & \ddots & & & \\ & & & \cdots & A_{n-1} & B_{n-1} & C_{n-1} \\ C_n & & & \cdots & & A_n & B_n \end{pmatrix}$ $Latex: \begin{pmatrix} \theta_0 \\ \theta_1 \\ \vdots \\ \theta_{n-1} \\ \theta_n \end{pmatrix}$ $Latex: = \begin{pmatrix} D_0 \\ D_1 \\ \vdots \\ D_{n-1} \\ D_n \end{pmatrix}$

bilden. Unter Zuhilfenahme der Sherman-Morrison-Formel kann auch dieses System mit linearer Zeitkomplexität gelöst werden.

Quellen

John D. Hobby: Smooth, Easy to Compute Interpolating Splines, 1985
ftp://db.stanford.edu/pub/cstr/reports/cs/tr/85/1047/
Martin Kerscher, Joachim Puls, Sigmung Stintzing: Numerik für Physiker, 2007
http://www.mathematik.uni-muenchen.de/~kerscher/vorlesungen/numeriksose08/doc/script.pdf
William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery:
Numerical Recipes in C++, Second Edition, Cambridge University Press

Downloads

Beispielprojekt [282559 Bytes]

Artikel als PDF [189103 Bytes]

Ihre Meinung

Falls Sie Fragen zu diesem Tutorial haben oder Ihre Erfahrung mit anderen Nutzern austauschen möchten, dann teilen Sie uns diese bitte in einem der unten vorhandenen Themen oder über einen neuen Beitrag mit. Hierzu können sie einfach einen Beitrag in einem zum Thema passenden Forum anlegen, welcher automatisch mit dieser Seite verknüpft wird.