Wir kommen nun auf das Problem der Division durch Null (vom Ende des 1. Abschnitts) zurück. Bei der folgenden Matrix versagt unser Algorithmus, obwohl sie invertierbar ist:
A = Matrix.fromString("0, 5, 2; 2, 6, 4; 2, 1, 1") luIter(A)
Wir können die Division durch Null hier verhindern, indem wir die ersten beiden Zeilen der Matrix vertauschen:
A1 = Matrix.fromString("2, 6, 4; 0, 5, 2; 2, 1, 1") luIter(A1) print(A1)
Die Vertauschung der der ersten beiden Zeilen lässt sich erreichen, indem man A von links mit der entsprechenden Permutationsmatrix multipliziert. Diese entsteht aus der Einheitsmatrix durch Vertauschung der der entsprechenden Zeilen:
In einer invertierbaren Matrix lassen sich die Zeilen immer so permutieren, dass unser Algortihmus anschließend erfolgreich ist. Es gilt nämlich:
Satz
Zu jeder invertierbaren Matrix A gibt es
so dass PA = LU.
Wir können das lineare Gleichungssystem auf beiden Seiten mit P multiplizieren, und dann mit der Zerlegung der der permutierten Matrix arbeiten:
Dieses Gleichungssystem lässt sich in zwei Schritten lösen:
In der nun folgenden Implementierung vertauschen wir nicht nur, um Divisionen durch Null zu verhindern, sondern bringen immer das dem Betrage nach größte Element an die Divisionsposition. Damit erhöht sich die numerische Stabilität beim Rechnen mit Gleitkommazahlen.
def lup(A): A = A.copy() # Kopie, damit übergebene Matrix unverändert bleibt n = A.height() P = Matrix.identity(n) # Ausgangsbasis für Permutationsmatrix for k in fromTo(0,n-2): # Suche das dem Betrag nach größte Element # der k-ten Spalte ab der k-ten Zeile: absMax, kMax = max([(abs(A[i,k]), i) for i in fromTo(k,n-1)]) assert absMax > 0 # Vertausche die Zeile des größten Elements mit der k-ten Zeile... P[k], P[kMax] = P[kMax], P[k] # ... in der Permutation A[k], A[kMax] = A[kMax], A[k] # ... und in der Matrix for i in fromTo(k+1,n-1): # Dividiere die restlichen Elemente # der k-ten Spalte durch das „Kopfelement“: A[i,k] /= A[k,k] # Subtrahiere in der Restmatrix überall das Produkt # aus den „Randelementen“ links und oberhalb: for j in fromTo(k+1,n-1): A[i,j] -= A[i,k] * A[k,j] # A in linke und rechte Dreiecksmatrix zerlegen L = Matrix.identity(n) U = Matrix.null(n) for i in range(n): for k in range(n): if k < i: L[i,k] = A[i,k] else: U[i,k] = A[i,k] return L, U, P
Ein Test mit der Matrix von oben:
A = Matrix.fromString("0, 5, 2; 2, 6, 4; 2, 1, 1")
L,U,P = lup(A)
print(L,U,P)
Und die Kontrolle:
print(P * A) print(L * U)
Die Funktion lup
ist in mathGUIde bereits eingebaut,
allerdings nicht als globale Funktion sondern als Matrix-Methode:
A = Matrix.fromString("0, 5, 2; 2, 6, 4; 2, 1, 1")
L,U,P = A.lup()
print(L,U,P)