Un poco de teoría
Como sabe, el determinante de una matriz cuadrada n * n es la suma de n! términos, cada uno de los cuales es un producto que contiene exactamente un elemento de matriz de cada columna y exactamente uno de cada fila. El signo de la siguiente pieza:
está determinada por la paridad de la sustitución:
\ begin {pmatrix} 1 & 2 &… & n \\ {i} _ {1} & {i} _ {2} &… & {i} _ {n} \ end {pmatrix}
Un método directo para calcular el determinante consiste en descomponerlo por los elementos de una fila o columna en la suma de los productos de los elementos de una fila o columna por sus complementos algebraicos. A su vez, el complemento algebraico del elemento de la matriz
Ahi esta
donde
- hay un elemento menor (i, j), es decir determinante obtenido del determinante original eliminando la i-ésima fila y la j-ésima columna.
Este método genera un proceso recursivo que permite calcular cualquier determinante. Pero el rendimiento de este algoritmo deja mucho que desear: ¡O (n!). Por lo tanto, se utiliza un cálculo directo de este tipo, excepto para cálculos simbólicos (y con determinantes de orden no demasiado elevado).
El método de Gauss resulta mucho más productivo. Su esencia se basa en las siguientes disposiciones:
1. Determinante de la matriz triangular superior \ begin {pmatrix} {a} _ {1,1} & {a} _ {1,2} &… & {a} _ {1, n} \\ 0 & {a} _ {2,2} &… & {a} _ {2, n} \\ 0 & 0 &… & ... \\ 0 & 0 &… & {a} _ {n, n} \\\ end { pmatrix} es igual al producto de sus elementos diagonales. Este hecho se sigue inmediatamente de la descomposición del determinante en los elementos de la primera fila o la primera columna. 2. Si en la matriz a los elementos de una fila suman los elementos de otra fila, multiplicados por el mismo número, entonces el valor del determinante no cambiará. 3. Si dos filas (o dos columnas) se intercambian en la matriz, entonces el valor del determinante cambiará su signo al opuesto. igual al producto de sus elementos diagonales. Este hecho se sigue inmediatamente de la descomposición del determinante en términos de los elementos de la primera fila o la primera columna.
Al elegir los coeficientes, podemos sumar la primera fila de la matriz con todas las demás y obtener ceros en la primera columna en todas las posiciones excepto en la primera. Para obtener cero en la segunda línea, debe agregar a la segunda línea la primera, multiplicada por
Para obtener cero en la tercera línea, agregue la primera línea multiplicada por
etc. Finalmente, la matriz se reducirá a una forma en la que todos los elementos
para n> 1 será igual a cero.
Si en la matriz el elemento
resulta ser igual a cero, entonces puede encontrar un elemento distinto de cero en la primera columna (suponga que está en el k-ésimo lugar) e intercambiar la primera y k-ésima fila. Con esta transformación, el determinante simplemente cambia de signo, lo que se puede tener en cuenta. Si no hay elementos distintos de cero en la primera columna, entonces el determinante es igual a cero.
Además, actuando de manera similar, puede obtener ceros en la segunda columna, luego en la tercera, etc. Es importante que al agregar las cadenas, los ceros obtenidos anteriormente no cambien. Si para cualquier cadena no es posible encontrar un elemento distinto de cero para el denominador, entonces el determinante es igual a cero y el proceso puede detenerse. La finalización normal del proceso gaussiano genera una matriz en la que todos los elementos ubicados debajo de la diagonal principal son iguales a cero. Como se mencionó anteriormente, el determinante de dicha matriz es igual al producto de los elementos diagonales.
Pasemos a la programación.
Trabajamos con datos de punto flotante. Representamos matrices como listas de cadenas. Primero, definamos dos tipos:
type Row = [Double]
type Matrix = [Row]
Recursividad simple
No más dudas, expandiremos el determinante por los elementos de la primera fila (es decir, cero). Necesitamos un programa para construir el menor, obtenido al eliminar la primera fila y la columna k.
-- k-
deln :: Matrix -> Int -> Matrix
deln matrix k = map (\ r -> (take (k) r)++(drop (k+1) r)) matrix
Y aquí está el menor:
-- k-
minor :: Matrix -> Int -> Double
minor matrix k = det $ deln (drop 1 matrix) k
Tenga en cuenta: el menor es un factor determinante. Llamamos a la función det, que aún no hemos implementado. Para implementar det, tendremos que formar la suma alterna de los productos del siguiente elemento de la primera fila por el determinante del siguiente menor. Para evitar expresiones engorrosas, creemos una función separada para formar el signo de suma:
sgn :: Int -> Double
sgn n = if n `rem` 2 == 0 then 1.0 else (-1.0)
Ahora podemos calcular el determinante:
--
det :: Matrix -> Double
det [[a,b],[c,d]] = a*d-b*c
det matrix = sum $ map (\c -> ((matrix !! 0)!!c)*(sgn c)*(minor matrix c)) [0..n]
where n = length matrix - 1
El código es muy simple y no requiere comentarios especiales. Para probar el desempeño de nuestras funciones, escribamos la función principal:
main = print $ det [[1,2,3],[4,5,6],[7,8,(-9)]]El valor de este determinante es 54, lo cual se puede verificar.
Método de Gauss
Necesitaremos algunas funciones de utilidad (que se pueden usar en otros lugares). El primero es el intercambio de dos filas en la matriz:
--
swap :: Matrix -> Int -> Int -> Matrix
swap matrix n1 n2 = map row [0..n]
where n=length matrix - 1
row k | k==n1 = matrix !! n2
| k==n2 = matrix !! n1
| otherwise = matrix !! k
Como puede ver en el código anterior, la función pasa línea por línea. En este caso, si se encuentra una línea con el número n1, la línea n2 se sustituye a la fuerza (y viceversa). El resto de las líneas permanecen en su lugar.
La siguiente función calcula la cadena r1 sumada con la cadena r2 multiplicada elemento por elemento por el número f:
-- r1+f*r2
comb :: Row -> Row -> Double -> Row
comb r1 r2 f = zipWith (\ x y -> x+f*y) r1 r2
Todo aquí es extremadamente transparente: las acciones se realizan en las filas de la matriz (es decir, en listas [dobles]). Pero la siguiente función realiza esta transformación en la matriz (y, naturalmente, obtiene una nueva matriz):
-- r1 r2, f
trans :: Matrix -> Int -> Int -> Double -> Matrix
trans matrix n1 n2 f = map row [0..n]
where n=length matrix - 1
row k | k==n1 = comb (matrix !! n1) (matrix !! n2) f
| otherwise = matrix !! k
La función getNz busca el número del primer elemento distinto de cero en la lista. Es necesario cuando el siguiente elemento diagonal es igual a cero.
--
getNz :: Row -> Int
getNz xs = if length tmp == 0 then (-1) else snd $ head tmp
where tmp=dropWhile (\ (x,k) -> (abs x) <= 1.0e-10) $ zip xs [0..]
Si todos los elementos de la lista son cero, la función devolverá -1.
La función de búsqueda comprueba si la matriz es adecuada para la siguiente transformación (debe tener un siguiente elemento diagonal distinto de cero). Si no es así, la matriz se transforma mediante permutación de filas.
--
search :: Matrix -> Int -> Matrix
search matrix k | (abs ((matrix !! k) !! k)) > 1.0e-10 = matrix
| nz < 0 = matrix --
| otherwise = swap matrix k p
where n = length matrix
lst = map (\ r -> r !! k) $ drop k matrix
nz = getNz lst
p = k + nz
Si no se puede encontrar el elemento principal (distinto de cero) (la matriz está degenerada), la función lo devolverá sin cambios. La función mkzero forma ceros en la siguiente columna de la matriz:
--
mkzero :: Matrix -> Int -> Int -> Matrix
mkzero matrix k p | p>n-1 = matrix
| otherwise = mkzero (trans matrix p k (-f)) k (p+1)
where n = length matrix
f = ((matrix !! p) !! k)/((matrix !! k) !! k)
La función triangular forma la forma triangular superior de la matriz:
--
triangle :: Matrix -> Int -> Matrix
triangle matrix k | k>=n = matrix
| (abs v) <= 1.0e-10 = [[0.0]] --
| otherwise = triangle (mkzero tmp k k1) k1
where n = length matrix
tmp = search matrix k
v = (tmp !! k) !! k --
k1 = k+1
Si en la siguiente etapa no fue posible encontrar el elemento principal, la función devuelve una matriz cero de primer orden. Ahora puede componer la función de desfile de reducir la matriz a la forma triangular superior:
--
gauss :: Matrix -> Matrix
gauss matrix = triangle matrix 0
Para calcular el determinante, necesitamos multiplicar los elementos diagonales. Para hacer esto, creemos una función separada:
--
proddiag :: Matrix -> Double
proddiag matrix = product $ map (\ (r,k) -> r !!k) $ zip matrix [0,1..]
Bueno, y el "arco" es el cálculo real del determinante:
--
det :: Matrix -> Double
det matrix = proddiag $ triangle matrix 0
Veamos cómo funciona esta función:
main = print $ det [[1,2,3],[4,5,6],[7,8,-9]]
[1 of 1] Compiling Main ( main.hs, main.o )
Linking a.out ...
54.0
¡Gracias a los que leyeron hasta el final!
El código se puede descargar aquí.