import numpy def backsubstitute(A, b): n,i = b.shape if i != 1: raise ValueError("b should be shape (n, 1), but is %s" % (b.shape,)) i,j = A.shape if i != n or j != n: raise ValueError("A should be shape (%d, %d), but is %s" % (n,n,A.shape)) x = b.copy() for i in range(n-1,-1,-1): x[i] = b[i] for j in range(n-1, i, -1): x[i] -= A[i,j] * x[j] x[i] /= A[i,i] return x def gauss_eliminate(A_in, b_in): # Kopien der Matrix und des Vektors anlegen A = A_in.copy() b = b_in.copy() # Gausselimination durchfuehren # .... return A, b def solve(A, b): Ares, bres = gauss_eliminate(A, b) x = backsubstitute(Ares, bres) return x N = 5 A1 = numpy.random.random((N,N)) A2 = A1.copy() A2[0,0] = 0.0 b2 = b1 = numpy.random.random((N,1)) import scipy.linalg print "----------" print "Solve eq 1" print "----------" print "A1 * x1 = b1" print "A1 = %s" % A1 print "b1 = %s" % b1 x1_scipy = scipy.linalg.solve(A1, b1) print "x1_scipy = %s" % x1_scipy x1_own = solve(A1, b1) print "x1_own = %s" % x1_own diff = abs(x1_own - x1_scipy) print "Maximal error = %s" % diff.max() print print "----------" print "Solve eq 2" print "----------" print "A2 * x2 = b2" print "A2 = %s" % A2 print "b2 = %s" % b2 x2_scipy = scipy.linalg.solve(A2,b2) print "x1_scipy = %s" % x2_scipy x2_own = solve(A2,b2) print "x2_own = %s" % x2_own diff = abs(x1_own - x1_scipy) print "Maximal error = %s" % diff.max()