Serie 10

Kuengjoe_S10/Kuengjoe_S10_Aufg3a.py

@@ -0,0 +1,58 @@
+import numpy as np
+from numpy.linalg import norm
+
+def Kuengjoe_S10_Aufg3a(A, b, x0, tol, opt):
+    dim = A.shape[0]
+
+    D_vec= np.diag(A)
+    D = np.diag(D_vec)
+    L = np.tril(A, -1)
+    U = np.triu(A, 1)
+
+    if  opt == 'Jacobi':
+        inv_D = np.diag(1.0 / D_vec)
+
+        B = -np.dot(inv_D, (L + U))
+        c = np.dot(inv_D, b)
+
+        del inv_D, L, U, D,
+
+    elif opt == 'Gauss-Seidel':
+        DL = D + L
+        inv_DL = np.linalg.inv(DL)
+
+        B = -np.dot(inv_DL, U)
+        c = np.dot(inv_DL, b)
+
+        del inv_DL, U, D, L,
+
+    else:
+        raise ValueError("Ungültige Option. Wählen Sie 'Jacobi' oder 'Gauss-Seidel'.")
+
+    norm_B = norm(B, np.inf)
+
+    x = x0.copy()
+    x1 = np.dot(B, x) + c
+    diff_x1_x0 = norm(x1 - x0, np.inf)
+
+
+    if norm_B < 1.0 and diff_x1_x0 > 0:
+        numerator = np.log((tol * (1 - norm_B)) / diff_x1_x0)
+        denominator = np.log(norm_B)
+        n2 = int(np.ceil(numerator / denominator))
+
+    else:
+        n2 = -1
+    n = 0
+    current_diff = tol +1
+
+    x = x1
+    n = 1
+
+    while current_diff > tol and n < 10000:
+        x_new = np.dot(B, x) + c
+        current_diff = norm(x_new - x, np.inf)
+        x = x_new
+        n += 1
+    xn = x
+    return xn, n, n2

Kuengjoe_S10/Kuengjoe_S10_Aufg3b.py

@@ -0,0 +1,77 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import time
+from Kuengjoe_S10_Aufg3a import Kuengjoe_S10_Aufg3a
+
+print("Datenerzeugung läuft (dim=3000)... bitte warten.")
+dim = 3000
+
+A = np.diag(np.diag(np.ones((dim, dim)) * 4000)) + np.ones((dim, dim))
+
+dum1 = np.arange(1, int(dim/2 + 1), dtype=np.float64).reshape((int(dim/2), 1))
+dum2 = np.arange(int(dim/2), 0, -1, dtype=np.float64).reshape((int(dim/2), 1))
+
+x_exact = np.append(dum1, dum2, axis=0)
+
+b = np.dot(A, x_exact)
+
+x0 = np.zeros((dim, 1))
+tol = 1e-4
+
+
+# A) Jacobi-Verfahren
+print("\nStarte Jacobi-Verfahren...")
+start_time = time.time()
+xn_jac, n_jac, n2_jac = Kuengjoe_S10_Aufg3a(A, b, x0, tol, 'Jacobi')
+time_jac = time.time() - start_time
+print(f"Jacobi: {time_jac:.4f} sek | Iterationen: {n_jac} (A-priori: {n2_jac})")
+
+# B) Gauss-Seidel-Verfahren
+print("\nStarte Gauss-Seidel-Verfahren...")
+start_time = time.time()
+xn_gs, n_gs, n2_gs = Kuengjoe_S10_Aufg3a(A, b, x0, tol, 'Gauss-Seidel')
+time_gs = time.time() - start_time
+print(f"Gauss-Seidel: {time_gs:.4f} sek | Iterationen: {n_gs} (A-priori: {n2_gs})")
+
+
+print("\nStarte Gauss-Verfahren (Numpy linalg.solve)...")
+start_time = time.time()
+xn_gauss = np.linalg.solve(A, b)
+time_gauss = time.time() - start_time
+print(f"Gauss (Numpy): {time_gauss:.4f} sek")
+
+
+"""
+KOMMENTAR ZU TEIL B:
+Ein manuelles Gauss-Verfahren (mit Python-Schleifen) wäre bei dim=3000 extrem langsam.
+Da ich hier aber np.linalg.solve (optimierter C-Code) benutze, ist Gauss hier
+schneller als Gauss-Seidel.
+
+A manual Gauss method (using Python loops) would be extremely slow for dim=3000.
+However, since I use np.linalg.solve (optimized C code) here, Gauss appears 
+faster than Gauss-Seidel in this test.
+"""
+
+
+err_jac = np.abs(xn_jac - x_exact)
+err_gs = np.abs(xn_gs - x_exact)
+err_gauss = np.abs(xn_gauss - x_exact)
+
+plt.figure(figsize=(10, 6))
+
+plt.plot(err_jac, label='Fehler Jacobi', color='blue', alpha=0.7)
+plt.plot(err_gs, label='Fehler Gauss-Seidel', color='red', alpha=0.7)
+plt.plot(err_gauss, label='Fehler Gauss (Exakt)', color='green', alpha=0.7)
+
+plt.title(f"Vergleich der absoluten Fehler (Dim={dim})")
+plt.xlabel("Index des Vektorelements")
+plt.ylabel("Absoluter Fehler |x_berechnet - x_exakt|")
+plt.yscale('log')
+plt.legend()
+plt.grid(True, which="both", ls="-", alpha=0.4)
+plt.show()
+
+print("\n--- Beobachtungen (Teil C) ---")
+print("1. Der Fehler des Gauss-Verfahrens (direkter Löser) liegt im Bereich der Maschinengenauigkeit (ca. 1e-15).")
+print("2. Die iterativen Verfahren (Jacobi, GS) stoppen, sobald die Fehlertoleranz (tol=1e-4) erreicht ist.")
+print("3. Gauss-Seidel konvergiert typischerweise schneller (weniger Iterationen) als Jacobi.")