Ausgleichende Gerade im Gauß-Helmert-Modell

Vorwort

Im Rahmen dieses Artikels soll die Anwendung des Gauß-Helmert-Modells anhand der Bestimmung ausgleichender Geraden vorgeführt werden. Grundlegendes zur Methode der kleinsten Quadrate finden Sie (hier). Die allgemeine Herleitung des Gauß-Helmert-Modells kann (hier) nachgelesen werden.

 

Modellierung der Gerade im Gauß-Helmert Modell

Es gibt mehrere Möglichkeiten eine Gerade in einem ebenen Koordinatensystem zu beschreiben. Mit dem Blick auf die anschließende Parameterschätzung wird im weiteren Verlauf zunächst eine Gerade in der Punkt-Steigung-Form und danach in der Normalenform besprochen.

Punkt-Steigungs-Form

Gerade, beschrieben in der Punkt-Steigungsform

Gerade, beschrieben in der Punkt-Steigungsform

 

Die Gerade wird in der euklidischen (ebenen) Geometrie als die kürzeste Verbindung zwischen zwei Punkten definiert. Aus ihrem mathematischen Zusammenhang

(1)   \begin{equation*} \frac{y-y_1}{x-x_1} = \frac{y_2-y_1}{x_2-x_1} \end{equation*}

kann über algebraische Umformung zunächst der Wertebereich von y in Abhängigkeit der übrigen Parameter angegeben werden.

(2)   \begin{equation*} y = \left(\frac{y_2-y_1}{x_2-x_1}\right)(x-x_1) + y_1 \end{equation*}

Betrachtet man nun den Winkel \alpha zwischen der Abszisse und der Gerade so kann \frac{y_2-y_1}{x_2-x_1} auch mit Hilfe von trigonometrischen Funktionen ausgedrückt werden

(3)   \begin{equation*} \frac{Gegenkathete}{Ankathete} = \frac{y_2-y_1}{x_2-x_1} = \frac{\Delta y}{\Delta x} = \tan{\alpha} \end{equation*}

Nach dieser Substitution entsteht die Punkt-Steigung-Form

(4)   \begin{equation*} y = \tan{\alpha}(x-x_1)+y_1 \end{equation*}

Da eine Gerade in jedem Punkt eine konstante Steigung besitzt ist auch der Winkel \alpha und damit auch der Ausdruck \tan{\alpha} konstant und kann zusammenfassend als Parameter a geführt werden so dass nun die folgende Form entsteht

(5)   \begin{equation*} y = a(x-x_1)+y_1 \end{equation*}

Betrachtet man diesen Term für einen Punkt im Ursprung der Abszisse so dass a(x-x_1) = 0, dann kann y_1 auch als der Schnitt mit der Ordinate betrachtet und als Parameter b geführt werden. In (5) wurde bereits gezeigt dass die Steigung der Gerade überall konstant ist, auch an der Stelle x - x_1 womit x_1 auch vernachlässigt werden kann. So kann nun eine lineare Funktion in Abhängigkeit der beiden Parameter a und b, der Steigung der Gerade und dem Schnitt mit der Ordinate angegeben werden

(6)   \begin{equation*} y = ax+b \end{equation*}

Überführt man diese Funktion in das Gauß-Helmert-Modell so entsteht unter Berücksichtung der Verbesserungen die endgültige Bedingungsfunktion. Vergleiche hierzu NEITZEL, F., PETROVIC, S. (2008)1

(7)   \begin{equation*} a\left(x_i + v_{x_i}\right) - \left(y_i + v_{y_i}\right) + b = 0 \end{equation*}

Diese Form genießt den Vorteil dass ihre Parameter ohne jegliche Umrechnung intuitiv nachvollziehbar sind. Eine ausführliche Herleitung sowie ein nummerisches Beispiel zu dieser Bedingungsfunktion kann der o.g. Publikation von NEITZEL, F., u. PETROVIC, S. entnommen werden. Dieser Artikel ist besonders empfehlenswert da unter anderem auch auf die Tücken und Folgen des falsch angewandten Gauß-Helmert-Modells eingegangen wird.

Hessesche Normalenform

Gerade, beschrieben in der Normalenform

Gerade, beschrieben in der Normalenform

 

Die Gerade kann in der Ebene auch durch zwei nicht identische Punkte \mathbf{P}_1 und \mathbf{P}_2 eindeutig beschrieben werden. Betrachtet man nun \mathbf{n} als denjenigen Vektor, der senkrecht auf dem Richtungsvektor \mathbf{P}_2\mathbf{P}_1 der Gerade steht, so kann durch Projektion eines beliebigen auf der Gerade liegenden Punktes auf \mathbf{n}, der kürzeste Abstand d der Gerade zum Ursprung berechnet werden.

(8)   \begin{equation*}  d=\frac{\mathbf{n} \cdot \mathbf{P}_i}{|\mathbf{n}|}=\mathbf{n}_o\mathbf{P}_i \end{equation*}

In der Hesseschen Normalenform wird die Lage der Gerade eindeutig mit dem normierten Normalenvektor \mathbf{n_0} und dem Abstand der Gerade zum Ursprung d beschrieben.

(9)   \begin{equation*}  G:\left[ \begin{array}{c} n_1^0 \\ n_2^0 \end{array} \right]\cdot \left[ \begin{array}{c} x_i \\ y_i \\ \end{array} \right]+d=0 \end{equation*}

Überführt man diese Funktion in das Gauß-Helmert-Modell so entsteht unter Berücksichtigung der Verbesserungen die endgültige Bedingungsfunktion

(10)   \begin{equation*}  n_1\left(x_i+v_{x_i}\right)+n_2\left(y_i+v_{y_i}\right)+d=0 \end{equation*}

Beschrieben in dieser Form sind die Parameter nicht so intuitiv wie in der Punkt-Steigung-Form. Allerdings hat diese Form keinerlei nummerische Probleme wenn die Steigung sehr groß ist, da hier die Steigungen in der Abszisse und Ordinate getrennt betrachtet werden.

Ausgeglichene Zielmarke bestehend aus konzentrischen Kreisen und Geraden

Ausgeglichene Zielmarke bestehend aus konzentrischen Kreisen und Geraden

 

In vielen Ausgleichungsfällen wie z. B. der Zentrumsbestimmung von Zielmarken können beliebig gelagerte (auch senkrechte) Geraden im Datenmatrial vorliegen, die meistens als eine stützende Geometrie in die Ausgleichung eingeführt werden. In diesen Fällen eignet sich die Normalenform, weil sie durch die getrennte Betrachtung der Steigungen in der Abszisse und Ordinate keine nummerischen Probleme verursacht. Auch können Zusatzbedingungen wie z. B. senkrecht aufeinander stehende Geraden oder Gerade durch einen erzwungenen Punkt (Mittelpunkt der Zielmarke) mittels einfacher Vektoralgebra in die Ausgleichung eingeführt werden.

Umrechnung aus der Normalenform

Berechnung von α und β

Berechnung von α und β

 

Im ersten Schritt wird Winkel \alpha, zwischen dem normierten Normalenvektor n^0 = n/\|n\| und der Y-Achse bestimmt

(11)   \begin{equation*} \alpha = \arccos \left(\frac{1}{\sqrt{n_x^2 + n_y^2}} \begin{bmatrix} n_x \\ n_y \\ \end{bmatrix} \cdot \begin{bmatrix} 0 \\ 1 \\ \end{bmatrix} \right) =\arccos \left ( n_y^0 \right ) \end{equation*}

Aus der Winkelsumme des ebenen Dreiecks und dem nun bekannten Winkel \alpha kann wiederum der Winkel \beta berechnen

(12)   \begin{equation*} \beta = 180^{\circ} - 90^{\circ} - \alpha = 90^{\circ} - \alpha \end{equation*}

Berechnung des Achsenabschnittes aus dem Sinussatz

Berechnung des Achsenabschnittes aus dem Sinussatz

 

Da nun alle Winkel des Dreiecks bestimmt sind, lässt sich aus dem Sinussatz auch der Schnitt mit der Ordinate im Nullpunkt der Abszisse angeben

(13)   \begin{equation*} b = \begin{cases} undef.& \; \alpha = \pm 90^{\circ} \\ \frac{d}{cos\left(\alpha\right)} & \; \alpha \neq \pm 90^{\circ} \\ \end{cases} \end{equation*}

wobei die folgenden Substitutionen verwendet werden

(14)   \begin{equation*} \frac{b}{\sin \left( 90^{\circ} \right)} = b \end{equation*}

(15)   \begin{equation*} \frac{d}{\sin \left(90^{\circ} - \alpha\right)} = \frac{d}{\cos \left(\alpha\right)} \end{equation*}

Berechnung des Steigungswinkels über die Winkelsumme im Kreis

Berechnung des Steigungswinkels über die Winkelsumme im Kreis

 

Die Steigung der Gerade (3) lässt sich aus der Winkelsumme berechnen. Aus

(16)   \begin{equation*} a = \tan \left( 360^{\circ} - 180^{\circ} - 2\beta - \alpha \right) \end{equation*}

folgt durch Substitution gemäß (12)

(17)   \begin{equation*} a = \tan \left( 360^{\circ} - 180^{\circ} - 2(90^{\circ}- \alpha)\right) \end{equation*}

und (11) die gesuchte Steigung der Gerade

(18)   \begin{equation*} a = \tan \left(\alpha\right) = \tan \left( \arccos \left( n_y^0 \right) \right) \end{equation*}

welche algebraisch vereinfacht auch wie folgt angegeben werden kann

(19)   \begin{equation*} a = \frac{ \sqrt{1 - \left(n_y^0\right)^2}}{n_y^0}, \forall n_y^0 \neq 0 \end{equation*}

Strenge Gauß-Helmert-Lösung der Normalenform

Die Jacobi- oder Modellmatrizen \mathbf{A} und \mathbf{B} werden an der Stelle \mathbf{x_0,v_0} differenziert und erhalten die folgende Form:

(20)   \begin{equation*}  \mathbf{A}(\mathbf{x,v})=\frac{\partial \mathbf{\Omega (x,v)}}{\partial \mathbf{x}}= \begin{bmatrix} \frac{\partial \psi _1}{\partial n_1^0} & \frac{\partial \psi _1}{\partial n_2^0} & \frac{\partial \psi _1}{\partial d} \\ \frac{\partial \psi _2}{\partial n_1^0} & \frac{\partial \psi _2}{\partial n_2^0} & \frac{\partial \psi _2}{\partial d} \\ \vdots & \vdots & \vdots \\ \frac{\partial \psi _n}{\partial n_1^0} & \frac{\partial \psi _n}{\partial n_2^0} & \frac{\partial \psi _n}{\partial d} \end{bmatrix} \end{equation*}

Wir erhalten bei drei Unbekannten für \mathbf{A} eine Matrix mit drei Spalten. Die Anzahl der Zeilen hängt von der Anzahl gegebener Punkte ab. Da Verbesserungen berücksichtigt wurden, erscheinen diese in der Modellmatrix \mathbf{A}

(21)   \begin{equation*}  \mathbf{A}(\mathbf{x,v})= \begin{bmatrix} x_1+v_{x_1} & y_1+v_{y_1} & 1 \\ x_2+v_{x_2} & y_2+v_{y_2} & 1 \\ \vdots & \vdots & \vdots \\ x_n+v_{x_n} & y_n+v_{y_n} & 1 \end{bmatrix}_{n\times 3} \end{equation*}

Die Modellmatrix \mathbf{B} wird als Bandmatrix mit zwei horizontal zusammengesetzten und nicht überlappenden Diagonalmatrizen definiert. Jedes Band beinhaltet nur einen Wert. Beide Bänder haben n > 2 Elemente und beginnen jeweils in der ersten Zeile.

(22)   \begin{equation*}  \mathbf{B}(\mathbf{x,v}) =\frac{\partial \mathbf{\Omega (x,v)}}{\partial \mathbf{v}} = \end{equation*}

    \begin{equation*} \begin{bmatrix} \frac{\partial \psi _1}{\partial v_{x_1}} & 0 & \cdots & 0 & \frac{\partial \psi _1}{\partial v_{y_1}} & 0 & \cdots & 0 \\ 0 & \frac{\partial \psi _2}{\partial v_{x_2}} & 0 & \vdots & 0 & \frac{\partial \psi _2}{\partial v_{y_2}} & \text{} & \vdots \\ \vdots & 0 & \ddots & 0 & \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & \frac{\partial \psi _n}{\partial v_{x_n}} & 0 & \cdots & 0 & \frac{\partial \psi _n}{\partial v_{y_n}} \end{bmatrix} \end{equation*}

Durch die Bildung der Ableitung entfallen die Verbesserungen. Dies ist ein Sonderfall und darf nicht verallgemeinert ohne eine explizite Überprüfung auf andere Modelle übertragen werden.

(23)   \begin{equation*}  \mathbf{B}(\mathbf{x,v})=\left[ \begin{array}{cccc} n_1^0 & 0 & \cdots & 0 \\ 0 & n_1^0 & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & n_1^0 \end{array} \right. \left. \begin{array}{cccc} n_2^0 & 0 & \cdots & 0 \\ 0 & n_2^0 & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & n_2^0 \end{array} \right]_{n\times (2n)} \end{equation*}

Explizit ausgeschrieben, lautet der Widerspruch des Gleichungssystems, \mathbf{w}_1

(24)   \begin{equation*}  \mathbf{w}_1= -\mathbf{Bv_0+\phi (x_0,v_0)} = \end{equation*}

    \begin{equation*} \begin{bmatrix} -n_1^0v_{x_1}-n_2^0v_{y_1}+n_1^0\left(x_1+v_{x_1}\right) +n_2^0\left(y_1+v_{y_1}\right)+d^0 \\ -n_1^0v_{x_2}-n_2^0v_{y_2}+n_1^0\left(x_2+v_{x_2}\right) +n_2^0\left(y_2+v_{y_2}\right)+d^0 \\ \vdots \\ -n_1^0v_{x_n}-n_2^0v_{y_n}+n_1^0\left(x_n+v_{x_n}\right) +n_2^0\left(y_n+v_{y_n}\right)+d^0 \end{bmatrix}_{n\times 1} \end{equation*}

Dieser Vektor kann nach der Auflösung der Klammern algebraisch so zusammengefasst werden, dass eine reduzierte Form entsteht, die keinerlei Verbesserungen beinhaltet.

(25)   \begin{equation*}  \mathbf{w}_1=\left[ \begin{array}{c} n_1^0x_1+n_2^0y_1+d^0 \\ n_1^0x_2+n_2^0y_2+d^0 \\ \vdots \\ n_1^0x_n+n_2^0y_n+d^0 \end{array} \right]_{n\times 1} \end{equation*}

An dieser Stelle muss noch einmal darauf hingewiesen werden, dass es im Allgemeinen nicht üblich ist, dass die Verbesserungen v_i in der strengen Lösung im Widerspruchsvektor entfallen. Auf diesen Sonderfall weisen bereits NEITZEL, PETROVIC (2008) am Beispiel der ausgleichenden Geraden hin.

Einführung der Zwangsbedingung

Beliebig verlängerbarer Normalenvektor

Beliebig verlängerbarer Normalenvektor

 

Da man den Normalenvektor der Bedingungsgleichung (10) mit einer beliebigen Zahl ungleich Null multiplizieren kann, gibt es unendlich viele solcher Gleichungen, die jeweils die gleiche Gerade beschreiben. Das Normalgleichungssystem bleibt in diesem Fall singulär. Geometrisch gesehen entspricht dies, einer Verlängerung des Normalenvektors mit gleichbleibender Orientierung.

Um diese Mehrdeutigkeit zu lösen fordern wir dass der Betrag des Normalenvektors der Gerade die Länge 1 annimmt und führen hierzu eine Zwangsbedingung (engl. constraint, siehe Lösung mit Nebenbedingungen ) ein

(26)   \begin{equation*}  \gamma (\mathbf{x})=\sqrt{\left(n_1^0\right){}^2+\left(n_2^0\right){}^2}=1 \end{equation*}

Diese wird nur an der Stelle \mathbf{x}_0 differenziert, weil hiermit eine Zwangsbedingung zwischen den Parametern beschreiben. Somit lautet unser Constraint-Vektor \mathbf{C}

(27)   \begin{equation*}  \mathbf{C}(\mathbf{x})=\left.\frac{\partial \gamma }{\partial \mathbf{x}}\right|_{\mathbf{x}_0} = \frac{1}{\gamma (\mathbf{x})}\left[ \begin{array}{cc} n_1^0 \ n_2^0 \end{array} \right]_{1\times 2} \end{equation*}

Der dazugehörige skalare Widerspruch w_2 beschreibt, wie groß die Abweichung von der Länge 1 ist

(28)   \begin{equation*}  w_2 = \gamma(\mathbf{x}) - 1 \end{equation*}

 

Quellcode

Nachfolgend sei der Octave Quellcode für die strenge Lösung der vier Parameter Ebene gelistet. Für die Bildung der Matrizen wurde auf Octave-Matrizenfunktionen zurückgegriffen, da die Bildung über Summen, die in Octave z. B. als for()-Schleifen realisiert werden können, wesentlich ineffizienter ist. Generell sollte man in Octave, soweit es möglich ist, auf Matrizenoperationen zurückgreifen, da diese intern über hocheffiziente, nummerisch stabile und teilweise mehrprozessoroptimierte BLAS2 und LAPACK Bibliotheken realisiert und nicht über die wesentlich langsamere M-Scriptsprache berechnet werden. Diese Behauptung kann natürlich nicht allgemein auf andere Programmiersprachen wie C/C++ übertragen werden, die auch mit for()-Schleifen sehr effizient eingesetzt werden können.

Quellcode als Octave Datei herunterladen:

Download hier: mbest_fit_line - v1.03 - m

%% Berechnung ausgleichender Geradenparameter im Gauss-Helmert-Modell
%
% v1.03
%
% Kisser Waldemar (2011)
% http://kisser.online/ausgleichung/ghm/gerade
 
clc;
clear all;
 
%% Datensatz waehlen
% Z.B. aus Datei laden: "Punkte = dlmread( 'C:/daten/linefit/daten.xy' );"
% oder
Punkte = ...
[ ...
     0.0 0.0
     1.0 1.0
     2.0 4.0
     3.0 9.0
];
 
%% Hilfsvariablen
 
% Anzahl der Beobachtungen , Dimension
[ n , d ] = size( Punkte );
 
% Initialisierung der Naeherungswerte
% nx, ny, d
X_dach = rand( 3 , 1 );
X_dach( 1 : 2 ) = X_dach( 1 : 2 ) ./ norm( X_dach( 1 : 2 ) );
X_dach( 3 ) = norm( mean( Punkte ) );
 
%% Schaetzung der exakten Minimierer X0,V0
 
% Vektor der Verbesserungen
v = zeros( d * n , 1 );
 
% Anzahl der Unbekannten
u = size( X_dach , 1 );
 
% Bedingungsvektor
c = 1;
C = zeros( c , u );
 
% Anzahl der Bedingungen
b = n * d;
 
% Legt den gewuenschten Konvergenzbetrag fest ab dem die
% Iterationen abgebrochen werden koennen.
EPSILON = 1e-12;
 
% Standardabweichung der Gewichtseinheit a priori
s0_prio = 1.0;
 
% Kofaktormatrix (hier gleichgenaue unkorrelierte Beobachtungen)
Cl = eye(2*n);
 
% Kovarianzmatrix
Ql = 1 / s0_prio^2 * Cl;
 
% Maximale Iterationsanzahl
maxit = 20;
 
% Beginn des Iterationsvorgangs
for iteration = 0 : maxit
 
    % Modellmatrix A
    A = [ Punkte + reshape( v , n , d ) , ones( n , 1 ) ];
 
	% Bildung des Bedingungsvektors
	C( 1 : u - 1 ) = X_dach( 1 : u - 1 ) ./ norm( X_dach(  1 : u - 1 ) );
 
	% Bildung der B-Modellmatrix
	% Da diese 2(n^2-n) Nullelemente hat, eignet sich der Sparse
	% Datentyp fuer duennbesetzte Matrizen.
	B = spdiags( repmat( X_dach( 1 : 2 ).' , n , 1 ) , [ 0 , n ] , n , 2 * n );
 
	% Bildung der Inversen zur B*Ql*B^T Matrix
	BQBT = B * Ql * B.';
 
	% Bildung der Normalgleichungsmatrix
	N = [ -A.' *( BQBT \ A ) ,  C.' ; C , zeros( c , c ) ];
 
	Ncond = rcond( N );
	Ndet = abs( det( N ) );
 
	if ( ( Ndet < EPSILON ) || ( Ncond < EPSILON ) )
        fprintf( 'Fehler: Normalgleichungsmatrix singul?r oder\n' );
        fprintf( 'schlecht konditioniert\n' );
        fprintf( 'Determinante    : %.3e\n' , Ndet );
        fprintf( 'Konditionierung : %.3e\n\n' , Ncond );
        return
	end 
 
	% Bilde den Widerspruchsvektor
	wm = [ Punkte , ones( n , 1 ) ] * X_dach;
 
	% Bilde den Widerspruchsskalar
	wc = 1 - C * X_dach;
 
	% Berechne die Differenzloesung
	dx = N \ [ A.' * ( BQBT \ wm ) ; wc ];
 
	% Berechne die Verbesserungen dieser Iteration
	v = Ql * B.' *( BQBT \ ( wm - A * dx( 1 : u ) ) );
 
	% Addiere die Differenzloesung
	X_dach = X_dach + dx( 1 : u ); 
 
	% Zweite Abbruchbedingung
    if max( abs( dx( 1 : u ) ) ) < EPSILON
        break
    end
end
 
%% Normalenvektor gemaess Hesseschen Normalenform
vz = 1.0;
if X_dach( end ) < 0.0
    X_dach = - X_dach;
    vz = -1.0;
end
 
%% Ausgabe
fprintf( '--- AUSGLEICHUNGSPROTOKOLL ---' );
fprintf( '\n\n-- Ausgleichungsvorgang --' );
fprintf( '\nModell\t\t\t\t: Gauss-Helmert-Modell' );
fprintf( '\nVersion\t\t\t\t: 1.03' );
fprintf( '\nForm\t\t\t\t: 3 Parameter Gerade' );
 
if iteration < maxit
 
    fprintf( '\nKonvergenz\t\t\t: Erfolgt' );
    fprintf( '\nKonvergenzgrenze\t: %.1e' , EPSILON );
    fprintf( '\nAnzahl Iterationen\t: %d' , iteration );
 
    fprintf( '\n\n-- Statistik --' );
 
    redundanz = n - u + c;
    vTPv =  v.' * ( Ql \ v );
    s0_post = sqrt( vTPv / redundanz );
    Qxdach = inv( - N );
 
    fprintf( '\nAnzahl Beobachtungen\t: %d' , n );    
    fprintf( '\nAnzahl Parameter\t\t: %d' , size( X_dach , 1 ) );
    fprintf( '\nAnzahl Bedingungen\t\t: %d' , c );
    fprintf( '\nGesamtredundanz\t\t\t: %d' , redundanz );
    fprintf( '\n\nvTPv\t: %12.8f' , vTPv );
    fprintf( '\ns0_prio\t: %12.8f' , s0_prio );
    fprintf( '\ns0_post\t: %12.8f' , s0_post ); 
 
    fprintf( '\n\n-- Parameter --' );
    fprintf( [ '\nnx\t\t: %12.8f\t' , char(177) , '%12.8f' ] , X_dach( 1 ) , s0_post * sqrt( Qxdach( 1 , 1 ) ) );
    fprintf( [ '\nny\t\t: %12.8f\t' , char(177) , '%12.8f' ] , X_dach( 2 ) , s0_post * sqrt( Qxdach( 2 , 2 ) ) );
    fprintf( [ '\nd\t\t: %12.8f\t' , char(177) , '%12.8f' ] , X_dach( 3 ) , s0_post * sqrt( Qxdach( 3 , 3 ) ) );
 
    %% Berechnung der Punkt-Steigung-Form
    if abs( X_dach( 2 ) ) > EPSILON
        a0 = tan( acos( X_dach( 2 ) ) );
        b0 = vz * X_dach( 3 ) / X_dach( 2 );
 
        %% Varianzfortpflanzung
        A2 = ...
        [ ...
            0 , - vz / X_dach( 2 )^2 / sqrt( 1 - X_dach( 2 )^2 ) , 0 ;
            0 , - vz * X_dach( 3 ) / X_dach( 2 )^2 , vz / X_dach( 2 ) 
        ];
 
        sab = s0_post * sqrt( diag( A2 * Qxdach( 1 : u , 1 : u ) * A2.' ) );
 
        fprintf( '\n\n-- 2D-Geradenparameter --' );
        fprintf( [ '\na\t\t: %12.8f\t' , char(177) , '%12.8f' ] , a0 , sab( 1 ) );
        fprintf( [ '\nb\t\t: %12.8f\t' , char(177) , '%12.8f' ] , b0 , sab( 2 ) );
 
    end
 
    fprintf( '\n' );
else
    fprintf( '\nKonvergenz: Nicht erfolgt\n' );    
end

 

Quellen und Fußnoten

  1. NEITZEL, F., PETROVIC, S. (2008): Total Least Squares (TLS) im Kontext der Ausgleichung nach kleinsten Quadraten am Beispiel der ausgleichenden Geraden. Zeitschrift für Geodäsie, Geoinformation und Landmanagement 3/2008 S. 141-148, Wißner-Verlag, Augsburg []
  2. Basic Linear Algebra Subprograms ist eine Softwarebibliothek für elementare Matrizenoperationen. Dank der starken Unterstützung von Hardwareentwicklern wie Intel oder AMD ist BLAS derzeit ein Quasi-Standard und bildet damit eine feste Teilmenge von erweiterten nummerischen Bibliotheken wie LAPACK/Intel MKL/GSL/CUDA SDK und vielen mehr. []