Método Jacobi

Método iterativo utilizado para resolver un sistema lineal de ecuaciones.

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 de la diagonal y se introduce 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 de matrices . El método recibe su nombre de Carl Gustav Jacob Jacobi .

Descripción

Sea un sistema cuadrado de n ecuaciones lineales, donde: ${\displaystyle A\mathbf {x} =\mathbf {b} }$ $${\displaystyle 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 se conocen y , y se desconoce, podemos utilizar el método de Jacobi para aproximar . El vector denota nuestra estimación inicial para (a menudo para ). Denotamos como la k- ésima aproximación o iteración de , y es la siguiente (o k +1) iteración de . ${\displaystyle A}$ ${\displaystyle \mathbf {b} }$ ${\displaystyle \mathbf {x} }$ ${\displaystyle \mathbf {x} }$ ${\displaystyle \mathbf {x} ^{(0)}}$ ${\displaystyle \mathbf {x} }$ ${\displaystyle \mathbf {x} _{i}^{(0)}=0}$ ${\displaystyle i=1,2,...,n}$ ${\displaystyle \mathbf {x} ^{(k)}}$ ${\displaystyle \mathbf {x} }$ ${\displaystyle \mathbf {x} ^{(k+1)}}$ ${\displaystyle \mathbf {x} }$

Fórmula basada en matrices

Luego 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 entonces iterativamente mediante $${\displaystyle 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}}.}$$

${\displaystyle \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 es la siguiente: El cálculo de requiere cada elemento excepto él mismo. A diferencia del método de 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 dos vectores de tamaño n . ${\displaystyle i}$ $${\displaystyle 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.}$$ ${\displaystyle x_{i}^{(k+1)}}$ ${\displaystyle \mathbf {x} ^{(k)}}$ ${\displaystyle x_{i}^{(k)}}$ ${\displaystyle x_{i}^{(k+1)}}$

Algoritmo

Entrada:  estimación inicial x (0) a la solución , matriz A (diagonal dominante), vector b 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 anteriork = 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:

${\displaystyle \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 otros términos:

${\displaystyle \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.

Nótese que el método de Jacobi no converge para cada matriz definida positiva simétrica . Por ejemplo, $${\displaystyle 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 con estimación inicial está dado por ${\displaystyle Ax=b}$ ${\displaystyle x^{(0)}}$

${\displaystyle 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 , descrita anteriormente, para estimar . Primero, reescribimos la ecuación en una forma más conveniente , donde y . A partir de los valores conocidos , determinamos como Además, se encuentra como Con y calculado, estimamos como : La siguiente iteración produce Este proceso se repite hasta la convergencia (es decir, hasta que sea pequeño). La solución después de 25 iteraciones es ${\displaystyle x^{(k+1)}=D^{-1}(b-(L+U)x^{(k)})}$ ${\displaystyle x}$ ${\displaystyle D^{-1}(b-(L+U)x^{(k)})=Tx^{(k)}+C}$ ${\displaystyle T=-D^{-1}(L+U)}$ ${\displaystyle C=D^{-1}b}$ $${\displaystyle 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}}.}$$ ${\displaystyle T=-D^{-1}(L+U)}$ $${\displaystyle 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}}.}$$ ${\displaystyle C}$ $${\displaystyle 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}}.}$$ ${\displaystyle T}$ ${\displaystyle C}$ ${\displaystyle x}$ ${\displaystyle x^{(1)}=Tx^{(0)}+C}$ $${\displaystyle 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}}.}$$ $${\displaystyle 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}}.}$$ ${\displaystyle \|Ax^{(n)}-b\|}$

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

Pregunta de ejemplo 2

Supongamos que nos dan el siguiente sistema lineal:

${\displaystyle {\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 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. $${\displaystyle {\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}}}$$

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

Ejemplo de Python

importar  numpy  como  npLÍMITE_ITERACIÓN  =  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  =  [ f " { A [ i , j ] } *x { j + 1 } " para j en el rango ( A . forma [ 1 ])]        imprimir ( f ' { " + " . join ( fila ) } = { b [ i ] } ' )imprimir ()x  =  np . ceros_similares ( b )para  it_count  en  el rango ( ITERATION_LIMIT ): si  it_count  !=  0 : imprimir ( f "Iteración { it_count } : { x } " ) x_new  =  np .zeros_like ( 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_nuevo [ i ]  ==  x_nuevo [ i - 1 ]: romper si  np .allclose ( x , x_new , atol = 1e-10 , rtol = 0. ) :    romper x  =  x_nuevoimprimir ( "Solución: " )imprimir ( x )error  =  np . punto ( A ,  x )  -  bimprimir ( "Error:" )imprimir ( error )

Método de Jacobi ponderado

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

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

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

${\displaystyle \mathbf {x} ^{(k+1)}=\omega D^{-1}\mathbf {b} +\left(I-\omega D^{-1}A\right)\mathbf {x} ^{(k)}}$.

Convergencia en el caso simétrico positivo definido

En el caso de que la matriz del sistema sea de tipo simétrico positivo definido se puede demostrar convergencia. ${\displaystyle A}$

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

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

donde es el valor propio máximo. ${\displaystyle \lambda _{\text{max}}}$

El radio espectral se puede minimizar para una elección particular de la siguiente manera, donde es el número de condición de la matriz . ${\displaystyle \omega =\omega _{\text{opt}}}$ $${\displaystyle \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)}}\,,}$$ ${\displaystyle \kappa }$

Véase también

• Método de Gauss-Seidel
• Relajación excesiva sucesiva
• Método iterativo § Sistemas lineales
• Propagación de creencias gaussianas
• División de matriz

Referencias

1. Saad, Yousef (2003). Métodos iterativos para sistemas lineales dispersos (2.ª ed.). SIAM . p. 414. ISBN 0898715342.

Enlaces externos

• Este artículo incorpora texto del artículo Jacobi_method en CFD-Wiki que está bajo la licencia GFDL .

• Black, Noel; Moore, Shirley y Weisstein, Eric W. "El método de Jacobi". MathWorld .
• Método Jacobi de www.math-linux.com