LU-Zerlegung - Step by step¶
Als Einleitung zum eigentlichen Praktikum setzen Sie hier die LU-Zerlegung Schritt für Schritt an einem Beispiel um. Gegeben ist das LGS \(A \vec{x}= \vec{b}\), mit \(A, \vec{b}\)
Erste Variante: ohne Pivotisierung¶
Die Matrix \(A\) ist vorgegeben. Wir erstellen gleich eine kleine Funktion, die diese Matrix zurückgibt, um später jederzeit Zugriff auf das Original zu haben.
[1]:
import numpy as np
def getA():
return np.array([[3.,6,3],[1,3,6],[6,3,3]])
[2]:
A = getA()
print(A)
[[3. 6. 3.]
[1. 3. 6.]
[6. 3. 3.]]
erste Zeilenoperation¶
Subtrahieren Sie geeignete Vielfache der ersten Zeile von der zweiten und dritten. Berechnen Sie jeweils die benötigten Faktoren aus den Matrixeinträgen!
[3]:
A = getA()
L = A[1:,0] / A[0,0] # Faktoren, diese ergeben nachher die Matrix L
for k in range(1,3):
A[k] -= L[k-1]*A[0] # Matrix R
print(A, "\n\n", L)
[[ 3. 6. 3.]
[ 0. 1. 5.]
[ 0. -9. -3.]]
[0.33333333 2. ]
zweite Zeilenoperation¶
[4]:
L = A[2:,1] / A[1,1]
A[2:3] = A[2:3] - L*A[1]
print(A, "\n\n", L)
[[ 3. 6. 3.]
[ 0. 1. 5.]
[ 0. 0. 42.]]
[-9.]
Damit haben wir die Zeilenstufenform von \(A\) - also die Matrix \(R\) - gefunden. Jetzt nochmal von vorn, aber diesmal erstellen Sie die Matrix \(L\) (ohne die Diagonale) im linken unteren Dreieck gleich mit:
[5]:
A = getA()
print(A)
[[3. 6. 3.]
[1. 3. 6.]
[6. 3. 3.]]
Subtrahieren Sie geeignete Vielfache der ersten Zeile von der zweiten und dritten und Speichern Sie die Faktoren am richtigen Ort
[6]:
A = getA()
# erste Zeilenoperation
A[1:,0] = A[1:,0] / A[0,0]
for k in range(1,3):
A[k,1:] -= A[k,0]*A[0,1:] # Matrix R
print(A, "\n\n")
# zweite Zeilenoperation
A[2:,1] = A[2:,1] / A[1,1]
A[2:,2:] -= A[2:,1]*A[1,2:]
#Ausgabe
print(A)
[[ 3. 6. 3. ]
[ 0.33333333 1. 5. ]
[ 2. -9. -3. ]]
[[ 3. 6. 3. ]
[ 0.33333333 1. 5. ]
[ 2. -9. 42. ]]
Wir extrahieren die Matrizen L und R aus dem Resultat
[7]:
L = np.tril(A, -1) + np.identity(3)
R = np.triu(A)
print("L=\n", L,"\n\nR=\n", R)
L=
[[ 1. 0. 0. ]
[ 0.33333333 1. 0. ]
[ 2. -9. 1. ]]
R=
[[ 3. 6. 3.]
[ 0. 1. 5.]
[ 0. 0. 42.]]
Jetzt testen wir, ob tatsächlich \(LR = A\) gilt:
[8]:
print(L@R - getA()) #sollte die Nullmatrix geben, evtl. bis auf Rundungsfehler ~1e-16
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
Zweite Variante, mit Pivotisierung¶
Jetzt nochmal von vorn, diesmal mit Pivotisierung
[9]:
A = getA()
print(A)
[[3. 6. 3.]
[1. 3. 6.]
[6. 3. 3.]]
Initialisieren Sie einen Indexvektor idx = [0,1,2]. Dieser Vektor enthält am Ende der LU-Zerlegung die (neue) Reihenfolge der Zeilen.
[10]:
idx = np.arange(3)
idxc = idx.copy()
print(idx)
[0 1 2]
Als Vorübung: bestimmen Sie die Zeile \(p\) mit dem betragsmässig grössten Element der ersten Spalte von \(A\). *Hinweis: np.argmax
[11]:
p = np.argmax(np.abs(A[:,0]))
print(p)
2
Tauschen Sie in der Matrix \(A\) die erste Zeile (Index \(0\)) mit der Zeile \(p\) und ebenso den ersten Eintrag von idx mit dem Eintrag \(p\)
[12]:
A = getA()
idx = np.arange(3)
A[[0,p]] = A[[p,0]]
idx[[0,p]]= idx[[p,0]]
print("A=\n", A, "\n\nidx=\n", idx)
A=
[[6. 3. 3.]
[1. 3. 6.]
[3. 6. 3.]]
idx=
[2 1 0]
Jetzt führen Sie die erste Zeilenoperation mit der aktuellen Matrix \(A\) aus und speichern Sie die Faktoren \(L\) in der unteren linken Dreiecksmatrix
[13]:
# erste Zeilenoperation
A[1:,0] = A[1:,0] / A[0,0]
for k in range(1,3):
A[k,1:] -= A[k,0]*A[0,1:] # Matrix R
print(A)
[[6. 3. 3. ]
[0.16666667 2.5 5.5 ]
[0.5 4.5 1.5 ]]
bestimmen Sie die Zeile \(p\) mit dem betragsmässig grössten Element der zweiten Spalte von \(A\), von der zweiten Zeile an
[14]:
p = 1 + np.argmax(np.abs(A[1:,1]))
print(p)
2
Tauschen Sie die entsprechenden Zeilen von \(A\) (das betrifft die Anteile \(L\) und \(R\) gleichermassen)
[15]:
A[[1,p]] = A[[p,1]]
idx[[1,p]]= idx[[p,1]]
print("A=\n", A, "\n\nidx=\n", idx)
A=
[[6. 3. 3. ]
[0.5 4.5 1.5 ]
[0.16666667 2.5 5.5 ]]
idx=
[2 0 1]
Jetzt führen Sie die zweite Zeilenoperation mit der aktuellen Matrix \(A\) aus und speichern Sie wiederum \(L\)
[16]:
# zweite Zeilenoperation
A[2:,1] = A[2:,1] / A[1,1]
A[2:,2:] -= A[2:,1]*A[1,2:]
[17]:
print(A)
[[6. 3. 3. ]
[0.5 4.5 1.5 ]
[0.16666667 0.55555556 4.66666667]]
Jetzt sind wir fertig. Wir extrahieren (testhalber) die Matrizen \(L, R\)
[18]:
L = np.tril(A, -1) + np.identity(3)
R = np.triu(A)
print("L=\n", L,"\n\nR=\n", R)
L=
[[1. 0. 0. ]
[0.5 1. 0. ]
[0.16666667 0.55555556 1. ]]
R=
[[6. 3. 3. ]
[0. 4.5 1.5 ]
[0. 0. 4.66666667]]
jetzt sollte \(L R = A\) sein:
[19]:
print(L@R-getA()[idx])
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]