Ausgleichskurve.pm

Die Intervall-Zerlegung und die Kurve

Ich betrachte die vier Polynome P(0), P(1), P(2) und P(3):

P(0)(X) = 2 X**3 - 3 X**2 + 1
P(1)(X) = X**3 - 2 X**2 + X
P(2)(X) = -2 X**3 + 3 X**2
P(3)(X) = X**3 - X**2

P(0)'(X) = 6 X**2 - 6 X
P(1)'(X) = 3 X**2 - 4 X + 1
P(2)'(X) = -6 X**3 + 6 X
P(3)'(X) = 3 X**2 - 2 X

P(0)(0) = 1
P(1)(0) = 0
P(2)(0) = 0
P(3)(0) = 0

P(0)'(0) = 0
P(1)'(0) = 1
P(2)'(0) = 0
P(3)'(0) = 0

P(0)(1) = 0
P(1)(1) = 0
P(2)(1) = 1
P(3)(1) = 0

P(0)'(1) = 0
P(1)'(1) = 0
P(2)'(1) = 0
P(3)'(1) = 1

Es sei eine wachsende Folge von Zahlen a(0), a(1), ... a(n) gegeben, d(i) := a(i+1) - a(i) für 0 <= i < n. Ich definiere die folgenden Funktionen f(i) für 0 <= i <= 2n + 1:

  f(0)    : x -> 1                                   für         x < a(0)
            x -> P(0)((x - a(0))/d(0))               für a(0) <= x < a(1)
            x -> 0                                   für a(1) <= x

  f(1)    : x -> x - a(0)                            für         x < a(0)
            x -> d(0) * P(1)(x - a(0))/d(0))         für a(0) <= x < a(1)
            x -> 0                                   für a(1) <= x

  Für 1 <= k < n-1:

  f(2k)   : x -> 0                                   für           x < a(k-1)
            x -> P(2)((x - a(k-1))/d(k-1))           für a(k-1) <= x < a(k)
            x -> P(0)((x - a(k))  /d(k)  )           für a(k)   <= x < a(k+1)
            x -> 0                                   für a(k+1) <= x

  f(2k+1) : x -> 0                                   für           x < a(k-1)
            x -> d(k-1) * P(3)((x - a(k-1))/d(k-1))  für a(k-1) <= x < a(k)
            x -> d(k)   * P(1)((x - a(k))  /d(k)  )  für a(k)   <= x < a(k+1)
            x -> 0                                   für a(k+1) <= x

  f(2n)   : x -> 0                                   für           x < a(n-1)
            x -> P(2)((x - a(n-1))/d(n-1))           für a(n-1) <= x < a(n)
            x -> 1                                   für a(n)   <= x

  f(2n+1) : x -> 0                                   für           x < a(n-1)
            x -> d(n-1) * P(3)((x - a(n-1))/d(n-1))  für a(n-1) <= x < a(n)
            x -> x - a(n)                            für a(n)   <= x

Die Funktionen f(i) sind stetig differenzierbar, und es gilt

f(2k)(a(k)) = 1
f(2k)(a(j)) = 0   für j != k
f'(2k)(a(j) = 0

f(2k+1)(a(j))  = 0
f'(2k+1)(a(k)) = 1
f'(2k+1)(a(j)) = 0 für j != k

Gesucht ist eine Funktion F = ∑b(j)*f(j), die in einem bestimmten Sinne optimal ist. Das Maß der Abweichung vom Optimum („Fehler”) ist eine Funktion E : F -> Q(F) - 2 * L(F) + C.

Q ist eine Funktion F -> H(F,F) mit einer bilinearen Abbildung H. von der hier nur der symmetrische Teil relevant ist. Ich kann H also als symmetrisch annehmen. In der Praxis ist H auch positiv definit.

L ist eine lineare Abbildung.

Ich betrachte die Abbildungen Q und L auf Linearkombinationen der Funktionen f(i), 0 <= i < 2n + 2. Die Abbildung Q kann auf diesem Bereich durch eine symmetrische Matrix A mit 2n + 2 Zeilen und Spalten dargestellt werden. Die Abbildung L entspricht einem Vektor B mit 2n + 2 Komponenten, ebenso natürlich die Linearkombinationen der Funktionen f(i).

Der „Fehler” von ∑b(j)*f(j) ist ∑b(j)*b(k)*A(j,k) - 2 * ∑b(j)*B(j) + C. Falls die Matrix A positiv-definit ist, ist der „Fehler” minimal, wenn ∑A(j,k)*b(k) = B(j)

In der Praxis ist die „Fehlerfunktion” E die Summe mehrerer Funktionen E(i) : F -> H(i)(F,F) - 2 * L(i)(F) + C(i). Dementsprechend ist A die Summme mehrerer positiv-semidefiniter, symmetrischer Matrizen A(i) und B die Summe mehrerer Vektoren B(i).

Der Konstruktor new legt die Zerlegung durch die Punkte a(i) fest und initialisiert die Matrix A und den Vektor B.

Die Funktion add_points fügt einen „Teilfehler” zu A und B hinzu. Der „Teilfehler” beschreibt, wie weit die Ausgleichsfunktion sich den Punkten annähert.

Die Funktion get_cbezier liefert die „Ausgleichskurve” ∑b(j)*f(j) als „Stützpunkte” einer Folge kubischer Bezier-Kurven.

Funktionen

Die Objekt-Funktionen geben das Objekt selbst zurück, wenn kein anderer Rückgabewert angegeben ist.

use Herbaer::Ausgleichskurve;
$ausgl = new Herbaer::Ausgleichskurve ($sections);
$ausgl -> verbose ($v);
$ausgl -> add_curvweight ($weight);
$ausgl -> add_points ($points, $weight);
$ausgl -> get_cbezier ();
$ausgl -> fix_coeff ($i, $val);
$ausgl -> fix_value_at_point ($x, $val);
$ausgl -> fix_derivation_at_point ($x, $val);
$ausgl -> solve ();
$ausgl -> show_matrix ($handle);
$ausgl -> _get_list ($start, $step, $end, $modus);
$ausgl -> get_values  ($start, $step, $end);
$ausgl -> get_polygon ($start, $step, $end);

Herbaer::Ausgleichskurve::symmatrix_invert ($symmat);
Herbaer::Ausgleichskurve::symmatrix_mult_vector ($symmat, $vect);
Herbaer::Ausgleichskurve::symmatrix_show ($symmat, $handle);

new Herbaer::Ausgleichskurve ($sections)

Legt die Datenstruktur zur Berechnung einer neuen Ausgleichskurve an. $sections (ARRAY-Ref) ist eine steigenden Folge [$a0, $a1, ..., $an] der Punkte a(i), $a0 < $a1 < ... < $an. Die Datenstruktur enthält die folgenden Felder:

sect

Die Intervall-Zerlegung in der Form [a0, ..., an].

Die Referenz $sections wird in diesem Feld direkt gespeichert. Die Werte werden nicht kopiert. Das ARRAY, auf das $sections verweist, sollte bis zum Ende der Berechnung der Ausgleichskurve nicht geändert werden.

numintervals

Die Anzahl n der Teilintervalle: die Anzahl der Elemente von @$sections minus 1. Der Wert wird nicht mehr geändert.

vect

Als mögliche Ausgleichskurven werden Linearkombinationen ∑b(j)*f(j), 0 <= j < 2n + 2 der Funktionen f(j) (s. oben) betrachtet. Gesucht sind die Werte b(j), die einem Ideal möglichst nahe kommen. Die Abweichung vom Ideal beschreibt eine Abbildung b -> A(i,j)*b(i)*b(j) - 2*B(i)*b(i) + C (Summen über gleiche Indizes) mit einer symmetrischen Matrix A, einem Vektor B und einer Konstanten C. In der Regel ist A positiv definit. Dann ist b = inv(A) * B.

Unter diesem Schlüssel wird der B als ARRAY-Referez gespeichert.

pointdata

Eine Liste der Form [ { "points" => [ [x, y], ...], "weights" => [w, ...] }, ... ], s. add_points. Die Liste unter dem Schlüssel points ist „Punktwolken”, denen sich die Ausgleichskurve annähern soll, die positiven Zahlen der Liste unter dem Schlüssel weights sind die „Gewichte” der Punkte.

fixvalues

Eine Liste der Form [ [x, y], ... ]

Die Ausgleichskurve soll durch die Punkte [x, y] verlaufen. Die Werte x sind in der Liste der Intervallzerlegungspunkte sect enthalten, s. fix_value_at_point.

fixderivs

Eine Liste der Form [ [x, s], ... ]

s ist die vorgegebene Steigung der Ausgleichskurve beim gegebenen x-Wert. Die Werte x sind in der Liste der Intervallzerlegungspunkte sect enthalten, s. fix_derivation_at_point.

symmatrix

Die symmetrische Matrix A wird hier in der Form [ a(0,0), a(1,0), a(1,1), ... ] als Array mit 2n*n + 5n + 3 Elementen gespeichert. n ist hier die Anzahl der Teilintervalle (numintervals). Die Matrix-Indizes laufen von 0 bis 2n+1. Das Matrix-Element a(i,j) (= a(j,i), i <= j) wird an der Array-Position j(j+1)/2 + i gespeichert.

coeff

Eine Liste [ b(0), ... b(2n+1) ] von Zahlen. Die Ausgleichskurve ist F = ∑b(j)*f(j).

Die Ausgleichkurve läuft durch die Punkte [a(i), b(2i)], 0 <= i <= n, b(2i+1) ist die Steigung der Ausgleichskurve an der Stelle a(i) (s. sect).

coeff wird von der Funktion solve berechnet.

cbezier

Ein Array [ [x(0), y(0)], ... , [x(3n), y(3n)] ] von Stützpunkten kubischer Bezier-Kurven auf den Teilintervallen der Zerlegung.

x(3i)   = a(i)                                     für 0 <= i <= n
y(3i)   = b(2i)
x(3i+1) = (2a(i) +  a(i+1)) / 3                    für 0 <= i < n
y(3i+1) = b(2i) + b(2i+1) * (a(i+1) - a(i)) / 3
x(3i+2) = ( a(i) + 2a(i+1)) / 3
y(3i+2) = b(2i+2) - b(2i+3) * (a(i+1) - a(i)) / 3

Die zusammengesetzten Bezier-Kurven sind die Ausgleichskurve (s. get_cbezier).

status

Die Zahl unter dem Schlüssel status gibt den Stand der Berechnungen an. Die möglichen Werte sind:

status 0: Anfangswert, die Datenstruktur ist neu angelegt.
status 1: Es sind schon Teil-Abweichungsfunktionen (s. add_points, add_curvweight) definiert. Die Teilintervall-Zerlegung kann nicht mehr geändert werden.
status 2: Es wurden schon ein Punkt festgelegt, durch den die Ausgleichskurve verlaufen soll, oder die Steigung der Ausgliechskurve an einer Stelle festgelegt (s. fix_value_at_point, fix_derivation_at_point). Die Abweichungsfunktion kann nicht mehr geändert werden.
status 3: Die Matrix A (s. symmatrix) ist invertiert, die Ausgleichskurve berechnet.

verbose

Der Wert regelt den Umfang der Meldungen nach STDERR im Laufe der Berechnungen. Der Wert 0 bedeutet keine Meldung, ein größerer Wert bedeutet mehr Meldungen.

$ausgl -> verbose ($v);

Wenn $v definiert ist, wird der Wert unter dem Schlüssel verbose gespeichert. Der Rückgabewert ist der zuvor gespeicherte Wert unter dem Schlüssel verbose.

$ausgl -> add_curvweight ($weight)

Die Ausgleichskurve soll meist nicht zu „krumm” sein. Diese Funktion addiert das Integral des Quadrats der zweiten Ableitung der Ausgleichsfunktion F zur Fehlerfunktion.

Wenn $w eine Zahl ist, wird das Integral mit $w multipliziert. $w sollte positiv sein.

Wenn $w eine ARRAY-Referenz ist, wird das Quadrat der Ableitung auf dem i-ten Teilintervall mit der i-ten Komponente des Arrays multipliziert. Wenn die Zahl der Array-Komponenten kleiner ist als die Zahl der Teilintervalle, werden die Array-Komponenten zyklisch wiederholt.

$ausgl -> add_points ($points, $weight)

$points ist von der Form [ [x(0), y(0]], ... [x(m), y(m)] ], $weight von der Form [ w(0), ..., w(p) ].

Für jeden Punkt [x(i), y(i)] wird das gewichtete Quadrat der Abweichung w(i mod (p+1)) * (F(x(i)) - y(i)) ** 2 zur Fehlerfunktion addiert.

Dem Array pointdata wird die Komponente { "point" => $points, "weights" => $weight } angehängt.

Die Punkte in $points müssen nach wachsender x-Komponente geordnet sein, x(i) <= x(j) für i < j. Diese Einschränkung erlaubt eine effizientere Implementierung.

Wenn $weight eine einfche (positive) Zahl ist, wird dieser Wert so behandelt wie [ $weight ]. Das Gewicht $weight gilt für alle Punkte.

Wenn $weight nicht definiert ist, wird [ 1 ] angenommen.

$ausgl -> get_cbezier ()

Berechnet die Stützpunkte für die Darstellung der Ausgleichskurve als kubischer Bezier-Kurven auf den Teilintervallen, s cbezier.

$ausgl -> fix_coeff ($i, $val)

Der Koeffizient b($i) der Ausgleichskurve ∑b(j)*f(j) wird auf den Wert $val festgelegt.

Diese Funktion wird normalerweise von den Funktionen fix_value_at_point und fix_derivation_at_point aufgerufen. Sie darf für jeden Index $i, 0 <= $i < 2n + 2, höchstens einmal aufgerufen werden. Nach dem Aufruf dieser Funktion kann die Abweichungsfunktion, also die Matrix A (symmatrix) und der Vektor B (vect), nicht mehr geändert werden.

$ausgl -> fix_value_at_point ($x, $val)

Der Wert der Ausgleichskurve an der Stelle $x wird auf $val festgesetzt. $x muss einer der Teilungspunkte a(i) der Teilintervall-Zerlegung (s. sect) sein, sonst wirkt diese Funktion nicht. Sie ruft fix_coeff (2i, $val) auf. In jedem Fall wird das Paar [$x, $val] der Liste fixvalues angehängt.

$ausgl -> fix_derivation_at_point ($x, $val)

Der Steigung der Ausgleichskurve an der Stelle $x wird auf $val festgesetzt. $x muss einer der Teilungspunkte a(i) der Teilintervall-Zerlegung (s. sect) sein, sonst wirkt diese Funktion nicht. Sie ruft fix_coeff (2i+1, $val) auf. In jedem Fall wird das Paar [$x, $val] der Liste fixderivs angehängt.

$ausgl -> solve ()

Invertiert die symmetrische Matrix (s. symmatrix) und bestimmt coeff.

Diese Funktion wird normalerweise bei Bedarf durch eine der Funktionen get_values, get_polygon oder get_cbezier aufgerufen.

$ausgl -> show_matrix ($handle)

$handle ist das Handle eines Ausgabe-Streams (einer geöffneten Datei). Diese Funktion gibt die symmetrische Matrix symmatrix in Textform an $handle aus. Sie diente zur Fehlersuche.

$ausgl -> _get_list ($start, $step, $end, $modus)

Diese Funktion liefert eine Liste von Werten der Ausgleichsfunktion an den Stellen x(i) = $start + i * $step, i >= 0, x(i) <= $end. $step muss eine positive Zahl sein.

$modus kann die Werte 0 oder 1 annehmen. Im Falle $modus = 0 hat das Ergebnis die Form [ F(x(0)), ... ], im Falle $modus = 1 die Form [ [x(0), F(x(0))], ... ].

Diese Funktion wird normalerweise von get_values oder get_polygon aufgerufen.

$ausgl -> get_values ($start, $step, $end)

Ergibt eine Liste von Werten der Ausgleichsfunktion F in der Form [ F(x(0)), F(x(1)), ... ], x(i) = $start + i * $step, i >= 0, x(i) <= $end. $step muss eine positive Zahl sein.

$ausgl -> get_polygon ($start, $step, $end)

Ergibt eine Liste von Punkten der Ausgleichsfunktion F in der Form [ [x(0), F(x(0))], [x(1), F(x(1))], ... ], x(i) = $start + i * $step, i >= 0, x(i) <= $end. $step muss eine positive Zahl sein.

Herbaer::Ausgleichskurve::symmatrix_invert ($symmat)

Invertiert die symmetrische nxn-Matrix $symmat „auf der Stelle”.

$symmat ist eine ARRAY-Referenz der Form [ a(0,0), a(1,0), a(1,1), ... a(n-1,n-1) ]. Das Matrix-Element a(i,j), i <= j, steht an der Array-Position j(j+1)/2 + i. Das Array hat n*(n+1)/2 Elemente.

Die Matrix wird nach dem „Austauschverfahren” invertiert. In jedem der n Schritte wird der Index i mit dem größten Betrag |a(i,i)| gesucht und dann, aunschaulich gesagt, „x(i) gegen -y(i) ausgetauscht”. Das Verfahren bricht ab, wenn alle verbleibenden Diagonal-Elemente der Matrix Null sind. Wenn das vor dem n-ten Schritt eintritt, ist die Matrix nicht invertierbar. In diesem Fall wird eine Meldung nach STDERR ausgegeben.

Herbaer::Ausgleichskurve::symmatrix_mult_vector ($symmat, $vect)

Multipliziert die symmetrische nxn-Matrix $symmat (s. symmatrix_invert) mit dem Vektor $vect der Länge n.

$symmat und $vect sind ARRAY-Referenzen. Das Ergebnis ist eine neue ARRAY-Referenz.

Herbaer::Ausgleichskurve::symmatrix_show ($symmat, $handle)

Gibt die symmetrische nxn-Matrix $symmat an das geöffnete Ausgabe-Handle $handle in Textform aus. Diese Funktion diente zur Fehlersuche.