Interpolation mit Splinefunktionen

[Dipl Phys Ernst Winter]

Interpolation mit Splinefunktionen
Die Interpolation einer Funktion y = f(x), die mit n + 1 Stützstellen x0, x1,... xn definiert ist, beginnt mit wachsendem n zwischen den Stützstellen zu oszillieren (Satz von Faber).
Aus dieser Erkenntnis entwickelten sich die Verfahren zur Interpolation mit Splinefunktionen. Der Name ist aus der englischen Bezeichnung eines Zeichengerätes mit einem biegsamen Stahlband hergeleitet.
Kubische Splinefunktionen
Am häufigsten werden kubische Splinefunktionen verwendet, die mit stetiger Steigung und Krümmung durch die Stützpunkte gehen.
Durch die Stützpunkte und die an den inneren Stützpunkten zu fordernde Gleichheit von Funktionswert, 1. und 2. Ableitung stehen 4n-2 Gleichungen zur Berechnung der 4n Koeffizienten bereit. Für die 2 fehlende sind 2 Randbedingungen hinzuzufügen.

Ich habe alles notwendige in der Unit Spline zusammengefasst, sie exportiert
- eine Funktionen zur Erzeugung von Splinefunktionen und
- Funktionen zur Berechnung der Funktionswerte, der Werten der 1. und 2. Ableitung und einer Stammfunktion..

Delphi-Quellcode:

function MakeSpline(Flag, n: integer; var P: array of TP; 
                    c0, cn: extended; var K: array of TK):boolean;

erzeugt das Koeffizientenarray K einer Splinefunktion aus n Stützstellen und zwei Randbedingungen. Flag, c0 und cn legen die Randbedingungen fest:
Flag Art der Randbedingungen
0 Natürliche Spline: c0, yn'' = cn
1 Spline mit s'(0) = c0, s''(n) = cn
2 Spline mit s'(0) = c0, s'(n) = cn
3 Spline mit s'(n) = cn, s''(0) = c(0)
4 Periodische Spline
In K gibt die Funktion die Koeffizienten des Ansatzes (1) zurück. Mit MakeSpline=true wird angezeigt, dass die Koeffizienten in K gültig sind. Bei MakeSpline=false ist ein Fehler aufgetreten.

Delphi-Quellcode:

function Spl(n: integer; x: extended;  var P: array of TP;
                                       var K: array of TK): extended; 
function SplAbl1(n: integer; x: extended; var P: array of TP; 
                                          var K: array of TK): extended;
function SplAbl2(n: integer; x: extended; var P: array of TP;  
                                          var K: array of TK): extended;
function SplStamm(n: integer; x: extended; var P: array of TP; 
                          var K: array of TK; var Stamm: array of extended): extended;

berechnen Werte von Splinefunktionen bzw. ihren 1. und 2. Ableitungen. Die Spline wird mit den Arrays P der Stützstellen und K der Koeffizienten übergeben. Das aufrufenden Programms hat dafür zu sorgen, dass x im Definitionsbereich [x0, xn] liegt, wobei geringe Bereichsüberschreitungen durch Rundungsfehler zulässig sind.
Man kann eine Splinefunktion analytisch Integrieren und gewinnt so die zugehörige Stammfunktion zur Berechnung bestimmter Integrale. Dies ist ein effektives Verfahren zur numerischen Integration.

[brechi]

Und was willst du mit dem Beitrag jetzt sagen? Wenn das in die Code-LIB soll dann solltest du schon den Quellcode anhängen, bzw. den Quellcode im Post mit beifügen. Ansonsten hilt der Beitrag nix und ist nur Werbung.

Zu Splines: Der Name kommt aus dem Schiffsbau.
Von welchen Splines sprichst du? Anscheinend nicht von Beziére-Splines (die interpolieren nicht, sondern approximieren nur). Also B-Splines oder NURBS etc.? (solltest du eventl. erwähnen).
Deine Programme hab ich jetzt aber nicht getestet.

[Mithrandir]

Zitat:

Zitat von brechi

Deine Programme hab ich jetzt aber nicht getestet.

Wäre vielleicht besser gewesen, dann hättest du vielleicht gemerkt, dass es sich um selbstextrahierende Archive handelt. Erst gucken, dann pupen. ;)

[Daniel]

Zitat:

Zitat von brechi

Und was willst du mit dem Beitrag jetzt sagen? Wenn das in die Code-LIB soll dann solltest du schon den Quellcode anhängen, bzw. den Quellcode im Post mit beifügen. Ansonsten hilt der Beitrag nix und ist nur Werbung.

Offenbar bist Du mit dem falschen Fuß aufgestanden. Die Sourcen sind dabei und für diesen ätzenden Tonfalls gibt's derzeit echt keinen Grund. :roll:

[jbaertle]

Wer sich für weitere (allerdings etwas angestaubte, aber funktionsfähige) numerische Routinen interessiert, findet unter efg's Delphi Math Functions doch tatsächlich einen Link zu einer frühen Version der Numerical Recipes in Pascal (nicht Delphi).

[mschaefer]

Sowas kann man immer mal wieder gebrauchen. Hatte mal Temperaturwerte zu interpolieren.

Ein praktisches Problem habe ich bei interpolationsroutinen aber meist gehabt. Die Beispiele waren immer mit festgelegten Arrays und nie mit dynamischen. Leider weiss man oft einfach noch nicht wie lange die Routinen Messwerte sammeln soll. Da mathematische Routinen Geduld und Knobelgeist erforden ist das sicher etwas für die CodeLib !

// Grüße Martin

[brechi]

Tja lag wohl an der EXE-warum ich des nicht geladen habe :)
Wie wärs mit nem ZIP-Archiv?

[Dipl Phys Ernst Winter]

"brechi"

Zitat:

Wie wärs mit nem ZIP-Archiv?

Ich halte ein selbstentpackendes Archiv für allgemeiner.
Für ein Zip-Archiv wird ein Entpacker gebraucht. Ob den wirklich jeder hat?

[Dipl Phys Ernst Winter]

"mschaefer"

Zitat:

Leider weiss man oft einfach noch nicht wie lange die Routinen Messwerte sammeln soll.

Immer erst schön sammeln und dann Interpolieren.

NP: Messwerte sollte man meistens Approximieren statt sie zu Interpolieren!

[brechi]

Du kannst es ja zusätzlich noch als ZIP-Anhängen, dann würd ich mir das auch noch anschauen, sobla ich zu Hause bin. Interessant ist es auf jeden Fall und wenn der Source dabei sein sollte auch passend für die CodeLib.

[Fridolin Walther]

Zitat:

Zitat von Dipl Phys Ernst Winter

"brechi"

Zitat:

Wie wärs mit nem ZIP-Archiv?

Für ein Zip-Archiv wird ein Entpacker gebraucht. Ob den wirklich jeder hat?

Ähm ... ja? Ab ME gehörts zum Betriebssystemumfang ;).

[Mithrandir]

Dann solltest du vielleicht mal die Augen aufmachen, vielleicht mal die Suche bemühen und vielleicht mal mit ihm Kontakt per PN aufnehmen. Aus dem Alter, Menschen öffentlich diffamieren zu müssen, sollten die meisten hier langsam draußen sein.

;)

[brechi]

Bevor Emil noch nen Herzinfarkt bekommt entschuldige ich mich hiermit öffetnlich bei Herrn Dipl Phys Ernst Winter. Und ich glaube Daiel meinte eher du solltest mal nach meinen sonstigen Beiträgen suchen bevor du auf meinen IQ schließt bzw. mir ne PN schreiben wenn dich das hier so stört.

[Mithrandir]

Zitat:

Zitat von EWeiss

Bitte ?
Ich spreche von brechi nicht von Ernst..
Also erst lesen dann schimpfen.

Easy. ;) Ich spreche ebenfalls von brechi. ;)

@brechi: Genau das meinte ich. ;)

[sx2008]

Dann mach ich hier mal weiter mit konstruktiver Kritik:
* man müsste versuchen, alle Aufrufe von ShowMessage() in Unit Spline komplett zu eliminieren.

Statt ShowMessage() sollte man bei Fehlern eine Exception schmeisen:

Delphi-Quellcode:

function MakeSpline(Flag, n: integer; var P: array of TP; c0, cn: extended;var K: array of TK):boolean;
...
begin   
...
  if K[i].c<=0 then begin
    raise EMath.Create('Stützstellen nicht aufsteigend geordnet');
  // für die Exception-Klasse EMath müsste man die Unit Math einbinden

Dann zu den Funktionen Spl() und Co.:
Muss man wirklich das Intervall n als Parameter übergeben?
Ich denke es wäre günstig, wenn die Funktionen selbst das Interval anhand von x bestimmen.
Und da die Sützpunkte aufsteigend nach x sortiert sind, sollte hier eine Binärsuche sehr schnell das richtige Interval finden.

Und noch etwas zum Dateianhang:
Also ich habe bei selbstenpackenden Exedateien immer etwas Bauchweh, weil damit (und sei es unabsichtlich) ein Virus gestartet werden könnte.
Eine ZIP-Datei wäre da sicherer und kompakter.
Und von einem Delphiprogrammierer kann man sicher erwarten, dass er einen ZIP-Entpacker am Start hat.

[Dipl Phys Ernst Winter]

"sx2008"

Zitat:

Statt ShowMessage() sollte man bei Fehlern eine Exception schmeisen:

Vielen Dank für den Hinweis. Ich werde das einarbeiten.

Zitat:

Muss man wirklich das Intervall n als Parameter übergeben?

Ja! Zu jedem offenen Array gehört die Angabe seiner Länge.
Hüte dich davon abzuweichen! Bis dahin sollte der Anwender gültige Daten eingeschrieben haben, was dahinter folgt ist völlig unbestimmt. Deshalb funktioniert dein Vorschlag nicht: hnter den gewünschten Stützstellen könnten weitere in aufsteigender Folge stehen!

Zitat:

Und noch etwas zum Dateianhang:

Danke! Da habe ich mich schon umgestellt!

[Medium]

Zitat:

Zitat von Dipl Phys Ernst Winter

Ja! Zu jedem offenen Array gehört die Angabe seiner Länge.

Eigentlich tragen dynamische Arrays ihre Länge schon immer mit sich, und es ist an sich auch üblich in Methoden ein Array vorauszusetzen dass komplett gültige Werte beinhaltet. Wenn der Benutzder der Methode ein "zu langes" Array zu Grunde legt ist es an ihm das entsprechende Teil-Array herauszukopieren. Eine gute lösung für alle "Ideologien" wären ein paar Überladungen der Methode, so dass man es sich aussuchen kann.

[Dipl Phys Ernst Winter]

"Medium" schreibt:

Zitat:

Eigentlich tragen dynamische Arrays ihre Länge schon immer mit sich,...

Hier geht es nicht um dynamische Arrays, sondern um die im Proceduraufruf ohne Indexvereinbarung übergebenen nullbasierten 'offenen Arrays' beliebiger Größe.

[Dipl Phys Ernst Winter]

"Medium" schreibt:

Zitat:

Eigentlich tragen dynamische Arrays ihre Länge schon immer mit sich,...

Hier geht es nicht um dynamische Arrays, sondern um die im Proceduraufruf ohne Indexvereinbarung übergebenen nullbasierten 'offenen Arrays' beliebiger Größe.

[Medium]

Diese offenen Arrays als Parameter tragen dennoch ihre Länge mit sich, d.h. sie lässt sich wie von Arrays gewohnt mit Length() bestimmen.

[Luckie]

Ich musste hier leider ein paar Beiträge ausblenden. Bitte klärt eure persönlichen Differenzen an einer anderen Stelle, aber bitte nicht hier im Thread.

[alzaimar]

Preisfrage: Wie generalisiert man kubische Splines, sodaß die zwingende Forderung (für alle Stützpunkt P_i (X_i, Y_i) gilt: X_[i+1]>X_i) wegfällt?

Oder anders ausgedrückt: Kubische Splines können nur Kurven zeichnen, die von links nach rechts gehen, aber wie zeichne ich eine Ausgleichskurve mit Hilfe kubischer Splines durch die Punkte
(-1,0) - (0,-1) - (1,0) - (0,1) - (-1,0) ?

:gruebel: :mrgreen:

[Dipl Phys Ernst Winter]

"alzaimar"

Zitat:

Preisfrage: Wie generalisiert man kubische Splines, sodaß die zwingende Forderung (für alle Stützpunkt P_i (X_i, Y_i) gilt: X_[i+1]>X_i) wegfällt?

Du mußt die Kurve parametrisieren in x(t) und y(t).
Sieh dir hierzu das mitgelieferte Beispiel Profile an.

[Medium]

Alternativ in dem du das Bezugssystem geeignet drehst, oder aber Bezier-Splines (oder einer ihrer Abwandlungen) nimmst.

Deine Punkte in dieser Reihenfolge können so nicht durch eine eine einfache Funktion von X abgebildet werden, die zur Bedingung hat X->Y eindeutig abzubilden. Wikipedia dazu: "Eine Funktion f ordnet jedem Element x einer Definitionsmenge D genau ein Element y einer Zielmenge Z zu."

[R2009]

Hi Ernst,

ich versteh nicht ganz was du zu erreichen versuchst. Ohne jetzt tiefer in den mathematischen Hintergrund eindringen zu wollen,
berechnest du da Schmiegparabeln oder was machst du da exakt. Den Ausdruck spline hab ich in keinem meiner Fachbücher gefunden.
Bitte lass dich etwas mehr über den mathematischen Hintergrund aus.

Viele Grüsse

[R2009]

Hi Ernst,

brauchst nicht zu antworten habs in Wikipedia gefunden. Splines sind so etwas ähnliches wie Schmiegparabeln nur mit Bedingungen an den Stützstellen.

Viel Grüsse

[brechi]

Um wirklich jede beliebige Freiformfläche/Linie zu zeichnen solltest du NURBs verwenden, die werden von Opengl schon ohne weiteres unterstützt.