Eigenwerte einer Matrix

Die Eigenwerte einer Matrix erhalten wir durch den Befehl eigenvalues  im Paket linalg .  Geben wir z. B. die Matrix

>	A:= 1/8*evalm(matrix(3,3,[43, -15, 3, 21, 15, -3, 30, -6, 14]));

ein, so ergeben sich die drei Eigenwerte

>	with(linalg): eigenvalues(A);

Warning, the protected names norm and trace have been redefined and unprotected

Selbstverstaendlich haetten wir auch direkt vorgehen koennen, indem wir uns das charakteristische Polynom verschaffen und die Nullstellen bestimmen. Dazu koennen wir die uns schon bekannten Befehle det und solve nutzen.

>	E:=array(identity,1..3,1..3); p:=det(A-lambda*E); EW:=[solve(p=0,lambda)];

Fuer das charakteristische Polynom ist ein eigener Befehl charpoly  im Paket linalg vorgesehen. Wir erhalten, wie oben schon berechnet,

>	charpoly(A,lambda);

Auch die Eigenvektoren lassen sich direkt durch das Kommando   eigenvectors   bestimmen

>	eigenvectors(A);

Der Befehl liefert eine Liste der Eigenwerte mit ihrer algebraischen Vielfachheit und einer Basis des zugehoerigen Eigenraums.

Nun  ueberpruefen wir noch die Aussage, dass die Summe der Diagonalelemente einer Matrix, die sogenannte Spur ( trace ), gerade die Summe der Eigenwerte ist.

>	sum(EW[j],j=1..3); trace(A);

Bemerkung: Durch das Laden des Pakets linalg hat sich die Bedeutung des Kommandos   trace  geaendert. Ohne das Paket ist es ein Befehl zur Fehlersuche (debuggen) bei Prozeduren. Dieser Befehl wird durch das Paket linalg ueberschrieben.

Verifizieren wir noch, dass das Produkt der Eigenwerte gleich der Determinante der Matrix ist:

>	product(EW[j],j=1..3); det(A);

Die geometrische Bedeutung des Eigenwerts als Streckfaktor in Richtung  der Eigenvektoren, koennen wir uns natuerlich  veranschaulichen, indem wir einen Eigenvektor herausgreifen und sein Bild unter der linearen Abbildung  betrachten.  

>

x:=matrix(3,1,[1,5/3,2]); y:=evalm(A&*x); with(plots): with(plottools): p1:=plot3d(0,0..5,0..5,color=gray,style=patchnogrid,axes=boxed): p2:=sphere([seq(x[j,1],j=1..3)],0.1,style=patchnogrid,color=blue): p3:=sphere([seq(y[j,1],j=1..3)],0.1,style=patchnogrid,color=red): p4:=line([0,0,0],[seq(y[j,1],j=1..3)],linestyle=1,color=black): display(p1,p2,p3,p4);

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

Wenn Sie das Bild mit der linken Maustaste anklicken, koennen Sie es drehen, um es sich besser dreidimensonal vorstellen zu koennen. Offensichtlich ist y (rot) der um den Faktor 3 gestreckte Vektor x (blau).

Im  Satz von Cayley-Hamilton wird festgestellt, dass wir die Nullmatrix erhalten, wenn wir in das charakteristische Polynom einsetzen. Versuchen wir dies: Zunaechst wandeln wir in eine Funktion.

>	p:=unapply(p,lambda);

Setzen wir nun die Matrix als Argument in das Polynom ein und werten den Ausdruck als Matrix aus, so folgt, wie zu erwarten,

>	evalm(p(A));

Wie anschliessend im Buch auf Seite 595 erlaeutert, koennen wir daraus die Inverse zu bestimmen, indem wir  aus   ausklammern. Mit

>	q:=unapply(simplify((24-p(lambda))/lambda),lambda);

vergleichen wir

>	evalm(1/24*q(A));

und

>	inverse(A);

>

Diagonalisierbarkeit

Inwieweit eine Matrix diagonalisierbar ist, kann mit Hilfe des Befehls jordan ermittelt werden, der die Jordannormalform einer Matrix anzeigt (zur Vertiefung siehe den optionalen Abschnitt 18.8). In unserem Beispiel mit drei verschiedenen Eigenwerten ist die Matrix diagonalisierbar und wir erhalten die Diagonalmatrix mit den Eigenwerten auf der Diagonalen.

>	J:=jordan(evalm(A));

Geben wir im Befehl noch ein weiteres Argument an, etwa 'S',  so wird der Variable S die Transformationsmatrix zugeordnet, d.h. es ist  

>	jordan(evalm(A),'S'): evalm(S); evalm(inverse(S)&*A&*S);

Betrachten wir noch ein anderes  Beispiel:

>	B:=matrix(3,3,[[2, 0, 0],[-2, 2, 3], [3, 0, -1]]); jordan(evalm(B));

In diesem Fall ist die Matrix  nicht diagonalisierbar. (Am unteren Jordanblock erkennt man, dass Eigenwert der Matrix ist mit algebraischer Vielfachheit 2 und geometrischer Vielfachheit 1)  

Fuer Matrizen ist auch die Exponentialfunktion  ( exponential ) im Paket linalg vorgesehen. Ist die Matrix diagonalisierbar, so wird fuer die Matrix   die explizite Darstellung mit den entsprechenden Eigenwerten und Eigenvektoren berechnet.

>	exponential(A,t);

>

Numerische Berechnung von Eigenwerten

Wir schreiben uns ein eigenes Programm zur Vektoriteration, um den groessten Eigenwert einer Matrix zu approximieren.

>

viter:=proc(A, v0, n) local j, xn, vn, vmax, lambda, lambdalist;   vn         := v0; lambdalist := []; for j to n do   xn  := evalm(A&*vn);   vmax:=norm(vn,infinity);   member(vmax,[seq(vn[i,1],i=1..rowdim(vn))],'jmax');   lambda     := evalf( xn[jmax,1]/vn[jmax,1], 3);   lambdalist := [op(lambdalist), lambda] ;   vn         := evalm(xn); od; op(lambdalist) end:

Bemerkung: Mit dem Befehl op lassen sich die Operanden eines Ausdrucks, d.h. in diesem Fall die Eintraege in der Liste, ausgeben.

Probieren wir nun das Programm aus mit der anfaenglichen Matrix A und dem Startvektor (1,1,1). Es ergibt sich

>	A:= 1/8*evalm(matrix(3,3,[43, -15, 3, 21, 15, -3, 30, -6, 14])); v0:=matrix(3,1,[[1],[1],[1]]); viter(A,v0,20);

Offensichlich konvergiert die Iteration. Deutlich ist in der Liste der Naeherungen an den Eigenwert die relativ langsame Konvergenzgeschwindigkeit zu sehen.

Das am haeufigsten angewendete Verfahren zur Berechnung von Eigenwerten diagonalisierbarer Matrizen ist das QR-Verfahren, fuer dessen Beschreibung wir auf die Literatur zur Numerischen Mathematik verweisen. Es beruht auf einer Zerlegung der Matrix in das Produkt  einer  unitaeren Matrix Q und eine obere Dreiecksmatrix R. Diese Methode ist in Maple unter dem Kommando Eigenvals implementiert.  

>	evalf( Eigenvals( evalm(A) ) );

>

Aufgaben

1. Bestimmen Sie die Eigenwerte und zugehoerige Eigenvektoren zu

>	A:= evalm(matrix(3,3,[1, -3, 3, 3, -5, 3, 6, -6, 4]));

indem Sie Nulllstellen des charakteristischen Polynoms bestimmen und die Eigenvektoren aus den entsprechenden linearen Gleichungssystemen berechnen.  Ueberpruefen Sie ihr Ergebnis durch die Befehle eigenvalues und eigenvectors.

>

2. a) Wir greifen Aufgabe 18.10 auf: Bestimmen Sie eine orthogonale Matrix S, sodass S^(-1)*A*S Diagonalmatrix ist mit

>	A:= evalm(matrix(3,3,[10, 8, 8, 8, 10, 8, 8, 8, 10]));

b) Pruefen Sie auch noch das  Resultat in Aufgabe 18.9 (a) zur Matrix

>	A:= evalm(matrix(2,2,[1, I, I, -1]));

und berechnen Sie .

>

3. Fuer welche reellen Zahlen ist die Matrix

>	A:= evalm(matrix(3,3,[1, 0, 0, 0, 3, a, 1, a, 3]));

diagionalisierbar?

>

>