5.4 Relaxation methods

Gauss–Seidel moves each component from \(x_i^{(k-1)}\) to the value \(\tilde x_i\) given by (5.3). When the iterates approach the solution steadily from one side, moving a little further than that, \(x_i^{(k)} = x_i^{(k-1)}+\omega (\tilde x_i-x_i^{(k-1)})\) with \(\omega >1\), speeds convergence. This is the successive over-relaxation or SOR method: \begin {equation} \label {eq:sor} x_i^{(k)} = (1-\omega )x_i^{(k-1)}+\frac {\omega }{a_{ii}}\left (b_i-\sum _{j<i}a_{ij}x_j^{(k)}-\sum _{j>i}a_{ij}x_j^{(k-1)}\right ). \end {equation} For \(\omega = 1\) it is Gauss–Seidel; \(\omega <1\) is under-relaxation, occasionally used to make a divergent iteration converge.

Theorem 5.4.1.

(i)
(Kahan) If \(a_{ii}\neq 0\) for all \(i\), SOR can converge only if \(0<\omega <2\).
(ii)
(Ostrowski–Reich) If \(A\) is positive definite and \(0<\omega <2\), SOR converges for every \(\mathbf x^{(0)}\).
(iii)
If \(A\) is positive definite and tridiagonal, then \(\rho (T_g) = \rho (T_j)^2<1\), and the choice of \(\omega \) that minimizes \(\rho (T_\omega )\) is \[\omega _{\text {opt}} = \frac {2}{1+\sqrt {1-\rho (T_j)^2}},\qquad \text {with}\qquad \rho (T_{\omega _{\text {opt}}}) = \omega _{\text {opt}}-1 .\]

We prove (i): the SOR iteration matrix is \(T_\omega = (D-\omega L)^{-1}\bigl ((1-\omega )D+\omega U\bigr )\), a product of triangular matrices, so \[\det T_\omega = \frac {\det \bigl ((1-\omega )D\bigr )}{\det D} = (1-\omega )^n .\] Since the determinant is the product of the eigenvalues, \(\rho (T_\omega )^n\geq |\det T_\omega | = |1-\omega |^n\), so \(\rho (T_\omega )\geq |1-\omega |\), which is at least \(1\) unless \(0<\omega <2\). Parts (ii) and (iii) are proved in the texts by Burden and Faires and by Atkinson.

Example 5.4.2. Solve \[4x_1+3x_2 = 24,\qquad 3x_1+4x_2-x_3 = 30,\qquad -x_2+4x_3 = -24\] (solution \((3,4,-5)^T\)) by Gauss–Seidel and by SOR, and find the optimal \(\omega \).

Solution. The matrix is positive definite and tridiagonal. Here \(T_j\) has eigenvalues \(0\) and \(\pm \sqrt {0.625}\), so \(\rho (T_j) = 0.7906\), \(\rho (T_g) = 0.625\) and \[\omega _{\text {opt}} = \frac {2}{1+\sqrt {1-0.625}} = 1.24 .\] From \(\mathbf x^{(0)} = (1,1,1)^T\):

\(k\) Gauss–Seidel SOR, \(\omega = 1.25\)
\(1\) \((5.2500, 3.8125, -5.0469)\) \((6.3125, 3.5195, -6.6501)\)
\(2\) \((3.1406, 3.8828, -5.0293)\) \((2.6223, 3.9585, -4.6004)\)
\(3\) \((3.0879, 3.9268, -5.0183)\) \((3.1333, 4.0103, -5.0967)\)
\(5\) \((3.0343, 3.9714, -5.0072)\) \((3.0037, 4.0029, -5.0057)\)
\(7\) \((3.0134, 3.9888, -5.0028)\) \((3.0001, 4.0003, -5.0003)\)

Starting from \(\mathbf 0\) and stopping when successive iterates differ by less than \(10^{-7}\), Jacobi needs \(79\) iterations, Gauss–Seidel \(35\) and SOR with \(\omega = 1.25\) only \(15\) (Figure 7).

  0−−−51015       −7     (0)
Ek11111234JGSlIω111112571ωG1ir00000000aaOot.0.2.4.6.85050o5tercuRge0p–Sroos, srt:arbiscaa = 3ti–Sωlt5onpe =e:io1.seid sn2re1ts4tol.2rai5in10tgheeter ldainefrtedoiosmn m xeaagna li=inn0esatr ω convergence

Figure 7: Iterative methods on the system of Example 5.4.2. Left: the error \(\|\mathbf x^{(k)}-\mathbf x\|_\infty \) falls geometrically, by the factor \(\rho (T)\) per step: \(0.79\) for Jacobi, \(0.625\) for Gauss–Seidel and about \(0.25\) for SOR. Right: the number of SOR iterations as a function of \(\omega \) has a sharp minimum near the value \(\omega _{\text {opt}}\) predicted by Theorem 5.4.1.

Questions on this section

Stuck on something here? Ask below and it stays attached to this topic.