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

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

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

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

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

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

где:

$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 :

$A=D+L+U,$

где

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

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

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

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

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

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

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

$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.$

Это снова можно записать аналитически в матрично-векторной форме без необходимости обращения матрицы. $(D+\omega L)$:

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

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

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

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

Скорость конвергенции

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

• параметр релаксации подходит: $\omega \in (0,2)$
• Якоби Итерационная матрица $C_{\text{Jac}}:=I-D^{-1}A$ имеет только действительные собственные значения
• Метод Якоби сходится: $\mu :=\rho (C_{\text{Jac}})<1$
• матричное разложение $A=D+L+U$ удовлетворяет свойству, которое $\operatorname {det} (\lambda D+zL+{\tfrac {1}{z}}U)=\operatorname {det} (\lambda D+L+U)$ для любого $z\in \mathbb {C} \setminus \{0\}$ и $\lambda \in \mathbb {C}$.

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

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

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

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

В частности, для $\omega =1$ ( Гаусса-Зейделя ) считается, что $\rho (C_{\omega })=\mu ^{2}=\rho (C_{\text{Jac}})^{2}$. Для оптимального $\omega$ мы получаем $\rho (C_{\omega })={\frac {1-{\sqrt {1-\mu ^{2}}}}{1+{\sqrt {1-\mu ^{2}}}}}={\frac {\mu ^{2}}{4}}+O(\mu ^{3})$, который показывает, что SOR примерно в четыре раза эффективнее, чем метод Гаусса – Зейделя.

Последнее предположение выполнено для трехдиагональных матриц, поскольку $Z(\lambda D+L+U)Z^{-1}=\lambda D+zL+{\tfrac {1}{z}}U$ для диагонали $Z$ с записями $Z_{ii}=z^{i-1}$ и $\operatorname {det} (\lambda D+L+U)=\operatorname {det} (Z(\lambda D+L+U)Z^{-1})$.

Алгоритм

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

Inputs: A, b, ω
Output: φ

Choose an initial guess φ to the solution
repeat until convergence
    for i from 1 until n do
        set σ to 0
        for j from 1 until n do
            if j ≠ i then
                set σ to σ + aij φj
            end if
        end (j-loop)
        set φi to (1 − ω)φi + ω(bi − σ) / aii
    end (i-loop)
    check if convergence is reached
end (repeat)

Примечание

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

Пример

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

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

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

Итерация	$x_{1}$	$x_{2}$	$x_{3}$	$x_{4}$
0 1	0.25	−2.78125	1.6289062	0.5152344
0 2	1.2490234	−2.2448974	1.9687712	0.9108547
0 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
from scipy import linalg

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 = linalg.norm(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 = linalg.norm(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 , в которой

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

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

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

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

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

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

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

Основная статья: экстраполяция Ричардсона

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

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

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

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

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

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

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

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

См. также

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

Ссылки

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

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

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

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