In: Advanced Math
Write a function or script that will solve linear systems of any size by Gaussian elimination with partial pivoting in Python.
# The 'gauss' function takes two matrices, 'a' and 'b', with 'a' square, and it return the determinant of 'a' and a matrix 'x' such that a*x = b. # If 'b' is the identity, then 'x' is the inverse of 'a'. import copy from fractions import Fraction def gauss(a, b): a = copy.deepcopy(a) b = copy.deepcopy(b) n = len(a) p = len(b[0]) det = 1 for i in range(n - 1): k = i for j in range(i + 1, n): if abs(a[j][i]) > abs(a[k][i]): k = j if k != i: a[i], a[k] = a[k], a[i] b[i], b[k] = b[k], b[i] det = -det for j in range(i + 1, n): t = a[j][i]/a[i][i] for k in range(i + 1, n): a[j][k] -= t*a[i][k] for k in range(p): b[j][k] -= t*b[i][k] for i in range(n - 1, -1, -1): for j in range(i + 1, n): t = a[i][j] for k in range(p): b[i][k] -= t*b[j][k] t = 1/a[i][i] det *= a[i][i] for j in range(p): b[i][j] *= t return det, b def zeromat(p, q): return [[0]*q for i in range(p)] def matmul(a, b): n, p = len(a), len(a[0]) p1, q = len(b), len(b[0]) if p != p1: raise ValueError("Incompatible dimensions") c = zeromat(n, q) for i in range(n): for j in range(q): c[i][j] = sum(a[i][k]*b[k][j] for k in range(p)) return c def mapmat(f, a): return [list(map(f, v)) for v in a] def ratmat(a): return mapmat(Fraction, a)