Algoritmos de resolución indirecta

Algoritmos de resolución indirecta#

Un esquema iterativo mejora la solución hasta alcanzar convergen­cia. Los más comunes son el método de Jacobi y el de Gauss­-Seidel. Ambos utilizan el mismo principio, pero el esquema difiere en la forma de realizar la sustitución.

Método de Jacobi#

Considere que la matriz \(\bf A\) se puede dividir en tres partes que corresponden a tres conjuntos de coeficientes. De modo que la ecuación \(\bf AX = B\) es

(60)#\[ \bf \left(L + D + U \right) X = B, \]

donde \(\bf L\) es una matriz triangular inferior con ceros en la diagonal, \(\bf D\) es la matriz diagonal que contiene los elementos de la diagonal de \(\bf A\), y \(\bf U\) es una matriz triangular superior con ceros en la diagonal.

El método de Jacobi para la \(r\)-­ésima iteración se escribe como

(61)#\[\begin{equation} \bf DX^{(r)} = -(L+U)X^{(r-1)}+B, \quad \left(r=0,1,2,\dots\right). \end{equation}\]

Nota

El método de Jacobi converge para matrices diagonalmente dominantes, en el sentido que es matemáticamente preciso.

Para determi­nar si el método converge, se tiene que la matriz de iteración

(62)#\[\begin{equation} \bf -D^{-1}\left(L+U \right), \end{equation}\]

además de un término aditivo, \(\bf D^{−1}B\), mapea un conjunto de \(\bf X\) en el siguiente.

La rapidez de convergencia de este método está dada por la rapidez de disminución del eigenvalor más lento. Es decir, el factor con el módulo más grande. El módu­lo de este factor está acotado entre 0 y 1, y se le llama radio espectral del operador de relajación \((\rho_{s})\). El número de iteraciones \((r)\) necesarias para reducir el error total por un factor \((10^{-p})\) se aproxima como

(63)#\[\begin{equation} r \approx \frac{p\, \ln(10)}{-\ln(\rho_{s})}. \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([
    [9, 0, 1, 1, 3],
    [2, 8, 2, 1, 0],
    [1, 2, 8, 1, 2],
    [2, 3, 1, 7, 1],
    [3, 1, 0, 2, 9]
], dtype=float)

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

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

# Inicialización de matrices D, L y U en ceros.
D = np.zeros((N, N))
L = np.zeros((N, N))
U = np.zeros((N, N))

# Se asigna la diagonal de A a un vector y luego a la matriz D.
K1 = np.diag(A)
D = np.diag(K1)

# Se eliminan los elementos de la diagonal de la matriz A para trabajar con los fuera de diagonal.
A_off = A - D

# Ciclo para extraer las matrices triangulares L y U
for k in range(N - 1):
    L[k:N, k] = A_off[k:N, k]
    U[k, k:N] = A_off[k, k:N]

# Cálculo del número de iteraciones y convergencia
# Matriz de iteración de Jacobi: F = -D^(-1) * (L + U)
F = -np.linalg.inv(D) @ (L + U)

# Cálculo de los valores propios y obtención del radio espectral (Rs).
eigenvalues = np.linalg.eigvals(F)
Rs = np.max(np.abs(eigenvalues))  # Módulo del valor propio de disminución más lento.

# Régimen de convergencia para estimar iteraciones.
p = 6
Ni = (p * np.log(10)) / (-np.log(Rs))  # Número de iteraciones en valor real.
Nm = int(np.fix(Ni)) + 1               # Número de iteraciones enteras.

# Proceso iterativo de Jacobi
# Matriz donde se guardan los resultados de cada iteración (N filas x Nm columnas).
X = np.zeros((N, Nm))

# Condiciones iniciales: vector de unos para la primera iteración (columna 0 en Python).
X[:, 0] = np.ones(N)

# Ciclo iterativo del método.
for k in range(1, Nm):
    X[:, k] = F @ X[:, k - 1] + np.linalg.inv(D) @ B

# Mostrar los resultados de todas las iteraciones en pantalla
print("Matriz de resultados por iteración (X):")
print(X)
print(f"\nSolución final obtenida (última iteración): {X[:, -1]}")

Matriz de resultados por iteración (X):
[[ 1.          0.22222222  0.79761905  0.35394621  0.64036255  0.4327246
   0.57593645  0.47557464  0.54565687  0.49671123  0.53090766  0.50702196
   0.5237076   0.51205207  0.52019395  0.5145065   0.51847942  0.51570417
   0.5176428   0.51628859  0.51723456  0.51657376  0.51703536  0.51671291
   0.51693815  0.51678081  0.51689072  0.51681395  0.51686758  0.51683011
   0.51685628  0.516838    0.51685077  0.51684185  0.51684808  0.51684373
   0.51684677  0.51684465  0.51684613]
 [ 1.         -0.25        0.31051587 -0.16889881  0.14958211 -0.07460709
   0.08229016 -0.02707329  0.04940606 -0.00399617  0.03331164  0.00725097
   0.02545534  0.01273875  0.0216218   0.01541661  0.01975119  0.01672331
   0.01883841  0.01736092  0.01839301  0.01767205  0.01817567  0.01782387
   0.01806962  0.01789796  0.01801787  0.01793411  0.01799262  0.01795175
   0.0179803   0.01796035  0.01797428  0.01796455  0.01797135  0.0179666
   0.01796992  0.0179676   0.01796922]
 [ 1.          0.25        1.08829365  0.5713459   0.94450054  0.68468054
   0.86578406  0.73905067  0.82750583  0.76569845  0.80887044  0.7787128
   0.79977929  0.78506353  0.79534313  0.7881624   0.79317844  0.78967453
   0.79212216  0.79041239  0.79160673  0.79077243  0.79135523  0.79094812
   0.7912325   0.79103385  0.79117262  0.79107568  0.7911434   0.7910961
   0.79112914  0.79110606  0.79112218  0.79111092  0.79111878  0.79111329
   0.79111713  0.79111444  0.79111632]
 [ 1.         -0.42857143  0.57936508 -0.04724112  0.42713057  0.10686841
   0.33314531  0.1755009   0.28564396  0.2086875   0.26243608  0.22488772
   0.25111619  0.23279441  0.24559294  0.23665263  0.24289781  0.23853529
   0.2415827   0.23945396  0.24094097  0.23990223  0.24062784  0.24012097
   0.24047504  0.24022771  0.24040048  0.24027979  0.24036409  0.2403052
   0.24034634  0.24031761  0.24033768  0.24032366  0.24033345  0.24032661
   0.24033139  0.24032805  0.24033038]
 [ 1.          0.          0.71560847  0.23754409  0.57794916  0.34167434
   0.50696627  0.39151221  0.47214973  0.41581505  0.45516528  0.42767703
   0.44687863  0.4334655   0.44283514  0.43629005  0.44086207  0.43766832
   0.43989929  0.43834087  0.43942949  0.43866904  0.43920025  0.43882918
   0.43908838  0.43890732  0.4390338   0.43894545  0.43900716  0.43896405
   0.43899417  0.43897313  0.43898783  0.43897756  0.43898473  0.43897972
   0.43898322  0.43898078  0.43898248]]

Solución final obtenida (última iteración): [0.51684613 0.01796922 0.79111632 0.24033038 0.43898248]

Método de Gauss-Seidel#

Las ecuaciones (60) se necesitan escalar mediante la división entre los ele­mentos de la diagonal, de manera que \(\bf D\) será igual a la matriz identidad \(\bf I\). Mientras que el método de Jacobi transfiere todos los términos al lado derecho, excepto los términos de la diagonal, iterando como sigue

(64)#\[\begin{equation} \bf X^{(r+1)} = \left(-L-U \right)X^{(r)} + B, \quad (r=0,1,2,\dots), \end{equation}\]

el método de Gauss­Seidel introduce \(x_{1}^{(r+1)},x_{2}^{(r+1)},\dots,\) en el lado derecho una vez calculados. Por lo tanto, las ecuaciones de iteración están dadas como

(65)#\[\begin{equation} \bf X^{(r+1)} = -LX^{(r+1)} - UX^{(r)} + B, \end{equation}\]

o bien

(66)#\[\begin{equation} \bf (I+L)X^{(r+1)} = -UX^{(r)} + B. \end{equation}\]

O se puede escribir las ecuaciones de la siguiente forma

(67)#\[\begin{equation} \begin{aligned} x_{1}^{(r+1)} & = - \left( a_{12}\,x_{2}^{(r)} + a_{13}\,x_{3}^{(r)} + a_{14}\,x_{4}^{(r)} + \cdots + a_{1n}\,x_{n}^{(r)} \right) + b_{1}, \\ x_{2}^{(r+1)} & = - \left( a_{21}\,x_{1}^{(r+1)} + a_{23}\,x_{3}^{(r)} + a_{24}\,x_{4}^{(r)} + \cdots + a_{2n}\,x_{n}^{(r)} \right) + b_{2}, \\ x_{3}^{(r+1)} & = - \left( a_{31}\,x_{1}^{(r+1)} + a_{32}\,x_{2}^{(r+1)} + a_{34}\,x_{4}^{(r)} + \cdots + a_{3n}\,x_{n}^{(r)} \right) + b_{3}, \\ \vdots \quad &= \quad \quad \vdots \\ x_{n}^{(r+1)} & = - \left( a_{n1}\,x_{1}^{(r+1)} + a_{n2}\,x_{2}^{(r+1)} + a_{n3}\,x_{3}^{(r+1)} + \cdots + a_{n(n-1)}\,x_{n-1}^{(r+1)} \right) + b_{n}. \\ \end{aligned} \end{equation}\]

Nota

Cuando el método de Jacobi converge, el método de Gauss­-Seidel converge más rápido.

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

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

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

# Ciclo iterativo para normalizar los pivotes de las ecuaciones
for k in range(N):
    B[k] = B[k] / A[k, k]
    A[k, :] = A[k, :] / A[k, k]

# Inicialización de matrices D, L y U.
D = np.eye(N)         # La matriz D es la matriz identidad.
L = np.zeros((N, N))  # Se inicializa L en ceros.
U = np.zeros((N, N))  # Se inicializa U en ceros.

# Se elimina la diagonal de la matriz A para trabajar con los términos fuera de ella.
A = A - D

# Ciclo para hacer la asignación de las matrices triangulares
for k in range(N - 1):
    L[k:N, k] = A[k:N, k]
    U[k, k:N] = A[k, k:N]

# Cálculo del número de iteraciones y convergencia
# Matriz de iteración de Gauss-Seidel: F = -(D + L)^(-1) * U
inv_DL = np.linalg.inv(D + L)
F = -inv_DL @ U

# Cálculo de los valores propios y obtención del radio espectral (Rs).
eigenvalues = np.linalg.eigvals(F)
Rs = np.max(np.abs(eigenvalues))  # Módulo del valor propio de disminución más lento.

# Régimen de convergencia para estimar el número de iteraciones.
p = 6
Ni = (p * np.log(10)) / (-np.log(Rs))  # Número de iteraciones en valor real.
Nm = int(np.fix(Ni)) + 1               # Número de iteraciones enteras.

# Proceso iterativo del método de Gauss-Seidel
# Matriz de N filas x Nm columnas donde se guardan los resultados de cada iteración.
X = np.zeros((N, Nm))

# Condiciones iniciales: vector de unos para iniciar el proceso iterativo (columna 0).
X[:, 0] = np.ones(N)

# Ciclo iterativo del método.
for k in range(1, Nm):
    X[:, k] = F @ X[:, k - 1] + inv_DL @ B

# Mostrar en pantalla el resultado de todas las iteraciones
print("Matriz de resultados por iteración (X):")
print(X)
print(f"\nSolución final obtenida (última iteración): {X[:, -1]}")

Matriz de resultados por iteración (X):
[[ 1.         -0.44444444 -0.08626578 -0.04226703 -0.04337483 -0.04514031
  -0.04547892 -0.04548471 -0.04547266]
 [ 1.          0.34920635  0.55394536  0.5885329   0.589716    0.58872375
   0.58846134  0.58844181  0.58844776]
 [ 1.         -0.11309524 -0.09867789 -0.07093506 -0.06657988 -0.06677967
  -0.06700145 -0.06704382 -0.06704466]
 [ 1.          0.680839    0.55654442  0.54594568  0.54876863  0.54984041
   0.5499718   0.54995936  0.54995142]
 [ 1.          0.96006551  0.90812575  0.89440651  0.89311666  0.89334839
   0.89344985  0.89346346  0.89346253]]

Solución final obtenida (última iteración): [-0.04547266  0.58844776 -0.06704466  0.54995142  0.89346253]
Contenido