|
| 1 | +import numpy as np |
| 2 | + |
| 3 | +matrix = np.array( |
| 4 | + [ |
| 5 | + [5.0, -5.0, -3.0, 4.0, -11.0], |
| 6 | + [1.0, -4.0, 6.0, -4.0, -10.0], |
| 7 | + [-2.0, -5.0, 4.0, -5.0, -12.0], |
| 8 | + [-3.0, -3.0, 5.0, -5.0, 8.0], |
| 9 | + ], |
| 10 | + dtype=float, |
| 11 | +) |
| 12 | + |
| 13 | + |
| 14 | +def solve_linear_system(matrix: np.ndarray) -> np.ndarray: |
| 15 | + """ |
| 16 | + Solve a linear system of equations using Gaussian elimination with partial pivoting |
| 17 | +
|
| 18 | + Args: |
| 19 | + - matrix: Coefficient matrix with the last column representing the constants. |
| 20 | +
|
| 21 | + Returns: |
| 22 | + - Solution vector. |
| 23 | +
|
| 24 | + Raises: |
| 25 | + - ValueError: If the matrix is not correct (i.e., singular). |
| 26 | +
|
| 27 | + https://courses.engr.illinois.edu/cs357/su2013/lect.htm Lecture 7 |
| 28 | +
|
| 29 | + Example: |
| 30 | + >>> A = np.array([[2, 1, -1], [-3, -1, 2], [-2, 1, 2]], dtype=float) |
| 31 | + >>> B = np.array([8, -11, -3], dtype=float) |
| 32 | + >>> solution = solve_linear_system(np.column_stack((A, B))) |
| 33 | + >>> np.allclose(solution, np.array([2., 3., -1.])) |
| 34 | + True |
| 35 | + >>> solve_linear_system(np.array([[0, 0], [0, 0]], dtype=float)) |
| 36 | + array([nan, nan]) |
| 37 | + """ |
| 38 | + ab = np.copy(matrix) |
| 39 | + num_of_rows = ab.shape[0] |
| 40 | + num_of_columns = ab.shape[1] - 1 |
| 41 | + x_lst: list[float] = [] |
| 42 | + |
| 43 | + # Lead element search |
| 44 | + for column_num in range(num_of_rows): |
| 45 | + for i in range(column_num, num_of_columns): |
| 46 | + if abs(ab[i][column_num]) > abs(ab[column_num][column_num]): |
| 47 | + ab[[column_num, i]] = ab[[i, column_num]] |
| 48 | + if ab[column_num, column_num] == 0.0: |
| 49 | + raise ValueError("Matrix is not correct") |
| 50 | + else: |
| 51 | + pass |
| 52 | + if column_num != 0: |
| 53 | + for i in range(column_num, num_of_rows): |
| 54 | + ab[i, :] -= ( |
| 55 | + ab[i, column_num - 1] |
| 56 | + / ab[column_num - 1, column_num - 1] |
| 57 | + * ab[column_num - 1, :] |
| 58 | + ) |
| 59 | + |
| 60 | + # Upper triangular matrix |
| 61 | + for column_num in range(num_of_rows): |
| 62 | + for i in range(column_num, num_of_columns): |
| 63 | + if abs(ab[i][column_num]) > abs(ab[column_num][column_num]): |
| 64 | + ab[[column_num, i]] = ab[[i, column_num]] |
| 65 | + if ab[column_num, column_num] == 0.0: |
| 66 | + raise ValueError("Matrix is not correct") |
| 67 | + else: |
| 68 | + pass |
| 69 | + if column_num != 0: |
| 70 | + for i in range(column_num, num_of_rows): |
| 71 | + ab[i, :] -= ( |
| 72 | + ab[i, column_num - 1] |
| 73 | + / ab[column_num - 1, column_num - 1] |
| 74 | + * ab[column_num - 1, :] |
| 75 | + ) |
| 76 | + |
| 77 | + # Find x vector (Back Substitution) |
| 78 | + for column_num in range(num_of_rows - 1, -1, -1): |
| 79 | + x = ab[column_num, -1] / ab[column_num, column_num] |
| 80 | + x_lst.insert(0, x) |
| 81 | + for i in range(column_num - 1, -1, -1): |
| 82 | + ab[i, -1] -= ab[i, column_num] * x |
| 83 | + |
| 84 | + # Return the solution vector |
| 85 | + return np.asarray(x_lst) |
| 86 | + |
| 87 | + |
| 88 | +if __name__ == "__main__": |
| 89 | + from doctest import testmod |
| 90 | + from pathlib import Path |
| 91 | + |
| 92 | + testmod() |
| 93 | + file_path = Path(__file__).parent / "matrix.txt" |
| 94 | + try: |
| 95 | + matrix = np.loadtxt(file_path) |
| 96 | + except FileNotFoundError: |
| 97 | + print(f"Error: {file_path} not found. Using default matrix instead.") |
| 98 | + |
| 99 | + # Example usage: |
| 100 | + print(f"Matrix:\n{matrix}") |
| 101 | + print(f"{solve_linear_system(matrix) = }") |
![pre-commit-ci[bot]](https://avatars.githubusercontent.com/in/68672?v=4&size=40){kind=avatar}
0 commit comments