I am trying to write a program that can do Gaussian Elimination without partial pivoting. This is what I have so far, I know I'm messing something up but I can't seem to figure out what. It keeps printing out the same matrix it started with and not the 0s with 1s in the diagonal. Any help appreciated!
import numpy as npimport math#A = np.array([[1,2,-1],[5,2,2],[-3,5,-1]])#b = np.array([[2,9,1]])A = np.array([[8,-2,-1,0,0],[-2,9,-4,-1,0],[-1,-3,7,-1,2],[0,-4,-2,12,-5],[0,0,-7,-3,15]])b = np.array([[5],[2],[0],[1],[5]])def forward_elim(A, b, n):#calculates the forward part of Gaussian elimination.for row in range(0, n-1):for i in range(row+1, n):factor = A[i,row] / A[row,row]for j in range(row, n):A[i,j] = A[i,j] - factor * A[row,j]b[i] = b[i] - factor * b[row]print("A = \n%s and b = %s" % (A,b))return A, bdef back_sub(A, b, n):#back substitution and returns resultx = np.zeros((n,1))x[n-1] = b[n-1] / A[n-1, n-1]for row in range(n-2, -1, -1):sums = b[row]for j in range(row+1, n):sums = sums - A[row,j] * x[j]x[row] = sums / A[row,row]return xdef gauss_elim(A, b):#performs gaussian elimination without pivotingn = A.shape[0]#checks for zero diagonal elementsif any(np.diag(A)==0):raise ZeroDivisionError(("Can't divide by 0")) #raise used to raise exceptionA, b = forward_elim(A, b, n)return back_sub(A, b, n)print(A)Best Answer
From the given code it is not clear whether you effectively call the function gauss_elim(A, b).
To help finding out was is going on, try stepping through a really simple example. Best with the help of a debugger, or by printing out intermediate values. Follow the same example writing the steps on a sheet of paper.
A simple example could be:
A = np.array([[2,4],[3,9]])b = np.array([[2,30]])One of the problems is that you declare b as a 1xn matrix, but later use it as just a vector. A solution can be to declare b as a vector:
b = np.array([2,30])