from numpy import * from matplotlib.pyplot import * 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 # TASK 2.1 def gauss_eliminate(A_in, b_in): A = A_in.copy() b = b_in.copy() 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)) # which row is currently handled for r in range(0, n): # which row to eliminate for i in range(r+1, n): l = A[i,r]/A[r,r] A[i,r] = 0 b[i] -= l*b[r] A[i,r+1:] -= l*A[r,r+1:] 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]) print x1 - solve(A1, b1) # fails: #print solve(A2, b2) # TASK 2.2 figure() # get the points xs = linspace(0, 2*pi, 4) ys = sin(xs) plot(xs, ys, 'o') # create the matrix A = zeros((4,4)) A[:, 0] = 1 A[:, 1] = xs A[:, 2] = xs**2 A[:, 3] = xs**3 # compute the coefficients a = solve(A, ys) def polyeval(a,x): p = zeros_like(x) xfac = 1.0 for av in a: p += av*xfac xfac *= x return p # plot the functions xs = linspace(0, 2*pi, 100) plot(xs, sin(xs), label='sin(x)') plot(xs, polyeval(a,xs), label='p(x)') # TASK 2.3.1 def gauss_eliminate_partial(A_in, b_in): A = A_in.copy() b = b_in.copy() 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)) # which row is currently handled for r in range(0, n): i = r + abs(A[r:,r]).argmax() # swap the rows A[[r,i],:] = A[[i,r],:] b[[r,i]] = b[[i,r]] # which row to eliminate for i in range(r+1, n): l = A[i,r]/A[r,r] A[i,r] = 0 b[i] -= l*b[r] A[i,r+1:] -= l*A[r,r+1:] return A, b def solve_partial(A, b): n,i = A.shape Ares, bres = gauss_eliminate_partial(A, b) x = backsubstitute(Ares, bres) return x print x2 - solve_partial(A2, b2) # TASK 2.3.2 def gauss_eliminate_complete(A_in, b_in): A = A_in.copy() b = b_in.copy() 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)) original_index = arange(0,n) # which row is currently handled for r in range(0, n): index = abs(A[r:,r:]).argmax() i,j = unravel_index(index, (n-r, n-r)) # row to swap i += r # column to swap j += r # swap the rows A[[r,i],:] = A[[i,r],:] b[[r,i]] = b[[i,r]] # swap the columns A[:,[r,j]] = A[:,[j,r]] original_index[[r,j]] = original_index[[j,r]] # which row to eliminate for i in range(r+1, n): l = A[i,r]/A[r,r] A[i,r] = 0 b[i] -= l*b[r] A[i,r+1:] -= l*A[r,r+1:] return A, b, original_index def solve_complete(A, b): n,i = A.shape Ares, bres, original_index = gauss_eliminate_totalpivot(A, b) x = backsubstitute(Ares, bres) x[original_index,:] = x[range(0,n),:] return x show()