Donnerstag, 24. Dezember 2009

Splines - das berechnete Kurvenlineal

Anlässlich meines aktuellen Projekts astropatterns war ich auf der Suche nach einem Algorithmus in C oder C++ zur Interpolation mit periodischen Splines. Splines sind stückweise kubische Funktionen, die eine gegebene Menge von Stützpunkten einer Funktion interpolieren. Sie haben dabei die Eigenschaft, dass sie unter allen hinreichend regulären interpolierenden Funktionen das Integral des Quadrats der zweiten Ableitung minimieren. Anschaulich bedeutet dies, dass sie versuchen, ihre Aufgabe des Interpolierens mit möglichst wenig Krümmung zu erledigen: Jede andere interpolierende Funktion krümmt sich gewissermassen mehr als der Spline. Das macht Splines für die Numerik attraktiver als beispielsweise das klassische Interpolationspolynom N-ter Ordnung für N+1 Stützstellen, das zwar auch an den Stützstellen die gewünschten Funktionswerte annimmt, aber bei unpassender Wahl der Stützstellen in den Bereichen zwischen ihnen stark oszillieren kann. Ein Spline ist sozusagen die rechnerische Form eines Kurvenlineals (das ist ein im Bürohandel erhältlicher biegsamer Stab aus einem speziellen Kunststoff, mit dessen Hilfe man interpolierende Funktionen mit dem Bleistift zeichnen kann).

Obwohl ich mich selbst bereits vor 23 Jahren im Studium mit diesen Gebilden befasste und seitdem sicher Tausende von Studenten im Rahmen einer Programmierübung für Praktische Mathematik einen Algorithmus zu deren Berechnung abgeliefert haben dürften, fand ich keinen freien C-, Java- oder C++-Code im Internet, mit dem periodische Splines, wie ich sie benötigte, berechnet werden konnten.

Wohl gibt es für die sogenannten "natürlichen Splines" und für Splines mit an den Rändern vorgegebener zweiter Ableitung eine 1989 von Peter und Nigel verfasste Spline-Funktion, die mit der mingw-Windows-Version des GCC 4.4.0 auf Anhieb und isoliert lauffähig ist (http://www.mech.uq.edu.au/staff/jacobs/nm_lib/cmathsrc/) und ein schönes schlankes Binary von nur 22.8 KB erzeugt. Da es reines C ist, kann ich auch den Tiny C Compiler von Fabrice Bellard verwenden und komme sogar auf ein nur 5.6 KB grosses Binary inclusive eines aufrufenden kleinen Testprogramms.

Leider sind die von Nigel und Peter abgedeckten Splines mit vorgegebenen zweiten Ableitungen an den Randstellen nicht das, was ich für meine Zwecke benötigte. Die Funktion, die ich interpolieren möchte, lebt auf dem astrologischen Häuserkreis. Es handelt sich also um eine periodische Funktion, und hierfür benötigt man den periodischen Spline, der in der Ausrechnung etwas komplizierter wird, da er nicht mehr auf ein lineares Gleichungssystem mit tridiagonaler Koeffizientenmatrix führt: Auch die Ecken oben rechts und unten links in der Matrix können einen von Null verschiedenen Wert enthalten.

Für Leute wie mich, die ihre Lehrbücher über Numerische Mathematik irgendwo im Dachboden in einer Kiste verwahrt haben, ist es angenehm, eine ausgezeichnete Praktikums- oder Seminarvorbereitung eines anonymen Regensburger Studenten im Web vorzufinden:
Das Papier lässt keine Fragen zur Theorie der Splines offen und leitet die zu lösenden linearen Gleichungssysteme für die verschiedenen Arten kubischer Splines her, wobei auch der periodische Spline ausführlich erörtert wird.

So blieb nur noch die Aufgabe, das lineare Gleichungssystem der periodischen Splines in seiner allgemeinen Form zu lösen. Hierbei half das Wiki zu Computational Fluid Dynamics, das einige sehr klar geschriebene Artikel über tridiagonale Matrizen und ihre Varianten enthält (zum Einstieg eignet sich der Artikel über den Thomas-Algorithmus).

In C++ würde man den Spline als ein Objekt betrachten, das bei der Konstruktion oder durch Aufruf einer speziellen Methode die Knoten vorgelegt bekommt und dann die Koeffizienten aller auftretenden kubischen Polynome errechnet und als globale Membervariablen aufbewahrt. Für die eigentliche Berechnung des Funktionswerts an einem beliebigen Zwischenwert bietet die Klasse dann eine spezielle Methode eval() an.

Genau so funktioniert meine Klasse PeriodicSpline.cpp, die ich hiermit der Allgemeinheit für die Berechnung von periodischen Splines zur Verfügung stelle. Sie ist unabhängig vom Projekt astropatterns in beliebigen Kontexten verwendbar.

Als Beispiel für die Verwendung möge ein Programm dienen, das den Array der Stützstellen aus einer Datei data.txt einliest, den periodischen Spline errechnet, an 100 Punkten auswertet und das Ergebnis in eine Datei curveData.txt schreibt.

#include <cstdlib>
#include <iostream>
#include <fstream>

#include "PeriodicSpline.hpp"


int main(int argc, char *argv[])
{

double xi, yi, t, period;
size_t i;
vector<double> x,y;
string msg;

// Parse the knots as pairs of doubles from an input file
ifstream dataFile( "data.txt" );
while ( dataFile >> xi >> yi ) {
x.push_back(xi);
y.push_back(yi);
}

period = x[x.size()-1] - x[0];
cout << "dimension: " << x.size() << endl;

for (i = 0; i < x.size(); i++ ) {
cout << x[i] << ":" << y[i] << endl;
}
cout << endl;

try {

// Compute the spline
PeriodicSpline s(x,y);


// Evaluate the spline values for 100 points
ofstream curveData( "curveData.txt" );
for (i=0;i<100;i++) {
t = x[0] + i * period / 100;
curveData << t << '\t'
<< s.eval(t) << endl;
}


} catch(string msg) {
cerr << msg;
exit(EXIT_FAILURE);
}
}



Anmerkung.- Die verwendeten Ein- und Ausgabestreams dieses Beispielprogramms machen den Code klar und lesbar, haben jedoch ihren Preis. Durch blosses Includieren von <iostream> und <fstream> erhöhte sich die Grösse meiner Executables von ca. 50 KB auf über 500 KB! Ich weiss nicht, ob dies an der mingw-Installation liegt oder auch auf anderen Plattformen, z.B. Linux auftritt. Das Phänomen ist bekannt und wird als unproblematisch angesehen. Mich stört es trotzdem. Ich kann mir nicht vorstellen, dass der gesamte, in das Executable gepackte Code wirklich verwendet wird und erwarte von einem intelligenten Linker, dass er nur die wirklich benötigten Codeteile in die ausführbare Datei packt.

Keine Kommentare :