טיוטה:שיטת התייעלות חוזרת (SOR)

באלגברה ליניארית נומרית, שיטת ההתייעלות החוזרת (SOR) היא וריאציה של שיטת גאוס-זיידל לפתרון מערכת משוואות ליניארית, מה שמוביל להתכנסות מהירה יותר. שיטה דומה יכולה לשמש עבור כל הליך איטרטיבי המתכנס באיטיות.

היא פותחה במקביל על ידי דייוויד מ. יאנג ג'וניור וסטנלי פ. פרנקל בשנת 1950 במטרה לפתור מערכות ליניאריות באופן אוטומטי במחשבים דיגיטליים. שיטות התייעלות חוזרת כבר היו בשימוש לפני עבודתם של יאנג ופרנקל. דוגמה לכך היא השיטה של לואיס פריי ריצ'רדסון, והשיטות שפיתחו RV Southwell. עם זאת, שיטות אלו תוכננו לחישוב על ידי מחשבים אנושיים, ודורשו מומחיות מסוימת כדי להבטיח התכנסות לפתרון, מה שהפך אותן ללא מתאימות לתכנות במחשבים דיגיטליים. היבטים אלה נדונים בתזה של דייוויד מ. יאנג ג'וניור.

נוסחה

נתונה מערכת משוואות מרובעת של 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:

$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 או האיטרציה ה-k של $\mathbf {x}$, ו-$\mathbf {x} ^{(k+1)}$ היא האיטרציה הבאה או k + 1 של $\mathbf {x}$. עם זאת, תוך ניצול הצורה של המטריצה משולשת של (D+ωL), ניתן לחשב את רכיבי x ^{( k +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)

הפניות

• תבנית:CFDWiki
• Abraham Berman, Robert J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, 1994, SIAM. מסת"ב 0-89871-321-8ISBN 0-89871-321-8.
• Black, Noel & Moore, Shirley. "Successive Overrelaxation Method". MathWorld.
• A. Hadjidimos, Successive overrelaxation (SOR) and related methods, Journal of Computational and Applied Mathematics 123 (2000), 177–199.
• Yousef Saad, Iterative Methods for Sparse Linear Systems, 1st edition, PWS, 1996.
• Netlib's copy of "Templates for the Solution of Linear Systems", by Barrett et al.
• Richard S. Varga 2002 Matrix Iterative Analysis, Second ed. (of 1962 Prentice Hall edition), Springer-Verlag.
• David M. Young Jr. Iterative Solution of Large Linear Systems, Academic Press, 1971. (reprinted by Dover, 2003)