Solved Task 2b

Schenk_Brandenberger_S8_Aufg2.py

@@ -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
+    # A = np.array([[1, 2, -1], [4, -2, 6], [3, 1, 0]])
+    # b = np.array([
+    #     [9],
+    #     [-4],
+    #     [9]
+    # ])
+
+    # Example from Task 1
     A = np.array([
         [1, -2, 3],
         [-5, 4, 1],
@@ -60,4 +75,22 @@ if __name__ == '__main__':
         [5]
     ])
 
-Serie8_Aufg2(A)
+    [Q,R]=Serie8_Aufg2(A)
+
+    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)