Solved Task 2a

Schenk_Brandenberger_S8_Aufg2.py

@@ -32,32 +32,43 @@ 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__':
-    A = np.array([
-        [1, -2, 3],
-        [-5, 4, 1],
-        [2, -1, 3]
-    ])
-    b = np.array([
-        [1],
-        [9],
-        [5]
-    ])
+    A = np.array([[1, 2, -1], [4, -2, 6], [3, 1, 0]])
 
-Serie8_Aufg2(A)
+    # Example from Task 1
+    # A = np.array([
+    #     [1, -2, 3],
+    #     [-5, 4, 1],
+    #     [2, -1, 3]
+    # ])
+    # b = np.array([
+    #     [1],
+    #     [9],
+    #     [5]
+    # ])
+
+    [Q,R]=Serie8_Aufg2(A)
+    print("\nQ:\n", Q)
+    print("\nR:\n", R)