Lineare Algebra

Dieser Abschnitt erläutert

Lösen linearer Gleichungssysteme

Es gibt eine ganze Reihe von Build-In-Funktionen zur Ausführung von linearer Algebra. Das Matrizenprodukt wurde bereits besprochen. Die Operation

    >A\b

liefert auf x, einem Vektor des Typs (n,1), die Lösung des Gleichungssystems Ax=b, wobei A eine Matrix des Typs (n,n) und b ein Vektor des Typs (n,1) ist. Wenn in

    >A\B

B eine Matrix des Typs (m,n) ist, wird das System simultan für alle A\B[i] mit i über alle Spalten gelöst. Falls die Determinate von A kleiner als der interne Genauigkeitswert epsilon - d.h. nahe Null - sein sollte erfolgt eine Fehlermeldung.

Es gibt auch eine etwas präzisere Version, die eine Restwert-Iteration benutzt. Damit erhält man gewöhnlich sehr gute Ergebnisse. Natürlich ist die Funktion langsamer.

    >xlgs(A,b)

Sie können der Funktion zusätzlich eine maximale Anzahl von Iterationen mitgeben.

    >inv(A)

berechnet die Inverse der Matrix A. Diese Funktion ist als

    >A\id(cols(A))

definiert. Es gibt auch tiefergehende Funktionen wie

    >lu(A)

für Matrizen des Typs (m,n). Diese und die folgenden Funktionen wollen wir für die mathematisch Interessierten hier im Detail erklären. Sie liefert mehrere Werte zurück (siehe Mehrfache Zuweisungen). Sie können die Rückgabewerte mit

    >{Res,ri,ci,det}=lu(A)

Variablen zuweisen. Wenn Sie nur

    >Res=lu(A)

verwenden sind die anderen Ausgaben verloren. Um die Ausgaben zu erklären, wollen wir mit Res beginnen. Res ist eine Matrix des Typs (m,n), die die LU-Zerlegung (engl. für Lower/Upper) von A enthält; d.h., L.U=A mit einer unteren Dreiecksmatrix L und einer oberen Dreiecksmatrix U. L besitzt nur Einsen in der Diagonale, die übergangen werden, so dass L und U in Res gespeichert werden können. det ist natürlich die Determinate von A. ri enthält die Indizes der Zeilen von Res, da der Algorithmus Zeilen vertauscht haben könnte. ci ist uninteressant, wenn A nichtsingulär ist. Falls A jedoch singulär ist enthält Res das Ergebnis des Gauss-Algorithmus und ci 1 und 0, je nachdem ob die Spalten eine Basis der Spalten von A bilden.

Ein Beispiel:

    >A=random(3,3);
    >{LU,r,c,d}=lu(A);
    >LU1=LU[r];
    >L=band(LU1,-2,-1)+id(3); R=band(LU1,0,2);
    >B=L.R,

ergibt die Matrix A[r]. Um A zu erhalten, muss man die inverse Permutation r1 berechnen.

    >{rs,r1}=sort(r); B[r1]

B ist dann wieder identisch zu A.

Mit der LU-Zerlegung von A können wir sehr schnell lineare Gleichungssysteme der Form A.x=b lösen. Diese Form ist äquivalent zu A[r].x=b[r], wobei LU[r] die Zerlegung von A[r] ist, so dass x mit

    >{LU,r}=lu(A);
    >lusolve(LU[r],b[r])

berechnet werden kann. Eine etwas exaktere Version ist

    >xlusolve(A,b)

falls A bereits die LU-Form besitzt. Das heißt, diese Funktion arbeitet mit A in der Form einer oberen Dreiecksmatrix, z.B.

    >A=random(10,10); A=band(A,0,10);
    >xlusolve(A,b)

Diese Funktion kann also zur exakten Auswertung arthmetischer Ausdrücke benutzt werden.

lu wird von verschiedenen Funktionen in util.e benutzt. Z.B. ist

    >kernel(A)

eine Basis des Kerns von A; d.h. der Vektoren x bei Ax=0.

    >image(A)

ist eine Basis der Vektoren Ax. Sie können kernel und image einen zusätzlichen Wert als Parameter eps=... mitgeben, der das interne epsilon in diesen Funktionen ersetzt. Diese Funktionen normieren die Matrix mit

    >norm(A)

das die maximale Zeilensumme von abs(A) zurückliefert.

Auch die Singulärwertzerlegung ist implementiert. Die Basisfunktion ist

    >{U,w,V}=svd(A)

Diese Funktion liefert drei Werte zurück. A muss eine reelle Matrix des Typs (m,n) sein, dann ist U ebenfalls eine Matrix des Typs (m,n), w ein Vektor des Typs (1,n) und V eine Matrix des Typs (m,n). Die Spalten von U und V sind orthogonal. Wir erhalten A=U.W.V', wobei W eine Diagonalmatrix des Typs (n,n) mit dem Vektor w auf der Hauptdiagonalen ist, d.h.

    >A=U.diag(size(V),0,w).V'

Diese Zerlegung ist unter vielen Umständen nützlich. Die Datei svd.e (wird beim Programmstart geladen) enthält die Applikationen svdkernel und svdimage, die die orthogonale Basis des Kerns und des Abbilds von A berechnen. svddet berechnet die Determinate von A und svdcondition eine Konditionalzahl.

    >fit(A,b)
    >svdsolve(A,b)

findet die Lösung von A.x=b für singuläre Matrizen A (auch nicht-rechteckige) durch Minimierung der Norm von A.x-b. Die Funktion svdsolve ist stabiler und sollte bevorzugt werden. U, w und V können dazu benutzt werden Lösungen von A.x=b mit

    >x=V.diag(size(V),0,1/w).U'

zu berechnen, wenn w Zahlen ungleich Null enthält. Dies ist eine ähnliche Vorgehensweise wie bei der obigen Funktion lu. Übrigens würde

>svdsolve(A,id(cols(A));

die sogenannte Pseudoinverse von A berechnen. Ebenfalls gibt es svddet, das stabiler als det sein könnte.

Eigenwerte

Die Grundfunktion zur Berechnung von Eigenwerten ist

    >charpoly(A)

die das charakteristische Polynom von A berechnet. Sie wird von

    >eigenvalues(A)

benutzt, um die Eigenwerte von A zu berechnen.

    >eigenspace(A,l)

berechnet eine Basis des Eigenwertraumes von l. Diese Funktion benutzt kernel und bricht ab, falls der Eigenwert nicht exakt genug bestimmt werden kann.

    >{l,x}=xeigenvalue(A,l)

verbessert den Eigenwert l, der ein einfacher Eigenwert sein muss. Die Funktion liefert den verbesserten Wert l und einen Eigenvektor zurück. Sie können der Funktion eine weiteren Parameter zur Verfügung stellen, der eine Approximation des Eigenvektors sein muss.

    >{l,X}=eigen(A)

liefert die Eigenwerte von A in l zurück und die Eigenvektoren in X. Es gibt eine verbesserte, aber langsamere Version eigen1, die öfter zum Erfolg führt als eigen. Es gibt auch die Routine svdeigen, die die Singulärwertzerlegung zur Bestimmung des Kerns benutzt.

    >jacobi(a)

benutzt die Methode von Jacobi zur Berechnung der Eigenwerte einer symmetrischen Matrix des Typs (n,n).

Ein spezielles Feature ist die Generierung von Toeplitz-Matrizen mit der Funktion toeplitz. Der Parameter dieser Funktion ist ein Vektor r mit ungerader Anzahl von Elementen 2n-1 und die Ausgabe eine Matrix R mit R(i,j)=r(n-i+j), d.h. die letzte Zeile von R entspricht den ersten n Elementen von r und alle Spalten darüber wurden um eins nach links verschoben von unten nach oben. Sie können ein lineares Gleichungssystem mit einer Toeplitz-Matrix mit

    >toeplitzsolve(r,b)

lösen, wobei b ein Vektor des Typs (n,q) ist. Das Ergebnis ergibt wieder toeplitz(r).x=b.

Das interne Epsilon

    >epsilon()

ist ein internes epsilon, das von vielen Funktionen und dem Operator ~= benutzt wird. ~= vergleicht zwei Werte und gibt 1 zurück, wenn die absolute Differenz kleiner ist als epsilon. Dieses epsilon kann mit

    >setepsilon(value)

verändert werden.