successive over-relaxation#Other applications of the method

Given a square system of n linear equations with unknown x:

:$A\backslash mathbf\; x\; =\; \backslash mathbf\; b$

where:

:

Then A can be decomposed into a diagonal component D, and strictly lower and upper triangular components L and U:

:$A=D+L+U,$

where

:

The system of linear equations may be rewritten as:

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

for a constant ω > 1, called the relaxation factor.

The method of successive over-relaxation is an iterative technique that solves the left hand side of this expression for x, using the previous value for x on the right hand side. Analytically, this may be written as:

:$\backslash mathbf\{x\}^\{(k+1)\}\; =\; (D+\backslash omega\; L)^\{-1\}\; \backslash big(\backslash omega\; \backslash mathbf\{b\}\; -\; [\backslash omega\; U\; +\; (\backslash omega-1)\; D\; ]\; \backslash mathbf\{x\}^\{(k)\}\backslash big)=L\_\{\backslash omega\}\; \backslash mathbf\{x\}^\{(k)\}+\backslash mathbf\{c\},$

where $\backslash mathbf\{x\}^\{(k)\}$ is the kth approximation or iteration of $\backslash mathbf\{x\}$ and $\backslash mathbf\{x\}^\{(k+1)\}$ is the next or k + 1 iteration of $\backslash mathbf\{x\}$.

However, by taking advantage of the triangular form of (D+ωL), the elements of x^(k+1) can be computed sequentially using forward substitution:

:

This can again be written analytically in matrix-vector form without the need of inverting the matrix $(D+\backslash omega\; L)$:{{cite book |last1=Törnig |first1=Willi |title=Numerische Mathematik für Ingenieure und Physiker |date=1979 |publisher=Springer Berlin, Heidelberg |isbn=978-3-642-96508-1 |page=180 |edition=1 |doi=10.1007/978-3-642-96508-1 |url=https://link.springer.com/book/10.1007/978-3-642-96508-1 |access-date=20 May 2024 |ref=numerical-methods-for-engineers-and-physicists |language=en}}

:$\backslash mathbf\{x\}^\{(k+1)\}\; =\; (1-\backslash omega)\backslash mathbf\{x\}^\{(k)\}\; +\; \backslash omega\; D^\{-1\}\; \backslash big(\; \backslash mathbf\{b\}\; -\; L\; \backslash mathbf\{x\}^\{(k+1)\}\; -\; U\; \backslash mathbf\{x\}^\{(k)\}\; \backslash big).$

[[File:Spectral Radius.svg|thumb|Spectral radius $\backslash rho(C\_\backslash omega)$ of the iteration matrix for the SOR method $C\_\backslash omega$.

The plot shows the dependence on the spectral radius of the Jacobi iteration matrix $\backslash mu\; :=\; \backslash rho(C\_\backslash text\{Jac\})$.]]

The choice of relaxation factor ω is not necessarily easy, and depends upon the properties of the coefficient matrix.

In 1947, Ostrowski proved that if $A$ is symmetric and positive-definite then for .

Thus, convergence of the iteration process follows, but we are generally interested in faster convergence rather than just convergence.

The convergence rate for the SOR method can be analytically derived.

One needs to assume the following{{Cite book|title=Iterative Solution of Large Sparse Systems of Equations {{!}} SpringerLink|volume = 95|last=Hackbusch|first=Wolfgang|year=2016|isbn=978-3-319-28481-1|location=|pages=|language=en-gb|chapter=4.6.2|doi=10.1007/978-3-319-28483-5|series = Applied Mathematical Sciences}}{{Cite book|title=Iterative Methods for Solving Linear Systems|series = Frontiers in Applied Mathematics|volume = 17|last=Greenbaum|first=Anne|year=1997|isbn=978-0-89871-396-1|language=en-gb|chapter=10.1|doi=10.1137/1.9781611970937}}

• the relaxation parameter is appropriate: $\backslash omega\; \backslash in\; (0,2)$
• Jacobi's iteration matrix $C\_\backslash text\{Jac\}:=\; I-D^\{-1\}A$ has only real eigenvalues
• Jacobi's method is convergent:
• the matrix decomposition $A=D+L+U$ satisfies the property that $\backslash operatorname\{det\}(\backslash lambda\; D\; +\; zL\; +\; \backslash tfrac\{1\}\{z\}U)\; =\; \backslash operatorname\{det\}(\backslash lambda\; D\; +\; L\; +\; U)$ for any $z\backslash in\backslash mathbb\{C\}\backslash setminus\backslash \{0\backslash \}$ and $\backslash lambda\backslash in\backslash mathbb\{C\}$.

Then the convergence rate can be expressed as

:$$

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

where the optimal relaxation parameter is given by

:$$

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

In particular, for $\backslash omega\; =\; 1$ (Gauss-Seidel) it holds that $\backslash rho(C\_\backslash omega)=\backslash mu^2=\backslash rho(C\_\backslash text\{Jac\})^2$.

For the optimal $\backslash omega$ we get $\backslash rho(C\_\backslash omega)=\backslash frac\{1-\backslash sqrt\{1-\backslash mu^2\}\}\{1+\backslash sqrt\{1-\backslash mu^2\}\}\; =\; \backslash frac\{\backslash mu^2\}\{4\}\; +\; O(\backslash mu^3)$, which shows SOR is roughly four times more efficient than Gauss–Seidel.

The last assumption is satisfied for tridiagonal matrices since $Z(\backslash lambda\; D\; +\; L\; +\; U)Z^\{-1\}=\backslash lambda\; D\; +\; zL\; +\; \backslash tfrac\{1\}\{z\}U$ for diagonal $Z$ with entries $Z\_\{ii\}=z^\{i-1\}$ and $\backslash operatorname\{det\}(\backslash lambda\; D\; +\; L\; +\; U)\; =\; \backslash operatorname\{det\}(Z(\backslash lambda\; D\; +\; L\; +\; U)Z^\{-1\})$.

Since elements can be overwritten as they are computed in this algorithm, only one storage vector is needed, and vector indexing is omitted. The algorithm goes as follows:

Inputs: {{mvar|A}}, {{mvar|b}}, {{mvar|ω}}

Output: {{mvar|φ}}

Choose an initial guess {{mvar|φ}} to the solution

repeat until convergence

for {{mvar|i}} from 1 until {{mvar|n}} do

set {{mvar|σ}} to 0

for {{mvar|j}} from 1 until {{mvar|n}} do

if {{mvar|j}} ≠ {{mvar|i}} then

set {{mvar|σ}} to {{math|σ + a_ij φ_j}}

end if

end ({{mvar|j}}-loop)

set {{math|φ_i}} to {{math|(1 − ω)φ_i + ω(b_i − σ) / a_ii}}

end ({{mvar|i}}-loop)

check if convergence is reached

end (repeat)

;Note: $(1-\backslash omega)\backslash phi\_i\; +\; \backslash frac\{\backslash omega\}\{a\_\{ii\}\}\; (b\_i\; -\; \backslash sigma)$ can also be written $\backslash phi\_i\; +\; \backslash omega\; \backslash left(\; \backslash frac\{b\_i\; -\; \backslash sigma\}\{a\_\{ii\}\}\; -\; \backslash phi\_i\backslash right)$, thus saving one multiplication in each iteration of the outer for-loop.

We are presented the linear system

:$$

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

To solve the equations, we choose a relaxation factor $\backslash omega\; =\; 0.5$ and an initial guess vector $\backslash phi\; =\; (0,\; 0,\; 0,\; 0)$. According to the successive over-relaxation algorithm, the following table is obtained, representing an exemplary iteration with approximations, which ideally, but not necessarily, finds the exact solution, {{nowrap|(3, −2, 2, 1)}}, in 38 steps.

A simple implementation of the algorithm in Common Lisp is offered below.

;; 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+))))))

A simple Python implementation of the pseudo-code provided above.

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

1. 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)

The version for symmetric matrices A, in which

:$U=L^T,\backslash ,$

is referred to as Symmetric Successive Over-Relaxation, or (SSOR), in which

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

and the iterative method is

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

The SOR and SSOR methods are credited to David M. Young Jr.

{{main|Richardson extrapolation}}

A similar technique can be used for any iterative method. If the original iteration had the form

:$x\_\{n+1\}=f(x\_n)$

then the modified version would use

:$x^\backslash mathrm\{SOR\}\_\{n+1\}=(1-\backslash omega)x^\{\backslash mathrm\{SOR\}\}\_n+\backslash omega\; f(x^\backslash mathrm\{SOR\}\_n).$

However, the formulation presented above, used for solving systems of linear equations, is not a special case of this formulation if {{mvar|x}} is considered to be the complete vector. If this formulation is used instead, the equation for calculating the next vector will look like

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

where $L\_*\; =\; L\; +\; D$. Values of $\backslash omega>1$ are used to speed up convergence of a slow-converging process, while values of are often used to help establish convergence of a diverging iterative process or speed up the convergence of an overshooting process.

There are various methods that adaptively set the relaxation parameter $\backslash omega$ based on the observed behavior of the converging process. Usually they help to reach a super-linear convergence for some problems but fail for the others.

• Jacobi method
• Gaussian Belief Propagation
• Matrix splitting

{{Reflist}}

• {{CFDWiki|name=Successive_over-relaxation_method_-_SOR}}
• Abraham Berman, Robert J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, 1994, SIAM. {{isbn|0-89871-321-8}}.
• {{MathWorld|urlname=SuccessiveOverrelaxationMethod|title=Successive Overrelaxation Method|author=Black, Noel|author2=Moore, Shirley|name-list-style=amp}}
• A. Hadjidimos, [http://www.sciencedirect.com/science/article/pii/S0377042700004039 Successive overrelaxation (SOR) and related methods], Journal of Computational and Applied Mathematics 123 (2000), 177–199.
• Yousef Saad, [http://www-users.cs.umn.edu/%7Esaad/books.html Iterative Methods for Sparse Linear Systems], 1st edition, PWS, 1996.
• [http://www.netlib.org/utk/papers/templates/node11.html 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)

• [https://web.archive.org/web/20040317151200/http://math.fullerton.edu/mathews/n2003/SORmethodMod.html Module for the SOR Method]
• [https://web.archive.org/web/20140310122757/http://michal.is/projects/tridiagonal-system-solver-sor-c/ Tridiagonal linear system solver] based on SOR, in C++

{{Numerical linear algebra}}

Category:Numerical linear algebra

Category:Articles with example pseudocode

Category:Articles with example Python (programming language) code

Category:Relaxation (iterative methods)