0

I am trying to calculate the inverse matrix using the Gauss-Jordan Method. For that, I need to find the solution X to A.X = I (A and X being N x N matrices, and I the identity matrix).

However, for every column vector of the solution matrix X I calculate in the first loop, I have to use the original matrix A, but I don't know why it keeps changing when I did a copy of it in the beginning.

def SolveGaussJordanInvMatrix(A):
 N = len(A[:,0])
 I = np.identity(N)
 X = np.zeros([N,N], float)
 A_orig = A.copy()

for m in range(N): 
 x = np.zeros(N, float)
 v = I[:,m] 
 A = A_orig 

 for p in range(N): # Gauss-Jordan Elimination 
  A[p,:] /= A[p,p]
  v[p] /= A[p,p]

  for i in range(p): # Cancel elements above the diagonal element
    v[i] -= v[p] * A[i,p]
    A[i,p:] -= A[p,p:]*A[i,p]

  for i in range(p+1, N): # Cancel elements below the diagonal element
    v[i] -= v[p] * A[i,p]
    A[i,p:] -= A[p,p:]*A[i,p]


  X[:,m] = v # Add column vector to the solution matrix

return X

A = np.array([[2, 1, 4, 1 ],
          [3, 4, -1, -1],
          [1, -4, 7, 5],
          [2, -2, 1, 3]], float)

SolveGaussJordanInvMatrix(A)

Does anyone know how turn A back to its original form after the Gauss-Elimination loop?

I'm getting

array([[ 228.1,    0. ,    0. ,    0. ],
   [-219.9,    1. ,    0. ,    0. ],
   [ -14.5,    0. ,    1. ,    0. ],
   [-176.3,    0. ,    0. ,    1. ]])

and expect

    [[ 1.36842105 -0.89473684 -1.05263158  1.        ]
   [-1.42105263  1.23684211  1.13157895 -1.        ]
   [ 0.42105263 -0.23684211 -0.13157895 -0.        ]
   [-2.          1.5         1.5        -1.        ]]
lviveiros
  • 1
  • 1
  • 1
    I tried using deepcopy instead of copy, but it still doesn't work. The matrix A_orig changes after the 1st loop. – lviveiros May 18 '22 at 18:56
  • Would you mind sharing the code outside of the function so we can see how you're calling it and what you're displaying after the call? – whege May 18 '22 at 19:11
  • 1
    I've just updated it :) – lviveiros May 18 '22 at 19:20
  • Well for one thing, you're only returning X and not A (or any version thereof). If you want to get the original matrix A from your function, you also need to have it return A_orig. – whege May 18 '22 at 19:25
  • Actually, I only need to get the solution matrix, which is the inverse of A. For every Gauss loop, I solve the system A.x = v (x beign a column vector and v a column of the I matrix). When the loop is done, v gives the solution to the system Ax = v. I need to solve a system Ax = v N times. The problem is that when I try to solve the system a 2nd time, the Gauss loop doesn't use the original matrix A. – lviveiros May 18 '22 at 19:31
  • Gotcha. What result are you getting back, and what are you expecting? – whege May 18 '22 at 19:42
  • I've just updated it so you can see. Just so you know, the Gauss loop turns the original matrix A into the identity to find the x vector (Ax = v) – lviveiros May 18 '22 at 19:53

0 Answers0