Eliminación Gaussiana Usando Pivote en Python

La eliminación gaussiana también se conoce como el método de reducción de filas. Es un algoritmo comúnmente utilizado para resolver problemas lineales. El algoritmo implica una serie de operaciones de fila en una matriz de coeficientes extraídos de las ecuaciones lineales hasta que la matriz se reduce a la forma escalonada.

Las siguientes tres operaciones se realizan cuando se transforma una matriz en una forma escalonada. Éstas incluyen:

• Multiplicar una fila por un escalar y sumarlo o restarlo de otra fila.
• Intercambio de filas.
• Multiplicar filas con un escalar.

Una matriz está en forma escalonada de filas si el primer elemento de cada fila, también conocido como la entrada inicial, es un cero; la entrada inicial de cada fila está una columna a la derecha de la entrada inicial de la fila anterior. Además, las filas con elementos cero deben estar debajo de los elementos que tienen filas distintas de cero.

Las operaciones de fila se clasifican en gran medida en dos. La eliminación directa implica reducir una matriz a forma escalonada para determinar si existe una solución viable y si es finita o infinita. Por otro lado, la sustitución hacia atrás reduce aún más la matriz a la forma escalonada reducida por filas.

Eliminación gaussiana con pivote en Python

Pivotar es el intercambio de filas y columnas para obtener el elemento de pivote adecuado. Un elemento de pivote adecuado debe ser distinto de cero y significativamente grande pero más pequeño en comparación con las otras entradas de la fila.

El pivote se clasifica en pivote parcial y pivote completo. Bajo el algoritmo de pivote parcial, el elemento más grande se considera el elemento pivote para minimizar los errores de redondeo.

Por otro lado, el pivoteo completo incluye el intercambio de filas y columnas para obtener el mejor elemento de pivote, aumentando así la precisión.

Un conjunto de ecuaciones se considera lineal si ninguna variable tiene un exponente de más de uno. La eliminación gaussiana implica una serie de pasos; el primer paso consiste en crear una matriz de coeficientes.

Una matriz de coeficientes es simplemente una matriz de coeficientes extraídos de un conjunto de ecuaciones lineales. El siguiente paso consiste en crear una matriz aumentada que luego se somete a una serie de operaciones que la reducen a una forma escalonada.

Sin embargo, en los casos en que el elemento pivote sea cero o de muy pequeña magnitud, tenemos que intercambiar la fila pivote con una fila inferior.

Al implementar el método de eliminación de Gauss, debemos tener en cuenta los índices de matriz. Los iterables de Python, como listas y matrices, a menudo comienzan en el índice 0 y terminan en el índice n-1.

Luego podemos leer el contenido del array aumentada, aplicar el método de eliminación, sustitución hacia atrás y finalmente mostrar la respuesta.

from numpy import array, zeros, fabs, linalg

a = array(
    [
        [0, 6, -1, 2, 2],
        [0, 3, 4, 1, 7],
        [5, 1, 0, 3, -1],
        [3, 1, 3, 0, 2],
        [4, 4, 1, -2, 1],
    ],
    float,
)
# the b matrix constant terms of the equations
b = array([5, 7, 2, 3, 4], float)

print("Solution by NumPy:")


print(linalg.solve(a, b))

n = len(b)
x = zeros(n, float)

# first loop specifys the fixed row
for k in range(n - 1):
    if fabs(a[k, k]) < 1.0e-12:

        for i in range(k + 1, n):
            if fabs(a[i, k]) > fabs(a[k, k]):
                a[[k, i]] = a[[i, k]]
                b[[k, i]] = b[[i, k]]
                break

    # applies the elimination below the fixed row

    for i in range(k + 1, n):
        if a[i, k] == 0:
            continue

        factor = a[k, k] / a[i, k]
        for j in range(k, n):
            a[i, j] = a[k, j] - a[i, j] * factor
            # we also calculate the b vector of each row
        b[i] = b[k] - b[i] * factor
print(a)
print(b)


x[n - 1] = b[n - 1] / a[n - 1, n - 1]
for i in range(n - 2, -1, -1):
    sum_ax = 0

    for j in range(i + 1, n):
        sum_ax += a[i, j] * x[j]

    x[i] = (b[i] - sum_ax) / a[i, i]

print("The solution of the system is:")
print(x)

Producción :

Solution by NumPy:
[ 0.38947368  0.49473684 -0.10877193  0.12982456  0.83157895]
[[ 5.00000000e+00  1.00000000e+00  0.00000000e+00  3.00000000e+00  -1.00000000e+00]
 [ 0.00000000e+00  3.00000000e+00  4.00000000e+00  1.00000000e+00  7.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  4.50000000e+00  0.00000000e+00  6.00000000e+00]
 [ 0.00000000e+00  4.44089210e-16  0.00000000e+00  3.52702703e+00  2.95945946e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00  2.11354647e+00]]
[2.         7.         4.5        2.91891892 1.75758075]
The solution of the system is:
[ 0.38947368  0.49473684 -0.10877193  0.12982456  0.83157895]

En el código anterior, hemos utilizado la función array y la función fabs proporcionadas por la biblioteca NumPy para crear una matriz y leer valores absolutos. También hemos utilizado Linalg; una subbiblioteca NumPy utilizada para realizar operaciones como el cálculo de valores propios y vectores y determinantes.

from numpy import array, zeros, fabs, linalg

La variable a representa los coeficientes constantes de las variables, mientras que b represental array de términos constantes correspondientes a las ecuaciones. Generalmente, la matriz y el vector se extraen de la representación matricial de cuatro ecuaciones lineales que pretendemos resolver.

a = array(
    [
        [0, 6, -1, 2, 2],
        [0, 3, 4, 1, 7],
        [5, 1, 0, 3, -1],
        [3, 1, 3, 0, 2],
        [4, 4, 1, -2, 1],
    ],
    float,
)
# the b matrix constant terms of the equations
b = array([5, 7, 2, 3, 4], float)

El objetivo de la eliminación es transformar la matriz para que tenga ceros en la diagonal principal.

Al realizar la eliminación gaussiana, la primera fila debe permanecer fija. Sin embargo, si el elemento pivote en la primera fila es igual a cero, debemos intercambiar la primera fila con cualquier fila subsiguiente sin cero como elemento principal.

Hemos usado el bucle k para especificar la fila fija en el código anterior. En este caso, k representa el índice de la fila fija que va de 1 a n-1.

Las siguientes declaraciones if verifican si el elemento pivote es cero y lo intercambian con la fila apropiada. La i utilizada representa el índice de las filas debajo de la fila fija.

for k in range(n - 1):
    if fabs(a[k, k]) < 1.0e-12:

        for i in range(k + 1, n):
            if fabs(a[i, k]) > fabs(a[k, k]):
                a[[k, i]] = a[[i, k]]
                b[[k, i]] = b[[i, k]]
                break

En el siguiente código, usamos un bucle para garantizar que la eliminación no se aplique en la primera fila.

Luego hemos multiplicado cada fila debajo de la fila fija por un factor. Un factor es igual al elemento diagonal correspondiente a la fila fija dividido por el primer elemento de la fila correspondiente.

Luego, restamos las dos filas de la fila fija para obtener el elemento en la segunda columna debajo de la diagonal principal. Las filas subsiguientes también se pueden reducir en el mismo procedimiento.

Como hemos mencionado, la primera fila siempre permanecerá sin cambios. También hemos calculado el factor y el vector b para cada fila como se muestra a continuación.

for i in range(k+1, n):
    # first row should remain intact
    if a[i, k] == 0:
        continue
    # calculating the factor
    factor = a[k, k]/a[i, k]
    for j in range(k, n):
        a[i, j] = a[k, j] - a[i, j]*factor
        # calculating the b vector
    b[i] = b[k] - b[i]*factor
print(a)
print(b)

El término sustitución hacia atrás se refiere al hecho de que los valores de x se pueden calcular desde la última ecuación hasta la primera ecuación. En este caso, aprovechamos los valores correspondientes del vector escalar b multiplicados por los coeficientes de x. El resultado se divide luego por el elemento correspondiente en la diagonal principal.

Hemos priorizado el cálculo de valores de x a partir de n y el resto de valores de x desde n-1 hasta el primero. Por otro lado, la suma debe comenzar desde j+1 ya que solo conocemos los valores de x de los cálculos anteriores. Además, el valor de la suma es el resultado de los valores de x multiplicados por los valores de a. Ahora podemos calcular el valor de x ya que tenemos el valor de la suma.

# calculating the value of x beginning from n-1
x[n - 1] = b[n - 1] / a[n - 1, n - 1]
for i in range(n - 2, -1, -1):
    sum_ax = 0
    # calculating the value of the sum
    for j in range(i + 1, n):
        sum_ax += a[i, j] * x[j]
    # calculating the value of x
    x[i] = (b[i] - sum_ax) / a[i, i]

print(x)

Producción :

[[ 2.00000000e+00  3.00000000e+00  4.00000000e+00  1.00000000e+00
   7.00000000e+00]
 [ 0.00000000e+00  7.00000000e+00 -1.00000000e+00  3.00000000e+00
   1.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00 -1.30000000e+01  2.00000000e+00
  -2.10000000e+01]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  2.81250000e+00
   5.81250000e+00]
 [ 0.00000000e+00  8.88178420e-16  0.00000000e+00 -4.44089210e-16
   4.95734797e+00]]
[  7.           5.         -14.           0.625        0.15371622]
The solution of the system is:
[0.02170543 0.79224806 1.05116279 0.15813953 0.03100775]

La siguiente función NumPy también debería dar como resultado la misma respuesta.

print("The numpy solution:")
print(linalg.solve(a, b))

Producción :

Solution by Nuumpy:
[0.02170543 0.79224806 1.05116279 0.15813953 0.03100775]