diff --git a/Schenk_Brandenberger_S10_Aufg3.py b/Schenk_Brandenberger_S10_Aufg3.py new file mode 100644 index 0000000..8651b46 --- /dev/null +++ b/Schenk_Brandenberger_S10_Aufg3.py @@ -0,0 +1,66 @@ +import numpy as np + + +def Single_Jacobi_Iteration(B, c, xn): + return B @ xn + c + + +def Single_Gauss_Seidel_Iteration(B, c, xn): + return B @ xn + c + + +def A_Priori_Estimation(B, x0, x1, tol): + return np.log((tol / np.linalg.norm((x1 - x0), np.inf) * (1 - np.linalg.norm(B, np.inf)))) / np.log( + np.linalg.norm(B, np.inf)) + + +def A_Posteriori_Estimation(B, xn, xn_minus_one): + return (np.linalg.norm(B, np.inf) / (1 - np.linalg.norm(B, np.inf))) * np.linalg.norm(xn - xn_minus_one) + + +def Schenk_Brandenberger_S10_Aufg3a(A, b, x0, tol, opt): + A = A.astype("float64") + b = b.astype("float64") + x0 = x0.astype("float64") + x1 = np.copy(x0) + L = np.tril(A, k=-1) + D = np.diag(np.diag(A)) + R = np.triu(A, k=1) + n = 1 + + if opt == 0: + B = -np.linalg.inv(D) @ (L + R) + c = np.linalg.inv(D) @ b + x1 = Single_Jacobi_Iteration(B, c, x0) + elif opt == 1: + B = -np.linalg.inv(D + L) @ R + c = np.linalg.inv(D + L) @ b + x1 = Single_Gauss_Seidel_Iteration(B, c, x0) + xn = np.copy(x1) + xn_minus_one = np.copy(x0) + # a-priori Abschätzung + n2 = A_Priori_Estimation(B, x0, x1, tol) + + while tol < A_Posteriori_Estimation(B, xn, xn_minus_one): + xn_minus_one = np.copy(xn) + if opt == 0: + xn = Single_Jacobi_Iteration(B, c, xn) + elif opt == 1: + xn = Single_Gauss_Seidel_Iteration(B, c, xn) + n = n + 1 + return (xn, n, n2) + + +A = np.array([[8, 5, 2], [5, 9, 1], [4, 2, 7]]) +b = np.array([[19], [5], [34]]) +x0 = np.array([[1], [-1], [3]]) +opt = 1 # 0 is Jacobi; 1 is Gauss-Seidel + +[xn, n, n2] = Schenk_Brandenberger_S10_Aufg3a(A, b, x0, 0.0001, opt) + +print("Solution vector xn:") +print(xn) +print("Number of iterations:") +print(n) +print("Used steps with a priori:") +print(n2) \ No newline at end of file diff --git a/Schenk_Brandenberger_S10_Aufg3b_3c.py b/Schenk_Brandenberger_S10_Aufg3b_3c.py new file mode 100644 index 0000000..19d5b44 --- /dev/null +++ b/Schenk_Brandenberger_S10_Aufg3b_3c.py @@ -0,0 +1,77 @@ +import numpy as np +import time +from Schenk_Brandenberger_S10_Aufg3 import * +from Gauss_Algorithm import * +import matplotlib.pyplot as plt + +#3b +dim = 3000 +A = np.diag( np.diag( np.ones( ( dim , dim ) )*4000 ) )+ np.ones( ( dim , dim ) ) +dum1 = np.arange( 1 , np.int32( dim /2+1) , dtype = np.float64 ).reshape( ( np.int32( dim / 2 ) , 1 ) ) +dum2 = np.arange( np.int32( dim / 2 ) ,0 , -1 , dtype=np.float64 ).reshape( ( np.int32( dim / 2 ) , 1 ) ) +x = np.append( dum1 , dum2 , axis=0) +b = A@x +x0 = np.zeros( ( dim , 1 ) ) +tol = 1e-4 + +startLinalgSolve = time.time() +x_linalg_solve = np.linalg.solve(A,b) +stoppLinalgSolve = time.time() + +startJacobi = time.time() +[x_jacobi, n, n2] = Schenk_Brandenberger_S10_Aufg3(A,b,x0,tol,0) +stoppJacobi = time.time() + +startGaussSeidel = time.time() +[x_gauss_seidel, n, n2] = Schenk_Brandenberger_S10_Aufg3(A,b,x0,tol,0) +stoppGaussSeidel = time.time() + +startGauss = time.time() +x_gauss = Witschi_Floian_S6_Aufg2(A,b)[2] +stoppGauss = time.time() + +print("***********Time Estimation***********") +print("Time for np.linalg.solve:") +print(stoppLinalgSolve-startLinalgSolve) +print("Time for Jacobi:") +print(stoppJacobi-startJacobi) +print("Time for Gauss-Seidel:") +print(stoppGaussSeidel-startGaussSeidel) +print("Time for Gauss:") +print(stoppGauss-startGauss) + + +#In this algorithm for Jacobi and Gauss-Seidel is the B calculated every single time +""" Time for np.linalg.solve: +0.6146464347839355 +Time for Jacobi: +319.44120621681213 +Time for Gauss-Seidel: +307.4256772994995 +Time for Gauss: +95.12803220748901 """ + +#In this algorithm for Jacobi and Gauss-Seidel is the B calculated one time +""" Time for np.linalg.solve: +0.5016989707946777 +Time for Jacobi: +18.02965211868286 +Time for Gauss-Seidel: +17.64187240600586 +Time for Gauss: +105.31756782531738 """ + +print(x_gauss) + +#3c +#x_axis = np.array(["x0","x1","x2"]) +#plt.plot(x_axis, x_linalg_solve) +x_axis = np.arange(dim) +plt.plot(x_axis, x_gauss_seidel-x) +plt.plot(x_axis, x_jacobi-x) +plt.plot(x_axis, x_gauss-x) +plt.plot(x_axis, x_linalg_solve-x) +plt.legend(["Gauss Seidel", "Jacobi", "Gauss", "Linalg"]) +plt.show() + +#Das Gauss-Verfahren ist genauer wie das Jacobi und Gauss-Seidel Verfahren \ No newline at end of file