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.
use Herbaer::Ausgleichskurve; my $sections = [0, 10, 20]; my $points = [ [1,3], [4,8], [7,9], [9,8], [11, 4], [14,2], [17,0], [19,-3] ]; my $ausgl =new
Herbaer::Ausgleichskurve ($sections); $ausgl ->add_points
($points); my $cbezier = $ausgl ->get_cbezier
();
Mit dem Modul Herbaer::Punktediagramm
kann so eine grafische Darstellung erzeugt werden:
use Herbaer::Ausgleichskurve; use Herbaer::Punktediagramm; my $sections = [0, 10, 20]; my $points = [ [1,3], [4,8], [7,9], [9,8], [11, 4], [14,2], [17,0], [19,-3] ]; new Herbaer::Punktediagramm (3000, 2000, "Diagramm") -> add_pointset ("punkte", $points) -> add_cbezier ( "ausgleich", new Herbaer::Ausgleichskurve ($sections) -> add_points ($points) -> get_cbezier () ) -> y_autorange () -> x_autoticks () -> scale (115, 2960, 1940, 115) -> write_xml ("diagramm.xml"); system ("xsltproc pd_svg.xslt diagramm.xml > diagramm.svg"); system ("convert diagramm.svg diagramm.png");
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);
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.
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
.
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.
$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.
Berechnet die Stützpunkte für die Darstellung der Ausgleichskurve als kubischer Bezier-Kurven auf den Teilintervallen, s cbezier
.
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.
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
auf. In jedem Fall wird das Paar fix_coeff
(2i, $val)[$x, $val]
der Liste fixvalues
angehängt.
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
auf. In jedem Fall wird das Paar fix_coeff
(2i+1, $val)[$x, $val]
der Liste fixderivs
angehängt.
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.
$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.
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.
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.
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.
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.
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.
Gibt die symmetrische nxn
-Matrix $symmat
an das geöffnete Ausgabe-Handle $handle
in Textform aus. Diese Funktion diente zur Fehlersuche.
Das Programm ist mit Perl Version 5.24.1 entwickelt.