Método de Jacobi

En álgebra lineal numérica , el método de Jacobi (también conocido como método de iteración de Jacobi ) es un algoritmo iterativo para determinar las soluciones de un sistema de ecuaciones lineales estrictamente diagonalmente dominante . Se resuelve cada elemento diagonal y se sustituye por un valor aproximado. El proceso se itera hasta que converge. Este algoritmo es una versión simplificada del método de transformación de Jacobi para la diagonalización de matrices . El método recibe su nombre de Carl Gustav Jacob Jacobi .

Descripción

Dejar $A\mathbf {x} =\mathbf {b}$ Sea un sistema cuadrado de n ecuaciones lineales, donde: $A={\begin{bmatrix}a_{11}&a_{12}&\cdots &a_{1n}\\a_{21}&a_{22}&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &a_{nn}\end{bmatrix}},\qquad \mathbf {x} ={\begin{bmatrix}x_{1}\\x_{2}\\\vdots \\x_{n}\end{bmatrix}},\qquad \mathbf {b} ={\begin{bmatrix}b_{1}\\b_{2}\\\vdots \\b_{n}\end{bmatrix}}.$

Cuando $A$ y $\mathbf {b}$ son conocidos y $\mathbf {x}$ Si se desconoce, podemos usar el método de Jacobi para aproximarlo. $\mathbf {x}$ . El vector $\mathbf {x} ^{(0)}$ denota nuestra suposición inicial para $\mathbf {x}$ (a menudo $\mathbf {x} _{i}^{(0)}=0$ para $i=1,2,...,n$ ). Denotamos $\mathbf {x} ^{(k)}$ como la k- ésima aproximación o iteración de $\mathbf {x}$ , y $\mathbf {x} ^{(k+1)}$ es la siguiente (o k +1) iteración de $\mathbf {x}$ .

Fórmula basada en matrices

Entonces A se puede descomponer en un componente diagonal D , una parte triangular inferior L y una parte triangular superior U : $A=D+L+U\qquad {\text{where}}\qquad D={\begin{bmatrix}a_{11}&0&\cdots &0\\0&a_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &a_{nn}\end{bmatrix}}{\text{ and }}L+U={\begin{bmatrix}0&a_{12}&\cdots &a_{1n}\\a_{21}&0&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &0\end{bmatrix}}.$ La solución se obtiene entonces de forma iterativa mediante

\mathbf {x} ^{(k+1)}=D^{-1}(\mathbf {b} -(L+U)\mathbf {x} ^{(k)}).

Fórmula basada en elementos

La fórmula basada en elementos para cada fila $i$ es así: $x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j\neq i}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\ldots ,n.$ El cálculo de $x_{i}^{(k+1)}$ requiere cada elemento en $\mathbf {x} ^{(k)}$ excepto ella misma. A diferencia del método de Gauss-Seidel , no podemos sobrescribir $x_{i}^{(k)}$ con $x_{i}^{(k+1)}$ , ya que ese valor será necesario para el resto del cálculo. La cantidad mínima de almacenamiento son dos vectores de tamaño n .

Algoritmo

Entrada: estimación inicial x ⁽⁰⁾ de la solución , matriz A (diagonal dominante) , vector del lado derecho b , criterio de convergencia. Salida: solución cuando se alcanza la convergencia. Comentarios: pseudocódigo basado en la fórmula basada en elementos anterior. k = 0 mientras no se alcance la convergencia hacer para i := 1 paso hasta n hacer σ = 0 para j := 1 paso hasta n hacer si j ≠ i entonces σ = σ + a _ij x _j^{( k )} fin fin x _i^{( k +1)} = ( b _i − σ ) / a _ii fin incremento k fin

Convergencia

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:

\rho (D^{-1}(L+U))<1.

Una condición suficiente (pero no necesaria) para que el método converja es que la matriz A sea estrictamente o irreduciblemente diagonalmente dominante . La dominancia diagonal estricta por filas significa que, para cada fila, el valor absoluto del término diagonal es mayor que la suma de los valores absolutos de los demás términos:

\left|a_{ii}\right|>\sum _{j\neq i}{\left|a_{ij}\right|}.

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 simétricas definidas positivas . Por ejemplo, $A={\begin{pmatrix}29&2&1\\2&6&1\\1&1&{\frac {1}{5}}\end{pmatrix}}\quad \Rightarrow \quad D^{-1}(L+U)={\begin{pmatrix}0&{\frac {2}{29}}&{\frac {1}{29}}\\{\frac {1}{3}}&0&{\frac {1}{6}}\\5&5&0\end{pmatrix}}\quad \Rightarrow \quad \rho (D^{-1}(L+U))\approx 1.0661\,.$

Ejemplos

Ejemplo de pregunta

Un sistema lineal de la forma $Ax=b$ con estimación inicial $x^{(0)}$ es dado por

A={\begin{bmatrix}2&1\\5&7\\\end{bmatrix}},\ b={\begin{bmatrix}11\\13\\\end{bmatrix}}\quad {\text{and}}\quad x^{(0)}={\begin{bmatrix}1\\1\\\end{bmatrix}}.

Usamos la ecuación $x^{(k+1)}=D^{-1}(b-(L+U)x^{(k)})$ , descrito anteriormente, para estimar $x$ Primero, reescribimos la ecuación en una forma más conveniente. $D^{-1}(b-(L+U)x^{(k)})=Tx^{(k)}+C$ , dónde $T=-D^{-1}(L+U)$ y $C=D^{-1}b$ A partir de los valores conocidos $D^{-1}={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}},\ L={\begin{bmatrix}0&0\\5&0\\\end{bmatrix}}\quad {\text{and}}\quad U={\begin{bmatrix}0&1\\0&0\\\end{bmatrix}}.$ determinamos $T=-D^{-1}(L+U)$ como $T={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}}\left\{{\begin{bmatrix}0&0\\-5&0\\\end{bmatrix}}+{\begin{bmatrix}0&-1\\0&0\\\end{bmatrix}}\right\}={\begin{bmatrix}0&-1/2\\-5/7&0\\\end{bmatrix}}.$ Más, $C$ se encuentra como $C={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}}{\begin{bmatrix}11\\13\\\end{bmatrix}}={\begin{bmatrix}11/2\\13/7\\\end{bmatrix}}.$ Con $T$ y $C$ Calculado, estimamos $x$ como $x^{(1)}=Tx^{(0)}+C$ : $x^{(1)}={\begin{bmatrix}0&-1/2\\-5/7&0\\\end{bmatrix}}{\begin{bmatrix}1\\1\\\end{bmatrix}}+{\begin{bmatrix}11/2\\13/7\\\end{bmatrix}}={\begin{bmatrix}5.0\\8/7\\\end{bmatrix}}\approx {\begin{bmatrix}5\\1.143\\\end{bmatrix}}.$ La siguiente iteración produce $x^{(2)}={\begin{bmatrix}0&-1/2\\-5/7&0\\\end{bmatrix}}{\begin{bmatrix}5.0\\8/7\\\end{bmatrix}}+{\begin{bmatrix}11/2\\13/7\\\end{bmatrix}}={\begin{bmatrix}69/14\\-12/7\\\end{bmatrix}}\approx {\begin{bmatrix}4.929\\-1.714\\\end{bmatrix}}.$ Este proceso se repite hasta la convergencia (es decir, hasta $\|Ax^{(n)}-b\|$ es pequeño). La solución después de 25 iteraciones es

x={\begin{bmatrix}7.111\\-3.222\end{bmatrix}}.

Ejemplo de pregunta 2

Supongamos que se nos da el siguiente sistema lineal:

{\begin{aligned}10x_{1}-x_{2}+2x_{3}&=6,\\-x_{1}+11x_{2}-x_{3}+3x_{4}&=25,\\2x_{1}-x_{2}+10x_{3}-x_{4}&=-11,\\3x_{2}-x_{3}+8x_{4}&=15.\end{aligned}}

Si elegimos $(0, 0, 0, 0)$ como aproximación inicial, entonces la primera solución aproximada viene dada por ${\begin{aligned}x_{1}&=(6+0-(2*0))/10=0.6,\\x_{2}&=(25+0+0-(3*0))/11=25/11=2.2727,\\x_{3}&=(-11-(2*0)+0+0)/10=-1.1,\\x_{4}&=(15-(3*0)+0)/8=1.875.\end{aligned}}$ Utilizando las aproximaciones obtenidas, se repite el procedimiento iterativo hasta alcanzar la precisión deseada. A continuación se muestran las soluciones aproximadas tras cinco iteraciones.

La solución exacta del sistema es $(1, 2, - 1, 1)$ .

Ejemplo de Python

import numpy as npLÍMITE_DE_ITERACIÓN = 1000# inicializar la matrizA = np . array ([[ 10. , - 1. , 2. , 0. ],[ - 1. , 11. , - 1. , 3. ],[ 2. , - 1. , 10. , - 1. ],[ 0.0 , 3. , - 1. , 8. ]])# Inicializar el vector RHSb = np.array ([ 6 . , 25. , -11 . , 15. ] )# imprime el sistemaimprimir ( "Sistema:" )para i en rango ( A . forma [ 0 ]):fila = [ f " { A [ i , j ] } *x { j + 1 } " para j en rango ( A . forma [ 1 ])]print ( f ' { " + " . join ( fila ) } = { b [ i ] } ' )imprimir ()x = np.zeros_like ( b )para it_count en range ( ITERATION_LIMIT ):Si it_count es distinto de 0 :print ( f "Iteración { it_count } : { x } " )x_nuevo = np.zeros_like ( x )para i en rango ( A . forma [ 0 ]):s1 = np.dot ( 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 ]:romperif np . allclose ( x , x_new , atol = 1e-10 , rtol = 0. ):romperx = x_nuevoimprimir ( "Solución: " )imprimir ( x )error = np.punto ( A , x ) - bimprimir ( "Error:" )imprimir ( error )

Método de Jacobi ponderado

La iteración de Jacobi ponderada utiliza un parámetro $\omega$ para calcular la iteración como

\mathbf {x} ^{(k+1)}=\omega D^{-1}(\mathbf {b} -(L+U)\mathbf {x} ^{(k)})+\left(1-\omega \right)\mathbf {x} ^{(k)}

con $\omega =2/3$ siendo la opción habitual. ^{[ 1 ]} De la relación $L+U=A-D$ , esto también puede expresarse como

{\begin{aligned}\mathbf {x} ^{(k+1)}&=\omega D^{-1}\mathbf {b} +\left(I-\omega D^{-1}A\right)\mathbf {x} ^{(k)}\\&=\mathbf {x} ^{(k)}+\omega D^{-1}\mathbf {r} ^{(k)},\end{aligned}}

dónde $\mathbf {r} ^{(k)}=\mathbf {b} -A\mathbf {x} ^{(k)}$ es el residuo algebraico en la iteración $k$ .

Convergencia en el caso simétrico definido positivo

Si la matriz del sistema $A$ Si es simétrica definida positiva , se puede demostrar la convergencia.

Dejar $C=C_{\omega }=I-\omega D^{-1}A$ sea la matriz de iteración. Entonces, la convergencia está garantizada para

\rho (C_{\omega })<1\quad \Longleftrightarrow \quad 0<\omega <{\frac {2}{\lambda _{\text{max}}(D^{-1}A)}}\,,

dónde $\lambda _{\text{max}}$ es el valor propio máximo.

El radio espectral se puede minimizar para una elección particular de $\omega =\omega _{\text{opt}}$ como sigue $\min _{\omega }\rho (C_{\omega })=\rho (C_{\omega _{\text{opt}}})=1-{\frac {2}{\kappa (D^{-1}A)+1}}\quad {\text{for}}\quad \omega _{\text{opt}}:={\frac {2}{\lambda _{\text{min}}(D^{-1}A)+\lambda _{\text{max}}(D^{-1}A)}}\,,$ dónde $\kappa$ es el número de condición de la matriz .