Método de Gauss-Seidel

La solución se obtiene iterativamente mediante $${\displaystyle \mathbf {L} \mathbf {x} ^{(k+1)}=\mathbf {b} -\mathbf {U} \mathbf {x} ^{(k)},}$$ donde la matriz${\displaystyle \mathbf {A} }$se descompone en un componente triangular inferior${\displaystyle \mathbf {L} }$y un componente triangular superior estricto${\displaystyle \mathbf {U} }$de tal manera que${\displaystyle \mathbf {A} =\mathbf {L} +\mathbf {U} }$. ^{[ 4 ]} Más específicamente, la descomposición de${\displaystyle \mathbf {A} }$en${\displaystyle \mathbf {L} }$y${\displaystyle \mathbf {U} }$está dado por:

$${\displaystyle \mathbf {A} =\underbrace {\begin{bmatrix}a_{11}&0&\cdots &0\\a_{21}&a_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &a_{nn}\end{bmatrix}} _{\textstyle \mathbf {L} }+\underbrace {\begin{bmatrix}0&a_{12}&\cdots &a_{1n}\\0&0&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &0\end{bmatrix}} _{\textstyle \mathbf {U} }.}$$

Sin embargo, aprovechando la forma triangular de${\displaystyle \mathbf {L} }$, los elementos de${\displaystyle \mathbf {x} ^{(k+1)}}$se puede calcular secuencialmente para cada fila${\displaystyle i}$utilizando sustitución hacia adelante : ^{[ 5 ]}$${\displaystyle x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum _{j=i+1}^{n}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\dots ,n.}$$

Observe que la fórmula utiliza dos sumatorias por iteración, que pueden expresarse como una sola sumatoria.${\displaystyle \sum _{j\neq i}a_{ij}x_{j}}$que utiliza la iteración calculada más recientemente de${\displaystyle x_{j}}$. El procedimiento generalmente continúa hasta que los cambios realizados por una iteración estén por debajo de cierta tolerancia, como un residuo suficientemente pequeño .

La fórmula elemento a elemento para el método de Gauss-Seidel está relacionada con la del método de Jacobi (iterativo) , con una diferencia importante:

En Gauss-Seidel, el cálculo de${\displaystyle \mathbf {x} ^{(k+1)}}$utiliza los elementos de${\displaystyle \mathbf {x} ^{(k+1)}}$que ya han sido calculados, y solo los elementos de${\displaystyle \mathbf {x} ^{(k)}}$que no se han calculado en el${\displaystyle (k+1)}$-ésima iteración. Esto significa que, a diferencia del método de Jacobi, solo se requiere un vector de almacenamiento, ya que los elementos se pueden sobrescribir a medida que se calculan, lo que puede ser ventajoso para problemas muy grandes.

Sin embargo, a diferencia del método de Jacobi, los cálculos para cada elemento suelen ser mucho más difíciles de implementar en paralelo , ya que pueden tener una ruta crítica muy larga , y por lo tanto, son más factibles para matrices dispersas . Además, los valores en cada iteración dependen del orden de las ecuaciones originales.

Gauss-Seidel es lo mismo que la sobre-relajación sucesiva con${\displaystyle \omega =1}$.

Un sistema lineal mostrado como${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$está dado por: $${\displaystyle \mathbf {A} ={\begin{bmatrix}16&3\\7&-11\\\end{bmatrix}}\quad {\text{and}}\quad \mathbf {b} ={\begin{bmatrix}11\\13\end{bmatrix}}.}$$

Utilice la ecuación $${\displaystyle \mathbf {x} ^{(k+1)}=\mathbf {L} ^{-1}(\mathbf {b} -\mathbf {U} \mathbf {x} ^{(k)})}$$ en la forma $${\displaystyle \mathbf {x} ^{(k+1)}=\mathbf {T} \mathbf {x} ^{(k)}+\mathbf {c} }$$ dónde:

${\displaystyle \mathbf {T} =-\mathbf {L} ^{-1}\mathbf {U} \quad {\text{and}}\quad \mathbf {c} =\mathbf {L} ^{-1}\mathbf {b} .}$

Descomponer${\displaystyle \mathbf {A} }$en la suma de un componente triangular inferior${\displaystyle \mathbf {L} }$y un componente triangular superior estricto${\displaystyle U}$: $${\displaystyle \mathbf {L} ={\begin{bmatrix}16&0\\7&-11\\\end{bmatrix}}\quad {\text{and}}\quad \mathbf {U} ={\begin{bmatrix}0&3\\0&0\end{bmatrix}}.}$$

Lo contrario de${\displaystyle \mathbf {L} }$es: $${\displaystyle \mathbf {L} ^{-1}={\begin{bmatrix}16&0\\7&-11\end{bmatrix}}^{-1}={\begin{bmatrix}0.0625&0.0000\\0.0398&-0.0909\\\end{bmatrix}}.}$$

Ahora encuentra: $${\displaystyle {\begin{aligned}\mathbf {T} &=-{\begin{bmatrix}0.0625&0.0000\\0.0398&-0.0909\end{bmatrix}}{\begin{bmatrix}0&3\\0&0\end{bmatrix}}={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1194\end{bmatrix}},\\[1ex]\mathbf {c} &={\begin{bmatrix}0.0625&0.0000\\0.0398&-0.0909\end{bmatrix}}{\begin{bmatrix}11\\13\end{bmatrix}}={\begin{bmatrix}0.6875\\-0.7439\end{bmatrix}}.\end{aligned}}}$$

Con${\displaystyle \mathbf {T} }$y${\displaystyle \mathbf {c} }$los vectores${\displaystyle \mathbf {x} }$se puede obtener de forma iterativa.

En primer lugar, elige${\displaystyle \mathbf {x} ^{(0)}}$, Por ejemplo$${\displaystyle \mathbf {x} ^{(0)}={\begin{bmatrix}1.0\\1.0\end{bmatrix}}.}$$Cuanto más se acerque la estimación a la solución final, menos iteraciones necesitará el algoritmo.

Luego calcula: $${\displaystyle {\begin{aligned}\mathbf {x} ^{(1)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}1.0\\1.0\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.5000\\-0.8636\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(2)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.5000\\-0.8636\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8494\\-0.6413\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(3)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8494\\-0.6413\\\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8077\\-0.6678\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(4)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8077\\-0.6678\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8127\\-0.6646\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(5)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8127\\-0.6646\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8121\\-0.6650\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(6)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8121\\-0.6650\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(7)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}.\end{aligned}}}$$

Como era de esperar, el algoritmo converge a la solución: $${\displaystyle \mathbf {x} =\mathbf {A} ^{-1}\mathbf {b} \approx {\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}}$$.

De hecho, la matriz A es estrictamente diagonalmente dominante, pero no definida positiva.

Otro sistema lineal mostrado como${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$está dado por:

$${\displaystyle \mathbf {A} ={\begin{bmatrix}2&3\\5&7\\\end{bmatrix}}\quad {\text{and}}\quad \mathbf {b} ={\begin{bmatrix}11\\13\\\end{bmatrix}}.}$$

Utilice la ecuación $${\displaystyle \mathbf {x} ^{(k+1)}=\mathbf {L} ^{-1}(\mathbf {b} -\mathbf {U} \mathbf {x} ^{(k)})}$$ en la forma $${\displaystyle \mathbf {x} ^{(k+1)}=\mathbf {T} \mathbf {x} ^{(k)}+\mathbf {c} }$$ dónde:

${\displaystyle \mathbf {T} =-\mathbf {L} ^{-1}\mathbf {U} \quad {\text{and}}\quad \mathbf {c} =\mathbf {L} ^{-1}\mathbf {b} .}$

Descomponer${\displaystyle \mathbf {A} }$en la suma de un componente triangular inferior${\displaystyle \mathbf {L} }$y un componente triangular superior estricto${\displaystyle \mathbf {U} }$: $${\displaystyle \mathbf {L} ={\begin{bmatrix}2&0\\5&7\\\end{bmatrix}}\quad {\text{and}}\quad \mathbf {U} ={\begin{bmatrix}0&3\\0&0\\\end{bmatrix}}.}$$

Lo contrario de${\displaystyle \mathbf {L} }$es: $${\displaystyle \mathbf {L} ^{-1}={\begin{bmatrix}2&0\\5&7\\\end{bmatrix}}^{-1}={\begin{bmatrix}0.500&0.000\\-0.357&0.143\\\end{bmatrix}}.}$$

Ahora encuentra: $${\displaystyle {\begin{aligned}\mathbf {T} &=-{\begin{bmatrix}0.500&0.000\\-0.357&0.143\\\end{bmatrix}}{\begin{bmatrix}0&3\\0&0\\\end{bmatrix}}={\begin{bmatrix}0.000&-1.500\\0.000&1.071\\\end{bmatrix}},\\[1ex]\mathbf {c} &={\begin{bmatrix}0.500&0.000\\-0.357&0.143\\\end{bmatrix}}{\begin{bmatrix}11\\13\\\end{bmatrix}}={\begin{bmatrix}5.500\\-2.071\\\end{bmatrix}}.\end{aligned}}}$$

Con${\displaystyle \mathbf {T} }$y${\displaystyle \mathbf {c} }$los vectores${\displaystyle \mathbf {x} }$se puede obtener de forma iterativa.

En primer lugar, tenemos que elegir${\displaystyle \mathbf {x} ^{(0)}}$, Por ejemplo$${\displaystyle \mathbf {x} ^{(0)}={\begin{bmatrix}1.1\\2.3\end{bmatrix}}}$$

Luego calcula: $${\displaystyle {\begin{aligned}\mathbf {x} ^{(1)}&={\begin{bmatrix}0&-1.500\\0&1.071\\\end{bmatrix}}{\begin{bmatrix}1.1\\2.3\\\end{bmatrix}}+{\begin{bmatrix}5.500\\-2.071\\\end{bmatrix}}={\begin{bmatrix}2.050\\0.393\\\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(2)}&={\begin{bmatrix}0&-1.500\\0&1.071\\\end{bmatrix}}{\begin{bmatrix}2.050\\0.393\\\end{bmatrix}}+{\begin{bmatrix}5.500\\-2.071\\\end{bmatrix}}={\begin{bmatrix}4.911\\-1.651\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(3)}&=\cdots .\end{aligned}}}$$

En una prueba de convergencia encontramos que el algoritmo diverge. De hecho, la matriz${\displaystyle \mathbf {A} }$no es ni diagonalmente dominante ni definida positiva. Entonces, converge a la solución exacta. $${\displaystyle \mathbf {x} =\mathbf {A} ^{-1}\mathbf {b} ={\begin{bmatrix}-38\\29\end{bmatrix}}}$$ No está garantizado y, en este caso, no ocurrirá.

Supongamos que se da${\displaystyle n}$ecuaciones y un punto de partida${\displaystyle \mathbf {x} _{0}}$En cualquier paso de una iteración de Gauss-Seidel, resuelva la primera ecuación para${\displaystyle x_{1}}$en términos de${\displaystyle x_{2},\dots ,x_{n}}$; luego resuelve la segunda ecuación para${\displaystyle x_{2}}$en términos de${\displaystyle x_{1}}$recién encontrado y el resto${\displaystyle x_{3},\dots ,x_{n}}$; y continuar a${\displaystyle x_{n}}$Luego, repita las iteraciones hasta que se logre la convergencia, o deténgase si la divergencia en las soluciones comienza a superar un nivel predefinido.

Consideremos un ejemplo: $${\displaystyle {\begin{array}{rrrrl}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{array}}}$$

Resolver para${\displaystyle x_{1},x_{2},x_{3}}$y${\displaystyle x_{4}}$da: $${\displaystyle {\begin{aligned}x_{1}&=x_{2}/10-x_{3}/5+3/5,\\x_{2}&=x_{1}/11+x_{3}/11-3x_{4}/11+25/11,\\x_{3}&=-x_{1}/5+x_{2}/10+x_{4}/10-11/10,\\x_{4}&=-3x_{2}/8+x_{3}/8+15/8.\end{aligned}}}$$

Supongamos que (0, 0, 0, 0) es la aproximación inicial, entonces la primera solución aproximada viene dada por: $${\displaystyle {\begin{aligned}x_{1}&=3/5=0.6,\\x_{2}&=(3/5)/11+25/11=3/55+25/11=2.3272,\\x_{3}&=-(3/5)/5+(2.3272)/10-11/10=-3/25+0.23272-1.1=-0.9873,\\x_{4}&=-3(2.3272)/8+(-0.9873)/8+15/8=0.8789.\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 cuatro iteraciones.

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

El siguiente procedimiento iterativo produce el vector solución de un sistema de ecuaciones lineales:

import numpy as npLÍMITE_DE_ITERACIÓN = 1000# Inicializar la matriz A = np.array ( [ [ 10.0 , - 1.0 , 2.0 , 0.0 ] , [ - 1.0 , 11.0 , - 1.0 , 3.0 ], [ 2.0 , - 1.0 , 10.0 , - 1.0 ], [ 0.0 , 3.0 , - 1.0 , 8.0 ], ] ) # Inicializar el vector RHS b = np.array ( [ 6.0 , 25.0 , - 11.0 , 15.0 ] )print ( "Sistema de ecuaciones:" ) for i in range ( A . shape [ 0 ]): row = " + " . join ( f " { A [ i , j ] : 3g } *x { j + 1 } " for j in range ( A . shape [ 1 ])) print ( f "[ { row } ] = [ { b [ i ] : 3g } ]" )x = np.zeros_like ( b ) para it_count en range ( 1 , ITERATION_LIMIT ) : x_new = np.zeros_like ( x ) print ( f " Iteración { it_count } : { x } " ) para i en range ( A.shape [ 0 ] ) : s1 = A [ i , : i ] @ x_new [: i ] s2 = A [ i , i + 1 :] @ x [ i + 1 : ] x_new [ i ] = ( b [ i ] - s1 - s2 ) / A [ i , i ] if np.allclose ( x , x_new , rtol = 1e - 8 ) : break x = x_newprint ( f "Solución: { x } " ) error = A @ x - b print ( f "Error: { error } " )

Produce el siguiente resultado:

Sistema de ecuaciones: [ 10*x1 + -1*x2 + 2*x3 + 0*x4] = [ 6] [ -1*x1 + 11*x2 + -1*x3 + 3*x4] = [ 25] [ 2*x1 + -1*x2 + 10*x3 + -1*x4] = [-11] [ 0*x1 + 3*x2 + -1*x3 + 8*x4] = [ 15] Iteración 1: [ 0. 0. 0. 0.] Iteración 2: [ 0.6 2.32727273 -0.98727273 0.87886364] Iteración 3: [ 1.03018182 2.03693802 -1.0144562 0.98434122] Iteración 4: [ 1.00658504 2.00355502 -1.00252738 0.99835095] Iteración 5: [ 1.00086098 2.00029825 -1.00030728 0.99984975] Iteración 6: [ 1.00009128 2.00002134 -1.00003115 0.9999881 ] Iteración 7: [ 1.00000836 2.00000117 -1.00000275 0.99999922] Iteración 8: [ 1.00000067 [ 2.00000002 -1.00000021 0.99999996] Iteración 9: [ 1.00000004 1.99999999 -1.00000001 1. ] Iteración 10: [ 1. 2. -1. 1.] Solución: [ 1. 2. -1. 1.] Error: [ 2.06480930e-08 -1.25551054e-08 3.61417563e-11 0.00000000e+00]

El siguiente código utiliza la fórmula $${\displaystyle x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j<i}a_{ij}x_{j}^{(k+1)}-\sum _{j>i}a_{ij}x_{j}^{(k)}\right),\quad {\begin{array}{l}i=1,2,\ldots ,n\\k=0,1,2,\ldots \end{array}}}$$

función x = gauss_seidel ( A, b, x, iters ) para i = 1 : iters para j = 1 : tamaño ( A , 1 ) x ( j ) = ( b ( j ) - suma ( A ( j ,:) '.* x ) + A ( j , j ) * x ( j )) / A ( j , j ); fin fin fin