Serie 11 and 12

Kuengjoe_S11/Kuengjoe_S10_Aufg3.py

@@ -0,0 +1,39 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Nov 29 14:43:17 2020
+
+Höhere Mathematik 1, Serie 11, Aufgabe 3 (Gerüst)
+
+@author: kuengjoe
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+detail = 1000                           # number of pixels in x and y direction
+maxit = 120                             # maximum n for iterations (influences how detailed the structures are shown when zooming in)
+x_min = -2.0                            # minimum value of x-interval
+x_max = 0.7                             # maximum value of x-interval
+y_min = -1.4                            # minimum value of y-interval
+y_max = 1.4                             # maximum value of y-interval
+
+a = np.linspace(x_min, x_max, detail, dtype=np.float64)  # define real axis [x_min, x_max]
+b = np.linspace(y_min, y_max, detail, dtype=np.float64)  # define imaginary axis [y_min, y_max]
+
+B = np.zeros((detail, detail))          # for color values n
+
+[x, y] = np.meshgrid(a, b)              # to create the complex plane with the axes defined by a and b
+C = np.array(x + y*1j, np.complex128)   # creating the plane
+Z = np.zeros(C.shape, np.complex128)    # initial conditions (first iteration), Z has same dimension as C
+for n in np.arange(1, maxit + 1):       # start iteration
+    Z = Z**2 + C                        # calculating Z
+    expl = np.where(np.abs(Z) > 2)      # finding exploded values (i.e. with an absolute value > 2)
+    Z[expl] = 0                         # removing from iteration
+    C[expl] = 0                         # removing from plane
+    B[expl] = n                         # saving color value n
+
+plt.figure(1)
+B = B/np.max(np.max(B))                 # dividing by max value for correct color
+# display image
+plt.imshow(B, extent=[x_min, x_max, y_min, y_max], origin='lower', interpolation='bilinear')
+plt.show()

Kuengjoe_S11/Serie11_Aufg3_Gerüst.py

@@ -0,0 +1,38 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Nov 29 14:43:17 2020
+
+Höhere Mathematik 1, Serie 11, Aufgabe 3 (Gerüst)
+
+@author: knaa
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+detail = 1000                           #number of pixels in x and y direction
+maxit = 120                             #maximum n for iterations (influences how detailed the structures are shown when zooming in)
+x_min = """???"""                       #minimim value of x-interval
+x_max = """???"""                       #maximum value of x-interval
+y_min = """???"""                       #minimim value of y-interval
+y_max = """???"""                       #maximum value of y-interval
+
+a = np.linspace("""???""", detail, dtype=np.float64)  #define real axis [x_min, x_max]
+b = np.linspace("""???""", detail, dtype=np.float64)  #define imaginary axis [y_min, y_max]
+
+B = np.zeros((detail, detail))          #for color values n 
+
+[x,y] = np.meshgrid("""???""")          #to create the complex plane with the axes defined by a and b
+C = np.array(x + y*1j, np.complex128)   #creating the plane
+Z = np.zeros("""???""", np.complex128)  #initial conditions (first iteration), Z has same dimension as C
+for n in np.arange(1, maxit + 1):       #start iteration
+  Z = """???"""                         #calculating Z
+  expl = np.where("""???""")            #finding exploded values (i.e. with an absolute value > 2)
+  Z[expl] = 0                           #removing from iteration
+  C[expl] = 0                           #removing from plane
+  B[expl] = """???"""                   #saving color value n
+
+plt.figure(1)
+B = B/np.max(np.max(B))                 #dividing by max value for correct color
+#display image
+plt.imshow(B, extent=[x_min, x_max, y_min, y_max], origin='lower', interpolation='bilinear')

Kuengjoe_S12/Kuengjoe_S12_Aufg4.py

@@ -0,0 +1,48 @@
+import numpy as np
+
+def Kuengjoe_S12_Aufg4(in_matrix: np.ndarray, iteration: int):
+    current_iteration = np.array(in_matrix, dtype=float)
+    dimension = current_iteration.shape[0]
+    acc_orthogonal_matrix = np.eye(dimension, dtype=float)
+
+    for iteration in range(iteration):
+        q_matrix, r_matrix = np.linalg.qr(current_iteration)
+        current_iteration = r_matrix @ q_matrix
+        acc_orthogonal_matrix = acc_orthogonal_matrix @ q_matrix
+    return current_iteration, acc_orthogonal_matrix
+
+
+if __name__ == "__main__":
+    #4a)
+    test_matrix = np.array([
+        [1, -2, 0],
+        [2, 0, 1],
+        [0, -2, 1]
+    ], dtype=float)
+
+    ak_1, pk_1 = Kuengjoe_S12_Aufg4(test_matrix, 1)
+    print("A1 =\n", np.round(ak_1, 6))
+    print("P1 =\n", np.round(pk_1, 6))
+
+    #4b)
+    symmetric_matrix = np.array([
+        [6, 1, 2, 1, 2],
+        [1, 5, 0, 2, -1],
+        [2, 0, 5, -1, 0],
+        [1, 2, -1, 6, 1],
+        [2, -1, 0, 1, 7]
+    ], dtype=float)
+
+    ak_100, pk_100 = Kuengjoe_S12_Aufg4(symmetric_matrix, 100)
+
+    orthogonality_residual = np.linalg.norm(pk_100.T @ pk_100 - np.eye(5))
+    print("||P^T P - I|| =", orthogonality_residual)
+
+    approx_eigenvalues_from_qr = np.diag(ak_100)
+    print("Eigenvalues approx (diag(Ak)) =", approx_eigenvalues_from_qr)
+
+    #4c)
+
+    eigenvalues_numpy, eigenvectors_numpy = np.linalg.eig(symmetric_matrix)
+    print(eigenvalues_numpy)
+

Kuengjoe_S12/Kuengjoe_S12_Aufg5.py

@@ -0,0 +1,68 @@
+
+import numpy as np
+
+def Kuengjoe_S12_Aufg5(input_matrix: np.ndarray,
+                       initial_vector: np.ndarray,
+                       tolerance: float = 1e-4,
+                       max_iterations: int = 1000):
+    current_normalized_vector = initial_vector.astype(float)
+    current_normalized_vector = current_normalized_vector / np.linalg.norm(current_normalized_vector, ord=2)
+
+    number_of_performed_iterations = 0
+    last_difference_norm = None
+
+    for iteration_index in range(max_iterations):
+        multiplied_vector = input_matrix @ current_normalized_vector
+        next_normalized_vector = multiplied_vector / np.linalg.norm(multiplied_vector, ord=2)
+
+        difference_norm = np.linalg.norm(next_normalized_vector - current_normalized_vector, ord=2)
+        number_of_performed_iterations = iteration_index + 1
+        last_difference_norm = difference_norm
+
+        if difference_norm < tolerance:
+            current_normalized_vector = next_normalized_vector
+            break
+
+        current_normalized_vector = next_normalized_vector
+
+    rayleigh_quotient_eigenvalue_estimate = float(
+        current_normalized_vector.T @ (input_matrix @ current_normalized_vector)
+    )
+
+    return rayleigh_quotient_eigenvalue_estimate, current_normalized_vector, number_of_performed_iterations, last_difference_norm
+
+
+if __name__ == "__main__":
+    matrix_a = np.array([
+        [1,  1, 0],
+        [3, -1, 2],
+        [2, -1, 3]
+    ], dtype=float)
+
+    initial_vector_x0 = np.array([1, 0, 0], dtype=float)
+
+    dominant_eigenvalue_estimate, dominant_eigenvector_estimate, iteration_count, final_difference_norm = (
+        Kuengjoe_S12_Aufg5(matrix_a, initial_vector_x0, tolerance=1e-4)
+    )
+
+    print("von-Mises result")
+    print("iterations =", iteration_count)
+    print("final ||x_{k+1}-x_k||2 =", final_difference_norm)
+    print("dominant eigenvalue (Rayleigh) =", dominant_eigenvalue_estimate)
+    print("dominant eigenvector =", dominant_eigenvector_estimate)
+
+    eigenvalues_numpy, eigenvectors_numpy = np.linalg.eig(matrix_a)
+    index_of_dominant_eigenvalue = int(np.argmax(np.abs(eigenvalues_numpy)))
+
+    dominant_eigenvalue_numpy = float(eigenvalues_numpy[index_of_dominant_eigenvalue])
+    dominant_eigenvector_numpy = eigenvectors_numpy[:, index_of_dominant_eigenvalue].astype(float)
+    dominant_eigenvector_numpy = dominant_eigenvector_numpy / np.linalg.norm(dominant_eigenvector_numpy, ord=2)
+
+    print("\nnp.linalg.eig check")
+    print("dominant eigenvalue (eig) =", dominant_eigenvalue_numpy)
+    print("dominant eigenvector (eig) =", dominant_eigenvector_numpy)
+
+    # confronto robusto (segno +/-)
+    difference_same_sign = np.linalg.norm(dominant_eigenvector_estimate - dominant_eigenvector_numpy, ord=2)
+    difference_opposite_sign = np.linalg.norm(dominant_eigenvector_estimate + dominant_eigenvector_numpy, ord=2)
+    print("\nvector agreement (min with +/- sign) =", min(difference_same_sign, difference_opposite_sign))