diff --git a/Kuengjoe_S11/Kuengjoe_S10_Aufg3.py b/Kuengjoe_S11/Kuengjoe_S10_Aufg3.py new file mode 100644 index 0000000..ced6523 --- /dev/null +++ b/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() diff --git a/Kuengjoe_S11/Serie11_Aufg3_Gerüst.py b/Kuengjoe_S11/Serie11_Aufg3_Gerüst.py new file mode 100644 index 0000000..3e41dd6 --- /dev/null +++ b/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') diff --git a/Kuengjoe_S12/Kuengjoe_S12_Aufg4.py b/Kuengjoe_S12/Kuengjoe_S12_Aufg4.py new file mode 100644 index 0000000..9b55227 --- /dev/null +++ b/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) + diff --git a/Kuengjoe_S12/Kuengjoe_S12_Aufg5.py b/Kuengjoe_S12/Kuengjoe_S12_Aufg5.py new file mode 100644 index 0000000..5b43a46 --- /dev/null +++ b/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)) \ No newline at end of file