def gauss(A,b,offset=(0,0)): ox,oy = offset for i in range(oy,len(A)): if A[i][ox] != 0: A[i], A[oy] = A[oy], A[i] b[i], b[oy] = b[oy], b[i] break pprint(A) for i in range(oy+1,len(A)): if A[i][ox] != 0: c = A[i][0]/A[oy][0] A[i] = [ (A[i][x] - A[oy][x]*c) for x in range(ox,len(A[i])) ] b[i] = b[i] - b[oy]*c if oy < len(A)-1: A,b = gauss(A,b,offset=(ox+1,oy+1)) b[oy] /= A[oy][ox] x = b[oy] for y in range(oy): b[y] -= A[y][ox]*x return A,b def vprint(A): print("Vector ----") for i in range(len(A)): print("%10.4f" % A[i]) print("-----") def pprint(A): print("Matrix ----") for i in range(len(A)): print( " ".join(["%10.4f" % (A[i][j]) for j in range(len(A[i])) ]) ) print("-----") A = [ [ 0, 0, 0,-1,-1, 1 ], [ -1, 1, 0, 0, 0, 0 ], [ 0, -1, 1, 0, 0, 0 ], [-1000, 0, 0, 1, 0, 0 ], [ 0, -1000, 0, 0, 1, 0 ], [ 0, 0, 0, 0, 0, 1 ] ] b = [ 0, 0, 0, 0, 0, 5 ] pprint(A) vprint(b) # _, x = gauss(A,b) from numpy.linalg import solve x = solve(A,b) vprint(x)