Created
September 12, 2025 16:15
-
-
Save carlosal1015/2ab6f3d237aafe32cab1673434758d95 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| """Ejemplo 8.6.12 | |
| An Introduction to Numerical Analysis | |
| https://math.science.cmu.ac.th/docs/qNA2556/ref_na/Katkinson.pdf | |
| """ | |
| import numpy as np | |
| import logging | |
| from time import perf_counter | |
| from tabulate import tabulate | |
| logging.basicConfig(level=logging.INFO) | |
| np.set_printoptions(precision=12) | |
| def ForwardGaussSeidel( | |
| systemMatrix, rightHandSide, startVector, maxIter, epsilon, exact_solution | |
| ): | |
| """ | |
| Performs the forward Gauss-Seidel iteration to solve a linear system. | |
| Returns: | |
| - The final solution vector. | |
| - The number of iterations performed. | |
| - A list of lists containing the data for the iteration table. | |
| """ | |
| diagonal = np.diag(systemMatrix) | |
| solution = startVector.copy() | |
| dim = solution.shape[0] | |
| table_data = [] | |
| # Iteration 0 | |
| error_vec = exact_solution - solution | |
| error_norm = np.linalg.norm(error_vec, np.inf) | |
| table_data.append([0, *solution, error_norm, ""]) | |
| prev_error_norm = error_norm | |
| for iter_count in range(maxIter): | |
| for i in range(dim): | |
| row = systemMatrix[i] | |
| # The trick of setting solution[i] to 0 is to exclude the diagonal element's contribution | |
| # from the dot product. This works for Gauss-Seidel as the solution vector is updated in place. | |
| solution[i] = 0.0 | |
| solution[i] = (rightHandSide[i] - row.dot(solution)) / diagonal[i] | |
| error_vec = exact_solution - solution | |
| error_norm = np.linalg.norm(error_vec, np.inf) | |
| ratio = error_norm / prev_error_norm if prev_error_norm != 0 else 0.0 | |
| table_data.append([iter_count + 1, *solution, error_norm, ratio]) | |
| prev_error_norm = error_norm | |
| residual = rightHandSide - systemMatrix.dot(solution) | |
| residual_norm = np.linalg.norm(residual, np.inf) | |
| if residual_norm < epsilon: | |
| return solution, iter_count + 1, table_data | |
| return solution, maxIter, table_data | |
| A = np.array( | |
| object=[ | |
| [10, 3, 1], | |
| [2, -10, 3], | |
| [1, 3, 10], | |
| ], | |
| dtype=float, | |
| ) | |
| b = np.array(object=[14, -5, 14], dtype=float) | |
| start_direct = perf_counter() | |
| # Set initial guess, max iterations and tolerance | |
| x0 = np.zeros_like(b) | |
| max_iter = 100 | |
| epsilon = 1e-9 | |
| # For verification, solve with numpy's direct solver | |
| x_exact = np.linalg.solve(a=A, b=b) | |
| logging.info(f"Exact solution using np.linalg.solve: {x_exact}") | |
| solution, iterations, table_data = ForwardGaussSeidel( | |
| systemMatrix=A, | |
| rightHandSide=b, | |
| startVector=x0, | |
| maxIter=max_iter, | |
| epsilon=epsilon, | |
| exact_solution=x_exact, | |
| ) | |
| logging.info( | |
| f"Tiempo total de ForwardGaussSeidel: {perf_counter() - start_direct:.4f} segundos" | |
| ) | |
| logging.info(f"Solution found: {solution}") | |
| if iterations < max_iter: | |
| logging.info(f"Converged in {iterations} iterations.") | |
| else: | |
| logging.warning(f"Did not converge in {iterations} iterations.") | |
| # Generate and print table | |
| dim = len(x0) | |
| headers = ["m", *[f"x{i+1}(m)" for i in range(dim)], "||e(m)||inf", "Ratio"] | |
| float_formats = (".0f", *[".6f"] * dim, ".2E", ".3f") | |
| print("\nNumerical results for the Gauss-Seidel method") | |
| print(tabulate(table_data, headers=headers, floatfmt=float_formats)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment