Finding inverse of a matrix using Gauss - Jordan Method

Question

I am trying to write a program that finds the inverse of an nxn matrix in Mathematica. I've been able to find matrices that don't find the inverse of the matrix but whose product has a determinant of 1, lol. Can you tell me where I went wrong?

m = Input["row"];
n = Input["column"];
mat = Table[Subscript[b, i, j], {i, m}, {j, n}];
mat2 = Table[Subscript[c, i, j], {i, m}, {j, n}];
For[i = 1, i <= m, i++,
  For[j = 1, j <= n, 
   j++, {Subscript[c, i, i] = 1, Subscript[c, i, j] = 0, 
    Subscript[c, m, n] = 1}]];
For[i = 1, i <= m, i++,
  For[j = 1, j <= n, j++,
   Subscript[b, i, j] = RandomInteger[{1, 3}]]];
Print[mat // MatrixForm]
For[k = 1, k <= n - 1, k++,
  For[l = k + 1, l <= n, l++,
   If[Subscript[b, l, k] == 0, continue,
    coe = Subscript[b, l, k]/Subscript[b, k, k];
    For[t = k, t <= n, t++,
     Subscript[b, l, t] = Subscript[b, l, t] - coe*Subscript[b, k, t];
     Subscript[c, l, t] = 
      Subscript[c, l, t] - coe*Subscript[c, k, t]]]]];
For[i = 1, i <= m, i++,
  coe2 = 1/Subscript[b, i, i];
  For[j = 1, j <= n, j++,
   Subscript[b, i, j] = Subscript[b, i, j]*coe2;
   Subscript[c, i, j] = Subscript[c, i, j]*coe2]];
For[k = n - 1, k >= 1, k--,
  For[l = n, l >= k + 1, l--,
   If[Subscript[b, k, l] == 0, continue,
    coe3 = Subscript[b, k, l]/Subscript[b, k, k];
    For[t = k, t <= n, t++,
     Subscript[b, t, l] = Subscript[b, t, l] - coe3*Subscript[b, t, k];
     Subscript[c, t, l] = 
      Subscript[c, t, l] - coe3*Subscript[c, t, k]]]]];
Print[mat // MatrixForm];
Print[mat2 // MatrixForm];
Quit[];

Do you have a description of the particular Gauss-Jordan algorithm you are trying to code? — Chris Degnen, Jan 06 '23 at 17:50

Chris Degnen · Answer 1 · 2023-01-06T16:01:56.963

This is not a complete answer. Just a work in progress.

For an nxn matrix you won't need m and n. Also, using Subscript makes the code difficult to read so I have replaced Subscript[b, i, j] with b[i, j] etc.

Clear[c, b, m, n, i, j, coe, coe2, coe3, k, l, t]

(* m = n = Input["n"]; *)
m = n = 5;

a = mat = Table[b[i, j], {i, m}, {j, n}];
mat2 = Table[c[i, j], {i, m}, {j, n}];
For[i = 1, i <= m, i++,
  For[j = 1, j <= n, j++,
   {c[i, i] = 1, c[i, j] = 0, c[m, n] = 1}]];
For[i = 1, i <= m, i++,
  For[j = 1, j <= n, j++,
   b[i, j] = RandomInteger[{1, 3}]]];
Print[mat // MatrixForm]

This is what you are aiming for: i.e. using RowReduce for the elimination, as per this post.

MatrixForm[inverse = RowReduce[Join[mat, mat2, 2]][[All, n + 1 ;;]]]
a.inverse == IdentityMatrix[n]

These matrices are not always invertible, e.g.

a = mat = {{1, 2, 3, 3, 3}, {1, 3, 1, 2, 1},
   {2, 3, 3, 2, 3}, {3, 2, 1, 2, 1}, {3, 1, 3, 1, 3}};

So far unable to correct your code.

For[k = 1, k <= n - 1, k++,
  For[l = k + 1, l <= n, l++,
   If[b[l, k] != 0,
    coe = b[l, k]/b[k, k];
    For[t = k, t <= n, t++,
     b[l, t] = b[l, t] - coe*b[k, t];
     c[l, t] = c[l, t] - coe*c[k, t]]]]];

For[i = 1, i <= m, i++,
  coe2 = 1/b[i, i];
  For[j = 1, j <= n, j++,
   b[i, j] = b[i, j]*coe2;
   c[i, j] = c[i, j]*coe2]];

For[k = n - 1, k >= 1, k--,
  For[l = n, l >= k + 1, l--,
   If[b[k, l] != 0,
    coe3 = b[k, l]/b[k, k];
    For[t = k, t <= n, t++,
     b[t, l] = b[t, l] - coe3*b[t, k];
     c[t, l] = c[t, l] - coe3*c[t, k]]]]];

IdentityMatrix[n] == mat

MatrixForm[inverse = mat2]
a.inverse == IdentityMatrix[n]

Thanks for your answer but I don't want to use wolfram language. — yasirbagci, Jan 07 '23 at 16:11
What algorithm are you trying to use? The above doesn't quite look like this one: https://en.wikipedia.org/wiki/Gaussian_elimination#Pseudocode — Chris Degnen, Jan 07 '23 at 18:02
Even if you don't want to use Wolfram the code can be easily transcribed once it's working. — Chris Degnen, Jan 07 '23 at 18:04

Finding inverse of a matrix using Gauss - Jordan Method

1 Answers1