En álgebra lineal numérica , el método de Jacobi es un algoritmo iterativo para determinar las soluciones de un sistema de ecuaciones lineales estrictamente diagonalmente dominante . Se resuelve cada elemento de la diagonal y se inserta un valor aproximado. Luego, el proceso se itera hasta que converge. Este algoritmo es una versión simplificada del método de transformación de Jacobi de diagonalización matricial . El método lleva el nombre de Carl Gustav Jacob Jacobi .
Dejar
ser un sistema cuadrado de n ecuaciones lineales, donde:
Entonces A se puede descomponer en un componente diagonal D , una parte triangular inferior L y una parte triangular superior U :
La solución se obtiene luego iterativamente a través de
dónde es la k- ésima aproximación o iteración de y es la siguiente o k + 1 iteración de. La fórmula basada en elementos es así:
El cálculo de requiere cada elemento en x ( k ) excepto él mismo. A diferencia del método Gauss-Seidel , no podemos sobrescribir con , ya que ese valor será necesario para el resto del cálculo. La cantidad mínima de almacenamiento es de dos vectores de tamaño n .
Entrada: conjetura inicial a la solución , matriz (diagonal dominante) , vector del lado derecho , criterio de convergencia Salida: solución cuando se alcanza la convergencia Comentarios: pseudocódigo basado en la fórmula basada en elementos anterior mientras no se alcanza la convergencia hacer para i: = 1 paso hasta n hacer para j: = 1 paso hasta n hacer si j ≠ i entonces fin fin final final
La condición de convergencia estándar (para cualquier método iterativo) es cuando el radio espectral de la matriz de iteración es menor que 1:
Una condición suficiente (pero no necesaria) para que el método converja es que la matriz A sea estricta o irreductiblemente diagonalmente dominante . El dominio diagonal estricto de filas significa que para cada fila, el valor absoluto del término diagonal es mayor que la suma de los valores absolutos de otros términos:
El método de Jacobi a veces converge incluso si no se cumplen estas condiciones.
Tenga en cuenta que el método de Jacobi no converge para todas las matrices definidas positivas simétricas . Por ejemplo
Ejemplo 1
Un sistema lineal de la forma con estimación inicial es dado por
Usamos la ecuación , descrito anteriormente, para estimar . Primero, reescribimos la ecuación en una forma más conveniente, dónde y . De los valores conocidos
determinamos como
Más, se encuentra como
Con y calculado, estimamos como :
La siguiente iteración produce
Este proceso se repite hasta la convergencia (es decir, hasta es pequeño). La solución después de 25 iteraciones es
Ejemplo 2
Supongamos que se nos da el siguiente sistema lineal:
Si elegimos (0, 0, 0, 0) como aproximación inicial, entonces la primera solución aproximada viene dada por
Utilizando las aproximaciones obtenidas, se repite el procedimiento iterativo hasta alcanzar la precisión deseada. Las siguientes son las soluciones aproximadas después de cinco iteraciones.
| | | |
---|
0,6 | 2.27272 | -1,1 | 1.875 |
1.04727 | 1.7159 | -0,80522 | 0.88522 |
0,93263 | 2.05330 | -1.0493 | 1.13088 |
1.01519 | 1.95369 | -0,9681 | 0,97384 |
0,98899 | 2.0114 | -1.0102 | 1.02135 |
La solución exacta del sistema es (1, 2, −1, 1) .
Ejemplo de Python
importar numpy como npITERATION_LIMIT = 1000# inicializar la matrizA = np . matriz ([[ 10. , - 1. , 2. , 0. ], [ - 1. , 11. , - 1. , 3. ], [ 2. , - 1. , 10. , - 1. ], [ 0.0 , 3. , - 1. , 8. ]])# inicializar el vector RHSb = np . matriz ([ 6. , 25. , - 11. , 15. ])# imprime el sistemaimprimir ( "Sistema:" )para i en rango ( A . forma [ 0 ]): fila = [ " {} * x {} " . formato ( A [ i , j ], j + 1 ) para j en rango ( A . forma [ 1 ])] print ( "+" . join ( fila ), "=" , b [ i ])imprimir ()x = np . ceros_como ( b )para it_count en rango ( ITERATION_LIMIT ): if it_count ! = 0 : print ( "Iteración {0} : {1} " . formato ( it_count , x )) x_new = np . ceros_como ( x ) para i en rango ( A . forma [ 0 ]): s1 = np . punto ( A [ i , : i ], x [: i ]) s2 = np . punto ( A [ i , i + 1 :], x [ i + 1 :]) x_nuevo [ i ] = ( b [ i ] - s1 - s2 ) / A [ i , i ] si x_new [ i ] == x_new [ i - 1 ]: rotura si np . allclose ( x , x_new , atol = 1e-10 , rtol = 0. ): rotura x = x_nuevoimprimir ( "Solución:" )imprimir ( x )error = np . punto ( A , x ) - bimprimir ( "Error:" )imprimir ( error )