Sobrerelajación sucesiva

Método de resolución de un sistema lineal de ecuaciones.

En álgebra lineal numérica , el método de sobrerelajación sucesiva ( SOR ) es una variante del método de Gauss-Seidel para resolver un sistema lineal de ecuaciones , lo que resulta en una convergencia más rápida. Se puede utilizar un método similar para cualquier proceso iterativo que converja lentamente .

Fue ideado simultáneamente por David M. Young Jr. y por Stanley P. Frankel en 1950 con el propósito de resolver automáticamente sistemas lineales en computadoras digitales. Los métodos de excesiva relajación se habían utilizado antes del trabajo de Young y Frankel. Un ejemplo es el método de Lewis Fry Richardson y los métodos desarrollados por RV Southwell . Sin embargo, estos métodos fueron diseñados para el cálculo por calculadoras humanas , lo que requería cierta experiencia para garantizar la convergencia a la solución, lo que los hacía inaplicables para la programación en computadoras digitales. Estos aspectos se discuten en la tesis de David M. Young Jr. ^[1]

Formulación

Dado un sistema cuadrado de n ecuaciones lineales con x desconocida :

${\displaystyle A\mathbf {x} =\mathbf {b} }$

dónde:

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

Entonces A se puede descomponer en una componente diagonal D y en las componentes triangulares estrictamente inferior y superior L y U :

${\displaystyle A=D+L+U,}$

dónde

${\displaystyle D={\begin{bmatrix}a_{11}&0&\cdots &0\\0&a_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &a_{nn}\end{bmatrix}},\quad L={\begin{bmatrix}0&0&\cdots &0\\a_{21}&0&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &0\end{bmatrix}},\quad U={\begin{bmatrix}0&a_{12}&\cdots &a_{1n}\\0&0&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &0\end{bmatrix}}.}$

El sistema de ecuaciones lineales se puede reescribir como:

${\displaystyle (D+\omega L)\mathbf {x} =\omega \mathbf {b} -[\omega U+(\omega -1)D]\mathbf {x} }$

para una constante ω > 1, llamado factor de relajación .

El método de sobrerelajación sucesiva es una técnica iterativa que resuelve el lado izquierdo de esta expresión para x , utilizando el valor anterior para x en el lado derecho. Analíticamente, esto puede escribirse como:

${\displaystyle \mathbf {x} ^{(k+1)}=(D+\omega L)^{-1}{\big (}\omega \mathbf {b} -[\omega U+(\omega -1)D]\mathbf {x} ^{(k)}{\big )}=L_{\omega }\mathbf {x} ^{(k)}+\mathbf {c} ,}$

donde es la k -ésima aproximación o iteración de y es la siguiente o k + 1 iteración de . Sin embargo, aprovechando la forma triangular de ( D + ωL ), los elementos de x ^{( k +1)} se pueden calcular secuencialmente mediante sustitución directa : ${\displaystyle \mathbf {x} ^{(k)}}$ ${\displaystyle \mathbf {x} }$ ${\displaystyle \mathbf {x} ^{(k+1)}}$ ${\displaystyle \mathbf {x} }$

${\displaystyle x_{i}^{(k+1)}=(1-\omega )x_{i}^{(k)}+{\frac {\omega }{a_{ii}}}\left(b_{i}-\sum _{j<i}a_{ij}x_{j}^{(k+1)}-\sum _{j>i}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\ldots ,n.}$

Esto se puede escribir nuevamente analíticamente en forma de matriz-vector sin la necesidad de invertir la matriz : ^[2] ${\displaystyle (D+\omega L)}$

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

Convergencia

Radio espectral de la matriz de iteración para el método SOR . El gráfico muestra la dependencia del radio espectral de la matriz de iteración de Jacobi . ${\displaystyle \rho (C_{\omega })}$ ${\displaystyle C_{\omega }}$ ${\displaystyle \mu :=\rho (C_{\text{Jac}})}$

La elección del factor de relajación ω no es necesariamente fácil y depende de las propiedades de la matriz de coeficientes. En 1947, Ostrowski demostró que si es simétrico y definido positivo entonces para . Por lo tanto, se produce la convergencia del proceso de iteración, pero generalmente estamos interesados ​​en una convergencia más rápida en lugar de simplemente convergencia. ${\displaystyle A}$ ${\displaystyle \rho (L_{\omega })<1}$ ${\displaystyle 0<\omega <2}$

Tasa de convergencia

La tasa de convergencia del método SOR se puede derivar analíticamente. Es necesario asumir lo siguiente ^{[3] [4]}

• el parámetro de relajación es apropiado: ${\displaystyle \omega \in (0,2)}$
• La matriz de iteración de Jacobi sólo tiene valores propios reales ${\displaystyle C_{\text{Jac}}:=I-D^{-1}A}$
• El método de Jacobi es convergente: ${\displaystyle \mu :=\rho (C_{\text{Jac}})<1}$
• la descomposición matricial satisface la propiedad de que para cualquiera y . ${\displaystyle A=D+L+U}$ ${\displaystyle \operatorname {det} (\lambda D+zL+{\tfrac {1}{z}}U)=\operatorname {det} (\lambda D+L+U)}$ ${\displaystyle z\in \mathbb {C} \setminus \{0\}}$ ${\displaystyle \lambda \in \mathbb {C} }$

Entonces la tasa de convergencia se puede expresar como

${\displaystyle \rho (C_{\omega })={\begin{cases}{\frac {1}{4}}\left(\omega \mu +{\sqrt {\omega ^{2}\mu ^{2}-4(\omega -1)}}\right)^{2}\,,&0<\omega \leq \omega _{\text{opt}}\\\omega -1\,,&\omega _{\text{opt}}<\omega <2\end{cases}}}$

donde el parámetro de relajación óptimo está dado por

${\displaystyle \omega _{\text{opt}}:=1+\left({\frac {\mu }{1+{\sqrt {1-\mu ^{2}}}}}\right)^{2}=1+{\frac {\mu ^{2}}{4}}+O(\mu ^{3})\,.}$

En particular, para ( Gauss-Seidel ) se cumple que . Para el óptimo obtenemos , lo que muestra que SOR es aproximadamente cuatro veces más eficiente que Gauss-Seidel. ${\displaystyle \omega =1}$ ${\displaystyle \rho (C_{\omega })=\mu ^{2}=\rho (C_{\text{Jac}})^{2}}$ ${\displaystyle \omega }$ ${\displaystyle \rho (C_{\omega })={\frac {1-{\sqrt {1-\mu ^{2}}}}{1+{\sqrt {1-\mu ^{2}}}}}={\frac {\mu ^{2}}{4}}+O(\mu ^{3})}$

El último supuesto se cumple para matrices tridiagonales ya que para diagonales con entradas y . ${\displaystyle Z(\lambda D+L+U)Z^{-1}=\lambda D+zL+{\tfrac {1}{z}}U}$ ${\displaystyle Z}$ ${\displaystyle Z_{ii}=z^{i-1}}$ ${\displaystyle \operatorname {det} (\lambda D+L+U)=\operatorname {det} (Z(\lambda D+L+U)Z^{-1})}$

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:

Entradas: A , b , ω
Salida: φElija una suposición inicial φ para la solución repita hasta la convergencia para  i  desde 1 hasta  n  establezca σ en 0 para  j  desde 1 hasta  n  haga  si  j ≠ i  luego establezca σ en σ + a ij  φ j  end if  end ( j -loop ) establezca φ i en (1 − ω ) φ i + ω ( b i − σ ) / a ii  end ( i -bucle) comprobar si se alcanza la convergenciaterminar (repetir)

Nota

${\displaystyle (1-\omega )\phi _{i}+{\frac {\omega }{a_{ii}}}(b_{i}-\sigma )}$También se puede escribir , ahorrando así una multiplicación en cada iteración del bucle for externo . ${\displaystyle \phi _{i}+\omega \left({\frac {b_{i}-\sigma }{a_{ii}}}-\phi _{i}\right)}$

Ejemplo

Se nos presenta el sistema lineal.

${\displaystyle {\begin{aligned}4x_{1}-x_{2}-6x_{3}+0x_{4}&=2,\\-5x_{1}-4x_{2}+10x_{3}+8x_{4}&=21,\\0x_{1}+9x_{2}+4x_{3}-2x_{4}&=-12,\\1x_{1}+0x_{2}-7x_{3}+5x_{4}&=-6.\end{aligned}}}$

Para resolver las ecuaciones, elegimos un factor de relajación y un vector de estimación inicial . Según el algoritmo de sobrerelajación sucesiva, se obtiene la siguiente tabla, que representa una iteración ejemplar con aproximaciones, que idealmente, pero no necesariamente, encuentra la solución exacta, (3, −2, 2, 1) , en 38 pasos. ${\displaystyle \omega =0.5}$ ${\displaystyle \phi =(0,0,0,0)}$

A continuación se ofrece una implementación simple del algoritmo en Common Lisp.

;; Establezca el formato de punto flotante predeterminado en "flotación larga" para ;; garantizar el funcionamiento correcto en una gama más amplia de números. ( setf *read-default-float-format* 'long-float )  ( defparameter +NÚMERO-MÁXIMO-DE-ITERACIONES+ 100 "El número de iteraciones más allá del cual el algoritmo debe dejar de  funcionar, independientemente de su solución actual. Un número mayor de  iteraciones puede proporcionar un resultado más preciso, pero impone requisitos de rendimiento más altos  ." )   ( declamar ( tipo ( entero 0 * ) +NÚMERO-MÁXIMO-DE-ITERACIONES+ ))     ( defun get-errors ( solución-calculada- solución-exacta ) "Para cada componente del vector SOLUCIÓN-COMPUTADA, recupera su  error con respecto al vector SOLUCIÓN-EXACTA esperado, devolviendo un  vector de valores de error.  ---  Mientras ambos ingresan los vectores deben ser iguales en tamaño, esta condición  no se verifica y el más corto de los dos determina el  número de elementos del vector de salida.  ---  La fórmula establecida es la siguiente:  Let resultVectorSize = min(computedSolution.length, exactitudSolution.length)  Let. VectorResultado = nuevo vector de TamañoVectorResultado  Para i de 0 a (TamañoVectorResultado - 1)  VectorResultado[i] = Soluciónexacta[i] - Solucióncalculada[i]  Devuelve VectorResultado" ( declarar ( tipo ( número de vector * ) solución-calculada )) ( declarar ( escriba ( número de vector * ) solución-exacta )) ( mapa ' ( número de vector * ) #' - solución-exacta solución -calculada ))                       ( defun is-convergent ( errores &key ( error-tolerance 0.001 )) "Comprueba si se alcanza la convergencia con respecto al  vector ERRORES que registra la discrepancia entre el  vector solución calculado y el exacto.  ---  La convergencia se cumple si y solo si cada componente de error absoluto  es menor o igual a la TOLERANCIA DE ERRORES, es decir:  Para todo e en ERRORES, se cumple: abs(e) <= errorTolerance." ( declarar ( escribir ( número de vector * ) errores )) ( declarar ( escribir número de error-tolerancia )) ( flet (( el error-es-aceptable ( error ) ( declarar ( escribir número de error )) ( <= ( abs error ) error-tolerancia ))) ( cada #' error-es -errores aceptables )))                              ( defun make-zero-vector ( tamaño ) "Crea y devuelve un vector de TAMAÑO con todos los elementos establecidos en 0." ( declara ( tipo ( entero 0 * ) tamaño )) ( tamaño de creación de matriz : elemento inicial 0.0 : tipo de elemento 'número ))               ( defun relajación excesiva sucesiva ( A b omega & clave ( phi ( make-zero-vector ( longitud b ))) ( verificación de convergencia #' ( lambda ( iteración phi ) ( declarar ( ignorar phi )) ( >= iteración + NÚMERO-MÁXIMO-DE-ITERACIONES+ )))) "Implementa el método de sobre-relajación sucesiva (SOR), aplicado sobre  las ecuaciones lineales definidas por la matriz A y el vector B del lado derecho  , empleando el factor de relajación OMEGA, devolviendo el  vector de solución calculado.  ---  El primer paso del algoritmo, la elección de una estimación inicial de PHI, está  representada por el parámetro de palabra clave opcional PHI, que por defecto  es un vector cero de la misma estructura que B. Si se proporciona, este  vector será modificado destructivamente. En cualquier caso, el vector PHI  constituye el valor del resultado de la función.  ---  La condición de terminación se implementa mediante CONVERGENCE-CHECK,  un predicado opcional  lambda(iteración phi) => booleano generalizado  que devuelve T, lo que significa lo inmediato. terminación, al lograr  la convergencia, o NIL, que indica operación continua, en caso contrario. En  su configuración predeterminada, CONVERGENCE-CHECK simplemente respeta la  ascensión de la iteración al ``+NÚMERO-MÁXIMO-DE-ITERACIONES+'',  ignorando la precisión lograda del vector PHI." ( declarar ( tipo ( número de matriz ( * * ) ) A )) ( declarar ( tipo ( número de vector * ) b )) ( declarar ( número de tipo omega )) ( declarar ( tipo ( número de vector * ) phi )) ( declarar ( tipo ( función (( entero 1 * ) ( vector número * )) * ) verificación de convergencia )) ( let (( n                                                         ( dimensión-matriz A 0 ))) ( declarar ( tipo ( entero 0 * ) n )) ( bucle para iteración de 1 por 1 hacer ( bucle para i desde 0 debajo de n por 1 hacer ( let (( rho 0 )) ( declarar ( escriba el número rho )) ( bucle para j desde 0 debajo de n por 1 hacer ( cuando ( /= j i ) ( let (( a[ij] ( aref A i j )) ( phi[j] ( aref phi j ))) ( incf rho ( * a[ij] phi[j] ))))) ( setf ( aref phi i ) ( + ( * ( - 1 omega ) ( aref phi i )) ( * ( / omega ( aref A i i )) ( - ( aref b i ) rho )))))) ( formato T "~&~d. solución = ~a" iteración phi ) ;; Compruebe si se alcanza la convergencia ( cuando ( funcall convergence-. comprobar iteración phi ) ( retorno )))) ( el ( número de vector * ) phi ))                                                                                                       ;; Invoque la función con los parámetros de ejemplo. ( let (( A ( make-array ( lista 4 4 ) : contenido-inicial ' (( 4 -1 -6 0 ) ( -5 -4 10 8 ) ( 0 9 4 -2 ) ( 1 0 -7 5 ) ))) ( b ( vector 2 21 -12 -6 )) ( omega 0.5 ) ( solución exacta ( vector 3 -2 2 1 ))) ( sobre-relajación sucesiva A b omega :verificación-de-convergencia #' ( lambda ( iteración phi ) ( declarar ( tipo ( entero 0 * ) iteración )) ( declarar ( tipo ( número de vector * ) phi )) ( let (( errores ( get-errors phi solución-exacta ))) ( declarar ( tipo ( vector número * ) errores )) ( formato T "~&~d. errores = ~a" errores de iteración ) ( o ( errores es-convergentes :tolerancia-error 0.0 ) ( >= iteración +NÚMERO-MÁXIMO-DE-ITERACIONES+ )) ))))                                                                                        

Una implementación simple en Python del pseudocódigo proporcionado anteriormente.

importar  numpy  como  np desde  scipy  importar  linalgdef  sor_solver ( A ,  b ,  omega ,  inicial_guess ,  convergence_criteria ): """  Esta es una implementación del pseudocódigo proporcionado en el artículo de Wikipedia.  Argumentos:  A: matriz numpy nxn.  b: vector numpy de n dimensiones.  omega: relajación factor.inicial_guess :  una suposición  de solución inicial para  que el solucionador  comience  . ] residual = nalg . norma ( A @ phi - b ) # Residual inicial mientras residual > criterios_convergencia : para i en el rango ( A. forma [ 0 ]): sigma = 0 para j en el rango ( A. forma [ 1 ] ): si j ! = i : sigma += A [ i , j ] * phi [ j ] phi [ i ] = ( 1 - omega ) * phi [ i ] + ( omega / A [ i , i ]) * ( b [ i ] - sigma ) residual = nalg . norma ( A @ phi - b ) paso += 1 print ( "Paso {} Residual: {:10.6g} " . formato ( paso , residual )) return phi                                                                      # Un caso de ejemplo que refleja el del artículo de Wikipedia residual_convergence  =  1e-8 omega  =  0.5  # Factor de relajaciónA  =  np . matriz ([[ 4 ,  - 1 ,  - 6 ,  0 ],  [ - 5 ,  - 4 ,  10 ,  8 ],  [ 0 ,  9 ,  4 ,  - 2 ],  [ 1 ,  0 ,  - 7 ,  5 ]])b  =  np . matriz ([ 2 ,  21 ,  - 12 ,  - 6 ])conjetura_inicial  =  np . ceros ( 4 )phi  =  sor_solver ( A ,  b ,  omega ,  conjetura_inicial ,  convergencia_residual ) print ( phi )

Sobrerelajación sucesiva simétrica

La versión para matrices simétricas A , en la que

${\displaystyle U=L^{T},\,}$

se conoce como sobrerelajación sucesiva simétrica , o ( SSOR ), en la que

${\displaystyle P=\left({\frac {D}{\omega }}+L\right){\frac {\omega }{2-\omega }}D^{-1}\left({\frac {D}{\omega }}+U\right),}$

y el método iterativo es

${\displaystyle \mathbf {x} ^{k+1}=\mathbf {x} ^{k}-\gamma ^{k}P^{-1}(A\mathbf {x} ^{k}-\mathbf {b} ),\ k\geq 0.}$

Los métodos SOR y SSOR se atribuyen a David M. Young Jr.

Otras aplicaciones del método

Se puede utilizar una técnica similar para cualquier método iterativo. Si la iteración original tuviera la forma

${\displaystyle x_{n+1}=f(x_{n})}$

entonces la versión modificada usaría

${\displaystyle x_{n+1}^{\mathrm {SOR} }=(1-\omega )x_{n}^{\mathrm {SOR} }+\omega f(x_{n}^{\mathrm {SOR} }).}$

Sin embargo, la formulación presentada anteriormente, utilizada para resolver sistemas de ecuaciones lineales, no es un caso especial de esta formulación si se considera que x es el vector completo. Si en su lugar se utiliza esta formulación, la ecuación para calcular el siguiente vector se verá así

${\displaystyle \mathbf {x} ^{(k+1)}=(1-\omega )\mathbf {x} ^{(k)}+\omega L_{*}^{-1}(\mathbf {b} -U\mathbf {x} ^{(k)}),}$

dónde . Los valores de se utilizan para acelerar la convergencia de un proceso de convergencia lenta, mientras que los valores de se utilizan a menudo para ayudar a establecer la convergencia de un proceso iterativo divergente o acelerar la convergencia de un proceso excesivo . ${\displaystyle L_{*}=L+D}$ ${\displaystyle \omega >1}$ ${\displaystyle \omega <1}$

Existen varios métodos que establecen de forma adaptativa el parámetro de relajación en función del comportamiento observado del proceso convergente. Por lo general, ayudan a alcanzar una convergencia superlineal para algunos problemas, pero fracasan en otros. ${\displaystyle \omega }$

Ver también

• método jacobi
• Propagación de la creencia gaussiana
• División de matrices

Notas

1. Young, David M. (1 de mayo de 1950), Métodos iterativos para resolver ecuaciones en diferencias parciales de tipo elíptico (PDF) , tesis doctoral, Universidad de Harvard , consultado el 15 de junio de 2009.
2. Törnig, Willi. Numerische Mathematik für Ingenieure und Physiker (1 ed.). Springer Berlín, Heidelberg. pag. 180.ISBN​ 978-3-642-96508-1. Consultado el 20 de mayo de 2024 .
3. Hackbusch, Wolfgang (2016). "4.6.2". Solución iterativa de grandes sistemas dispersos de ecuaciones | Enlace Springer . Ciencias Matemáticas Aplicadas. vol. 95. doi :10.1007/978-3-319-28483-5. ISBN 978-3-319-28481-1.
4. Greenbaum, Anne (1997). "10.1". Métodos iterativos para resolver sistemas lineales . Fronteras en Matemática Aplicada. vol. 17. doi :10.1137/1.9781611970937. ISBN 978-0-89871-396-1.

Referencias

• Este artículo incorpora texto del artículo Successive_over-relaxation_method_-_SOR en CFD-Wiki que se encuentra bajo la licencia GFDL .

• Abraham Berman, Robert J. Plemmons , Matrices no negativas en las ciencias matemáticas , 1994, SIAM. ISBN 0-89871-321-8 . 
• Negro, Noel y Moore, Shirley. "Método de sobrerelajación sucesiva". MundoMatemático .
• A. Hadjidimos, Sobrerelajación sucesiva (SOR) y métodos relacionados , Journal of Computational and Applied Mathematics 123 (2000), 177–199.
• Yousef Saad , Métodos iterativos para sistemas lineales dispersos , 1.ª edición, PWS, 1996.
• Copia de Netlib de "Plantillas para la solución de sistemas lineales", de Barrett et al.
• Richard S. Varga 2002 Análisis iterativo matricial , segunda ed. (de la edición de Prentice Hall de 1962), Springer-Verlag.
• David M. Young Jr. Solución iterativa de grandes sistemas lineales , Academic Press, 1971. (reimpreso por Dover, 2003)

enlaces externos

• Módulo para el Método SOR
• Solucionador de sistemas lineales tridiagonales basado en SOR, en C++