Método iterativo utilizado para resolver un sistema lineal de ecuaciones.
En álgebra lineal numérica , el método de Gauss-Seidel , también conocido como método de Liebmann o método de desplazamiento sucesivo , es un método iterativo utilizado para resolver un sistema de ecuaciones lineales . Lleva el nombre de los matemáticos alemanes Carl Friedrich Gauss y Philipp Ludwig von Seidel , y es similar al método de Jacobi . Aunque se puede aplicar a cualquier matriz con elementos distintos de cero en las diagonales, la convergencia solo está garantizada si la matriz es estrictamente diagonalmente dominante , [1] o simétrica y definida positiva . Sólo se mencionó en una carta privada de Gauss a su alumno Gerling en 1823. [2] Seidel no entregó una publicación antes de 1874. [3]
Descripción
Sea un sistema cuadrado de n ecuaciones lineales, donde:![{\textstyle A\mathbf {x} =\mathbf {b} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\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{bmatriz}}.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Cuando y se conocen y se desconocen, podemos utilizar el método de Gauss-Seidel para aproximar . El vector denota nuestra estimación inicial de (a menudo, de ). Lo denotamos como la k -ésima aproximación o iteración de , y es la siguiente (o k +1) iteración de .![{\displaystyle A}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {b} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} ^{(0)}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} _{i}^{(0)}=0}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle i=1,2,...,n}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} ^{(k)}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} ^{(k+1)}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Fórmula basada en matrices
La solución se obtiene iterativamente mediante
![{\displaystyle L_{*}\mathbf {x} ^{(k+1)}=\mathbf {b} -U\mathbf {x} ^{(k)},}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
triangular inferiortriangular estrictamente superior[4]![{\displaystyle A}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle L_{*}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle U}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle A=L_{*}+U}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle A}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle L_{*}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle U}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle 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 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 U}.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Por qué funciona la fórmula matricial
El sistema de ecuaciones lineales se puede reescribir como:
![{\displaystyle {\begin{alignedat}{1}A\mathbf {x} &=\mathbf {b} \\(L_{*}+U)\mathbf {x} &=\mathbf {b} \\L_ {*}\mathbf {x} +U\mathbf {x} &=\mathbf {b} \\L_{*}\mathbf {x} &=\mathbf {b} -U\mathbf {x} \end{ alineado}}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
El método Gauss-Seidel ahora resuelve el lado izquierdo de esta expresión para , utilizando el valor anterior para en el lado derecho. Analíticamente, esto puede escribirse como:![{\displaystyle \mathbf {x} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} ^{(k+1)}=L_{*}^{-1}\left(\mathbf {b} -U\mathbf {x} ^{(k)}\right) .}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Fórmula basada en elementos
Sin embargo, aprovechando la forma triangular de , los elementos de se pueden calcular secuencialmente para cada fila mediante sustitución directa : [5]![{\displaystyle L_{*}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} ^{(k+1)}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle i}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\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,\puntos,n.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Observe que la fórmula utiliza dos sumas por iteración que se pueden expresar como una suma que utiliza la iteración calculada más recientemente de . 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 .![{\displaystyle \sum _ {j\neq i}a_ {ij}x_ {j}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle x_{j}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Discusión
La fórmula de elementos del método de Gauss-Seidel es similar a la del método de Jacobi .
El cálculo de utiliza los elementos que ya se han calculado, y solo los elementos que no se han calculado en la ( k +1) -ésima iteración. Esto significa que, a diferencia del método de Jacobi, sólo se requiere un vector de almacenamiento ya que los elementos se pueden sobrescribir a medida que se calculan, lo que puede resultar ventajoso para problemas muy grandes.![{\displaystyle \mathbf {x} ^{(k+1)}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} ^{(k+1)}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} ^{(k)}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Sin embargo, a diferencia del método de Jacobi, los cálculos para cada elemento son generalmente 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 una relajación excesiva sucesiva con .![{\displaystyle \omega =1}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Convergencia
Las propiedades de convergencia del método de Gauss-Seidel dependen de la matriz A. Es decir, se sabe que el procedimiento converge si:
El método de Gauss-Seidel a veces converge incluso si no se cumplen estas condiciones.
Golub y Van Loan dan un teorema para un algoritmo que se divide en dos partes. Supongamos que no es singular. Sea el radio espectral de . Entonces las iteraciones definidas por convergen para cualquier vector inicial si es no singular y . [8]![{\displaystyle A}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle A=MN}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle r=\rho (M^{-1}N)}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle M^{-1}N}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle x^{(k)}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle Mx^{(k+1)}=Nx^{(k)}+b}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle x=A^{-1}b}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle x^{(0)}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle M}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle r<1}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Algoritmo
Dado que los elementos se pueden sobrescribir a medida que se calculan en este algoritmo, solo se necesita un vector de almacenamiento y se omite la indexación de vectores. El algoritmo es el siguiente:
algoritmo El método de Gauss-Seidel es entradas: A , b salida: φ Elija una suposición inicial φ para la solución repita hasta la convergencia para i desde 1 hasta n haga σ ← 0 para j desde 1 hasta n haga si j ≠ i entonces σ ← σ + a ij φ j final si final ( j -bucle) φ i ← ( b i − σ ) / a ii final ( i -bucle) comprobar si se alcanza la convergencia terminar (repetir)
Ejemplos
Un ejemplo para la versión matricial.
Un sistema lineal mostrado como está dado por:![{\displaystyle A\mathbf {x} =\mathbf {b} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle A={\begin{bmatrix}16&3\\7&-11\\\end{bmatrix}}\quad {\text{and}}\quad b={\begin{bmatrix}11\\13\end {bmatriz}}.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
queremos usar la ecuacion
![{\displaystyle \mathbf {x} ^{(k+1)}=L_{*}^{-1}(\mathbf {b} -U\mathbf {x} ^{(k)})}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} ^{(k+1)}=T\mathbf {x} ^{(k)}+C}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle T=-L_{*}^{-1}U\quad {\text{and}}\quad C=L_{*}^{-1}\mathbf {b} .}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Debemos descomponer en la suma de una componente triangular inferior y una componente triangular superior estricta :![{\displaystyle A}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle L_{*}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle U}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle L_{*}={\begin{bmatrix}16&0\\7&-11\\\end{bmatrix}}\quad {\text{and}}\quad U={\begin{bmatrix}0&3\\ 0&0\end{bmatriz}}.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
La inversa de es:![{\displaystyle L_{*}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle L_{*}^{-1}={\begin{bmatrix}16&0\\7&-11\end{bmatrix}}^{-1}={\begin{bmatrix}0,0625&0,0000\\0,0398 &-0.0909\\\end{bmatriz}}.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Ahora podemos encontrar:
![{\displaystyle {\begin{aligned}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]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}}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Ahora tenemos y y podemos usarlos para obtener los vectores de forma iterativa.![{\displaystyle T}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle C}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
En primer lugar, tenemos que elegir : sólo podemos adivinar. Cuanto mejor sea la suposición, más rápido funcionará el algoritmo.![{\displaystyle \mathbf {x} ^{(0)}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Elegimos un punto de partida:
![{\displaystyle x^{(0)}={\begin{bmatrix}1.0\\1.0\end{bmatrix}}.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Entonces podemos calcular:
![{\displaystyle {\begin{aligned}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]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]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]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]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]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] 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}}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Como era de esperar, el algoritmo converge a la solución exacta:
![{\displaystyle \mathbf {x} =A^{-1}\mathbf {b} \approx {\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
De hecho, la matriz A es estrictamente diagonalmente dominante (pero no definida positiva).
Otro ejemplo para la versión matricial.
Otro sistema lineal mostrado está dado por:![{\displaystyle A\mathbf {x} =\mathbf {b} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle A={\begin{bmatrix}2&3\\5&7\\\end{bmatrix}}\quad {\text{and}}\quad b={\begin{bmatrix}11\\13\\\end {bmatriz}}.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
queremos usar la ecuacion
![{\displaystyle \mathbf {x} ^{(k+1)}=L_{*}^{-1}(\mathbf {b} -U\mathbf {x} ^{(k)})}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} ^{(k+1)}=T\mathbf {x} ^{(k)}+C}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle T=-L_{*}^{-1}U\quad {\text{and}}\quad C=L_{*}^{-1}\mathbf {b} .}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Debemos descomponer en la suma de una componente triangular inferior y una componente triangular superior estricta :![{\displaystyle A}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle L_{*}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle U}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle L_{*}={\begin{bmatrix}2&0\\5&7\\\end{bmatrix}}\quad {\text{y}}\quad U={\begin{bmatrix}0&3\\0&0\ \\end{bmatriz}}.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
La inversa de es:![{\displaystyle L_{*}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle L_{*}^{-1}={\begin{bmatrix}2&0\\5&7\\\end{bmatrix}}^{-1}={\begin{bmatrix}0,500&0,000\\- 0,357 y 0,143\\\end{bmatriz}}.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Ahora podemos encontrar:
![{\displaystyle {\begin{aligned}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]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{alineado}}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Ahora tenemos y y podemos usarlos para obtener los vectores de forma iterativa.![{\displaystyle T}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle C}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle \mathbf {x} }](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
En primer lugar, tenemos que elegir : sólo podemos adivinar. Cuanto mejor sea la suposición, más rápido ejecutará el algoritmo.![{\displaystyle \mathbf {x} ^{(0)}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Suponemos:
![{\displaystyle x^{(0)}={\begin{bmatrix}1.1\\2.3\end{bmatrix}}.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Entonces podemos calcular:
![{\displaystyle {\begin{aligned}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]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]x^{( 3)}&=\cdots .\end{alineado}}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Si probamos la convergencia encontraremos que el algoritmo diverge. De hecho, la matriz A no es diagonalmente dominante ni definida positiva. Entonces, la convergencia a la solución exacta.
![{\displaystyle \mathbf {x} =A^{-1}\mathbf {b} ={\begin{bmatrix}-38\\29\end{bmatrix}}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Un ejemplo para la versión de la ecuación.
Supongamos que se dan k ecuaciones donde x n son vectores de estas ecuaciones y el punto de partida x 0 . De la primera ecuación resuelve x 1 en términos de Para las siguientes ecuaciones sustituye los valores anteriores de x s.![{\displaystyle x_{n+1},x_{n+2},\dots,x_{n}.}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Para que quede claro, 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}}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Resolviendo para y da:![{\ Displaystyle x_ {1}, x_ {2}, x_ {3}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\displaystyle x_{4}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
![{\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}}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Supongamos que elegimos (0, 0, 0, 0) como 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}}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
Utilizando las aproximaciones obtenidas, se repite el procedimiento iterativo hasta alcanzar la precisión deseada. Las siguientes son las soluciones aproximadas después de cuatro iteraciones.
La solución exacta del sistema es (1, 2, −1, 1) .
Un ejemplo usando Python y NumPy
El siguiente procedimiento numérico simplemente se itera para producir el vector solución.
importar numpy como npITERACIÓN_LIMIT = 1000# inicializamos la matriz A = np . matriz ( [ [ 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 . matriz ([ 6.0 , 25.0 , - 11.0 , 15.0 ]) print ( "Sistema de ecuaciones:" ) para i en el rango ( A . forma [ 0 ]): fila = [ f " { A [ i , j ] : 3g } *x { j + 1 } " para j en el rango ( A. forma [ 1 ])] print ( "[ {0} ] = [ {1:3g} ] " . formato ( " + " . join ( fila ), b [ i ] )) x = np . zeros_like ( b , np . float_ ) para it_count en el rango ( 1 , ITERATION_LIMIT ): x_new = np . zeros_like ( x , dtype = np . float_ ) print ( f "Iteración { it_count } : { x } " ) para i en el rango ( A . forma [ 0 ]): s1 = np . punto ( A [ i ,: i ], x_new [: i ] ) s2 = np . punto ( A [ i , i + 1 :], x [ i + 1 :]) x_new [ i ] = ( b [ i ] - s1 - s2 ) / A [ i , i ] si np . allclose ( x , x_new , rtol = 1e-8 ): romper x = x_new print ( f "Solución: { x } " ) error = np . punto ( A , x ) - b imprimir ( f "Error: { error } " )
Produce la salida:
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.006585 04 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.0000 0002 -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]
Programa para resolver arbitrario no. de ecuaciones usando Matlab
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{matriz}}}](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)
función x = gauss_seidel ( A, b, x, iteros ) para i = 1 : iteros 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
Ver también
Notas
- ^ Sauer, Timoteo (2006). Análisis numérico (2ª ed.). Pearson Education, Inc. pág. 109.ISBN 978-0-321-78367-7.
- ^
Gauss 1903, pag. 279; enlace directo.
- ^ Seidel, Ludwig (1874). "Über ein Verfahren, die Gleichungen, auf welche die Methode der kleinsten Quadrate führt, sowie lineäre Gleichungen überhaupt, durch sucesive Annäherung aufzulösen" [Sobre un proceso para resolver por aproximaciones sucesivas las ecuaciones a las que conduce el método de mínimos cuadrados y el método lineal ecuaciones en general]. Abhandlungen der Mathematisch-Physikalischen Klasse der Königlich Bayerischen Akademie der Wissenschaften (en alemán). 11 (3): 81–108.
- ^ Préstamo Golub y Van 1996, pág. 511.
- ^ Golub & Van Loan 1996, ecuación (10.1.3)
- ^ Golub & Van Loan 1996, Thm 10.1.2.
- ^ Bagnara, Roberto (marzo de 1995). "Una prueba unificada de la convergencia de los métodos de Jacobi y Gauss-Seidel". Revisión SIAM . 37 (1): 93–97. CiteSeerX 10.1.1.26.5207 . doi :10.1137/1037008. JSTOR 2132758.
- ^ Golub & Van Loan 1996, Thm 10.1.2
Referencias
- Gauss, Carl Friedrich (1903), Werke (en alemán), vol. 9, Gotinga: Köninglichen Gesellschaft der Wissenschaften.
- Golub, Gene H .; Van Loan, Charles F. (1996), Computaciones matriciales (3.ª ed.), Baltimore: Johns Hopkins, ISBN 978-0-8018-5414-9.
- Negro, Noel y Moore, Shirley. "Método Gauss-Seidel". MundoMatemático .
Este artículo incorpora texto del artículo Gauss-Seidel_method en CFD-Wiki que se encuentra bajo la licencia GFDL .
enlaces externos
- "Método Seidel", Enciclopedia de Matemáticas , EMS Press , 2001 [1994]
- Gauss-Seidel de www.math-linux.com
- Gauss-Seidel del Instituto de Métodos Numéricos Holísticos
- Iteración de Gauss Siedel de www.geocities.com
- El método Gauss-Seidel
- bickson
- código matlab
- ejemplo de código C