Relajación excesiva sucesiva - Successive over-relaxation

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 resulta en 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 relajación excesiva 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 ser computados 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.

Formulación

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

{\ 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 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 se puede reescribir 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 sobre-relajació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 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_ {w} \ 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, al aprovechar la forma triangular de ( D + ωL ), los elementos de x ⁽^k⁺¹⁾ se pueden calcular secuencialmente usando 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 )} \ derecha), \ quad i = 1,2, \ ldots, n.}

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 positivo-definido, entonces para . Por lo tanto, sigue la convergencia del proceso de iteración, pero generalmente estamos interesados en una convergencia más rápida en lugar de solo una convergencia. ${\ Displaystyle A}$ ${\ Displaystyle \ rho (L _ {\ omega}) <1}$ ${\ Displaystyle 0 <\ omega <2}$

Tasa de convergencia

La tasa de convergencia para el método SOR se puede derivar analíticamente. Uno necesita asumir lo siguiente

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}}: = ID ^ {- 1} A}$
El método de Jacobi es convergente: ${\ Displaystyle \ mu: = \ rho (C _ {\ text {Jac}}) <1}$
existe una solución única: . ${\ Displaystyle \ det A \ neq 0}$

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} \ ,.}

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:

Inputs:  $A$ ,  $b$ ,  $ω$ 
Output:  $\phi$

Choose an initial guess  $\phi$  to the solution
repeat until convergence
    for  $i$  from 1 until  $n$  do
         $\sigma \leftarrow 0$ 
        for  $j$  from 1 until  $n$  do
            if  $j$  ≠  $i$  then
                 $\sigma \leftarrow \sigma +a_{ij}\phi _{j}$ 
            end if
        end ( $j$ -loop)
         $\phi _{i}\leftarrow (1-\omega )\phi _{i}+{\frac {\omega }{a_{ii}}}(b_{i}-\sigma )$ 
    end ( $i$ -loop)
    check if convergence is reached
end (repeat)

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 {alineado} 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 {alineado}}}

Para resolver las ecuaciones, elegimos un factor de relajación y un vector de aproximación inicial . Según 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)}$

Iteración	${\ Displaystyle x_ {1}}$	${\ Displaystyle x_ {2}}$	${\ Displaystyle x_ {3}}$	${\ Displaystyle x_ {4}}$
1	0,25	-2,78125	1.6289062	0.5152344
2	1.2490234	-2,2448974	1.9687712	0.9108547
3	2.070478	-1.6696789	1.5904881	0,76172125
...	...	...	...	...
37	2.9999998	-2,0	2.0	1.0
38	3,0	-2,0	2.0	1.0

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

;; Set the default floating-point format to "long-float" in order to
;; ensure correct operation on a wider range of numbers.
(setf *read-default-float-format* 'long-float)

(defparameter +MAXIMUM-NUMBER-OF-ITERATIONS+ 100
  "The number of iterations beyond which the algorithm should cease its
   operation, regardless of its current solution. A higher number of
   iterations might provide a more accurate result, but imposes higher
   performance requirements.")

(declaim (type (integer 0 *) +MAXIMUM-NUMBER-OF-ITERATIONS+))

(defun get-errors (computed-solution exact-solution)
  "For each component of the COMPUTED-SOLUTION vector, retrieves its
   error with respect to the expected EXACT-SOLUTION vector, returning a
   vector of error values.
   ---
   While both input vectors should be equal in size, this condition is
   not checked and the shortest of the twain determines the output
   vector's number of elements.
   ---
   The established formula is the following:
     Let resultVectorSize = min(computedSolution.length, exactSolution.length)
     Let resultVector     = new vector of resultVectorSize
     For i from 0 to (resultVectorSize - 1)
       resultVector[i] = exactSolution[i] - computedSolution[i]
     Return resultVector"
  (declare (type (vector number *) computed-solution))
  (declare (type (vector number *) exact-solution))
  (map '(vector number *) #'- exact-solution computed-solution))

(defun is-convergent (errors &key (error-tolerance 0.001))
  "Checks whether the convergence is reached with respect to the
   ERRORS vector which registers the discrepancy betwixt the computed
   and the exact solution vector.
   ---
   The convergence is fulfilled if and only if each absolute error
   component is less than or equal to the ERROR-TOLERANCE, that is:
   For all e in ERRORS, it holds: 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)
  "Creates and returns a vector of the SIZE with all elements set to 0."
  (declare (type (integer 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+))))
  "Implements the successive over-relaxation (SOR) method, applied upon
   the linear equations defined by the matrix A and the right-hand side
   vector B, employing the relaxation factor OMEGA, returning the
   calculated solution vector.
   ---
   The first algorithm step, the choice of an initial guess PHI, is
   represented by the optional keyword parameter PHI, which defaults
   to a zero-vector of the same structure as B. If supplied, this
   vector will be destructively modified. In any case, the PHI vector
   constitutes the function's result value.
   ---
   The terminating condition is implemented by the CONVERGENCE-CHECK,
   an optional predicate
     lambda(iteration phi) => generalized-boolean
   which returns T, signifying the immediate termination, upon achieving
   convergence, or NIL, signaling continuant operation, otherwise. In
   its default configuration, the CONVERGENCE-CHECK simply abides the
   iteration's ascension to the ``+MAXIMUM-NUMBER-OF-ITERATIONS+'',
   ignoring the achieved accuracy of the vector PHI."
  (declare (type (array  number (* *)) A))
  (declare (type (vector number *)     b))
  (declare (type number                omega))
  (declare (type (vector number *)     phi))
  (declare (type (function ((integer 1 *)
                            (vector number *))
                           *)
                 convergence-check))
  (let ((n (array-dimension A 0)))
    (declare (type (integer 0 *) n))
    (loop for iteration from 1 by 1 do
      (loop for i from 0 below n by 1 do
        (let ((rho 0))
          (declare (type number rho))
          (loop for j from 0 below n by 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)
      ;; Check if convergence is reached.
      (when (funcall convergence-check iteration phi)
        (return))))
  (the (vector number *) phi))

;; Summon the function with the exemplary parameters.
(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)))
  (successive-over-relaxation
    A b omega
    :convergence-check
    #'(lambda (iteration phi)
        (declare (type (integer 0 *)     iteration))
        (declare (type (vector number *) phi))
        (let ((errors (get-errors phi exact-solution)))
          (declare (type (vector number *) errors))
          (format T "~&~d. errors   = ~a" iteration errors)
          (or (is-convergent errors :error-tolerance 0.0)
              (>= iteration +MAXIMUM-NUMBER-OF-ITERATIONS+))))))

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

import numpy as np

def sor_solver(A, b, omega, initial_guess, convergence_criteria):
    """
    This is an implementation of the pseudo-code provided in the Wikipedia article.
    Arguments:
        A: nxn numpy matrix.
        b: n dimensional numpy vector.
        omega: relaxation factor.
        initial_guess: An initial solution guess for the solver to start with.
        convergence_criteria: The maximum discrepancy acceptable to regard the current solution as fitting.
    Returns:
        phi: solution vector of dimension n.
    """
    step = 0
    phi = initial_guess[:]
    residual = np.linalg.norm(np.matmul(A, phi) - b)  # Initial residual
    while residual > convergence_criteria:
        for i in range(A.shape[0]):
            sigma = 0
            for j in range(A.shape[1]):
                if j != i:
                    sigma += A[i, j] * phi[j]
            phi[i] = (1 - omega) * phi[i] + (omega / A[i, i]) * (b[i] - sigma)
        residual = np.linalg.norm(np.matmul(A, phi) - b)
        step += 1
        print("Step {} Residual: {:10.6g}".format(step, residual))
    return phi

# An example case that mirrors the one in the Wikipedia article
residual_convergence = 1e-8
omega = 0.5  # Relaxation factor

A = np.array([[4, -1, -6, 0],
              [-5, -4, 10, 8],
              [0, 9, 4, -2],
              [1, 0, -7, 5]])

b = np.array([2, 21, -12, -6])

initial_guess = np.zeros(4)

phi = sor_solver(A, b, omega, initial_guess, residual_convergence)
print(phi)

Relajación excesiva sucesiva simétrica

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

{\ Displaystyle U = L ^ {T}, \,}

se conoce como sobre-relajación sucesiva simétrica , o ( SSOR ), en el 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 acreditan 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 usa esta formulación en su lugar, la ecuación para calcular el siguiente vector se verá como

{\ 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 a menudo se utilizan 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 basándose en el comportamiento observado del proceso convergente. Por lo general, ayudan a alcanzar una convergencia súper lineal para algunos problemas, pero fallan para otros. ${\ Displaystyle \ omega}$

Ver también

Notas

Referencias

Este artículo incorpora texto del artículo Successive_over-relajación_method _-_ SOR en CFD-Wiki que está 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 , primera edición, PWS, 1996.
La copia de Netlib de "Plantillas para la solución de sistemas lineales", de Barrett et al.
Richard S. Varga 2002 Matrix Iterative Analysis , 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 sistema lineal tridiagonal basado en SOR, en C ++

Languages

In other projects