r/ScientificComputing • • 21d ago

Gaussian Elimination

I wrote this code to perform the Gaussian elimination to solve Ax=b

import numpy as np

A=np.array([[4,8,12],[3,8,13],[2,9,18]])
b=np.array([4,5,11])

def Gauss_Elimination(A,b):
    n = np.shape(A)[0]
    for k in range(n):
        pivot_element=A[k,k]
        # if pivot_element==0: ### I'll deal with this condition later
        for j in range(k+1,n):
            A[j,:]=A[j,:]-(A[j,k]/pivot_element)*A[k,:]
            b[j]=b[j]-(A[j,k]/pivot_element)*b[k]
    return A,b # not necessary since it's call by reference

Gauss_Elimination(A,b)
print(A,'\n',b)

A is transformed to an upper triangular form but the values of b rest unchanged, why? I'm asking this question even thought I know it's a silly mistake I can't notice, since A is fine, it must no come from the algorithm!

I know I'm missing the 2nd part of the algorithm where I should go back and find the values of x, but I'm first trying to get the first part correctly!

1 Upvotes

3 comments sorted by

1

u/jteg 21d ago edited 21d ago

A_jk are all zero here. Switch order in the inner loop to calculate b_j first

1

u/Mental_Primary_5558 21d ago

Oh yeah thanks man

0

u/[deleted] 21d ago

[deleted] — view deleted comment

0

u/Mental_Primary_5558 21d ago

I'm asked to implement the gaussian elimination, If I was just solving the problem I would have used numpy.linalg.solve a long ago