Последовательное чрезмерное расслабление - Successive over-relaxation

В численной линейной алгебре метод последовательной сверхрелаксации ( SOR ) является вариантом метода Гаусса – Зейделя для решения линейной системы уравнений , что приводит к более быстрой сходимости. Подобный метод можно использовать для любого медленно сходящегося итеративного процесса .

Он был разработан одновременно Дэвидом М. Янгом-младшим и Стэнли П. Франкелем в 1950 году с целью автоматического решения линейных систем на цифровых компьютерах. Методы чрезмерной релаксации использовались до работ Янга и Франкеля. Примером может служить метод Льюиса Фрая Ричардсона и методы, разработанные Р. В. Саутвеллом . Однако эти методы были разработаны для вычислений с помощью калькуляторов , требующих некоторого опыта для обеспечения сходимости к решению, которое сделало их неприменимыми для программирования на цифровых компьютерах. Эти аспекты обсуждаются в диссертации Дэвида М. Янга-младшего.

Формулировка

Дана квадратная система из n линейных уравнений с неизвестным x :

${\ Displaystyle А \ mathbf {x} = \ mathbf {b}}$

где:

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

Тогда A можно разложить на диагональную компоненту D и строго нижнюю и верхнюю треугольные компоненты L и U :

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

где

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

Систему линейных уравнений можно переписать как:

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

для постоянной ω > 1, называемой фактором релаксации .

Метод последовательной чрезмерной релаксации - это итерационный метод, который решает левую часть этого выражения для x , используя предыдущее значение для x в правой части. Аналитически это можно записать как:

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

где - k- е приближение или итерация, а - следующая или k + 1 итерация . Однако, пользуясь преимуществом треугольной формы ( D + ωL ), элементы x ^{( k +1)} могут быть вычислены последовательно с использованием прямой подстановки : ${\ Displaystyle \ mathbf {х} ^ {(к)}}$${\ displaystyle \ mathbf {x}}$${\ Displaystyle \ mathbf {х} ^ {(к + 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.}$

Конвергенция

Спектральный радиус итерационной матрицы для метода SOR . На графике показана зависимость итерационной матрицы Якоби от спектрального радиуса . ${\ displaystyle \ rho (C _ {\ omega})}$${\ displaystyle C _ {\ omega}}$${\ displaystyle \ mu: = \ rho (C _ {\ text {Jac}})}$

Выбор коэффициента релаксации ω не обязательно прост и зависит от свойств матрицы коэффициентов. В 1947 году Островский доказал , что если есть симметричная и положительно определенная , то для . Таким образом, следует сходимость итерационного процесса, но нас обычно интересует более быстрая сходимость, а не просто сходимость. ${\ displaystyle A}$${\ displaystyle \ rho (L _ {\ omega}) <1}$${\ displaystyle 0 <\ omega <2}$

Скорость сходимости

Скорость сходимости для метода SOR может быть получена аналитически. Необходимо предположить следующее

• параметр релаксации подходящий: ${\ displaystyle \ omega \ in (0,2)}$
• Матрица итераций Якоби имеет только действительные собственные значения ${\ displaystyle C _ {\ text {Jac}}: = ID ^ {- 1} A}$
• Метод Якоби сходится: ${\ displaystyle \ mu: = \ rho (C _ {\ text {Jac}}) <1}$
• уникальное решение существует: . ${\ Displaystyle \ Det А \ neq 0}$

Тогда скорость сходимости можно выразить как

${\ displaystyle \ rho (C _ {\ omega}) = {\ begin {case} {\ 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 {case}}}$

где оптимальный параметр релаксации определяется выражением

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

Алгоритм

Поскольку элементы могут быть перезаписаны по мере их вычисления в этом алгоритме, требуется только один вектор хранения, а индексирование вектора опускается. Алгоритм следующий:

Inputs: A, b, ω
Output: ${\displaystyle \phi }$

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

Примечание

${\ displaystyle (1- \ omega) \ phi _ {i} + {\ frac {\ omega} {a_ {ii}}} (b_ {i} - \ sigma)}$ также может быть записан , что позволяет сэкономить одно умножение на каждой итерации внешнего цикла for . ${\ displaystyle \ phi _ {i} + \ omega \ left ({\ frac {b_ {i} - \ sigma} {a_ {ii}}} - \ phi _ {i} \ right)}$

Пример

Нам представлена ​​линейная система

${\ displaystyle {\ begin {align} 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 {выравнивается}}}$

Чтобы решить уравнения, мы выбираем коэффициент релаксации и вектор начального предположения . В соответствии с алгоритмом последовательной избыточной релаксации получается следующая таблица, представляющая примерную итерацию с приближениями, которая в идеале, но не обязательно, находит точное решение (3, -2, 2, 1) за 38 шагов. ${\ displaystyle \ omega = 0,5}$${\ Displaystyle \ phi = (0,0,0,0)}$

Итерация	${\ 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

Ниже предлагается простая реализация алгоритма на 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+))))))

Простая реализация псевдокода Python, представленного выше.

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)

Симметричное последовательное чрезмерное расслабление

Версия для симметричных матриц A , в которой

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

называется симметричной последовательной чрезмерной релаксацией , или ( SSOR ), в которой

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

и итерационный метод

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

Методы SOR и SSOR приписываются Дэвиду М. Янгу-младшему .

Другие применения метода

Подобный прием можно использовать для любого итеративного метода. Если бы исходная итерация имела вид

${\ Displaystyle х_ {п + 1} = е (х_ {п})}$

тогда измененная версия будет использовать

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

Однако представленная выше формулировка, используемая для решения систем линейных уравнений, не является частным случаем этой формулировки, если x считается полным вектором. Если вместо этого использовать эту формулировку, уравнение для вычисления следующего вектора будет выглядеть так:

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

где . Значения используются для ускорения сходимости медленно сходящегося процесса, в то время как значения часто используются, чтобы помочь установить сходимость расходящегося итеративного процесса или ускорить сходимость процесса превышения . ${\ Displaystyle L _ {*} = L + D}$${\ displaystyle \ omega> 1}$${\ displaystyle \ omega <1}$

Существуют различные методы, которые адаптивно устанавливают параметр релаксации на основе наблюдаемого поведения процесса схождения. Обычно они помогают достичь суперлинейной сходимости для одних проблем, но не помогают для других. ${\ displaystyle \ omega}$

Смотрите также

• Метод Якоби
• Распространение гауссовской веры
• Расщепление матрицы

Заметки

Рекомендации

• Эта статья включает текст из статьи Successive_over-Relax_method _-_ SOR на CFD-Wiki, которая находится под лицензией GFDL .

• Абрахам Берман, Роберт Дж. Племмонс , Неотрицательные матрицы в математических науках , 1994, SIAM. ISBN   0-89871-321-8 .
• Блэк, Ноэль и Мур, Ширли. «Метод последовательной сверхрелаксации» . MathWorld .
• А. Хаджидимос, Последовательная избыточная релаксация (SOR) и родственные методы , Журнал вычислительной и прикладной математики 123 (2000), 177–199.
• Юсеф Саад , Итерационные методы для разреженных линейных систем , 1-е издание, PWS, 1996.
• Копия Netlib "Шаблоны для решения линейных систем" Барретта и др.
• Ричард С. Варга 2002 Матричный итерационный анализ , второе изд. (из издания Prentice Hall 1962 г.), Springer-Verlag.
• Дэвид М. Янг-младший. Итерационное решение больших линейных систем , Academic Press, 1971. (перепечатано Dover, 2003)

Внешние ссылки

• Модуль для метода SOR
• Решатель трехдиагональной линейной системы на основе SOR в C ++