Jacobi y Gauss-Seidel

El método de Gauss-Seidel es un método iterativo para resolver sistemas de ecuaciones lineales. El método consiste en despejar una variable de una ecuación y sustituirla en otra ecuación. Este proceso se repite hasta que se obtiene una solución aproximada.

Para deducir el método de Gauss-Seidel, se parte de un sistema de ecuaciones lineales de la forma:

\[Ax = b\] donde la matriz \(A\) es una matriz cuadrada de \(n \times n\), \(x\) es un vector de \(n\) elementos y \(b\) es un vector de \(n\) elementos. De esta forma la matriz \(A\) se puede descomponer en la suma de dos matrices \(D\) y \(M\) de la siguiente forma:

\[A = L+D+U\] donde \(D\) es la matriz de los elementos de la diagonal y \(L\) es la matriz triangular estrictamente inferior con \(a_{ij}=l_{ij}\) si \(i>j\) y \(l_{ii}=0\) y \(U\) es la matriz estrictamente triangular superior, es decir \(a_{ij}=u_{ij}\) si \(i<j\) y \(u_{ii}=0\) . De esta forma, el sistema de ecuaciones se puede reescribir como:

\[Dx = b - (L+U)x\] De esta forma podemos crear un método iterativo de la forma \[x^{n+1} = D^{-1}(b - (L+U)x^n)\]

El método de Jacobi es una generalización del método del punto fijo. Los críterios de convergencia de este método es

Pseudocódigo

Dados \(A_{n\times n}\) y \(b_{n\times 1}\), max_iteraciones, tolerancia y \(x_0\):

Para k=0,1,2,…,max_iteraciones:

Para i=1,2,...,n:

\[x_i^{k+1} = \frac{1}{a_{ii}}(b_i - \sum_{j=1, j \neq i}^{n} a_{ij}x_j^k)\]

si \(||x^{k+1} - x^k||\) < tolerancia$ para todo \(i=1,2,...,n\): Imprima la solución y detenga el algoritmo sino \[x^k = x^{k+1}\]

Ejercicio observe el siguietne código si es necesario,comente cada línea de código, compruebe que corresponde al pseudocódigo anterior (en caso contrario realice los respectivos cambios) y ejecute el código.

import numpy as np

def jacobi(A, b, x0, tol, max_iter):
    n = len(b)
    x = x0.copy()
    iter_count = 0
    
    while iter_count < max_iter:
        iter_count += 1
        new_x = np.zeros(n)
        
        for i in range(n):
            sum_ = 0
            for j in range(n):
                if j != i:
                    sum_ += A[i][j] * x[j]
            new_x[i] = (b[i] - sum_) / A[i][i]
        
        # Comprobar convergencia
        if np.all(np.abs(x - new_x) < tol):
            return new_x
        
        x = new_x.copy()
    
    raise Exception("El método de Jacobi no convergió después de {} iteraciones.".format(max_iter))

# Matriz de coeficientes A y vector de términos constantes b
A = np.array([[3, 1, -1],
              [1, 4, 1],
              [2, -1, 3]])

b = np.array([1, 5, 3])

# Estimación inicial
x0 = np.array([0, 0, 0])

# Tolerancia y número máximo de iteraciones
tolerancia = 1e-6
max_iteraciones = 1000

# Llamada a la función para resolver el sistema
solucion = jacobi(A, b, x0, tolerancia, max_iteraciones)

# Imprimir la solución
print("La solución aproximada es:", solucion)
La solución aproximada es: [0.38297877 0.89361697 1.04255282]

Metodo de Gauss-Seidel

Pseudocódigo Dados \(A_{n\times n}\) y \(b_{n\times 1}\), max_iteraciones, tolerancia y \(x_0\):

Para k=0,1,2,…,max_iteraciones:

Para i=1,2,…,n:

\[x_i^{k+1} = \frac{1}{a_{ii}}(b_i - \sum_{j=1}^{i-1} a_{ij}x_j^{k+1} - \sum_{j=i+1}^{n} a_{ij}x_j^k)\]

si \(||x^{k+1} - x^k|| < tolerancia\) para todo \(i=1,2,...,n\): Imprima la solución y detenga el algoritmo

sino \[x^k = x^{k+1}\]

Ejercicio observe el siguietne código si es necesario,comente cada línea de código, compruebe que corresponde al pseudocódigo anterior (en caso contrario realice los respectivos cambios) y ejecute el código.

import numpy as np
def gaussseidel(A, b, x0, tol, max_iter):
    n = len(b)
    x = x0.copy()
    iter_count = 0
    
    while iter_count < max_iter:
        iter_count += 1
        new_x = np.zeros(n)
        
        for i in range(n):
            sum_ = 0
            for j in range(n):
                if j != i:
                    if j < i:
                        sum_ += A[i][j] * new_x[j]
                    else:
                        sum_ += A[i][j] * x[j]
            new_x[i] = (b[i] - sum_) / A[i][i]
        
        # Comprobar convergencia
        if np.all(np.abs(x - new_x) < tol):
            return new_x
        
        x = new_x.copy()
    
    raise Exception("El método de Gauss-Seidel no convergió después de {} iteraciones.".format(max_iter))


# Matriz de coeficientes A y vector de términos constantes b

A = np.array([[3, 1, -1],
              [1, 4, 1],
              [2, -1, 3]])

b = np.array([1, 5, 3])

# Estimación inicial
x0 = np.array([0, 0, 0])

# Tolerancia y número máximo de iteraciones
tolerancia = 1e-6
max_iteraciones = 1000

# Llamada a la función para resolver el sistema
solucion = gaussseidel(A, b, x0, tolerancia, max_iteraciones)

print("La solución aproximada es:", solucion)
La solución aproximada es: [0.38297865 0.89361694 1.04255321]

Note que el método de Gauss-Seidel se puede escribir de forma matricial como:

\[Ax = b\] \[(L+D)x^{k+1} = b- Ux^k\] \[Dx^{k+1} = b-Lx^{k+1} Ux^k \] \[x^{k+1} = D^{-1}(b-Lx^{k+1} - Ux^k)\]

Note que los métodos de Jacobi y Gauus-Seidel son muy similares, la diferencia radica en que en el método de Gauss-Seidel se utiliza la solución de la iteración actual para calcular la siguiente iteración. Además, se puede escribir de forma matricial como:

\[x^{k+1}=Bx^k+c\]

donde para el método de Jacobi \[B_{jacobi}=D^{-1}(L+U)\] y para el método de Gauss-Seidel \[B_{gauss-seidel}=(D-L)^{-1}U\]

La convergencia de los dos métodos se estudia a partir del radio espectral de la matriz \(B\), es decir los autovalores de \(B\).

Teorema

Un método iterativo de la forma \(x^{k+1}=Bx^k+c\) converge si y solo si \(\rho(B)<1\), donde \(\rho(B)\) es el radio espectral de \(B\)., en otras palabras, el método converge si \[max_{\lambda \in \sigma(B)}|\lambda|<1,\] donde \(\lambda\) es un autovalor de \(B\) y \(\sigma(B)\) es el conjunto de autovalores de \(B\).

Convergencia de los métodos

Mostrar que que un método converge no es simple, pero podemos establecer algunos criterios para determinar si un método converge o no.

  • Lo métodos iterativos de Jacobi y Gauss-Seidel convergen si la matriz \(A\) es diagonal dominante, es decir que el valor absoluto de la diagonal debe ser mayor que la suma de los valores absolutos de los elementos de la fila. \[|a_{ii}| > \sum_{j=1, j \neq i}^{n} |a_{ij}|\]

  • Si la matriz \(A\) es simétrica y definida positiva, entonces el método de Gauss-Seidel converge.

Note un sistema puede converger usando Jacobi o Gauss-Seidel pero podría no satisfacer alguno de estos criterios.

Ejercicio

Sea el siguiente sistema

\[\begin{bmatrix} 8 & 2 & 1 \\ 5 & 4 & 1 \\ 4 & 3 & 4 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 1 \\ 5 \\ 3 \end{bmatrix}\] + Muestre que la matriz \(A\) no es diagonal dominante. + Muestre que la matriz \(A\) no es simétrica y definida positiva.

  • Encuentre la matriz \(B\) para el método de Jacobi y Gauss-Seidel, y encuentre el radio espectral de \(B\), para cada caso.

  • ¿El método de Jacobi converge? ¿El método de Gauss-Seidel converge?

Ejercicio

Sea el siguiente sistema

\[Ax=b\]

donde

\[A=\begin{bmatrix} 5&1&1&0&0&0&0&0\\1&5&1&1&0&0&0&0\\1&1&5&1&1&0&0&0\\0&1&1&5&1&1&0&0\\ 0&0&1&1&5&1&1&0\\0&0&0&1&1&5&1&1\\ 0&0&0&0&1&1&5&1\\ 0&0&0&0&0&1&1&5\\ \end{bmatrix}\]