Algoritmos de resolución directa

Contenido

Algoritmos de resolución directa#

Suponga que se tienen \(n\) ecuaciones lineales con \(n\) incógnitas que se pueden escribir de la siguiente forma

(8)#\[\begin{split} \begin{aligned} a_{11}x_{1} + a_{12}x_{2} + \cdots + a_{1n}x_{n} &= b_{1}, \\ a_{21}x_{1} + a_{22}x_{2} + \cdots + a_{2n}x_{n} &= b_{2}, \\ \vdots \quad \quad \quad \vdots \quad \quad \quad \quad &= \vdots, \\ a_{n1}x_{1} + a_{n2}x_{2} + \cdots + a_{nn}x_{n} &= b_{n}, \\ \end{aligned} \end{split}\]

donde el conjunto de coeficientes \(a_{i,j}\left(i,j=1,2,\dots n \right)\) se representan con \(\bf A\) y los elementos \(x_{j}\) y \(b_{j}\left(j=1,2,\dots,n \right)\) se representan con los vectores \(\bf X\) y \(\bf B\), respectivamente. De modo que, utilizando álgebra matricial, el conjunto de ecuaciones se puede representar de la siguiente forma

(9)#\[\begin{equation} {\bf A}{\bf X} = {\bf B}. \end{equation}\]

Recordemos que estas ecuaciones tienen una solución única, si y sólo si el \(\text{det}(A) \neq 0\).

Hay dos tipos de métodos para resolver el sistema de ecuaciones (8)

Los métodos directos. Se basa en la eliminación de elementos. Por ejemplo, el método de Gauss, el método de Gauss-Jordan y el cálculo de la matriz inversa. También se tienen los métodos con base en la factorización: descomposición triangular inferior y superior, factorización de Doolittle-Crout, el método de Cholesky y la factorización QR.
Los métodos indirectos, o métodos iterativos. calculan una sucesión de aproximaciones que se espera converjan a un valor suficientemente cercano a la solución verdadera. Entre estos destacan el método de Jacobi y el método de Gauss-Seidel.

La ventaja de los métodos directos es que la cantidad de cálculos es fija. Caso contrario, para los métodos iterativos, cuyos cálculos deben continuar indefinidamente hasta que las soluciones tengan una convergencia apropiada.

Eliminación gaussiana#

Suponga un sistema de ecuaciones de la forma (8) donde todos los divisores son no nulos. Entonces, la primera ecuación se almacena para ser usada posteriormente y la variable \(x_{1}\) se elimina de las restantes \(n-1\) ecuaciones al restar un múltiplo apropiado de la primera ecuación de cada una de las otras ecuaciones. Para esto, la siguientes notaciones son utilizadas para hacer referencia a los coeficientes originales

(10)#\[\begin{equation} a_{ij}^{(1)} = a_{ij}, \quad \left(i,j=1,2,\dots,n \right), \end{equation}\]

(11)#\[\begin{equation} b_{i}^{(1)} = b_{i}, \quad \left(i=1,2,\dots,n \right), \end{equation}\]

mientras que esta es para los nuevos coeficientes que se encuentran al usar los multiplicadores

(12)#\[\begin{equation} m_{i1} = \frac{a_{i1}^{(1)}}{a_{11}^{(1)}}, \quad \left(i,j=2,3,\dots,n \right), \end{equation}\]

con ellos conformamos los nuevos elementos para la matriz \(\bf A\)

(13)#\[\begin{equation} a_{ij}^{(2)} = a_{ij}^{(1)} - m_{i1}\, a_{1j}^{(1)}, \quad \left(i=2,3,\dots,n \right), \quad \left(j=1,2,\dots,n \right), \end{equation}\]

mientras que para el vector \(\bf B\)

(14)#\[\begin{equation} b_{i}^{(2)} = b_{i}^{(1)} - m_{i1}\, b_{i}^{(1)}, \quad \left(i=1,2,\dots,n \right). \end{equation}\]

De esto se observa que los elementos de la primera columna \(j = 1\), tienen los siguientes valores

(15)#\[\begin{equation} a_{i1}^{(2)} = a_{i1}^{(1)} - \frac{a_{i1}^{(1)}}{a_{11}^{(1)}}a_{11}^{(1)} = 0 \end{equation}\]

por lo que la primera variable \(x_{1}\) se ha eliminado de las últimas \(n−1\) ecuaciones. Si se repite este procedimiento \(n −1\) veces, la ecuación restante tendrá una sola incógnita y se puede resolver en forma directa.

En cada etapa del proceso, cuando la variable \(x_{k}\) se va a eliminar, los multiplicadores son

(16)#\[\begin{equation} m_{ik} = \frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}, \quad \left(i=k+1,k+2,k+3,\dots,n \right), \end{equation}\]

y los nuevos elementos formados son, para la matriz \(\bf A\)

(17)#\[\begin{equation} a_{ij}^{(k+1)} = a_{ij}^{(k)} - m_{ik}\, a_{kj}^{(k)}, \quad \left(i=k+1,k+2,\dots,n \right), (j=k,k+1,k+2,\dots,n) \end{equation}\]

y para el vector \(\bf B\)

(18)#\[\begin{equation} b_{i}^{(k+1)} = b_{i}^{(k)} - m_{ik}\, b_{k}^{(k)}, \quad \left(i=k+1,k+2,\dots,n \right). \end{equation}\]

El resultado de este proceso de eliminación es un conjunto triangular de ecuaciones de la forma

\[\begin{split} \begin{aligned} a_{11}^{(1)}x_{1} + a_{12}^{(1)}x_{2} + \cdots + a_{1n}^{(1)}x_{n} &= b_{1}^{(1)}, \\ a_{22}^{(2)}x_{2} + \cdots + a_{2n}^{(2)}x_{n} &= b_{2}^{(2)}, \\ \vdots \quad \quad \quad \vdots \quad \quad \quad \quad &= \vdots, \\ a_{nn}^{(n)}x_{n} &= b_{n}^{(n)}, \\ \end{aligned} \end{split}\]

Con un proceso de sustitución hacia atrás o regresiva, la última ecuación tiene la solución

(19)#\[\begin{equation} x_{n} = \frac{b_{n}^{(n)}}{a_{nn}^{(n)}}. \end{equation}\]

import numpy as np
import matplotlib.pyplot as plt

# Suprimir advertencias menores de integración durante los bucles de optimización
import warnings
warnings.filterwarnings("ignore")

# Configurar el estilo de graficación para una mejor visualización
plt.style.use('ggplot')

# Matriz de coeficientes del sistema de ecuaciones lineales.
A = np.array([
    [5, 4, 8, 1, 6],
    [9, 5, 1, 6, 4],
    [4, 1, 9, 8, 6],
    [4, 5, 9, 7, 6],
    [6, 2, 9, 7, 1]
], dtype=float)

# Vector de términos independientes (entradas) del sistema.
B = np.array([9, 1, 2, 9, 1], dtype=float)

# Número de incógnitas del sistema de ecuaciones (calculado como el rango de A).
N = np.linalg.matrix_rank(A)

# Vector para almacenar la solución final.
x = np.zeros(N)

# Proceso de eliminación gaussiana
for k in range(N - 1):
    for m in range(k + 1, N):
        # Cálculo del factor multiplicador para anular los elementos debajo de la diagonal.
        MT = -A[m, k] / A[k, k]
        
        # Modificación de la fila m de la matriz A.
        A[m, :] = A[m, :] + MT * A[k, :]
        
        # Modificación del elemento m del vector B.
        B[m] = B[m] + MT * B[k]

# Proceso de sustitución regresiva
# Se calcula primero la última incógnita.
x[N - 1] = B[N - 1] / A[N - 1, N - 1]

# Se resuelven las demás incógnitas desde la penúltima hasta la primera.
for k in range(N - 2, -1, -1):
    suma = 0
    # Sumatoria de los términos con incógnitas ya despejadas.
    for m in range(k + 1, N):
        suma += A[k, m] * x[m]
    
    # Cálculo de la incógnita actual.
    x[k] = (B[k] - suma) / A[k, k]

# Mostrar la solución obtenida
print("Vector solución x:", x)

Vector solución x: [-0.69044297  1.63071155  0.5417684  -0.47715382  0.34539588]

Eliminación de Gauss-Jordan#

Para aplicar este método, se procede de una manera similar a la eliminación gaussiana, pero en cada etapa se elimina la variable \(x_{k}\), no sólo de las ecuaciones sucesivas, sino también de todas las precedentes. La eliminación de Jordan necesita aproximadamente \(\frac{n^{3}}{2}\) operaciones, en contraste con las \(\frac{n^{3}}{3}\) de la eliminación gaussiana.

Importante

No es recomendable para el uso general.

# Matriz de coeficientes del sistema de ecuaciones lineales.
A = np.array([
    [8, 4, 6, 6, 9],
    [5, 3, 2, 6, 3],
    [5, 9, 5, 1, 2],
    [9, 3, 4, 4, 9],
    [8, 9, 8, 7, 9]
], dtype=float)

# Vector de términos independientes (entradas) del sistema.
B = np.array([3, 3, 3, 5, 7], dtype=float)

# Número de incógnitas del sistema de ecuaciones (calculado como el rango de A).
N = np.linalg.matrix_rank(A)

# Proceso de eliminación Gauss-Jordan
for k in range(N):
    # Normalización del pivote: se divide el elemento k de B y la fila k de A 
    # por el elemento pivote A[k, k] para que la diagonal principal sea 1.
    B[k] = B[k] / A[k, k]
    A[k, :] = A[k, :] / A[k, k]
    
    # Eliminación: hacer ceros en todos los demás elementos de la columna k.
    for m in range(N):
        if m != k:
            # Factor multiplicador para anular el elemento A[m, k]
            MT = -A[m, k]
            
            # Modificación de la fila m de la matriz A.
            A[m, :] = A[m, :] + MT * A[k, :]
            
            # Modificación del elemento m del vector B.
            B[m] = B[m] + MT * B[k]

# NOTA: En el método de Gauss-Jordan, la matriz A se convierte en la matriz identidad 
# y la solución queda almacenada directamente en el vector B.
print("Vector solución x:", B)

Vector solución x: [-0.81301939  1.76315789 -2.52770083  0.23961219  1.79778393]

Inversa de una matriz#

En el caso de un conjunto de ecuaciones con diferentes entradas de \(\bf B\), se puede encontrar la solución calculando la inversa de la matriz \(\bf A\), para lo cual se requiere que ésta sea no singular. Ésta se determina al resolver \(\bf X\) de la ecuación \(\bf AX = I\), donde la matriz \(\bf A\) tiene los coeficientes del lado izquierdo y el lado derecho, \(\bf I\), representa a la matriz identidad

(20)#\[\begin{equation} I = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \cdots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}. \end{equation}\]

Resolviendo cada uno de estos vectores columna de \(\bf I\) se obtienen las columnas de la matriz inversa \(\bf A^{-1}\)

(21)#\[\begin{equation} A^{-1} = \frac{\text{adj}(A)}{\text{det}(A)}. \end{equation}\]

Truco

El determinante de una matriz se puede calcular como producto del proceso de la eliminación gaussiana.

# Matriz de coeficientes del sistema de ecuaciones lineales.
A = np.array([
    [9, 2, 9, 2, 5],
    [6, 1, 5, 2, 8],
    [6, 7, 8, 3, 7],
    [1, 1, 7, 1, 5],
    [9, 7, 2, 3, 5]
], dtype=float)

# Vector de términos independientes (entradas) del sistema.
B = np.array([4, 1, 6, 1, 5], dtype=float)

# Número de incógnitas del sistema de ecuaciones (calculado como el rango de A).
N = np.linalg.matrix_rank(A)

# Matriz unitaria (identidad de tamaño N x N).
Id = np.eye(N)

# Proceso de eliminación Gauss-Jordan para calcular la inversa
for k in range(N):
    # Normalización del pivote en la matriz Identidad (Id) y en A.
    Id[k, :] = Id[k, :] / A[k, k]
    A[k, :] = A[k, :] / A[k, k]
    
    # Eliminación: hacer ceros en todos los demás elementos de la columna k.
    for m in range(N):
        if m != k:
            # Factor multiplicador para anular el elemento A[m, k]
            MT = -A[m, k]
            
            # Modificación de la fila m de la matriz A.
            A[m, :] = A[m, :] + MT * A[k, :]
            
            # Modificación de la fila m de la matriz Id (acumulando las operaciones inversas).
            Id[m, :] = Id[m, :] + MT * Id[k, :]

# NOTA: En este punto, la matriz A se convierte en la matriz identidad 
# y la matriz 'Id' almacena la matriz inversa de A.

# Cálculo del vector solución multiplicando la matriz inversa (Id) por el vector B.
x = np.dot(Id, B)

print("Matriz Inversa (Id):")
print(Id)
print("\nVector solución x:")
print(x)

Matriz Inversa (Id):
[[-1.98019802e-02 -7.92079208e-01 -1.83168317e+00  2.21782178e+00
   1.63366337e+00]
 [-3.16831683e-01 -1.67326733e+00 -3.30693069e+00  4.48514851e+00
   3.13861386e+00]
 [ 9.90099010e-02 -3.96039604e-02  1.58415842e-01 -8.91089109e-02
  -1.68316832e-01]
 [ 1.27722772e+00  8.08910891e+00  1.76435644e+01 -2.30495050e+01
  -1.58712871e+01]
 [-3.26732673e-01 -1.06930693e+00 -2.72277228e+00  3.59405941e+00
   2.45544554e+00]]

Vector solución x:
[-1.47524752 -2.6039604   0.37623762 16.65346535 -2.84158416]