Relajación excesiva sucesiva

Método de resolución de un sistema de ecuaciones lineales

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

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 sobre-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 , requiriendo cierta experiencia para asegurar la convergencia a la solución, lo que los hizo 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}}.}$

Luego A se puede descomponer en un componente diagonal D y 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 puede reescribirse como:

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

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

El método de relajación sucesiva es una técnica iterativa que resuelve el lado izquierdo de esta expresión para x , utilizando el valor anterior de x en el lado derecho. Analíticamente, esto se puede escribir 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 utilizando sustitución hacia adelante : ${\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 puede escribirse analíticamente en forma de matriz-vector sin 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 del 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étrica y definida positiva, entonces para . Por lo tanto, se deduce la convergencia del proceso de iteración, pero en general nos interesa una convergencia más rápida en lugar de una simple 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 solo 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 cualquier 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 viene 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 demuestra 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 satisface 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 funciona de la siguiente manera:

Entradas: A , b , ω
Salida: φElija una estimación inicial φ para la solución repita hasta la convergencia para  i  desde 1 hasta  n  haga establezca σ en 0 para  j  desde 1 hasta  n  haga  si  j ≠ i  entonces establezca σ en σ + a ij  φ j  fin si  fin ( j -loop) establezca φ i en (1 − ω ) φ i + ω ( b i − σ ) / a ii  end ( i -bucle) comprobar si se alcanza la convergenciafin (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 aproximación inicial . De acuerdo con el algoritmo de sobre-relajació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 "long-float" para ;; garantizar el funcionamiento correcto en un rango más amplio de números. ( setf *read-default-float-format* 'long-float )  ( defparameter +MAXIMUM-NUMBER-OF-ITERATIONS+ 100 "El número de iteraciones más allá del cual el algoritmo debe dejar de  funcionar, independientemente de su solución actual. Un mayor número de  iteraciones puede proporcionar un resultado más preciso, pero impone mayores  requisitos de rendimiento." )   ( declaim ( tipo ( entero 0 * ) + NÚMERO MÁXIMO DE ITERACIONES + ))     ( defun get-errors ( calculated-solution exact-solution ) "Para cada componente del vector COMPUTED-SOLUTION, recupera su  error con respecto al vector EXACT-SOLUTION esperado, devolviendo un  vector de valores de error.  ---  Si bien ambos vectores de entrada 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, exactSolution.length)  Let resultVector = new vector of resultVectorSize  For i from 0 to (resultVectorSize - 1)  resultVector[i] = exactSolution[i] - calculatedSolution[i]  Return resultVector" ( declare ( type ( vector number * ) calculated-solution )) ( declare ( type ( vector number * ) exact-solution )) ( map ' ( vector number * ) #' - exact-solution calculated-solution ))                       ( defun is-convergent ( errors &key ( error-tolerance 0.001 )) "Comprueba si se alcanza la convergencia con respecto al  vector ERRORS 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 ERROR-TOLERANCE, es decir:  Para todo e en ERRORS, se cumple: abs(e) <= errorTolerance." ( declare ( type ( vector number * ) errors )) ( declare ( type number error-tolerance )) ( flet (( error-is-acceptable ( error ) ( declare ( type number error )) ( <= ( abs error ) error-tolerance ))) ( every #' error-is-acceptable errors )))                              ( defun make-zero-vector ( size ) "Crea y devuelve un vector del TAMAÑO con todos los elementos establecidos en 0." ( declare ( type ( entire 0 * ) size )) ( make-array size :initial-element 0.0 :element-type 'number ))               ( defun successive-over-relaxation ( A b omega &key ( phi ( make-zero-vector ( length b ))) ( convergence-check #' ( lambda ( iteration phi ) ( declare ( ignore phi )) ( >= iteration +MAXIMUM-NUMBER-OF-ITERATIONS+ )))) "Implementa el método de sobre-relajación sucesiva (SOR), aplicado sobre  las ecuaciones lineales definidas por la matriz A y el vector del lado derecho  B, empleando el factor de relajación OMEGA, devolviendo el  vector solución calculado.  ---  El primer paso del algoritmo, la elección de una estimación inicial 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 es implementada por CONVERGENCE-CHECK,  un predicado opcional  lambda(iteration phi) => booleano generalizado  que devuelve T, lo que significa la terminación inmediata, al lograr  la convergencia, o NIL, lo que indica una operación continua, en caso contrario. En  su configuración predeterminada, CONVERGENCE-CHECK simplemente respeta la  ascensión de la iteración al ``+MÁXIMO-NÚMERO-DE-ITERACIONES+'',  ignorando la precisión lograda del vector PHI. ( declare ( type ( array number ( * * )) A )) ( declare ( type ( vector number * ) b )) ( declare ( type number omega )) ( declare ( type ( vector number * ) phi )) ( declare ( type ( function (( entero 1 * ) ( vector number * )) * ) convergence-check )) ( let (( n                                                         ( array-dimension A 0 ))) ( declare ( type ( whole 0 * ) n )) ( loop for iteration from 1 for 1 do ( loop for i from 0 below n for 1 do ( let (( rho 0 )) ( declare ( type number rho )) ( loop for j from 0 below n for 1 do ( when ( /= 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 )))))) ( format T "~&~d. solution = ~a" iteration phi ) ;; Comprueba si se alcanza la convergencia. ( cuando ( funcall convergence-check iteration phi ) ( return )))) ( el ( vector number * ) phi ))                                                                                                       ;; Invoca la función con los parámetros de ejemplo. ( let (( A ( make-array ( list 4 4 ) :initial-contents ' (( 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 ) ( exact-solution ( vector 3 -2 2 1 ))) ( sucesiva-sobre-relajación A b omega :convergence-check #' ( lambda ( iteración phi ) ( declare ( type ( entero 0 * ) iteración )) ( declare ( type ( vector número * ) phi )) ( let (( errors ( get-errors phi exact-solution ))) ( declare ( type ( vector número * ) errores )) ( format T "~&~d. errores = ~a" iteración errores ) ( or ( errores convergentes : tolerancia a errores 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 ,  initial_guess ,  convergence_criteria ): """  Esta es una implementación del pseudocódigo proporcionado en el artículo de Wikipedia.  Argumentos:  A: matriz numpy de nxn.  b: vector numpy de n dimensiones.  omega: factor de relajación.  initial_guess: una estimación de solución inicial con la que el solucionador debe comenzar.  convergence_criteria: la discrepancia máxima aceptable para considerar que la solución actual es adecuada.  Devuelve:  phi: vector de solución de dimensión n.  """ step = 0 phi = initial_guess [:] residual = linalg . norm ( A @ phi - b ) # Residuo inicial while residual > convergence_criteria : para i en rango ( A . shape [ 0 ]): sigma = 0 para j en rango ( A . shape [ 1 ]): si j != i : sigma += A [ i , j ] * phi [ j ] phi [ i ] = ( 1 - omega ) * phi [ i ] + ( omega / A [ i , i ]) * ( b [ i ] - sigma ) residual = linalg . norm ( A @ phi - b ) paso += 1 print ( "Paso {} Residual: {:10.6g} " . format ( 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 ,  estimación_inicial ,  convergencia_residual ) print ( phi )

Relajación excesiva sucesiva simétrica

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

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

se denomina Sobre-relajació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 se utiliza esta formulación en su lugar, la ecuación para calcular el siguiente vector será similar a

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

donde . 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 de sobreimpulso . ${\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 fallan para otros. ${\displaystyle \omega }$

Véase también

• Método Jacobi
• Propagación de creencias gaussianas
• División de matriz

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. Recuperado el 20 de mayo de 2024 .
3. Hackbusch, Wolfgang (2016). "4.6.2". Solución iterativa de grandes sistemas dispersos de ecuaciones | SpringerLink . Applied Mathematical Sciences. 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 . Frontiers in Applied Mathematics. Vol. 17. doi :10.1137/1.9781611970937. ISBN 978-0-89871-396-1.

Referencias

• Este artículo incorpora texto del artículo Método_de_relajación_sucesiva_-_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 . 
• Black, Noel y Moore, Shirley. "Método de sobrerelajación sucesiva". MathWorld .
• 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 Matrix Iterative Analysis , segunda edición (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++