Merge remote-tracking branch 'origin/master'

Schenk_Brandenberger_S8_Aufg2.py

@@ -15,10 +15,10 @@ Example: A = np.array([[1,2,-1],[4,-2,6],[3,1,0]])
 @author: knaa
 """
 import numpy as np
+import timeit
 
+# Aufgabe 2a
 def Serie8_Aufg2(A):
-    unknown = 0
-
     
     A = np.copy(A)                       #necessary to prevent changes in the original matrix A_in
     A = A.astype('float64')              #change to float
@@ -32,23 +32,38 @@ def Serie8_Aufg2(A):
     R = A
 
     for j in np.arange(0,n-1):
-        a = np.copy(unknown).reshape(n-j,1)
-        e = np.eye(unknown)[:,0].reshape(n-j,1)
+        # erzeuge Nullen in R in der j-ten Spalte unterhalb der Diagonalen:
+        #a = (Q @ A)[j:,j:][:,0]
+        a = np.copy((Q @ A)[j:,j:][:,0]).reshape(n-j,1)
+        e = np.eye(n-j)[:,0].reshape(n-j,1)
         length_a = np.linalg.norm(a)
-        if a[0] >= 0: sig = unknown
-        else: sig = unknown
-        v = unknown
-        u = unknown
-        H = unknown
-        Qi = np.eye(n)
-        Qi[j:,j:] = unknown
-        R = unknown
-        Q = unknown
+        if a[0] >= 0: sig = 1
+        else: sig = -1
+        v = a + sig * length_a * e          # vj := aj + sign(a1j) · |aj| · ej
+        u = 1 / np.linalg.norm(v) * v       # uj := 1/|vj|*vj
+        ut = u.T
+        ua = u @ ut
+        ub = 2 * ua
+        x = np.eye(n)
+        H = np.eye(n-j) - (2 * (u @ u.T))   # Hj := In − 2u1u1T bestimme die (n − j + 1) × (n − j + 1) Householder-Matrix Hj
+        Qj = np.eye(n)
+        Qj[j:,j:] = H                       # erweitere Hi durch einen Ii−1 Block links oben zur n × n Matrix Qi
+        R = Qj @ R                          # R := Qj · R
+        Q = Q @ Qj.T                        # Q := Q · QjT
         
     return(Q,R)
 
 
 if __name__ == '__main__':
+    # Beispiel aus Skript
+    # A = np.array([[1, 2, -1], [4, -2, 6], [3, 1, 0]])
+    # b = np.array([
+    #     [9],
+    #     [-4],
+    #     [9]
+    # ])
+
+    # Beispiel aus Aufgabe 1
     A = np.array([
         [1, -2, 3],
         [-5, 4, 1],
@@ -60,4 +75,46 @@ if __name__ == '__main__':
         [5]
     ])
 
-Serie8_Aufg2(A)
+    [Q,R]=Serie8_Aufg2(A)
+
+    # Aufgabe 2b
+    n = len(b) - 1
+    QTb = Q.T @ b
+    result = [0 for i in range(n+1)]
+    row = n
+    while row >= 0:
+        value = QTb[row][0]
+        column = n
+        while column > row:
+            value -= R[row][column] * result[column]
+            column -= 1
+        value = value / R[row,row]
+        result[row] = value
+        row -= 1
+
+    print("\nQ:\n", Q)
+    print("\nR:\n", R)
+    print("\Result:\n", result)
+
+    # Aufgabe 2c
+    t1 = timeit.repeat("Serie8_Aufg2(A)", "from __main__ import Serie8_Aufg2, A", number=100)
+    t2 = timeit.repeat("np.linalg.qr(A)", "from __main__ import np, A", number=100)
+    avg_t1 = np.average(t1) / 100
+    avg_t2 = np.average(t2) / 100
+
+    print("Geschwindigkeit mit 3x3 Matrix:")
+    print("Benötigte Zeit mit eigener Funktion:", avg_t1)
+    print("Benötigte Zeit mit Numpy:", avg_t2)
+
+    # Aufgabe 2d
+    Test = np.random.rand(100,100)
+    t1 = timeit.repeat("Serie8_Aufg2(Test)", "from __main__ import Serie8_Aufg2, Test", number=100)
+    t2 = timeit.repeat("np.linalg.qr(Test)", "from __main__ import np, Test", number=100)
+    avg_t1 = np.average(t1) / 100
+    avg_t2 = np.average(t2) / 100
+
+    print("Geschwindigkeit mit 100x100 Matrix:")
+    print("Benötigte Zeit mit eigener Funktion:", avg_t1)
+    print("Benötigte Zeit mit Numpy:", avg_t2)
+
+    # Die von Numpy zur Verfügung gestellt Funktion ist wesentlich effizienter.

Schenk_Brandenberger_S8_Aufg3.py

@@ -0,0 +1,79 @@
+import numpy as np
+
+
+def switchRows(matrix, row1, row2):
+    matrix[[row1, row2]] = matrix[[row2, row1]]
+    return matrix
+
+def Schenk_Brandenberger_S6_Aufg2(A, b):
+
+    def calculateRow(A, b, row, column):
+        b[row] = [b[row][0] - (A[row][column] / A[column][column]) * b[column][0]]
+        A[row] = [(A[row][i] - (A[row][column] / A[column][column]) * A[column][i]) for i in range(len(A[row]))]
+        return A, b
+
+    # Erstelle obere Dreiecksmatrix
+    countRowSwitch = 0
+    columnsToEdit = []
+    for row in range(1, len(A)):
+        columnsToEdit.append(row - 1)
+        for column in columnsToEdit:
+            if(A[row-1][column] == 0 and (row - 1 == column)):
+                rowToSwitch = row
+                if(row == 1):
+                    while(A[rowToSwitch][column] == 0):
+                        if(len(A) > rowToSwitch + 1):
+                            rowToSwitch += 1
+                        else:
+                            return "Matrix ist nicht regulär!"
+                A = switchRows(A, row - 1, rowToSwitch)
+                b = switchRows(b, row - 1, rowToSwitch)
+                countRowSwitch += 1
+            else:
+                A, b = calculateRow(A, b, row, column)
+    #print("\nObere Dreiecksmatrix A:\n", A, "\nb:\n", b)
+
+    # Rückwärtseinsetzen
+    columnsToEdit = []
+    for row in range((len(A) - 2), -1, -1):
+        columnsToEdit.append(row + 1)
+        for column in columnsToEdit:
+            A, b = calculateRow(A, b, row, column)
+        row -= 1
+    #print("\nA:\n", A, "\nb:\n", b)
+
+    det = 1
+    result = []
+    for i in range(len(A)):
+        result.append(b[i][0] / A[i][i])
+        det *= A[i][i]
+    if(countRowSwitch % 2 == 1):
+        det *= (-1)
+    return result, det
+
+
+
+if __name__ == '__main__':
+    # Zahlen von Serie 7 Aufgabe 1c
+    # A = np.array([[20000.0, 30000.0, 10000.0],
+    #                [10000.0, 17000.0, 6000.0],
+    #                [2000.0, 3000.0, 2000.0]])
+    #
+    # b = np.array([[5720000.0],
+    #                [3300000.0],
+    #                [836000.0]])
+
+    # Aufgabe 3c
+    A = np.array([[20000.0 - 100.0, 30000.0 - 100.0, 10000.0 - 100.0],
+                   [10000.0 - 100.0, 17000.0 - 100.0, 6000.0 - 100.0],
+                   [2000.0 - 100.0, 3000.0 - 100.0, 2000.0 - 100.0]])
+
+    b = np.array([[5720000.0 + 100000.0],
+                   [3300000.0 + 100000.0],
+                   [836000.0 + 100000.0]])
+
+    result, det = Schenk_Brandenberger_S6_Aufg2(A, b)
+    print("Ergebnis:")
+    for i in range(len(result)):
+        print("x" + str(i) + ": " + str(result[i]))
+    print("Determinante: " + str(det))