El algoritmo que voy a explicar, es el más usado y el que mejor rendimiento ofrece, en comparación al algoritmo recursivo, que se usa para hallar la determinante de una matriz nxn. Este algoritmo se centra básicamente en convertir una matriz nxn, en una matriz triangular superior, usando el método de eliminación de Gauss. En donde el determinante de toda matriz triangular superior, es la multiplicación de todos los elementos de su diagonal. Para ello debemos tener en cuenta ciertas reglas o propiedades de las matrices:

  1. El intercambio de una fila o columna de una matriz nxn, produce el cambio de signo del determinante.
  2. Una matriz nxn formada con una fila o columna de ceros tiene un determinante cero.
  3. Si a una fila o columna de una matriz nxn se le suma o resta, otra paralela a ella multiplicada por un número cualquiera, su determinante no varía.

Para este primer acercamiento al algoritmo, haremos uso de la 3ra propiedad. Tomemos en cuenta la siguiente matriz, como ejemplo a desarrollar:

Algoritmo para hallar el determinante de una matriz nxn en Free Pascal - 01.

La idea principal es convertir los números 1, 5, -8 de color rojo de las filas 2,3 y 4 en ceros, para ello debemos obtener un número cualquiera, que al multiplicarlo y restarlo por una de las filas, nos permita convertir los números 1, 5, -8 en ceros, esto lo podemos hacer del siguiente modo:

f2 = f2 - (f1*1/3)

f3 = f3 - (f1*5/3)

f4 = f4 - (f1*-8/3)

el número cualquiera se obtiene dividiendo los números 1, 5, -8 al cuál denominaremos como valores auxiliares, con el primer elemento de la diagonal de la matriz, que en este caso es el 3, a los elementos de una diagonal se les conoce como pivotes. Después de hacer las operaciones indicadas obtendremos la siguiente matriz:

Algoritmo para hallar el determinante de una matriz nxn en Free Pascal - 01.

Procedemos de modo similar en donde esta vez se hacen las divisiones con el pivote 2/3 y los valores auxiliares -2/3, 5/3. Las operaciones serían entonce las siguientes:

f3 = f3 - (f2*-2:3/2:3)

f4 = f4 - (f2*5:3/2:3)

Con lo cuál obtenemos la siguiente matriz:

Algoritmo para hallar el determinante de una matriz nxn en Free Pascal - 01.

Siguiendo el mismo mecanismo la última operación a realizar sería la siguiente:

f4 = f4 - (f3*16/4)

Con lo cuál obtenemos la siguiente matriz triangular:

Algoritmo para hallar el determinante de una matriz nxn en Free Pascal - 01.

En donde el determinante se obtiene multiplicando los pivotes de la matriz triangular.

Det = 2 · (2/3) · 4 · 1 = 8

Ahora lo que tenemos que hacer es convertir todas estas operaciones, en un algoritmo que se pueda resolver con el uso de una computadora. Primero se debe empezar por un bucle que recorra todos los pivotes, el bucle entonces será algo como esto:

Para d=1 hasta columnas hacer
pivote = M[d,d]
...
FinPara

Aquí M es la matriz, d se refiere a la diagonal, obviamente la cantidad de pivotes que tiene la matriz será la cantidad de columnas o filas que tiene la matriz nxn, ya que esta es una matriz cuadrada.

Lo siguiente es anidar un segundo bucle para obtener los valores auxiliares, del siguiente modo:

Para d=1 hasta columnas hacer
 pivote = M[d,d]
  Para f = d+1 hasta filas hacer
   aux = M[f,d]
   ...
  FinPara
FinPara

Los valores auxiliares, empiezan a partir de la segunda fila, o de la siguiente fila del pivote anteriormente obtenido, es la razón por la cual f, empieza con d+1.

Ahora se hace un tercer bucle anidado a los dos anteriores, que recorre todas las columnas para operar los elementos correspondientes a la fila a operar.

Para d=1 hasta columnas hacer
   pivote = M[d,d]
   Para f = d+1 hasta filas hacer
      aux = M[f,d]
      Para c = 1 hasta columnas hacer
         ...
      FinPara
   FinPara
FinPara

Y finalmente la operación, que permite poner a cero, los elementos por debajo de la diagonal de la matriz:

Para d=1 hasta columnas hacer

   pivote = M[d,d]

   Para f = d+1 hasta filas hacer

      aux = M[f,d]

      Para c=1 hasta columnas hacer

         M[f,c] = M[f,c] - (M[d,c]*(pivote/aux))

      FinPara

   FinPara

FinPara

Para hallar la determinante sólo es necesario hacer un bucle que recorra la diagonal y los multiplique:

Det = 1

Para d = 1 hasta columnas hacer Det = Det * M[d,d]

La implementación de este primer acercamiento del algoritmo en es el siguiente:

Algo a tomar en cuenta en la implementación del algoritmo, el programa usa arreglos dinámicos con lo que el índice de las matrices empieza siempre en 0. Este algoritmo funcioná bien siempre y cuando los pivotes sean diferentes a cero. En la siguiente sección de este apunte se explica como mejorar este algoritmo para detectar un pivote nulo y proceder a intercambiarlo con otra fila, cuando sea necesario.

En el primer acercamiento del algoritmo para hallar la determinante de una matriz nxn, no se vio la posibilidad de que los pivotes, puedan tener un cero.

El programa que usaremos a continuación, es el mismo que el anterior, pero para hacernos la vida un poco más fácil le añadí dos procedimientos, Llenar y Mostrar, Llenar nos permite llenar la matriz de una manera más cómoda, y Mostrar se creo para no estar repitiendo código, cada vez que queremos mostrar la matriz. Veamos ahora el primer programa de prueba:

en este programa se esta usando la siguiente matriz:

Algoritmo para hallar el determinante de una matriz nxn en Free Pascal - 02.

El determinante de esta matriz es 293, pero debido al cero en el primer pivote de la diagonal de la matriz, el programa fallará y nos producirá un error por una división con cero.

f2 = f2 - (f1*1/0)

Ahora veamos otro caso, en donde la matriz aparentemente no tiene pivotes con ceros en su diagonal, como por ejemplo esta matriz que tiene como determinante 50/3.

Algoritmo para hallar el determinante de una matriz nxn en Free Pascal - 01.

El problema aparecerá con el 4/3, ya que si hacemos la operación correspondiente

f3 = f3 - (f2*1/3)

y reemplazando valores

M[2,2] = 4/3 - (4*1/3)

haremos que el segundo pivote se convierta de 4/3 a un 0.

Ahora haremos los cambios necesarios para mejorar el algoritmo, obviamente la solución es añadir el código que nos permita intercambiar filas para evitar tener un pivote con cero. Primero debemos detectar si el pivote tiene un cero, para ello dentro del primer bucle y antes de ingresar al bucle que recorre las filas de la matriz, se verifica con una estructura condicional Si / entonces, si el pivote tiene un cero.

Para d=1 hasta columnas hacer

    Si M[d,d] = 0 entonces

       ...

    FinSi

    pivote = M[d,d]

    Para f = d+1 hasta filas hacer

       aux = M[f,d]

       Para c = 1 hasta columnas hacer

          M[f,c] = M[f,c] - (M[d,c] * (pivote/aux))

       FinPara

    FinPara

FinPara

Debido a que vamos a intercambiar las filas entonces el determinante de la matriz, cambiará de signo. Para controlar esto añadimos una nueva variable z, que se multiplicará por -1, cada vez que se realice un intercambio de filas. También necesitamos una variable p, que nos indicará si realmente hubo intercambios, este debe empezar con un valor 0 que nos indicá que no hay intercambios, si p es distinto de cero entonces nos indica que se pudo intercambiar la fila, y p tendrá el valor de la fila del intercambio, pero en caso tenga un valor -1 este nos indicará que no se pudo intercambiar la fila, o mejor dicho que no se pudo encontrar una fila con la cual se pueda intercambiar la fila con el pivote 0.

z=1

Para d=1 hasta columnas hacer

    p=0

    Si M[d,d]=0 entonces

       p=-1

       ...

    FinSi

    pivote = M[d,d]

    Para f = d+1 hasta filas hacer

       aux = M[f,d]

       Para c = 1 hasta columnas hacer

          M[f,c] = M[f,c] - (M[d,c] * (pivote/aux))

       FinPara

    FinPara

FinPara

Ahora debemos implementar el bucle que recorra las filas disponibles para hacer un intercambio, esto se hace con un bucle Mientras, que recorre las filas que hay debajo del pivote, para encontrar una con la cual intercambiarla.

z=1

Para d=1 hasta columnas hacer

    p=0

    Si M[d,d] = 0 entonces

       p = -1

       f = d

       Mientras (f <= filas) y (p=-1)

           Si M[f,d] <> 0 entonces

              p=f FilaIntercambio(d,p)

              z = z * -1 //por cada intercambio se debe multiplicar

            FinSi

            f = f + 1

       FinMientras

    FinSi

    pivote = M[d,d]

    Para f = d + 1 hasta filas hacer

       aux = M[f,d]

       Para c = 1 hasta columnas hacer

          M[f,c] = M[f,c] - (M[d,c] * (pivote/aux))

       FinPara

    FinPara

FinPara

 

En este algoritmo se hizo un cambio importante, se añadió una estructura Si / entonces, antes de hacer las operaciones para convertir la matriz en triangular superior. Esta estructura Si / entonces comprueba el valor de p, si este es -1 nos indica que hay un pivote 0 y que no se pudo encontrar una fila para intercambiarla, con lo cual ya no se deben hacer más operaciones para convertir la matriz en triangular superior.

A continuación el programa, con la implementación del algoritmo completo, que nos permite hallar el determinante de la siguiente matriz.

Algoritmo para hallar el determinante de una matriz nxn en Free Pascal - 01.

Como punto final de este artículo, debo mencionar que es muy importante considerar la precisión del número flotante que se va a usar, para implementar este algoritmo. En la implementación hemos estado usando un número flotante de 80 bits (extended), que tiene una precisión mucho mayor que el tipo de dato double (64 bits). Si implementamos esto con otro lenguaje de programación que no soporte un numero flotante de 80 bits, o mejor dicho si sólo usamos double. La siguiente línea de código:

M[f,c] := M[f,c] - (M[d,c] * (aux / pivote))

con los siguientes valores

M[f,c] := 4/3 - (4 * (1/3))

M[f,c] := 1.333 ... 3 - (4 * (0.333 ... 3))

M[f,c] := 1.333 ... 3 - (1.333 ... 2)

M[f,c] := 0.000 ... 1

No producirá exactamente un cero (0.000...1), haciendo que el algoritmo fallé. Para ello se puede usar un artificio, que consiste en convertir el número obtenido M[f,c] a una cadena de caracteres indicándole una cantidad menor de dígitos decimales, y de nuevo volver a convertirlo a un valor numérico, para luego asignarlo a M[f,c], en freepascal sería algo así:

M[f,c] := M[f,c] - (M[d,c] * (aux / pivote));

Str(M[f,c]:0:10, cad);

Val(cad, M[f,c], error);

Para el algoritmo implementado en Java, se hace uso de la clase String para dar formato y precisión en los decimales, con 10 cifras significativas en la parte decimal y luego se convierte al tipo de dato double.