import numpy def backsubstitute(A, b): n = len(b) 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 # Sample solutions A1 = array([[ 1. , 0.2 , 0.04 , 0.008 ], [ 1. , 0.53333333, 0.28444444, 0.1517037 ], [ 1. , 0.86666667, 0.75111111, 0.65096296], [ 1. , 1.2 , 1.44 , 1.728 ]]) x1 = array([ 1.27785932, -0.17219435, -8.75422788, 7.22564732]) b1 = array([ 0.95105652, -0.20791169, -0.74314483, 0.95105652]) A2 = array([[ 0. , 0.2 , 0.04 , 0.008 ], [ 1. , 0.53333333, 0.28444444, 0.1517037 ], [ 1. , 0.86666667, 0.75111111, 0.65096296], [ 1. , 1.2 , 1.44 , 1.728 ]]) x2 = array([ -0.85418404, 8.06213979, -18.74818115, 11.0694755 ]) b2 = array([ 0.95105652, -0.20791169, -0.74314483, 0.95105652])